arbitrarySqrt.frink

Download or view 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, units of measure, 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] :=
{
   if (n conforms dimensionless)
      units = 1
   else
   {
      [n, units] = factorUnits[n]
      units = units^(1/2)
   }
   
   // Handle negative numbers
   if n < 0
      return sqrt[-n, digits] * i * units

   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 * units
         } 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 units
   }
   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] :=
{
   if (n conforms dimensionless)
      units = 1
   else
   {
      [n, units] = factorUnits[n]
      units = units^(1/2)
   }
   
   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 units
   }
   return result units
}


Download or view 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 20145 days, 11 hours, 50 minutes ago.