pi.frink

View or download pi.frink in plain text format


// Program to calculate pi using binary splitting.
// See http://numbers.computation.free.fr/Constants/Algorithms/splitting.html

use karatsuba.frink

class Pi
{
   class var digitsPerIteration = 14.1816474627254776555

   class var largestDigits = undef

   class var cachePi = undef

   /** This is the main public method to get the value of pi to a certain
       number of digits, calculating it if need be.  If pi has already been
       calculated to a sufficient number of digits, this returns it from the
       cache.
   */

   class getPi[digits = getPrecision[]] :=
   {
      origPrec = getPrecision[]
      try
      {
         setPrecision[digits]
         if (largestDigits != undef) and (largestDigits >= digits)
            return 1. * cachePi
         else
            return 1. * calcPi[digits]
      }
      finally
         setPrecision[origPrec]
   }

   /** This is the main public method to get the value of 2 * pi to a certain
       number of digits, calculating it if need be.  If pi has already been
       calculated to a sufficient number of digits, this returns it from the
       cache.
   */

   class get2Pi[digits = getPrecision[]] :=
   {
      origPrec = getPrecision[]
      try
      {
         setPrecision[digits]
         if (largestDigits != undef) and (largestDigits >= digits)
            return 2 * cachePi
         else
            return 2 * calcPi[digits]
      }
      finally
         setPrecision[origPrec]
   }      

   /** This is an internal method that calculates the digits of pi if
       necessary. */

   class calcPi[digits] :=
   {
      oldPrec = getPrecision[]

      // Find number of terms to calculate
      k = max[floor[digits/digitsPerIteration], 1]

      try
      {
         setPrecision[digits+5]

         sum = 0

         p = p[0,k]
         q = q[0,k]

         sqC = karatsubaSqrt[640320, digits+8, true]

         piprime = p * 53360 / (q + 13591409 * p)
         piRational = piprime * sqC

         // Convert to floating-point
         setPrecision[digits]
         pi = 1. *  piRational

         largestDigits = digits
         cachePi = pi
         return pi
      }
      finally
         setPrecision[oldPrec]
   }

   /** Internal method for binary splitting. */
   class q[a,b] :=
   {
      if (b-a) == 1
         return (-1)^b * g[a,b] * (13591409 + 545140134 b)

      m = (a+b) div 2
      return q[a,m] p[m,b] + q[m,b] g[a,m]
   }

   /** Internal method for binary splitting. */
   class p[a,b] :=
   {
      if (b-a) == 1
         return 10939058860032000 b^3

      m = (a+b) div 2
      return p[a,m] p[m,b]
   }

   /** Internal method for binary splitting. */
   class g[a,b] :=
   {
      if (b-a) == 1
         return (6b-5)(2b-1)(6b-1)

      m = (a+b) div 2
      return g[a,m] g[m,b]
   }
}


View or download pi.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 17651 days, 5 hours, 3 minutes ago.