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 17592 days, 15 hours, 29 minutes ago.