karatsuba.frink

View or download karatsuba.frink in plain text format


/** This library calculates nearest square roots of integers using the
    Karatsuba square root recursive algorithm.  It also contains
    functions that build upon the integer routines to calculate floating-point
    (and rational) square roots of integers and floating-point numbers to an
    arbitrary number of digits. 

    For reference to the Karatsuba square root algorithm see:
       https://hal.inria.fr/inria-00072854/PDF/RR-3805.pdf
*/



/** Perform Karatsuba square root, returning a floating-point square root of an
    integer or calculating the floating-point square root of a floating-point
    number (by first converting it to a giant integer.)

    This works using a rather simple trick.  If we need a lot of
    digits in the result, we "shift left" the number by the number of bits
    we need (by multiplying by 2^m bits, (see how m is calculated below))
    and round it to an integer.  We then do an integer square root of this big
    integer, and divide the result by 2^(m/2) bits which turns the result into
    the correct floating-point (or actually rational) square root of the
    number.

    (Note that this trick also elegantly works to *reduce* the number of digits
     we have to calculate by "shifting right" when the requested number of
     digits is small.  For example, karatsubaFloat[10^1000, 3] would
     normally take a 1000-digit number and use the full force of Frink to
     calculate 3 digits of its square root.)
   
    This algorithm is able to calculate a limited number of digits, which
    makes it suitable for providing a good initial guess (with a large number
    of correct digits) for an operation like Newton iteration for larger
    numbers.  Newton iteration approximately doubles the number of correct
    digits per iteration if given a good guess, so an optimal algorithm *may*
    use this algorithm and Newton iteration.  Note that as of Java 1.8,
    the integer operations from Frink were incorporated into Java, making
    integer operations much faster than floating-point for large numbers.
   
    If rational==false, (the default,) this returns a floating-point number.
    If rational==true,  this returns a rational number (which may be faster)
*/

karatsubaSqrt[n, digits=getPrecision[], rational=false] :=
{
   oldPrec = getPrecision[]

   try
   {
      m = ceil[digits log[10]/log[2] - approxLog2[n]]
      setPrecision[digits+3]

      // Make the number of bits even
      if m mod 2 == 1
         m = m + 1

      shift = 2^m
      
      [num, remainder] = sqrtRem[round[n * shift]]
      denom = 2^(m/2)

      // println["$num / $denom, remainder = $remainder"]
      
      result = num/denom
      
      /* This adjusts the result using the remainder and a Taylor series
         approximation:

         The Karatsuba square root algorithm only returns integers, so it
         estimates the nearest integer to the result (called s in the paper)
         as:
         result = floor[sqrt[n]]

         And the remainder (called r in the paper) is:
         remainder = n - result^2
      
         thus, the correct value of sqrt[n] is:
         sqrt[n] = sqrt[result^2 + remainder]

         However, this is not easily turned into a number without doing
         another whole square root algorithm, so we correct it by calculating
         a Taylor series for the corrections in the remainder.  That is, the
         Taylor series of the above:  sqrt[result^2 + remainder]

         Due to the choice of bits that we chose in the initial calculation,
         1 term in the Taylor series seems to be just right at enough.
         Without using the remainder, we have to calculate *twice* as many
         digits in the recursive algorithm above to get correct answers.  By
         using the remainder, we can get all digits correct by calculating a
         couple more digits than requested, rather than doubling the number of
         digits.

         The following terms are the first terms of the Taylor series.  

         First term:  r / (2 s)
      */

      radj = remainder / (2 result shift)
      result = result + radj
      //println["Remainder adjustment is $radj"]

      /* Second-order Taylor series remainder adjustment.  This is generally
         not necessary and may only affect the last digit or so.
      
         Second term:  -r^2 (8 s^3)
      */

      
      /* radj2 = - remainder^2 / (8 (result shift)^3)
      //  println["Remainder adjustment (second order) is $radj2"]
      result = result + radj2 */

      
      setPrecision[digits]
      
      if rational == false
         return 1. * result   // Return floating-point number
      else
         return result        // Return rational number
   }
   finally
   {
      setPrecision[oldPrec]
   }
}

/* Karatsuba recursive square root.

   For explication of algorithm, see:
     https://hal.inria.fr/inria-00072854/PDF/RR-3805.pdf

   Thanks to Jeremy Roach for much hard work at turning my broken
   implementation into a working implementation.
*/

sqrtRem[n] :=
{
   // Fallback case when numbers are smallish (TODO: tune this threshold)
   if n < 2^256
   {
      s = intSqrt[n]
      return [s, n - s^2]
   }

   [w, z] = divideAndRemainder[bitLength[n], 4]
   w = ((4 w + z) + ((4-z) mod 4)) / 4
   b = 2^w

   if z == 1 || z == 2
   {
      fixup = true
      n = 4 n
   } else
   {
      fixup = false
   }

   a0 = bitAnd[n, b-1]
   n = shiftRight[n, w]
   a1 = bitAnd[n, b-1]
   n = shiftRight[n, w]
   a2 = bitAnd[n, b-1]
   n = shiftRight[n, w]
   a3 = bitAnd[n, b-1]

   [s1, r1] = sqrtRem[a3*b + a2]
   [q, u] = divideAndRemainder[r1*b + a1, 2*s1]
   s = s1*b + q
   r = u*b + a0 - q^2
   if r < 0
   {
      r = r + 2*s - 1
      s = s - 1
   }
   if fixup
      return [s/2, r/4]
   
   return [s, r]
}

/* Integer square root: returns the greatest integer less than or equal to the 
   square root of n.  Only call this with integers.

   This is Exercise 5.7 in Bressoud with my own modifications for better
   initial guess.
*/

intSqrt[n] :=
{
   a = 2^((bitLength[n]+1) div 2)
   b = a - 1

   while b<a
   {
      a = b
      b = (a^2 + n) div (2*a)
   }

   return a
}


View or download karatsuba.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 17591 days, 17 hours, 47 minutes ago.