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 20218 days, 0 hours, 8 minutes ago.