we compute an approximation of Riemann zeta(2) = pi^2/6 using the Euler product formula over primes, then solve for pi:
? for(n=0, 10, z=1; forprime(p=1, 10^n, z*=(1-(1.0/p^2))); print("10^"n" ",sqrt(6/z)))
pi estimated from the product of all primes less than 10^n:
10^0 2.4494897427831780981972840747058913920
10^1 3.0935921676911454192536940842087145469
10^2 3.1387371950717209555231899287777680519
10^3 3.1413931847879282123372947308786365411
10^4 3.1415772351603451399115522305263190950
10^5 3.1415913934202453316750296375573376397
10^6 3.1415925471279879791648160810012226847
10^7 3.1415926443848380990491814377585237930
10^8 3.1415926527790129584805149679019884804
10^9 3.1415926535173431626391708456483210697
10^10 3.1415926535832448165410657025850891080
cpu time = 7min, 50,161 ms, real time = 7min, 50,326 ms.
each power of 10 yields slightly more than one more correct digit. asymptotically, probably just one digit, so this is not an efficient method of calculating pi, though it is nice to have gotten the digits above without too much effort.
previously, estimating e with primes.
No comments :
Post a Comment