[ynwwppvz] estimating e by counting primes

let n = 2^22 = 4194304.  there are p = primepi(n) = 295947 primes less than n.

using the prime number theorem n/log(n) ~= primepi(n), then e ~= n^(p/n) ~= 2.93 .

using n/(log(n) - 1) ~= primepi(n), then e ~= n^(p/(n+p)) ~= 2.732 .

(the same n was used previously computing e with primorial.  using the values of primes as we did then yielded slightly more precision (2.717) than using just their count as we do now.)

the typical way of computing fractional powers is with log and exp, but if you have such functions, then e = exp(1).

we chose n to be a power of 2 to permit easily computing the fractional power ^(p/n) in the first approximation without needing log and exp: the 1/n power in the can be computed by repeated square roots.  n^(p/n) = (n^p)^(1/n).  compute n^p first by augmenting machine precision floating point with an exponent that can become much larger than provided by the machine type.  the only operations one needs are multiplication and square root.  example implementation in Haskell.  there are some bits (pun!) of symmetry between exponentiation by squaring for the numerator of p/n followed by repeated square root for the denominator.

for the second approximation n^(p/(n+p)), one could have attempted to choose a different n such that n+p is a power of two.  but 3916608 + primepi(3916608) = 2^22-1 and 3916609 + primepi(3916609) = 2^22+1 (because 3916609 is prime), so it would have required some fudging.

better is to approximate the fractional exponent p/(n+p) as round(2^64 * p/(n+p))/2^64 .  (64 bits should be enough for everyone.)  for n = 2^22, p = 295947, this works out to to 1215802539408625636 / 2^64.

or, use Newton's method to compute the (n+p)th root of x^p.

comparing with methods that don't involve primes, but inputs of a similar scale:

e ~= (1+1/n)^n : for n=2^22, e ~= 2.718281504 .

sum(n=0, 10, 1.0/n!) ~= 2.71828180 , choosing 10 because 10! is the largest factorial less than 2^22 .