// Functions for performing calculations to arbitrary precision.

// 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
         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[]
   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
   
   setPrecision[origPrec]
   return retval
}

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