gauss-legendre.md

20 million digits of pi in under a minute with Julia

I recently discovered a relatively obscure algorithm for calculating the digits of pi: https://en.wikipedia.org/wiki/Gauss–Legendre_algorithm. Well, at least obscure compared to Chudnovsky's. Wikipedia notes that it is "memory-intensive" but is it really? Let's compare to the MPFR pi function:

function gauss_legendre(prec)
    setprecision(BigFloat, prec, base=10)
    GC.enable(false)
    println(precision(BigFloat, base=2))
    x::BigFloat = a::BigFloat = BigFloat(1, precision=precision(BigFloat, base=2))
    b::BigFloat = BigFloat(a / sqrt(BigFloat(2, precision=precision(BigFloat, base=2))))
    t::BigFloat = BigFloat(a / 4, precision=precision(BigFloat, base=2))
    y::BigFloat = a
    for _=0:ceil(BigInt, log2(prec))
        y, a = a::BigFloat, (a::BigFloat + b::BigFloat) / 2
        b = sqrt(b::BigFloat*y::BigFloat)
        t = t::BigFloat-(x::BigFloat * (y::BigFloat-a::BigFloat)^2)
        x = x::BigFloat*2
    end
    (a + b)^2 / (4 * t)
end
function mpfr_pi(ndigits::Int64)
    setprecision(ndigits, base=10)
    BigFloat(π)
end
gauss_legendre(1)
@time begin
    my_pi = gauss_legendre(20000000)
end
mpfr_pi(1)
@time begin
    mpfrs_pi = mpfr_pi(20000000)
end
println(Float64((my_pi::BigFloat-mpfrs_pi::BigFloat)::BigFloat*20000000))

20 million digits is a fair amount! Let's see how they run:

 59.159281 seconds (477.69 k allocations: 18.049 GiB)
 42.238366 seconds (390.11 k allocations: 12.369 GiB)
-0.0

All benchmarks shown are run on my 2020 MBP 13' M1. That last number is the "error" (comparing our implementation to MPFR). Only ~17 seconds slower, and with about 6 more gigs of memory allocated. However--my algorithm is written in pure Julia, whereas MPFR is in C. Perhaps this is the new holy grail of pi-computation algorithms?

Oh, and I mentioned Chudnovsky's algorithm:

 32.878370 seconds (1.48 M allocations: 25.777 GiB, 0.01% compilation time)

That was for only 100k digits. Perhaps I'm missing something but why has no one set a world record with Gauss-Legendre? If anyone has a super powerful computer and wants to try this out, please post the results below. I wanna see how far you can push this.

This is really cool! that said none of your ::BigFloat annotations are doing anything here.

Can probably squeeze some more performance

function gauss_legendre2(prec)
    half = BigFloat(0.5)
    setprecision(BigFloat, prec, base = 10)
    prec₂ = precision(BigFloat, base = 2)
    a = y = BigFloat(1, prec₂)
    b = a / sqrt(BigFloat(2, prec₂))
    t = a / 4
    for i in Iterators.map(x -> 2^x, 0:ceil(BigInt, log2(prec)))
        y = a
        a += b
        a *= half
        t -= i * (y - a)^2
        b *= y
        b = √b
    end
    (a + b)^2 / 4t
end

Running the same 20 million benchmark,

122.178489 seconds (910 allocations: 2.297 GiB) 
 54.441351 seconds (936 allocations: 1.702 GiB) # gauss_legendre2