ArbitraryPrecision.frink

View or download ArbitraryPrecision.frink in plain text format


// Functions for performing calculations to arbitrary precision.

use pi.frink

// See http://en.literateprograms.org/Arbitrary-precision_elementary_mathematical_functions_(Python)

arbitraryExp[x, digits=getPrecision[], debug=false] :=
{
   if debug
      println["in arbitraryExp[$x]"]
   origPrec = getPrecision[]
   setPrecision[digits+3]

   try
   {
      if x < 0
         s = 1.0 / arbitraryExp[abs[x],digits]
      else
      {
         xpower = 1.
         ns = 0.
         s = 1
         n = 0.
         factorial = 1.
         while s != ns
         {
            s = ns
            term = xpower / factorial
            if debug
               println["term is $term"]
            ns = s + term
            xpower = xpower * x
            n = n + 1.
            factorial = factorial * n
            if debug
               println["s is $s"]
         }
      }
      setPrecision[digits]
      retval = 1.0 s
      if debug
         println["arbitraryExp returning $retval"]
      return retval
   }
   finally
      setPrecision[origPrec]
}


// Natural log to arbitrary precision.
arbitraryLn[x, digits=getPrecision[], debug=false] :=
{
   if debug
      println["in arbitraryLn[$x]"]
   
   if x <= 0
      return "Error:  Arbitrary logs of negative numbers not yet implemented."
   
   origPrec = getPrecision[]
   try
   {
      setPrecision[digits+3]

      eps = 10.^-(digits+1)

      // A good initial estimate is needed
      if (x < 10^308)
         r2 = ln[x]
      else                         // TODO:  Store ln[2] somewhere.
         r2 = approxLog2[x] * ln[2]

      //   println["Epsilon is $eps"]

      // Use Newton's method
      do
      {
         r = r2
  //       println["r is $r"]
         r2 = r + x/arbitraryExp[r,digits+2,debug] - 1.
      } while abs[r2-r] > eps
      setPrecision[digits]
      retval = 1. r

      if debug
         println["arbitraryLn about to return $retval"]
      return retval
   }
   finally
      setPrecision[origPrec]
}

// Arbitrary-precision power  x^p
// This uses the relationship that x^p = exp[p * ln[x]]
arbitraryPow[x, p, digits = getPrecision[], debug=false ] :=
{
   // TODO:  Make this work much faster for integer and rational powers.
   
   prec = getPrecision[]

   try
   {
      workingdigits = digits + 2
      if digits <= 12
         workingdigits = digits + 4
      
      setPrecision[workingdigits]

      ret  = arbitraryExp[p * arbitraryLn[x, workingdigits, debug],
                          workingdigits, debug]

      setPrecision[digits]
      return 1. * ret
   }
   finally
      setPrecision[prec]
}

// Arbitrary log to the base 10.
arbitraryLog[x, digits=getPrecision[]] :=
{
   origPrec = getPrecision[]
   try
   {
      setPrecision[digits+2]
      // TODO:  Store ln[10] somewhere.
      retval = arbitraryLn[x, digits+2] / arbitraryLn[10, digits+2]
      setPrecision[digits]
      return 1. retval
   }
   finally
      setPrecision[origPrec]
}

