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] :=
{
   // 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"]

   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 18261 days, 16 hours, 0 minutes ago.