sqrtWayne.frink

View or download sqrtWayne.frink in plain text format


/** This is an attempt to convert 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.

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

    Thanks to Jeremy Roach for the research.
*/

sqrt[n, digits] :=
{
   // Handle negative numbers
   if n < 0
      return sqrt[-n, digits] * i
   
   // 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

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

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


sqrtHalley[n, digits] :=
{
   // Intial estimate.  Find smallest power of 2 >= sqrt[n]
   // TODO:  Make this work with arbitrary precision?  Or faster?
   initial = 2.^ceil[ln[n] / ln[2] / 2]
//   println["initial is $initial"]

   x = initial
   precision = 1

//   other = sqrt[n, digits]

   origPrecision = getPrecision[]
   try
   {
      setPrecision[15]
      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]
}


View or download sqrtWayne.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 18370 days, 13 hours, 58 minutes ago.