Stoneham's series

[Albert Chan]

(06-07-2024 07:55 PM)Johnh Wrote:  Jovial JRPN15 incremented to 2,664,645 and got to a Pi estimate of 3.141607465.

Jovial wins the test, much faster, doing over 1.3 million loops in a minute!, though not much extra convergence is achieved with this formula for all the extra cycles.

The problem is (1 - 1/b^2), when b goes big ... it is rounded to 1.
Doing product directly over-estimated what the algorithm should offer.

(1 - a) * (1 - b) = 1 - a - b + a*b = 1 - (a + (1-a)*b)

Code:

function p(n) -- 1 - product((1-1/b^2), b = 3 .. n step 2)
    local a = 0
    for b=3,n,2 do a = a + (1-a)/(b*b) end
    return a
end

lua> p(2664645)
0.21460168922870704
lua> (1-_) * 4
3.141593243085172

Or, we do product backwards, for more accuracy.

Code:

function q(n) -- 1 - product((1-1/b^2), b = n-(n+1)%2 .. 3 step -2)
    local a = 0
    for b=n-(n+1)%2,3,-2 do a = a + (1-a)/(b*b) end
    return a
end

lua> q(2664645)
0.2146016892287097
lua> (1-_) * 4
3.141593243085161


For reference, this is without rounding error (40 digits shown)

3.141593243085161166771739970151562977338...