### Pi Day

Sunday, March 14, 2010 by gradbot

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]

Filed under , , having