Having noticed it was Pi day before Google alerted me to the fact, I'm feeling pretty dorky. Following some links I found the current record holder for the most digits of pi computed. Seeing his algorithm is recursive I thought "Hey I could write that in F# with bigint." So here I am on Pi day computing pi.

There is no square root method built in so I had to make my own. It takes about 30s on my i7 to compute the first 100,000 digits. I confirmed the results against the record holders data.

105790 // digits computed

"70150789337728658035712790913767420805655493624646" // digits 99951..100000

let A = 13591409I

let B = 545140134I

let C = 640320I

let a n = if n % 2I = 0I then (A + B * n) else -(A + B * n)

let p n = (2I*n - 1I) * (6I*n - 5I) * (6I*n - 1I)

let q n = n*n*n * C*C*C / 24I

let div2floor x =

if x % 2I = 0I then

x / 2I

else

(x - 1I) / 2I

let rec pqt n1 n2 =

if n1 + 1I = n2 then

(p n2, q n2, (a n2) * (p n2))

else

let m = div2floor (n1 + n2)

let P1, Q1, T1 = pqt n1 m

let P2, Q2, T2 = pqt m n2

(P1*P2, Q1 * Q2, T1*Q2 + P1*T2)

let sqrtC length =

let rec sqrt x y =

let e = y - x * x

let p = e / (2I * x)

if p = 0I then

x

else

sqrt (x + p) y

let guess = 800I * bigint.Pow(10I, length / 2)

let accuracy = bigint.Pow(10I, 2*(length / 2))

sqrt guess (C * accuracy)

let pi n =

let p, q, t = pqt 0I n

let dividend = q

let divisor = 12I * (t + q * A)

let sign =

if divisor < 0I then

-divisor

else

divisor

let length = int (bigint.Log10(sign))

C * (sqrtC length) * dividend / divisor

let test = string (pi 8000I)

printfn "%A" test.Length

printfn "%A" test.[99951..100000]