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 17592 days, 6 hours, 22 minutes ago.