home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e031 / 3.ddi / MATHZIP2 / STARTUP / PRIMEPI.M < prev    next >
Encoding:
Text File  |  1991-09-23  |  4.2 KB  |  112 lines
  1.  
  2. Begin["System`"]
  3.  
  4. Unprotect[ PrimePi]
  5.  
  6. Clear[ PrimePi]
  7.  
  8. (*
  9.       PrimePi[x] gives the number of primes less than or equal x
  10.       (pi(x) in number theory terms). The prime number theorem states 
  11.       that this is aymptotically  LogIntegral[x]. The Riemann Hypothesis 
  12.       asserts that the absolute error is  < constant * Sqrt[x] Log[x] 
  13.       (good luck!).
  14. *)  
  15.  
  16. PrimePi::usage =
  17.     "PrimePi[x] gives the number of primes less than or equal to x."
  18.  
  19. Begin["`Private`"]
  20.  
  21. rlist = N[{0.6079271018540267, 0.2079768431451769,
  22.     0.05132991127342167, 0.01004570146280482,
  23.     0.001638254320440967, 0.0002295647814440843, 
  24.     0.00002822956919290008, 3.093984476471857*10^-6,
  25.     3.05888207019471*10^-7, 2.7543707437655*10^-8,
  26.     2.276904083198473*10^-9, 1.739516287122248*10^-10, 
  27.     1.235235408505764*10^-11, 8.19313909898923*10^-13, 
  28.     5.098031245110584*10^-14, 2.987150519283226*10^-15, 
  29.     1.653792071899792*10^-16, 8.67732064659652*10^-18, 
  30.     4.326646002346259*10^-19, 2.05515783148352*10^-20,
  31.     9.3204459029594*10^-22, 4.043995605352677*10^-23,
  32.     1.681813017415686*10^-24, 6.71557301275289*10^-26,
  33.     2.578780075325957*10^-27, 9.536908633654*10^-29, 
  34.     3.401366603549346*10^-30, 1.171389011057329*10^-31, 
  35.     3.899987198591185*10^-33, 1.256662542353454*10^-34, 
  36.     3.922984004097956*10^-36, 1.18762211075385*10^-37, 
  37.     3.489798672758064*10^-39, 9.96222804536054*10^-41, 
  38.     2.765026559568875*10^-42, 7.467278517408549*10^-44, 
  39.     1.963637886250443*10^-45, 5.031482118517906*10^-47, 
  40.     1.257043527310022*10^-48, 3.064043597819572*10^-50, 
  41.     7.291002017418843*10^-52, 1.694620650307293*10^-53, 
  42.     3.849327599400236*10^-55, 8.54964291189126*10^-57, 
  43.     1.857700188262819*10^-58, 3.950685655568404*10^-60, 
  44.     8.22686917863953*10^-62, 1.678224181406505*10^-63, 
  45.     3.355050425135873*10^-65, 6.575898833266313*10^-67,
  46.     1.264109733422975*10^-68, 2.384230636263747*10^-70, 
  47.     4.413670099171053*10^-72, 8.02210271797208*10^-74, 
  48.     1.432044782712371*10^-75, 2.511558132945803*10^-77, 
  49.     4.328939841334718*10^-79, 7.335005081928624*10^-81, 
  50.     1.222149654558633*10^-82, 2.002967489415538*10^-84, 
  51.     3.229724519347816*10^-86, 5.125213207081604*10^-88, 
  52.     8.00612796268731*10^-90, 1.231411283323488*10^-91, 
  53.     1.865333068229662*10^-93, 2.783440069672361*10^-95, 
  54.     4.09238237020218*10^-97, 5.929706289004022*10^-99, 
  55.     8.46922973434727*10^-101, 1.192605819734615*10^-102, 
  56.     1.656068386856241*10^-104, 2.268149218109434*10^-106, 
  57.     3.064491343664464*10^-108, 4.085242295242986*10^-110, 
  58.     5.374363197297437*10^-112, 6.978484068512949*10^-114, 
  59.     8.94526546140975*10^-116, 1.132125970625494*10^-117, 
  60.     1.414930711565271*10^-119, 1.746555097088381*10^-121, 
  61.     2.129620603064632*10^-123, 2.565426365976133*10^-125, 
  62.     3.053635680215457*10^-127, 3.592003422022151*10^-129, 
  63.     4.176170068510183*10^-131, 4.799546455156376*10^-133, 
  64.     5.453309487956776*10^-135, 6.126522797678712*10^-137, 
  65.     6.806388160531835*10^-139, 7.478624028238683*10^-141, 
  66.     8.12795752374691*10^-143, 8.73870669495473*10^-145, 
  67.     9.29542162025477*10^-147, 9.78354697469097*10^-149, 
  68.     1.019006554704655*10^-150, 1.050408232388696*10^-152, 
  69.     1.071731218081781*10^-154, 1.082444066575727*10^-156, 
  70.     1.082333624369158*10^-158, 1.071510288125467*10^-160}];
  71.  
  72. RiemannR[x_] := Block[{logx = Log[N[x]]}, 1+Round[rlist . (logx^Range[100])]]
  73.  
  74. PrimePi[x_] := PrimePi[Floor[x]] /; !IntegerQ[x] && IntegerQ[Floor[x]]
  75.  
  76. PrimePi[n_Integer] := Which[n < 2, 0, n < 3, 1, n < 5, 2, True, 3] /; n < 7
  77.  
  78. roundout[x_] := Sign[x] Ceiling[Abs[x]];
  79.  
  80. PrimePi[n_Integer] :=   
  81.         Block[{s = 1/Log[N[n]], t, diff, a, b = RiemannR[N[n]], sign},
  82.             diff = roundout[(Prime[b] - n) s];
  83.             sign = Sign[diff];
  84.             a = b - diff;
  85.             diff = Prime[a] - n; 
  86.             While[diff sign > 0,
  87.                 b = a;
  88.                 a -= roundout[diff s];
  89.                 diff = Prime[a] - n
  90.                 ];
  91.             If[diff == 0, Return[a]];
  92.             t = a;
  93.             If[a > b, diff = b; b = a; a = diff];
  94.             (* a and b are now lower and upper bounds *)
  95.             (* search with mixed bisection/Newton *)
  96.             While[True,
  97.                 (* Prime[a] < n < Prime[b] *)
  98.                 If[b == a+1, Return[a]];
  99.                 t = t - Round[diff s];
  100.                 If[Not[a < t < b], t = Round[(a+b)/2]];
  101.                 diff = Prime[t] - n;
  102.                 Which[    diff == 0, Return[t],
  103.                     diff  > 0, b = t,
  104.                     True, a = t
  105.                     ]
  106.                 ]
  107.             ] /; n > 6
  108. Protect[PrimePi]
  109.  
  110. End[] (* "`Private`" *)
  111. End[]; (* "System`" *)
  112.