
Begin["System`"]

Unprotect[ PrimePi]

Clear[ PrimePi]

(*
      PrimePi[x] gives the number of primes less than or equal x
      (pi(x) in number theory terms). The prime number theorem states 
      that this is aymptotically  LogIntegral[x]. The Riemann Hypothesis 
      asserts that the absolute error is  < constant * Sqrt[x] Log[x] 
      (good luck!).
*)  

PrimePi::usage =
	"PrimePi[x] gives the number of primes less than or equal to x."

Begin["`Private`"]

rlist = N[{0.6079271018540267, 0.2079768431451769,
	0.05132991127342167, 0.01004570146280482,
	0.001638254320440967, 0.0002295647814440843, 
	0.00002822956919290008, 3.093984476471857*10^-6,
	3.05888207019471*10^-7, 2.7543707437655*10^-8,
	2.276904083198473*10^-9, 1.739516287122248*10^-10, 
	1.235235408505764*10^-11, 8.19313909898923*10^-13, 
	5.098031245110584*10^-14, 2.987150519283226*10^-15, 
	1.653792071899792*10^-16, 8.67732064659652*10^-18, 
	4.326646002346259*10^-19, 2.05515783148352*10^-20,
	9.3204459029594*10^-22, 4.043995605352677*10^-23,
	1.681813017415686*10^-24, 6.71557301275289*10^-26,
	2.578780075325957*10^-27, 9.536908633654*10^-29, 
	3.401366603549346*10^-30, 1.171389011057329*10^-31, 
	3.899987198591185*10^-33, 1.256662542353454*10^-34, 
	3.922984004097956*10^-36, 1.18762211075385*10^-37, 
	3.489798672758064*10^-39, 9.96222804536054*10^-41, 
	2.765026559568875*10^-42, 7.467278517408549*10^-44, 
	1.963637886250443*10^-45, 5.031482118517906*10^-47, 
	1.257043527310022*10^-48, 3.064043597819572*10^-50, 
	7.291002017418843*10^-52, 1.694620650307293*10^-53, 
	3.849327599400236*10^-55, 8.54964291189126*10^-57, 
	1.857700188262819*10^-58, 3.950685655568404*10^-60, 
	8.22686917863953*10^-62, 1.678224181406505*10^-63, 
	3.355050425135873*10^-65, 6.575898833266313*10^-67,
	1.264109733422975*10^-68, 2.384230636263747*10^-70, 
	4.413670099171053*10^-72, 8.02210271797208*10^-74, 
	1.432044782712371*10^-75, 2.511558132945803*10^-77, 
	4.328939841334718*10^-79, 7.335005081928624*10^-81, 
	1.222149654558633*10^-82, 2.002967489415538*10^-84, 
	3.229724519347816*10^-86, 5.125213207081604*10^-88, 
	8.00612796268731*10^-90, 1.231411283323488*10^-91, 
	1.865333068229662*10^-93, 2.783440069672361*10^-95, 
	4.09238237020218*10^-97, 5.929706289004022*10^-99, 
	8.46922973434727*10^-101, 1.192605819734615*10^-102, 
	1.656068386856241*10^-104, 2.268149218109434*10^-106, 
	3.064491343664464*10^-108, 4.085242295242986*10^-110, 
	5.374363197297437*10^-112, 6.978484068512949*10^-114, 
	8.94526546140975*10^-116, 1.132125970625494*10^-117, 
	1.414930711565271*10^-119, 1.746555097088381*10^-121, 
	2.129620603064632*10^-123, 2.565426365976133*10^-125, 
	3.053635680215457*10^-127, 3.592003422022151*10^-129, 
	4.176170068510183*10^-131, 4.799546455156376*10^-133, 
	5.453309487956776*10^-135, 6.126522797678712*10^-137, 
	6.806388160531835*10^-139, 7.478624028238683*10^-141, 
	8.12795752374691*10^-143, 8.73870669495473*10^-145, 
	9.29542162025477*10^-147, 9.78354697469097*10^-149, 
	1.019006554704655*10^-150, 1.050408232388696*10^-152, 
	1.071731218081781*10^-154, 1.082444066575727*10^-156, 
	1.082333624369158*10^-158, 1.071510288125467*10^-160}];

RiemannR[x_] := Block[{logx = Log[N[x]]}, 1+Round[rlist . (logx^Range[100])]]

PrimePi[x_] := PrimePi[Floor[x]] /; !IntegerQ[x] && IntegerQ[Floor[x]]

PrimePi[n_Integer] := Which[n < 2, 0, n < 3, 1, n < 5, 2, True, 3] /; n < 7

roundout[x_] := Sign[x] Ceiling[Abs[x]];

PrimePi[n_Integer] :=   
		Block[{s = 1/Log[N[n]], t, diff, a, b = RiemannR[N[n]], sign},
			diff = roundout[(Prime[b] - n) s];
			sign = Sign[diff];
			a = b - diff;
			diff = Prime[a] - n; 
			While[diff sign > 0,
				b = a;
				a -= roundout[diff s];
				diff = Prime[a] - n
				];
			If[diff == 0, Return[a]];
			t = a;
			If[a > b, diff = b; b = a; a = diff];
			(* a and b are now lower and upper bounds *)
			(* search with mixed bisection/Newton *)
			While[True,
				(* Prime[a] < n < Prime[b] *)
				If[b == a+1, Return[a]];
				t = t - Round[diff s];
				If[Not[a < t < b], t = Round[(a+b)/2]];
				diff = Prime[t] - n;
				Which[	diff == 0, Return[t],
					diff  > 0, b = t,
					True, a = t
					]
				]
			] /; n > 6
Protect[PrimePi]

End[] (* "`Private`" *)
End[]; (* "System`" *)
