(* :Name: NumberTheory`Ramanujan` *)

(* :Title: Ramanujan's tau-Dirichlet Series *)

(* :Author: Jerry B. Keiper *)

(* :Summary:
This package evaluates Ramanujan's tau function, Ramanujan's tau-Dirichlet
series, and a related function that is real along the critical line.
*)

(* :Context: NumberTheory`Ramanujan` *) 

(* :Package Version: 2.0 *)

(* :Copyright: Copyright 1990-1999,  Wolfram Research, Inc.
*)

(* :History:
	Written by Jerry B. Keiper, February 1991.
*)

(* :Keywords: Ramanujan, tau *)

(* :Source:
	G. H. Hardy, Ramanujan: Twelve Lectures on Subjects Suggested by
	His Life and Work, Chelsea, New York, 1959, MR 21 #4881.

	Robert Spira, Calculation of the Ramanujan tau-Dirichlet Series,
	Mathematics of Computation, vol. 27, num. 122, Apr. 1973, pp. 379-385.

	Hiroyuki Yoshida, On Calculations of Zeros of L-Functions Related
	with Ramanujan's Discriminant Function on the Critical Line, 
	J. Ramanujan Math. Soc. 3 (1) 1988, pp. 87-95.
*)

(* :Mathematica Version: 2.0 *)

(* :Limitation:
	Evaluation of the tau function itself is only practical for small
	integers (up to several thousand).
*)

(* :Discussion:
	RamanujanTauZ[t] is a real function for real t and is analogous to
	RiemannSiegelZ in the theory of Riemann's zeta function.
	A conjecture due to Ramanujan is that the nontrivial zeros of
	RamanujanTauZ are all real.

	The number of zeros in the critical strip from t == 0 to t == T
	is given by N[T] == (RamanujanTauTheta[T] + 
	        Im[Log[RamanujanTauDirichletSeries[6 + I T]]])/Pi
	where the logarithm is defined by analytic continuation from
	t == 10.  In fact, continuity from 12 + I T is sufficient, because
	RamanujanTauDirichletSeries[s] remains in the right half plane
	for Re[s] > 12.

	There are several ways to evaluate the Ramanujan tau-Dirichlet
	series.  Two of the methods (numerical integration or a sum
	involving incomplete Gamma functions) are not very practical,
	since precision loss is a major problem: (empirically) about
	t/2 digits are lost where t is the argument of RamanujanTauZ
	or the imaginary part of the argument of
	RamanujanTauDirichletSeries.  The method used here is iterated
	Abel summation of the Dirichlet series itself.  Since this requires
	knowledge of RamanujanTau[n] for rather large n, this method is
	only practical for t up to about 1000.  Precision loss is also
	a problem here, but not nearly so much as in the other methods.
*)

BeginPackage["NumberTheory`Ramanujan`"]

RamanujanTau::usage =
"RamanujanTau[n] gives the coefficient of x^n in the series expansion of
x Product[1 - x^k, {k, 1, Infinity}]^24."

RamanujanTauDirichletSeries::usage =
"RamanujanTauDirichletSeries[s] gives the value of the Dirichlet series
Sum[RamanujanTau[n] n^(-s), {n, 1, Infinity}]."

RamanujanTauGeneratingFunction::usage =
"RamanujanTauGeneratingFunction[z] evaluates the generating function of
RamanujanTau[n], i.e., Sum[RamanujanTau[n] z^n, {n, 1, Infinity}]."

RamanujanTauTheta::usage =
"RamanujanTauTheta[t] is analogous to RiemannSiegelTheta[t] in the theory
of the Riemann zeta function.  In particular, Exp[I RamanujanTauTheta[t]] *
RamanujanTauDirichletSeries[6 + I t] is RamanujanTauZ[t] and is a real
valued function if t is real."

RamanujanTauZ::usage =
"RamanujanTauZ[t] evaluates the function
Gamma[6 + I t] RamanujanTauDirichletSeries[6 + I t] (2 Pi)^(-I t)
Sqrt[Sinh[Pi t]/(Pi t Product[k^2 + t^2, {k, 5}])]."

Unprotect[RamanujanTauZ, RamanujanTauDirichletSeries, RamanujanTauTheta,
	RamanujanTauGeneratingFunction, RamanujanTau];

Begin["NumberTheory`Ramanujan`Private`"]

$RT = {{1}, {1}, {1}, {1}};

RamanujanTau[n_Integer] := 0 /; n <= 0

nextlistextension[j_, n_] :=
    Module[{i, rev = Reverse[$RT[[j]]]},
	Table[Take[$RT[[j]], i] . Take[rev, -{i, 1}],
		{i, Length[$RT[[j+1]]] + 1, Min[Length[$RT[[j]]], n]}]
	]

RamanujanTau[n_Integer] :=
    Module[{tmp, k, k0, k1},
	If[n > Length[$RT[[4]]],
	    If[n > Length[$RT[[1]]],
		k0 = Floor[Sqrt[2.0 (Length[$RT[[1]]]-1)]] + 1;
		k1 = Round[Sqrt[2. n]] + 1;
		tmp = Flatten[Table[Append[Table[0, {k-1}], (-1)^k (2k + 1)],
				{k, k0, k1}]];
		AbortProtect[$RT[[1]] = Join[$RT[[1]], tmp]]
		];
	    If[n > Length[$RT[[2]]],
		tmp = nextlistextension[1, n];
		AbortProtect[$RT[[2]] = Join[$RT[[2]], tmp]]
		];
	    If[n > Length[$RT[[3]]],
		tmp = nextlistextension[2, n];
		AbortProtect[$RT[[3]] = Join[$RT[[3]], tmp]]
		];
	    tmp = nextlistextension[3, n];
	    AbortProtect[$RT[[4]] = Join[$RT[[4]], tmp]]
	    ];
	$RT[[4, n]]
	] /; n > 0

RamanujanTauDirichletSeries[s_] :=
    RamanujanTauDirichletSeries[12-s] Gamma[12-s] (2Pi)^(2s-12)/
	Gamma[s] /; NumberQ[s] && Re[s] < 6 && Precision[s] < Infinity

as1[n_] = {};

as2[0, i_] = 0;

aslist = {};

    (* The only reason for the function as2 is that we want to cache
	values of the iterated partial sums of RamanujanTau.  The basic
	idea is rather simple, but the implementation is complicated, since
	the caching needs to be both extensible and sparse. *)

as2[n_, i_] :=
    Module[{tmp, j, k, np, n0},
	If[Length[as1[n]] < i,
	    n0 = If[n === 50, 0, 50 Round[n/Sqrt[5000.]]];
	    np = {n};
	    While[n0 > 0 && as1[n0] == {},
		PrependTo[np, n0];
		n0 = If[n0 === 50, 0, 50 Round[n0/Sqrt[5000.]]]
		];
	    If[0 == (j = Length[as1[n]]),
		aslist = Join[aslist, Table[RamanujanTau[k], {k, n0+1, n}]]];
	    While[j++ < i,
		tmp = as2[n0, j];
		AbortProtect[
	    	    Do[aslist[[k]] = (tmp += aslist[[k]]), {k, n0+1, n}];
		    Map[AppendTo[as1[#], aslist[[#]]]&, np]
		    ];
		While[n0 > 0 && Length[as1[n0]] < j,
		    PrependTo[np, n0];
		    n0 = If[n0 === 50, 0, 50 Round[n0/Sqrt[5000.]]]
		    ]
		]
	    ];
	as1[n][[i]]
	];

RamanujanTauDirichletSeries[ss_] :=
    Module[{r, oldr, n, tmp, i, u = {}, s, prec},
	(* iterated Abel summation on the Dirichlet series *)
	tmp = Log[8, Abs[Im[ss]] + 8.];
	prec = Precision[ss] + Floor[tmp] + 10;
	s = SetPrecision[ss, prec];
	n = prec (Abs[Im[s]] + 10) tmp/200;
	n = 50 Round[2.^(Ceiling[Log[2, n n]]/2)];
	r = Sum[RamanujanTau[i] i^(-s), {i, n, 1, -1}];
	oldr = 0;
	s = SetPrecision[ss, 2 prec];
	While[r != oldr,
	    oldr = r;
	    tmp = 0;
	    AppendTo[u, (n + 1 + Length[u])^(-s)];
	    Do[u[[i]] -= u[[i+1]], {i, Length[u]-1, 1, -1}];
	    r -= as2[n, Length[u]] First[u]
	    ];
	If[MachineNumberQ[Re[ss]+Im[ss]],
		N[r],
	    (* else *)
		SetPrecision[r, Min[Precision[ss], Precision[r]]]]
	] /; NumberQ[ss] && Precision[ss] < Infinity

RamanujanTauGeneratingFunction[z_] :=
    Module[{g},
	g /; NumberQ[g = ramg[-Log[z]]]
	] /; NumberQ[z] && Precision[z] < Infinity

RamanujanTauTheta[t_] :=
    Module[{prec = Accuracy[t]},
	If[ prec < 13 Log[10, Abs[Im[t]]] - 10,
	    N[(11*Pi)/4-15472974061/(156*t^13)+1742885234309/(360360*t^11)- 
		295081381/(1188*t^9)+23237999/(1680*t^7)-1115101/(1260*t^5)+ 
		26999/(360*t^3)-181/(12*t)+t(Log[t/(2Pi)] - 1), prec],
	  (* else *)
	    prec=N[(LogGamma[6+I t]-LogGamma[6-I t])/(2I) - t Log[2Pi], prec];
	    If[Head[t] === Real, Re[prec], prec]
	    ]
	] /; NumberQ[t] && Precision[t] < Infinity

RamanujanTauZ[t_] :=
    Module[{z, pi = N[Pi, Accuracy[t]], k, norm},
	z = RamanujanTauDirichletSeries[6 + I t] Gamma[6 + I t] (2pi)^(-I t);
	If[Head[t] =!= Complex, z = Re[z]];
	norm = 1/Product[k^2 + t^2, {k, 5}];
	If[t != 0.0, norm *= Sinh[pi t]/(pi t)];
	z Sqrt[norm]
	] /; NumberQ[t] && Precision[t] < Infinity

ramg[0.] = 0

ramg[y_] :=
    Module[{pi2 = N[2Pi, Precision[y]], ypr},
	ypr = pi2/y;
	If[Re[y] > Re[ypr], ramg0[y], (ypr)^12 ramg0[pi2 ypr]]
	]

ramg0[y_] :=
    Module[{x = Exp[-y], k},
	x NProduct[1 - x^k, {k, 1, Infinity}, Method -> SequenceLimit,
		VerifyConvergence -> False, WorkingPrecision -> Precision[x],
		NProductFactors -> 1 + Floor[(5 - Precision[x])/Log[10, x]],
		NProductExtraFactors -> 11, WynnDegree -> Infinity]^24
	] /; Re[y] > .1

End[ ]	(* "NumberTheory`Ramanujan`Private`" *)

Protect[RamanujanTauZ, RamanujanTauDirichletSeries, RamanujanTauTheta,
	RamanujanTauGeneratingFunction, RamanujanTau];

EndPackage[ ]	(* "NumberTheory`Ramanujan`" *)
