Little math problem(s) November 2020

[Albert Chan]

(11-21-2020 10:45 AM)grsbanks Wrote:  

(11-21-2020 10:10 AM)Paul Dale Wrote:  \(\pi=\sqrt{6\zeta(2)}\)

I used that here

http://techy.horwits.com/2020/09/using-r...ue-of.html

Instead of random number simulation, we can also do brute force.
Probability of a,b co-prime to each other (1 ≤ a,b ≤ n) = (number of co-primes pairs) / n^2

Number of co-prime pairs = sum(sum(gcd(x,y)==1, y=1 .. n), x=1 .. n)

We can simplify with symmetry: gcd(x,y) = gcd(y,x), gcd(x,x) = x:

XCAS> sum(sum(gcd(x,y)==1, y=1 .. x-1), x=1 .. n) * 2 + 1             (*)

+1 is due to under-counted gcd(1,1) = 1 (1 is co-prime to all integers)
We could also include the diagonal, and remove the over-counted gcd(1,1).

XCAS> sum(sum(gcd(x,y)==1, y=1 .. x), x=1 .. n) * 2 - 1

Inside sum is Euler totient function. (ref: "C++ for Mathematicians", page 67)

XCAS> P(n) := (sum(Phi(x), x=1 .. n) * 2. - 1) / n^2
XCAS> estimated_Pi(n) := sqrt(6./P(n))

XCAS> estimated_Pi(200)      → 3.13220921695, error = 0.2987%
XCAS> estimated_Pi(1e3)      → 3.14041534038, error = 0.0375%       (**)
XCAS> estimated_Pi(1e4)      → 3.14153423902, error = 0.0019%
XCAS> estimated_Pi(1e5)      → 3.14158477584, error = 0.0003%       // on my laptop, this take 2 seconds

(*) XCas sum(f(x), x=a .. b) had weird meaning if a > b+1
If a = b+1, result is 0
If a > b+1, result is -sum(f(x), x=b+1 .. a-1)       // see XCas sum limit confusion

(**) the blog wrongly stated for n=1000, error = 0.06%