// This method computes sine of a number to an arbitrary precision.
// This method is actually a dispatcher function which conditions the values
// and tries to dispatch to the appropriate method which will be most likely
// to converge rapidly.
arbitrarySin[x, digits=getPrecision[], debug=false] :=
{
   origPrec = getPrecision[]
   try
   {
      // If x is something like 10^50, we actually need to work with
      // 50+digits at this point to get a meaningful result.
      extradigits = max[0, ceil[approxLog2[abs[x]]/ 3.219]]   // Approx. log 10

      if debug
         println["Extradigits is " + extradigits]

      setPrecision[digits+extradigits+4]

      if debug
         println["Dividing by pi to " + (digits + extradigits + 4) + " digits"]
      
      pi = Pi.getPi[digits+extradigits+4]
      
      // Break up one period of a sinewave into octants, each with width pi/4
      octant = floor[(x / (pi/4)) mod 8]

      // Adjust x into [0, 2 pi]
      x = x mod (2 pi)

      if debug
         println["Octant is $octant"]

      if debug
         println["Adjusted value is $x"]
      
      if octant == 0
         val = arbitrarySinTaylor[x, digits]
      else
         if octant == 1
            val = arbitraryCosTaylor[-(x - pi/2), digits]
         else
            if octant == 2
               val = arbitraryCosTaylor[x - pi/2, digits]
            else
               if octant == 3 or octant == 4
                  val = -arbitrarySinTaylor[x-pi, digits]
               else
                  if octant == 5
                     val = -arbitraryCosTaylor[-(x - 3/2 pi), digits]
                  else
                     if octant == 6
                        val = -arbitraryCosTaylor[x - 3/2 pi, digits]
                     else
                        val = arbitrarySinTaylor[x - 2 pi, digits]

      setPrecision[digits]
      return 1. * val            
   }
   finally
      setPrecision[origPrec]
}

/* This method computes cosine of a number to an arbitrary precision.
   This method actually just calls arbitrarySin[x + pi/2]
*/

arbitraryCos[x, digits=getPrecision[]] :=
{
   origPrec = getPrecision[]

   // If x is something like 10^50, we actually need to work with
   // 50+digits at this point to get a meaningful result.
   extradigits = max[0, ceil[approxLog2[abs[x]]/ 3.219]]   // Approx. log 10
   
   if debug
      println["Extradigits is " + extradigits]

   if debug
      println["Dividing by pi to " + (digits + extradigits + 4) + " digits"]
      
   pi = Pi.getPi[digits+extradigits+4]
      
   try
   {
      setPrecision[digits+extradigits+4]
      pi = Pi.getPi[digits+extradigits+4]
      arg = x+pi/2
      return arbitrarySin[arg, digits]
   }
   finally
      setPrecision[origPrec]
}

// Arbitrary-precision sine
arbitrarySinTaylor[x, digits=getPrecision[], returnInterval = false] :=
{
   origPrec = getPrecision[]
   eps = 10.^-(digits+3)
   terms = new array

   try
   {
      setPrecision[digits+5]
      x = x / radian               // Factor out radians if we use them
      pi = Pi.getPi[digits+5]
      x = x mod (2 pi)
      if x > pi
         x = x - 2 pi
      num = x
      sum = x
      term = 3
      denom = 1
      factor = -x*x
      terms.push[sum]
      do
      {
         prevSum = sum
         num = num * factor
         denom = denom * (term-1) * term
         part = num/denom
         sum = sum + part
         term = term + 2
         terms.push[part]
      } while prevSum != sum

//      println["terms for sin is $term"]
      sum = sum[reverse[terms]]
      setPrecision[digits]
      return 1. * sum
   }
   finally
      setPrecision[origPrec]
}

// Cosine for arbitrary digits.  We could write this in terms of the sine
// function (cos[x] = sin[x + pi/2]) but it's faster and more accurate
// (especially around pi/2) to write it as the Taylor series expansion.
arbitraryCosTaylor[x, digits=getPrecision[]] :=
{
   origPrec = getPrecision[]
   eps = 10.^-(digits+4)
   terms = new array

   try
   {
      setPrecision[digits+4]
      x = x / radian               // Factor out radians if we use them
      pi = Pi.getPi[digits+4]
      x = x mod (2 pi)
//      println["Effective x is $x"]
      num = 1
      sum = 1
      term = 2
      denom = 1
      factor = -x*x
      terms.push[sum]
      do
      {
         prevSum = sum
         num = num * factor
         denom = denom * (term-1) * term
         part = num/denom
         sum = sum + part
//         println["$term $part $sum"]
         term = term + 2
         terms.push[part]
      } while prevSum != sum

//      println["terms for cos is $term"]
      sum = sum[reverse[terms]]
      setPrecision[digits]
      return 1. * sum
   }
   finally
      setPrecision[origPrec]
}


