
(* :Name: NumberTheory`Ramanujan` *)

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

(* :Author: Jerry B. Keiper *)

(* :Summary:
	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: Mathematica 2.0 *)

(* :Copyright: Copyright 1990,  Wolfram Research, Inc.
	Permission is hereby granted to modify and/or make copies of
	this file for any purpose other than direct profit, or as part
	of a commercial product, provided this copyright notice is left
	intact.  Sale, other than for the cost of media, is prohibited.

	Permission is hereby granted to reproduce part or all of
	this file, provided that the source is acknowledged.
*)

(* :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}])]."

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`" *)
