/** Functions for performing calculations to arbitrary precision.
    A good reference is Ronald W. Potter,
    "Arbitrary Precision Calculation of Selected Higher
    Functions."  References to "Potter" and equation numbers are from this book.

    https://www.lulu.com/shop/ronald-w-potter/arbitrary-precision-calculation-of-selected-higher-functions/paperback/product-1nwy7p7q.html?q=potter+functions&page=1&pageSize=4

    Also see Henrik Vestermark's excellent site which is a treasure trove:
    http://www.hvks.com/

    especially

    http://www.hvks.com/Numerical/papers.html
*/

use pi2.frink


/** This is a new arbitrary precision exponentiation function with a binary
    splitting implementation.  It is orders of magnitude faster than the
    previous implementation.

    http://www.hvks.com/

    especially

    http://www.hvks.com/Numerical/papers.html
    http://www.hvks.com/Numerical/arbitrary_precision.html

    and specifically outlined in the paper:
    http://www.hvks.com/Numerical/Downloads/HVE%20Fast%20Exp()%20calculation%20for%20arbitrary%20precision.pdf
*/
arbitraryExp[x, digits=getPrecision[], debug=false] :=
{
   if x == 0
      return 1
   
   origPrec = getPrecision[]
   try
   {
      setPrecision[18]
      prec = digits + 2 + ceil[log[digits]]
      c1 = 1
      c2 = 2
      v = x
      xp = 1

      // e^-x == 1/e^x
      if realSignum[v] == -1
         v = negate[v]

      // Automatically calculate optimal reduction factor as a power of two
      r = 8 * ceil[ln[2] * ln[prec]]
      rexp = floor[log[v,2]]
      r = r + rexp + 1       // r += v.exponent() + 1
      r = max[0, r]
      // Adjust the precision
      // println["rexp is $rexp"]
      prec = prec + floor[log[digits]] * r
      // println["Final precision is $prec"]

      // Calculate needed Taylor terms
      k = xstirlingApprox[prec, -(r-rexp)]
      if k < 2
         k = 2 // Minimum 2 terms otherwise it can't split

      setPrecision[prec+1]
      
      // println["r is $r"]
      // println["v was    $v"]
      v = v * 2^-r  //v.adjustExponent(-r);
      // println["v is now $v"]

      [xp, p, q] = binarySplittingExp[v, 1, 0, k]
      // println["Out of binary splitting."]
      // println["xp=$xp"]
      // println["p=$p"]
      // println["q=$q"]

      // Adjust and calculate exp[x]
      pp = q
      p = p + pp
      p = p / q
      p = 1. p

      // println["Before adjust, p = $p"]

      // Reverse argument reduction
      // Brent enhancement avoids loss of significant digits when x is small.
      if r > 0
      {
         p = p - 1
         while r > 0
         {
            p =  p * (p+2)
            r = r - 1
         }
         p = p + 1
      }

      if realSignum[x] == -1
         p = 1/p

      setPrecision[digits]
      return 1. p
   }
   finally
   {
      setPrecision[origPrec]
   }
}


/** Performs a step of binarySplitting and returns [xp, p, q]
    This is based on Henrik Vestermark's algorithm,

    http://www.hvks.com/

    especially

    http://www.hvks.com/Numerical/papers.html
    http://www.hvks.com/Numerical/arbitrary_precision.html

and specifically outlined in the paper:
    http://www.hvks.com/Numerical/Downloads/HVE%20Fast%20Exp()%20calculation%20for%20arbitrary%20precision.pdf
*/
binarySplittingExp[x, xp, a, b] :=
{
//   println["In binarySplittingExp a = $a   b=$b"]
   diff = b-a
   if diff == 1
   {
      p = xp * x
      return [p, p, b]
   }

   if diff == 2
   {
      xp = xp * x
      p = (x+b) xp
      xp = xp * x
      q = b (b-1)
      return [xp, p, q]
   }

   mid = (a+b) div 2
   [xp, p, q]   = binarySplittingExp[x, xp, a, mid]  // Interval a...mid
   [xp, pp, qq] = binarySplittingExp[x, xp, mid, b]  // Interval a...mid
   return [xp, p*qq + pp, q * qq]
}

