/** This provides arbitrary-precision square root routines that use both
    Newton's algorithm (which approximately doubles the number of correct
    digits with each iteration) and Halley's algorithm (which approximately
    triples the number of correct digits with each iteration, but has more
    overhead.)

    The initial algorithm is a conversion of the algorithm of Prof. Wayne at
    Princeton for arbitrary precision sqrt (which only allows integer
    arguments) to a more generalized algorithm that allows floating-point
    arguments, negative numbers, etc.

    These algorithms work orders of magnitude faster with Java 1.8 and later
    combined with Frink: The Next Generation
    (see https://frinklang.org/experimental.html )

    https://introcs.cs.princeton.edu/java/92symbolic/ArbitraryPrecisionSqrt.java.html

    Thanks to Jeremy Roach for the research.
*/

/** This function returns the square root of n to the specified number of
    digits using Newton's method. */
sqrt[n, digits] :=
{
   // Handle negative numbers
   if n < 0
      return sqrt[-n, digits] * i

   origPrecision = getPrecision[]
   try
   {
      setPrecision[20]
      // Intial estimate.  Find smallest power of 2 >= sqrt[n]
      // TODO:  Make this return exact integers or rational numbers when possible.
      if n > 10^308 or n < 10^-308
         initial = 2.^ceil[approxLog2[n]/2]
      else
         initial = 2.^ceil[log[n, 2] / 2]

      //   println["initial is $initial"]

      t = initial
      precision = 1

      while (precision <= 2 * digits)
      {
         setPrecision[1 + 2 * precision]
         t = (n/t + t) / 2.
         precision = precision * 2
      }

      setPrecision[digits]
      return 1. t
   }
   finally
      setPrecision[origPrecision]
}


/** This uses Halley's algorithm to find the square root of a number to the
    specified number of digits.  Halley's algorithm approximately triples the
    number of correct digits with each iteration.  It switches to Newton's
    algorithm for the final iterations, as this is found to be more efficient.
*/
sqrtHalley[n, digits] :=
{
   // Handle negative numbers
   if n < 0
      return sqrtHalley[-n, digits] * i

   origPrecision = getPrecision[]
   try
   {
      setPrecision[15]
      // Intial estimate.  Find smallest power of 2 >= sqrt[n]
      // TODO:  Make this work with arbitrary precision?  Or faster?
      if n > 10^308 or n < 10^-308
         initial = 2.^ceil[approxLog2[n]/2]
      else
         initial = 2.^ceil[log[n, 2] / 2]

      //   println["initial is $initial"]

      x = initial
      precision = 1

      minError = 10.^-digits
      maxPrecision = ceil[digits * 1.05]
      do
      {
         if precision * 4 < maxPrecision 
         {
            newPrecision = 1 + 3 * precision
            if newPrecision > maxPrecision
               newPrecision = maxPrecision
            println["Using Halley for $newPrecision"]
            setPrecision[newPrecision]
            x = (x^3  + 3 n x) / (3x^2 + n)
            precision = precision * 3
         } else
         {
            newPrecision = 1 + 2 * precision
            if newPrecision > maxPrecision
               newPrecision = maxPrecision
            precision = newPrecision
            println["Using Newton for $newPrecision"]
            setPrecision[newPrecision]
            x = (n/x + x) / 2.
            precision = precision * 2
         }
         
         calcError = undef
         if precision > digits
         {
            println["Calculating error..."]
            calcError = abs[x^2 - n]
            println["calcError is " + formatSci[calcError,1,3]]
         }

//         println["Error: " + formatSci[other - x, 1, 3] + "\t" + (calcError != undef ? formatSci[calcError, 1, 3] : "") + "\tPrecision = " + newPrecision]
      } while (calcError == undef) or (calcError > minError)

//      println["Done.  calcError is " + formatSci[calcError,1,3]]
      setPrecision[digits]
      return 1. x
   }
   finally
      setPrecision[origPrecision]
}