home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMBERTH.PAK / RAMANUJA.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  8.3 KB  |  253 lines
  1.  
  2. (* :Name: NumberTheory`Ramanujan` *)
  3.  
  4. (* :Title: Ramanujan's tau-Dirichlet Series *)
  5.  
  6. (* :Author: Jerry B. Keiper *)
  7.  
  8. (* :Summary:
  9.     Evaluates Ramanujan's tau function, Ramanujan's tau-Dirichlet series,
  10.     and a related function that is real along the critical line.
  11. *)
  12.  
  13. (* :Context: NumberTheory`Ramanujan` *) 
  14.  
  15. (* :Package Version: Mathematica 2.0 *)
  16.  
  17. (* :Copyright: Copyright 1990,  Wolfram Research, Inc.
  18.     Permission is hereby granted to modify and/or make copies of
  19.     this file for any purpose other than direct profit, or as part
  20.     of a commercial product, provided this copyright notice is left
  21.     intact.  Sale, other than for the cost of media, is prohibited.
  22.  
  23.     Permission is hereby granted to reproduce part or all of
  24.     this file, provided that the source is acknowledged.
  25. *)
  26.  
  27. (* :History:
  28.     Written by Jerry B. Keiper, February, 1991.
  29. *)
  30.  
  31. (* :Keywords: Ramanujan, tau *)
  32.  
  33. (* :Source:
  34.     G. H. Hardy, Ramanujan: Twelve Lectures on Subjects Suggested by
  35.     His Life and Work, Chelsea, New York, 1959, MR 21 #4881
  36.  
  37.     Robert Spira, Calculation of the Ramanujan tau-Dirichlet Series,
  38.     Mathematics of Computation, vol. 27, num. 122, Apr. 1973, pp. 379-385
  39.  
  40.     Hiroyuki Yoshida, On Calculations of Zeros of L-Functions Related
  41.     with Ramanujan's Discriminant Function on the Critical Line, 
  42.     J. Ramanujan Math. Soc. 3 (1) 1988, pp. 87-95
  43. *)
  44.  
  45. (* :Mathematica Version: 2.0 *)
  46.  
  47. (* :Limitation:
  48.     Evaluation of the tau function itself is only practical for small
  49.     integers (up to several thousand).
  50. *)
  51.  
  52. (* :Discussion:
  53.     RamanujanTauZ[t] is a real function for real t and is analogous to
  54.     RiemannSiegelZ in the theory of Riemann's zeta function.
  55.     A conjecture due to Ramanujan is that the nontrivial zeros of
  56.     RamanujanTauZ[ ] are all real.
  57.  
  58.     The number of zeros in the critical strip from t == 0 to t == T
  59.     is given by N[T] == (RamanujanTauTheta[T] + 
  60.             Im[Log[RamanujanTauDirichletSeries[6 + I T]]])/Pi
  61.     where the logarithm is defined by analytic continuation from
  62.     t == 10.  In fact, continuity from 12 + I T is sufficient, because
  63.     RamanujanTauDirichletSeries[s] remains in the right half plane
  64.     for Re[s] > 12.
  65.  
  66.     There are several ways to evaluate the Ramanujan tau-Dirichlet
  67.     series.  Two of the methods (numerical integration or a sum
  68.     involving incomplete Gamma functions) are not very practical,
  69.     since precision loss is a major problem: (empirically) about
  70.     t/2 digits are lost where t is the argument of RamanujanTauZ[ ]
  71.     or the imaginary part of the argument of
  72.     RamanujanTauDirichletSeries[ ].  The method used here is iterated
  73.     Abel summation of the Dirichlet series itself.  Since this requires
  74.     knowledge of RamanujanTau[n] for rather large n this method is
  75.     only practical for t up to about 1000.  Precision loss is also
  76.     a problem here, but not nearly so much as in the other methods.
  77. *)
  78.  
  79. BeginPackage["NumberTheory`Ramanujan`"]
  80.  
  81. RamanujanTau::usage =
  82. "RamanujanTau[n] gives the coefficient of x^n in the series expansion of
  83. x Product[1 - x^k, {k, 1, Infinity}]^24."
  84.  
  85. RamanujanTauDirichletSeries::usage =
  86. "RamanujanTauDirichletSeries[s] gives the value of the Dirichlet series
  87. Sum[RamanujanTau[n] n^(-s), {n, 1, Infinity}]."
  88.  
  89. RamanujanTauGeneratingFunction::usage =
  90. "RamanujanTauGeneratingFunction[z] evaluates the generating function of
  91. RamanujanTau[n], i.e. Sum[RamanujanTau[n] z^n, {n, 1, Infinity}]."
  92.  
  93. RamanujanTauTheta::usage =
  94. "RamanujanTauTheta[t] is analogous to RiemannSiegelTheta[t] in the theory
  95. of the Riemann zeta function.  In particular, Exp[I RamanujanTauTheta[t]] *
  96. RamanujanTauDirichletSeries[6 + I t] is RamanujanTauZ[t] and is a real
  97. valued function if t is real."
  98.  
  99. RamanujanTauZ::usage =
  100. "RamanujanTauZ[t] evaluates the function
  101. Gamma[6 + I t] RamanujanTauDirichletSeries[6 + I t] (2 Pi)^(-I t)
  102. Sqrt[Sinh[Pi t]/(Pi t Product[k^2 + t^2, {k, 5}])]."
  103.  
  104. Begin["NumberTheory`Ramanujan`Private`"]
  105.  
  106. $RT = {{1}, {1}, {1}, {1}};
  107.  
  108. RamanujanTau[n_Integer] := 0 /; n <= 0
  109.  
  110. nextlistextension[j_, n_] :=
  111.     Module[{i, rev = Reverse[$RT[[j]]]},
  112.     Table[Take[$RT[[j]], i] . Take[rev, -{i, 1}],
  113.         {i, Length[$RT[[j+1]]] + 1, Min[Length[$RT[[j]]], n]}]
  114.     ]
  115.  
  116. RamanujanTau[n_Integer] :=
  117.     Module[{tmp, k, k0, k1},
  118.     If[n > Length[$RT[[4]]],
  119.         If[n > Length[$RT[[1]]],
  120.         k0 = Floor[Sqrt[2.0 (Length[$RT[[1]]]-1)]] + 1;
  121.         k1 = Round[Sqrt[2. n]] + 1;
  122.         tmp = Flatten[Table[Append[Table[0, {k-1}], (-1)^k (2k + 1)],
  123.                 {k, k0, k1}]];
  124.         AbortProtect[$RT[[1]] = Join[$RT[[1]], tmp]]
  125.         ];
  126.         If[n > Length[$RT[[2]]],
  127.         tmp = nextlistextension[1, n];
  128.         AbortProtect[$RT[[2]] = Join[$RT[[2]], tmp]]
  129.         ];
  130.         If[n > Length[$RT[[3]]],
  131.         tmp = nextlistextension[2, n];
  132.         AbortProtect[$RT[[3]] = Join[$RT[[3]], tmp]]
  133.         ];
  134.         tmp = nextlistextension[3, n];
  135.         AbortProtect[$RT[[4]] = Join[$RT[[4]], tmp]]
  136.         ];
  137.     $RT[[4, n]]
  138.     ] /; n > 0
  139.  
  140. RamanujanTauDirichletSeries[s_] :=
  141.     RamanujanTauDirichletSeries[12-s] Gamma[12-s] (2Pi)^(2s-12)/
  142.     Gamma[s] /; NumberQ[s] && Re[s] < 6 && Precision[s] < Infinity
  143.  
  144. as1[n_] = {};
  145.  
  146. as2[0, i_] = 0;
  147.  
  148. aslist = {};
  149.  
  150.     (* The only reason for the function as2[ ] is that we want to cache
  151.     values of the iterated partial sums of RamanujanTau[ ].  The basic
  152.     idea is rather simple, but the implementation is complicated, since
  153.     the caching needs to be both extensible and sparse. *)
  154.  
  155. as2[n_, i_] :=
  156.     Module[{tmp, j, k, np, n0},
  157.     If[Length[as1[n]] < i,
  158.         n0 = If[n === 50, 0, 50 Round[n/Sqrt[5000.]]];
  159.         np = {n};
  160.         While[n0 > 0 && as1[n0] == {},
  161.         PrependTo[np, n0];
  162.         n0 = If[n0 === 50, 0, 50 Round[n0/Sqrt[5000.]]]
  163.         ];
  164.         If[0 == (j = Length[as1[n]]),
  165.         aslist = Join[aslist, Table[RamanujanTau[k], {k, n0+1, n}]]];
  166.         While[j++ < i,
  167.         tmp = as2[n0, j];
  168.         AbortProtect[
  169.                 Do[aslist[[k]] = (tmp += aslist[[k]]), {k, n0+1, n}];
  170.             Map[AppendTo[as1[#], aslist[[#]]]&, np]
  171.             ];
  172.         While[n0 > 0 && Length[as1[n0]] < j,
  173.             PrependTo[np, n0];
  174.             n0 = If[n0 === 50, 0, 50 Round[n0/Sqrt[5000.]]]
  175.             ]
  176.         ]
  177.         ];
  178.     as1[n][[i]]
  179.     ];
  180.  
  181. RamanujanTauDirichletSeries[ss_] :=
  182.     Module[{r, oldr, n, tmp, i, u = {}, s, prec},
  183.     (* iterated Abel summation on the Dirichlet series *)
  184.     tmp = Log[8, Abs[Im[ss]] + 8.];
  185.     prec = Precision[ss] + Floor[tmp] + 10;
  186.     s = SetPrecision[ss, prec];
  187.     n = prec (Abs[Im[s]] + 10) tmp/200;
  188.     n = 50 Round[2.^(Ceiling[Log[2, n n]]/2)];
  189.     r = Sum[RamanujanTau[i] i^(-s), {i, n, 1, -1}];
  190.     oldr = 0;
  191.     s = SetPrecision[ss, 2 prec];
  192.     While[r != oldr,
  193.         oldr = r;
  194.         tmp = 0;
  195.         AppendTo[u, (n + 1 + Length[u])^(-s)];
  196.         Do[u[[i]] -= u[[i+1]], {i, Length[u]-1, 1, -1}];
  197.         r -= as2[n, Length[u]] First[u]
  198.         ];
  199.     If[MachineNumberQ[Re[ss]+Im[ss]],
  200.         N[r],
  201.         (* else *)
  202.         SetPrecision[r, Min[Precision[ss], Precision[r]]]]
  203.     ] /; NumberQ[ss] && Precision[ss] < Infinity
  204.  
  205. RamanujanTauGeneratingFunction[z_] :=
  206.     Module[{g},
  207.     g /; NumberQ[g = ramg[-Log[z]]]
  208.     ] /; NumberQ[z] && Precision[z] < Infinity
  209.  
  210. RamanujanTauTheta[t_] :=
  211.     Module[{prec = Accuracy[t]},
  212.     If[ prec < 13 Log[10, Abs[Im[t]]] - 10,
  213.         N[(11*Pi)/4-15472974061/(156*t^13)+1742885234309/(360360*t^11)- 
  214.         295081381/(1188*t^9)+23237999/(1680*t^7)-1115101/(1260*t^5)+ 
  215.         26999/(360*t^3)-181/(12*t)+t(Log[t/(2Pi)] - 1), prec],
  216.       (* else *)
  217.         prec=N[(LogGamma[6+I t]-LogGamma[6-I t])/(2I) - t Log[2Pi], prec];
  218.         If[Head[t] === Real, Re[prec], prec]
  219.         ]
  220.     ] /; NumberQ[t] && Precision[t] < Infinity
  221.  
  222. RamanujanTauZ[t_] :=
  223.     Module[{z, pi = N[Pi, Accuracy[t]], k, norm},
  224.     z = RamanujanTauDirichletSeries[6 + I t] Gamma[6 + I t] (2pi)^(-I t);
  225.     If[Head[t] =!= Complex, z = Re[z]];
  226.     norm = 1/Product[k^2 + t^2, {k, 5}];
  227.     If[t != 0.0, norm *= Sinh[pi t]/(pi t)];
  228.     z Sqrt[norm]
  229.     ] /; NumberQ[t] && Precision[t] < Infinity
  230.  
  231. ramg[0.] = 0
  232.  
  233. ramg[y_] :=
  234.     Module[{pi2 = N[2Pi, Precision[y]], ypr},
  235.     ypr = pi2/y;
  236.     If[Re[y] > Re[ypr], ramg0[y], (ypr)^12 ramg0[pi2 ypr]]
  237.     ]
  238.  
  239. ramg0[y_] :=
  240.     Module[{x = Exp[-y], k},
  241.     x NProduct[1 - x^k, {k, 1, Infinity}, Method -> SequenceLimit,
  242.         VerifyConvergence -> False, WorkingPrecision -> Precision[x],
  243.         NProductFactors -> 1 + Floor[(5 - Precision[x])/Log[10, x]],
  244.         NProductExtraFactors -> 11, WynnDegree -> Infinity]^24
  245.     ] /; Re[y] > .1
  246.  
  247. End[ ]    (* "NumberTheory`Ramanujan`Private`" *)
  248.  
  249. Protect[RamanujanTauZ, RamanujanTauDirichletSeries, RamanujanTauTheta,
  250.     RamanujanTauGeneratingFunction, RamanujanTau];
  251.  
  252. EndPackage[ ]    (* "NumberTheory`Ramanujan`" *)
  253.