/** Stirling approximation helper function for calculating e^x. */
xstirlingApprox[digits, xexpo] :=
{
   // println["xexpo is $xexpo"]
   test = (digits + 1) * ln[10]
   // x^n/k!<10^-p, where p is the precision of the number
   // x^n~2^x’exponent
   // Stirling approximation of k!~Sqrt(2*pi*k)(k/e)^k.
   // Taking ln on both sides you get:
   //
   //   -k*log(2^xexpo) + k*(log((k)-1)+0.5*log(2*pi*m)=test=>
   //   -k*xexpo*log(2) + k*(log((k)-1)+0.5*log(2*pi*m)=test
   // Use the Newton method to find in less than 4-5 iteration
   xold = 5
   xnew = 0
   NEWTONLOOP:
   while true
   {
      f = -xold * ln[2] * xexpo + xold * (ln[xold]-1) + 0.5 ln[2 pi xold]
      // println["f is $f"]
      f1 = 0.5 / xold + ln[xold] - xexpo * ln[2]
      xnew = xold - (f - test) / f1
      if ceil[xnew] == ceil[xold]
         break NEWTONLOOP
      xold = xnew
   }

   // println["xStirlingApprox returning " + ceil[xnew]]
   return ceil[xnew]
}

/** Natural log to arbitrary precision.  This uses a cubic convergence
    algorithm (that is, the number of correct digits in the result
    approximately triple with each iteration) with adaptive precision using
    equation 3.47 in Potter.  It is significantly faster than the previous
    algorithm that did not have adaptive precision and used a quadratic
    Newton's method algorithm.
*/
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."

   if x == 1
      return 0

   extraDigits = 5 + digitLength[x]
   
   origPrec = getPrecision[]
   try
   {
      setPrecision[10]
      eps = 10.^-(digits+1)
      
      // A good initial estimate is needed.

      if (x > 0.999) and (x < 1.001)
      {
         // Use Taylor series around 1 because Math.log returns such a bad
         // value.
         y = arbitraryLnTaylor[x, digits*2, debug]
         prec = 15
      } else
      {
         if (x < 10^290) and (x > 10^-308)  // If within the range of a double
         {
            y = ln[x]
            prec = 15
         } else                         // TODO:  Store ln[2] somewhere.
         {
            y = approxLog2[x] * ln[2]
            prec = 1   // Is this reasonable?  approxLog2 has a lot of latitude.
         }
      }

      if debug
         println["Epsilon is $eps"]

      // Use Newton's method
      do
      {
         setPrecision[max[prec,15]+extraDigits]
         y2 = y
         if debug
            println["About to do arbitraryExp[" + (-y) + "]"]
         le = arbitraryExp[-y, max[prec,15]+extraDigits, debug]
         if debug
            println["Out of arbitraryExp, value is $le"]
         zn = 1 - x le
         y = y - zn(1 + zn/2)
         if debug
            println["prec is $prec, y is $y"]
         prec = prec * 3
         if (prec > digits + extraDigits)
            prec = digits + extraDigits
      } while (prec < digits) or (abs[y2-y] > eps)
      setPrecision[digits]
      retval = 1. y

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

/** This is an arbitrary-precision natural logarithm that uses a Taylor series
    for arctanh.  It is theoretically valid for any real number x>0 but should
    probably only be used closely around x=1 where the binary splitting version
    does not converge well.  This is probably because the IEEE-754 unit
    returns a poor value for ln(x) around x=1?!  For example,

    log(1.000001) gives, with IEEE_754,    9.999994999180668e-07
    whereas the actual value is            9.99999500000333333083333533333e-07

    The expression is:
    ln[x] = 2 arctanh(r) = 2 (r + 1/3 r^3 + 1/5 r^5 + 1/7 r^7 ...)
    where
    r = (x-1)/(x+1)

    See:
http://www.hvks.com/Numerical/Downloads/HVE%20Fast%20Log()%20calculation%20for%20arbitrary%20precision.pdf
*/
arbitraryLnTaylor[x, digits=getPrecision[], debug=false] :=
{
   terms = 10
   origPrec = getPrecision[]
   workingPrecision = digits + 2
   try
   {
      setPrecision[workingPrecision]
      z = (x-1)/(x+1)
      z2 = z*z
      term = z
      sum = z
      denom = 3
      for i = 1 to terms
      {
         term = term * z2
         sum = sum + term/denom
         denom = denom+2
      }

      setPrecision[digits]
      return 2. sum
   }
   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[], debug=false] :=
{
   origPrec = getPrecision[]
   try
   {
      setPrecision[digits+2]
      // TODO:  Store ln[10] somewhere.
      retval = arbitraryLn[x, digits+2,debug] / arbitraryLn[10, digits+2, debug]
      setPrecision[digits]
      return 1. retval
   }
   finally
      setPrecision[origPrec]
}

// Arbitrary log to the specified base.
arbitraryLogBase[x, base, digits=getPrecision[], debug=false] :=
{
   origPrec = getPrecision[]
   try
   {
      setPrecision[digits+2]
      retval = arbitraryLn[x, digits+2,debug] / arbitraryLn[base, digits+2, debug]
      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.
//
//  We could also try using Newton's method to invert arctan[x] which
//  has a simple series expansion,
//  arctan[x] = sum[(-1)^k x^(2k+1) / (2k + 1),  {k, 0, infinity}]
//  but this only converges for abs[x] <= 1, x != +/- i
// 
//  See:
//   Fast Algorithms for High-Precision Computation of Elementary Functions,
//   Richard P. Brent, 2006
//   https://pdfs.semanticscholar.org/bf5a/ce09214f071251bfae3a09a91100e77d7ff6.pdf
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]
}

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

/** This is a somewhat naive implementation of e^x which has a Maclaurin
    series of

    exp[x] = 1 + x + x^2/2! + x^3/3! + ...

    It is now done much faster with a binary splitting implementation.
    See BinarySplittingExp.frink for another sample.
*/
arbitraryExpOld[x, digits=getPrecision[], debug=false] :=
{
   if debug
      println["in arbitraryExpOld[$x, $digits]"]
   origPrec = getPrecision[]
   setPrecision[digits+3]

   var s
   
   try
   {
      if x < 0
      {
         s = 1.0 / arbitraryExpOld[abs[x],digits,debug]
      } 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["arbitraryExpOld[$x, $digits] returning $retval"]
      return retval
   }
   finally
      setPrecision[origPrec]
}


/**
    This calculates the gamma function.  The gamma function is closely related
    to the factorial, but generalizes to the real numbers.
    At the integers, gamma[x] = (x-1)!   or, conversely, x! = gamma[x+1]

    This algorithm is based on Henrik Vestermark's algorithm for calculating
    using integration by parts.

   See:
   http://www.hvks.com/Numerical/papers.html

   and specifically:

   http://www.hvks.com/Numerical/Downloads/HVE%20Fast%20Gamma,%20Beta,%20Error,%20Zeta%20function%20for%20arbitrary%20precision.pdf
*/
gamma[x, digits=getPrecision[], debug=false] :=
{
   precision = digits + digitLength[x] + 8

   if isInteger[x]
      if x <= 0
      {
         println["Error:  gamma function is not defined for negative integers."]
         return undef
      } else
         return (x-1)!

   // TODO:  Implement half-integer case here

   origPrec = getPrecision[]      
   try
   {
      setPrecision[precision]

      fpip = x div 1  // Called fp2 in paper, integer part
      fp   = x mod 1  // Called fp in paper, fractional part
      if debug
      {
         println["fpip is $fpip"]
         println["fp   is $fp"]
      }
      
      if (x < 0)
      {
         fpip = fpip - 1
         // if T(x<0) then use Euler reflection formula T(x)=pi/(T(1-x)*sin(pi*x))
         
         fp = gamma[1 - x]
         pi = Pi.getPi[precision]
         fp = fp * arbitrarySin[pi * x, precision]
         fp = pi / fp
         setPrecision[digits]
         return 1. * fp
      }

      // Use the Integration by parts method
      // First integral as the sum of the Taylor series exp(-u) near u==0
      // integral[0,M]u^(x-1)*exp(-u)du=M^x*sum([n=0,M](-1)^n*M^n/(n+x)*n1)
      // Step 1 choose M>(P+ln(P))*ln(10). Due to the alternating sign working
      // precision needs to be 2P
      M = ceil[(precision + ln[precision]) * ln[10]]

      if debug
         println["M is $M"]

      // Step 2
      // calculate how many shifts is needed to bring x within [1-2]. 
      shifts = 0

      if x < 1
         shifts = 1   // x<1 set shifts to 1;
      else
         if x > 2 
            shifts = -floor[fpip] + 1
         
      fp = x + shifts

      // Step 3 calculate the series sum
      nmax = ceil[precision * ln[10] / 0.2785]
      powprec = ceil[nmax + log[M]]

      fp2 = 1
      fp3 = 1
      sum = 0
      for n = 0 to nmax
      {
         fp2 = fp2 * (fp + n)
         sum = sum + (fp3 / fp2)
         fp3 = fp3 * M
      }

      // Step 4 finalize the gamma value
      sum = sum * arbitraryExp[-M, precision]
      fp = arbitraryPow[M, fp, precision] * sum

      if debug
         println["shifts is $shifts"]
      
      // Step 5 readjust for any shifted T[x]
      if shifts < 0
      {
         fp2 = 1
         while shifts < 0
         {
            fp2 = fp2 * (x + shifts)
            shifts = shifts + 1
         }

         fp = fp * fp2
      } else
      {
         if shifts > 0
         {
            // Max 1 shift
            fp2 = x
            fp = fp / fp2
         }
      }

      setPrecision[digits]
      return 1. * fp
   }
   finally
   {
      setPrecision[origPrec]
   }
}

/** The beta function is closely tied to Gamma.   It is defined as:
    beta[a, b] = gamma[a] gamma[b] / gamma[a+b]
*/
beta[a, b, digits=getPrecision[], debug=false] :=
{
   origPrecision = getPrecision[]

   try
   {
      workingPrec = digits + 4
      setPrecision[workingPrec]
      beta = gamma[a, workingPrec, debug] gamma[b, workingPrec, debug] / gamma[a+b, workingPrec, debug]
      setPrecision[digits]
      return 1. * beta
   }
   finally
      setPrecision[origPrecision]
}

/** Calculates the error function erf, which is widely used in statistics.
    For calculating the error function erf to arbitrary precision, see:

   See:
   http://www.hvks.com/Numerical/papers.html

   and specifically:

   http://www.hvks.com/Numerical/Downloads/HVE%20Fast%20Gamma,%20Beta,%20Error,%20Zeta%20function%20for%20arbitrary%20precision.pdf

   Which references:

   The functions erf and erfc computed with arbitrary precision and explicit
   error bounds, Journal: Information and Computation, : 2012, ISSN: 0890-5401

    https://core.ac.uk/display/82728431
*/
arbitraryErf[x, digits=getPrecision[], debug=false] :=
{
   if x == 0
      return 0
   
   ld = log[digits]
   workprec = digits + 20 + digitLength[x] + ceil[ld]
   bitprecision = ceil[(digits + ld) * log[10,2]]
   extra = 9

   origPrec = getPrecision[]
   try
   {
      setPrecision[workprec]
      xsq = x*x   //   x^2
      y = 2 xsq   // 2 x^2
      yexpo = floor[log[y,2]]
      eps = bitprecision + extra
      
      // Compute N using low-precision arithmetic.
      // N/(ex^2)*log2(N/(ex^2)) >= a, where
      // a=(targetprecision+3+max(0,x.exponent())-x^2*log2(e))/)ex^2)
      e1xsq = e * xsq
      xexpo = floor[log[abs[x],2]]
      a = (eps + max[0, xexpo] - xsq * log[e,2]) / e1xsq
      if a >= 2.0
         nmax = ceil[2 a e1xsq / log[a,2]]
      else
         if (a >= 0) // [0..2]
            nmax = ceil[e1xsq * 2.0^(1/4) * 2.0^(a/2)]
         else
         {
            // a<0
            nmax = ceil[e1xsq * 2.0^a]
            if nmax < 2 xsq
               nmax = ceil[2 xsq]
         }
         
      eps = negate[eps]
      eps = eps + min[0, xexpo- 1]

      // Compute L. Optimal with L~sqrt(N)
      lmax = ceil[sqrt[nmax]]

      ypl = y^lmax
      S = new array[[lmax+1], 0]
      fpacc = 2 abs[x] * arbitraryExp[-xsq] / sqrt[Pi.getPi[workprec], workprec]
      if fpacc == 0
         return x.realSignum[]  // Return either +1 or -1

      // Need to do a least kmin loops before we can consider exiting the loop
      kmin = max[y, 0.9 nmax]
      i=0
      KLOOP:
      for k=1 to nmax
      {
         S@i = S@i + fpacc
         i = i + 1
         if i == lmax
         {
            i = 0
            fpacc = fpacc * ypl
         }

         fpacc = fpacc / (2 k + 1)
         if ( k >= kmin and i==0 and (floor[log[fpacc,2]] < eps - yexpo * i))
               break KLOOP
       }

       res = S@(lmax-1)
       for i = lmax-2 to 0 step -1
          res = res y + S@i

       if x < 0
          res = negate[res]

       setPrecision[digits]
       return 1. res
   }
   finally
   {
      setPrecision[origPrec]
   }
}


/** Arbitrary-precision erfc function.  erfc[x] = 1 - erf[x] but for large x
    erf[x] may be very close to 1 and lose precision.  It should be much
    faster to call this function for large x.

    Note that this function is very limited in the number of correct digits it
    can produce for smaller x.  For example, at erfc[5], it can only produce
    9 correct digits.  It will print an error and return with an undef value if
    too many digits are requested.  The approximate number of achievable decimal
    digits is:

    exp[2 ln[x]] log[2, 10]


    THINK ABOUT:  Add a dispatcher function to only call this for larger values
    of x.

    THINK ABOUT:  If we can't produce enough digits, call 1-arbitraryErf[x]
    instead.   Look at how statistics.frink implements its dispatchers.

    This is implemented from eq. 3 and algorithm 5 of:

    The functions erf and erfc computed with arbitrary precision and explicit
    error bounds, Journal: Information and Computation, : 2012, ISSN: 0890-5401

    https://core.ac.uk/display/82728431

    In his paper
    http://www.hvks.com/Numerical/Downloads/HVE%20Fast%20Gamma,%20Beta,%20Error,%20Zeta%20function%20for%20arbitrary%20precision.pdf

    Vestermark recommends just using 1-arbitraryErf[x] instead.
*/
arbitraryErfc[x, digits=getPrecision[], debug=false] :=
{
   if x == 0
      return 0

   origPrec = getPrecision[]
   
   if x < 0
   {
      try
      {
         setPrecision[digits]
         return 2-arbitraryErfc[-x, digits, debug]
      }
      finally
         setPrecision[origPrec]
   }
   
   // Recipe 5
   tprime = ceil[digits log[10,2]]  // tprime is *bits* of precision, not digits
   if debug
      println["tprime is $tprime"]

   // We can peform the pass/(no pass) criterion at lower precision
   xsq1 = x^2
   ex2 = e xsq1
   ad = (-tprime - 3) / ex2
   if ad < -log[e,2]/e
   {
      // Frink solver to calculate number of possibly correct digits:
      // Does not always solve specified term.  Fix this.
      // https://frinklang.org/fsp/solve2.fsp?equations=xsq1+%3D+x%5E2%0D%0Aex2+%3D+e+xsq1%0D%0Aad+%3D+%28-tprime+-+3%29+%2F+ex2%0D%0Aad+%3D+-log%5Be%2C2%5D%2Fe%0D%0Adigits+%3D+tprime+%2F+log%5B10%2C2%5D%0D%0A&solveFor=&f=m&ev=on&sel_ad=S&val_ad=27%2F312500000000000+s&sel_digits=S&val_digits=&sel_e=L&val_e=2.71828182845904523536&sel_ex2=S&val_ex2=&sel_tprime=S&val_tprime=&sel_x=L&val_x=5.1&sel_xsq1=S&val_xsq1=&resultAs=&solved=on&showorig=on
      maybeDigits = floor[(xsq1 - ln[8])/ln[10]]
      println["arbitraryErfc cannot be used for x=$x and $digits digits.  It may be accurate up to $maybeDigits digits for this value of x."]
      return undef
   }

   N0 = ceil[ex2 ad / log[-ad, 2]]
   if N0 <= xsq1
      N = N0
   else
   {
      n1part = xsq1/ex2
      if n1part log[n1part, 2] <= ad
         N = ceil[xsq1]
      else
      {
         println["arbitraryErfc cannot be used for $x and $digits digits (final test)"]
         return undef
      }
   }

   // Recipe 6
   workprec = ceil[(tprime + 9 + ceil[log[N,2]]) / log[10,2]] + digitLength[x]
   if debug
      println["workprec is $workprec"]

   try
   {
      setPrecision[workprec]
      E = floor[log[x,2]]
      xsq = x*x   //   x^2
      
      G = ceil[xsq log[e,2]]
      acc = xsq
      y = 1 / (2 acc)
      acc = arbitraryExp[-acc, workprec]
      tmp = x * sqrt[Pi.getPi[workprec], workprec]
      acc = acc / tmp

      F = floor[log[abs[y],2]]

      if debug
         println["N is $N"]
      
      L = ceil[sqrt[N]] + 1  // Not sure from the paper what L should be!
//      L = floor[xsq - 1/2]  
//      L = 2 N
//      L = 2N-1
      if debug
         println["L is $L"]
      
      z = arbitraryPow[y, L, workprec]
      S = new array[[L+1], 0]

      i = 0
      k = 0
      do
      {
         S@i = S@i + (-1)^k acc
         k = k+1
         if i == L - 1
         {
            i = 0
            acc = acc * z
         } else
            i = i + 1

         acc = acc * (2 k - 1)
      } until k == N or floor[log[acc,2]] < -tprime - 3 - F i - G - E

      // Now evaluate by Horner's rule
      R = S@(L-1)
      for i = L-2 to 0 step -1
         R = R y + S@i

      setPrecision[digits]
      return 1. R
   }
   finally
   {
      setPrecision[origPrec]
   }
}

/** Arbitrary-precision phi function.  This is the cumulative distribution
    function for the normal/Gaussian distribution with a mean of 0 and sd of 1.

    Calculate the probability (from 0 to 1) of a score *less than or equal to* z
    in a normally-distributed distribution with the given mean and standard
    deviation.
*/
arbitraryPhi[z, digits=getPrecision[]] :=
{
   origPrec = getPrecision[]
   workPrec = digits + 3
   
   try
   {
      setPrecision[workPrec]
      phi = 1/2 (1 + arbitraryErf[z / sqrt[2, workPrec], workPrec])
      setPrecision[digits]
      return 1.0 phi
   }
   finally
   {
      setPrecision[origPrec]
   }
}

/** Arbitrary-precision phi function.  This is the cumulative distribution
    function for the normal/Gaussian distribution.

    Calculate the probability (from 0 to 1) of a score *less than or equal to* x
    in a normally-distributed distribution with the given mean and standard
    deviation.
    For example, to find the percentile score of an IQ of 130, given the
    assumption that IQ has mean of 100 and sd of 15, use:
        arbitraryPhi[130, 100, 15]
    which gives 0.977, or the 97.7th percentile.
*/
arbitraryPhi[x, mean, sd, digits=getPrecision[]] :=
{
   origPrec = getPrecision[]
   workPrec = digits + 3
   
   try
   {
      setPrecision[workPrec]
      z = (x - mean) / sd            // Convert to "z-score"
      phi = arbitraryPhi[z, digits]
      setPrecision[digits]
      return 1.0 phi
   }
   finally
   {
      setPrecision[origPrec]
   }
}


"ArbitraryPrecision included OK"

/*

digits = 1
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[2]
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[]
*/