// Tangent for arbitrary digits.  This is written in terms of
// sin[x]/cos[x] but it seems to behave badly around pi/2 where cos[x] goes to
// zero.
//
// TODO:  Make this a series expansion with the tangent numbers.  This might
// be more efficient.
//  See:  http://mathworld.wolfram.com/TangentNumber.html
//  also  TangentNumbers.frink
//            which calculate these numbers directly and efficiently.
arbitraryTan[x, digits=getPrecision[], debug=false] :=
{
   // If x is something like 10^50, we actually need to work with
   // 50+digits at this point to get a meaningful result.
   extradigits = max[0, ceil[approxLog2[abs[x]]/ 3.219]]   // Approx. log 10
   
   if debug
      println["Extradigits is " + extradigits]

   origPrec = getPrecision[]
   try
   {
      setPrecision[digits+extradigits+4]
      retval = arbitrarySin[x, digits+4] / arbitraryCos[x, digits+4]
      setPrecision[digits]
      return 1. * retval
   }
   finally
      setPrecision[origPrec]
}

// Polylogarithm.  See:
// https://en.wikipedia.org/wiki/Polylogarithm
polylog[s, z, digits = getPrecision[]] :=
{
//   if x <= 0
//      return "Error:  Arbitrary logs of negative numbers not yet implemented."
   
   origPrec = getPrecision[]
   try
   {
      setPrecision[digits+3]

      eps = 10.^-(digits+1)

      sum = 1. * z
      term = 0
      k = 2
      
      do
      {
         term = z^k / k^s
         sum = sum + term
         k = k + 1
//         println[sum]
      } while abs[term] > eps
      setPrecision[digits]
      retval = 1. sum
      
      return retval
   }
   finally
      setPrecision[origPrec]
}

// Binary logarithm (that is, logarithm to base 2.)
binaryLog[x, digits = getPrecision[]] :=
{
   origPrec = getPrecision[]
   try
   {
      setPrecision[digits+3]
      x = 1. x
      y = 0
      b = .5
      while x < 1
      {
         x = 2 x
         y = y - 1
      }

      while x >= 2
      {
         x = x / 2
         y = y + 1
      }

      setPrecision[15]
      epsilon = y * 10.^-(digits+3)
      setPrecision[digits+3]
      
      println["Epsilon is $epsilon"]
      do
      {
         x = x^2
         if x >= 2
         {
            x = x/2
            y = y + b
         }

         b = b/2

         //println["$y $x $b"]
         
      } while b > epsilon

      setPrecision[digits]
      return 1. * y
   }
   finally
      setPrecision[origPrec]
}



/*
digits = 400
setPrecision[digits]
pi = Pi.getPi[digits]
num = pi/4
collapseIntervals[false]
inum = new interval[num, num, num]

println[arbitrarySin[num,digits]]
println[arbitrarySin[-num,digits]]
println[sin[num]]
println[sin[inum]]
println[]

println[arbitraryCos[num,digits]]
//println[arbitraryCosAround2Pi[num,digits]]
println[arbitraryCos[-num,digits]]
println[cos[num]]
println[cos[inum]]
println[]

println[arbitraryTan[num,digits]]
println[arbitraryTan[-num,digits]]
println[tan[num]]
println[tan[-num]]
println[tan[inum]]

setPrecision[100]
g = new graphics
ps = new polyline
pc = new polyline
for x=-20 pi to 20 pi step .1
{
   ps.addPoint[x,-arbitrarySin[x]]
   pc.addPoint[x,-arbitraryCos[x]]
}

g.add[ps]
g.color[0,0,1]
g.add[pc]
g.show[]
*/


View or download ArbitraryPrecision.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 17360 days, 18 hours, 57 minutes ago.