### 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 = 13591409Ilet B = 545140134Ilet C = 640320Ilet 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 / 24Ilet div2floor x =    if x % 2I = 0I then        x / 2I    else        (x - 1I) / 2Ilet 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 / divisorlet test = string (pi 8000I)printfn "%A" test.Lengthprintfn "%A" test.[99951..100000]`

Filed under , , having