
(*********************************************************************

        Adapted from
        Roman E. Maeder: Programming in Mathematica,
        Second Edition, Addison-Wesley, 1991.

 *********************************************************************)


(* the function is called "MyPrimePi" to avoid naming conflicts
   with the built-in function PrimePi *)

MyPrimePi::usage = "MyPrimePi[x] returns the number of primes <= x."

Attributes[MyPrimePi] = {Listable}

Begin["`Private`"]

MyPrimePi[x_] := 0 /; x < 2

MyPrimePi[x_] :=
	Module[{li, n0, n1, m, nx = N[x]},
		li = LogIntegral[nx];
		n0 = Floor[li - LogIntegral[Sqrt[nx]]];
		n1 = Ceiling[li];
		While[ n1 - n0 > 1,
			m = Floor[(n0+n1)/2];                 (* midpoint *)
			If[ Prime[m] <= nx, n0 = m, n1 = m ]
		];
		n0
	] /; x >= 2

End[]
