arbitrarySqrt.frink

View or download arbitrarySqrt.frink in plain text format


/** This provides arbitrary-precision square root routines that uses
    Newton's algorithm (which approximately doubles the number of correct
    digits with each iteration).

    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[14]
      // Intial estimate.  Find smallest power of 2 >= sqrt[n]

      if n > 10^308 or n < 10^-308
         initial = 2.^ceil[approxLog2[n]/2]
      else
      {
         // The built-in routines will give adequate results for 14 digits or
         // less. 
         if digits <= 14
         {
            result = sqrt[n]
            setPrecision[digits]
            return 1. * result
         } 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 is like the arbitrary-precision square root above but it gives exact
    integer (or rational number?) results when possible.
*/

sqrtExact[n, digits] :=
{
   result = sqrt[n, digits]
   if isInteger[n]
   {
      multiplier = 1
      if n < 0    // Handle negative numbers
      {
         result = -i result
         n = -n
         multiplier = i
      }
      rounded = round[result]
      if rounded^2 == n
         return rounded multiplier
   }
   return result
}


View or download arbitrarySqrt.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 18888 days, 15 hours, 32 minutes ago.