// Functions for performing calculations to arbitrary precision.

use root.frink
use pi.frink

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

arbitraryExp[x, digits=getPrecision[]] :=
{
   origPrec = getPrecision[]
   setPrecision[digits+3]

   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
         //println["term is $term"]
         ns = s + term
         xpower = xpower * x
         n = n + 1.
         factorial = factorial * n
         //println["s is $s"]
      }
   }
   setPrecision[digits]
   retval = 1.0 s

   setPrecision[origPrec]
   return retval
}


// Natural log to arbitrary precision.
arbitraryLn[x, 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)

      // 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] - 1.
      } while abs[r2-r] > eps
      setPrecision[digits]
      retval = 1. r
      
      return retval
   }
   finally
      setPrecision[origPrec]
}

// 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[]] :=
{
   origPrec = getPrecision[]
   try
   {
      setPrecision[digits+4]
      pi = Pi.getPi[digits+4]
      // Break up one period of a sinewave into octants, each with width pi/4
      octant = floor[(x / (pi/4)) mod 8]
//      println["Octant is $octant"]
      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 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.
arbitraryCos[x, digits=getPrecision[]] :=
{
   origPrec = getPrecision[]
   try
   {
      setPrecision[digits+4]
      pi = Pi.getPi[digits+4]
      arg = x+pi/2
      return arbitrarySin[x+pi/2, 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.
arbitraryTan[x, digits=getPrecision[]] :=
{
   origPrec = getPrecision[]
   try
   {
      setPrecision[digits]
      return arbitrarySin[x, digits+4] / arbitraryCos[x, digits+4]
   }
   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[]
*/