root.frink

View or download root.frink in plain text format



// Functions for taking roots of a number using Newton's method to arbitrary
// precision using the root[] function.
//
// This function also performs integer and rational exponents to arbitrary
// precision using its pow[] function.
//
// These methods shouldn't be used if you're already working with imprecise
// floating-point numbers.
//
// Newton's method converges quadratically if you have a good guess.
// It can be horribly slow if you don't start with a good guess.
//
// This uses a few tricks:
//  * Uses floating-point math (when possible) for a good initial guess.
//  * Dynamically increases the working precision as estimates get better.


// Convenience function to do square root to current working precision.
sqrt[n] := sqrt[n, getPrecision[]]

// Convenience method for sqrt.
sqrt[n, digits, debug=false] := root[n, 2, digits, debug]

// Convenience function to do power to current working precision.
pow[n, p] := pow[n, p, getPrecision[]]

// Arbitrary-precision power
pow[n, p, digits, debug=false] :=
{
   if debug
      println["Entering pow $n $p $digits"]
   
   if p == 0
      return 1
   else
   {
      origPrec = getPrecision[]
      try
      {
         setPrecision[digits+4]
         return root[n, 1/p, digits, debug]
      }
      finally
         setPrecision[origPrec]
   }
}

// Convenience function to do root to current working precision.
root[n, p] := root[n, p, getPrecision[]]

// General arbitrary-precision root finder
// n is the number, p is the root (e.g. 2 for square root, 3 for cube root)
root[n, p, digits, debug=false] :=
{
   if debug
      println["in root[$n, $p]"]
   if p == 1
      return n
   if p == 0
      return 1
   if p < 0
   {
      origPrec = getPrecision[]
      try
      {
         setPrecision[digits+4]
         return 1/root[n, -p]
      }
      finally
         setPrecision[origPrec]
   }
   
   alter = false
   if n<0
   {
      if p == 2
      {
         alter=true
         factor = i
         n = -n
      } else
        if p mod 2 == 1         // Negative base, odd power
        {
           alter = true
           factor = -1
           n = -n               // Negative base, even power
        } else
        {
           println["Arbitrary-precision root cannot produce complex numbers.  Arguments were root[$n, $p, $digits]"]
           return undef
        }
   }

   origPrec = getPrecision[]
   try
   {
      // Handle exact rational numbers
      if isRational[p]
      {
         prec = getPrecision[]
         setPrecision[digits+3]

         // TODO:  This needs to use arbitary-precision powers!
         // We can't apparently replace it with call to pow[] because
         // the stack never terminates.
         retval = root[n, numerator[p], digits, debug]^denominator[p]
         if alter
            retval = retval * factor
         setPrecision[prec]
         return retval
      }
      
      prec = getPrecision[]
      err = 10.^-(ceil[digits]+1)

      // Initial guess
      setPrecision[5]
      if (n<1e+308 and n>1e-323)
      {
         x = n^(1/p)               // Use floating-point if we can for guess
         err = err * x / 10.
         scale = approxLog2[x] / 3.219   // Approx. log 10
      } else
      {
         x = 2.^(approxLog2[n]/p)  // Dumb guess; could use introot function below
         err = err * x
         scale = approxLog2[x] / 3.219  // Approx. log 10
      }

      if scale == 0
         scale = 1
      
      if (debug)
      {
         println["First guess: $x"]
         println["Err is: $err"]
         println["Scale is: $scale"]
      }

      newWorkingPrecision = ceil[min[30,digits+3]]
      if newWorkingPrecision < 30
         newWorkingPrecision = 30
      workingPrecision = 15
      diff = abs[x]

      scaledErr = abs[err*scale]
      
      while (workingPrecision < digits+p) || (abs[diff] > scaledErr)
      {
         workingPrecision = newWorkingPrecision
         if debug
            println["precision is $workingPrecision"]
         setPrecision[workingPrecision]
         oldx = x
         if p==2
         {
            x = (x + n / x) / 2
            diff = oldx - x
         } else
         {
            // TODO:  This needs to use arbitary-precision powers!
            // We can't apparently replace it with call to pow[] because
            // the stack never terminates.
            errterm = (x^p - n)/ (p x^(p-1))
            x = x - errterm
            diff = errterm
         }

         if debug
         {
            println["x is $x"]
            println["diff is $diff"]
         }
         // This slightly more than doubles  working digits.
         setPrecision[10]
         if diff == 0
            goodDigits = workingPrecision * 2
         else
         {
            if debug
               println["approxLog2 is " + approxLog2[abs[diff]]]
            
            goodDigits = -approxLog2[abs[diff]]/3.219+scale
         }

         if debug
            println["Good digits: $goodDigits"]
         if (goodDigits < 30)
            goodDigits = 30
         newWorkingPrecision = min[ceil[digits+p+1], ceil[goodDigits*2.1]] 
      }

      if (debug)
         println["Final diff is $diff"]

      if alter
      {
         setPrecision[digits+p+1]
         retval = factor * x
      } else
      retval = x
      
      return retval
   }
   finally
      setPrecision[origPrec]
}

// Inverse square root.  1/sqrt[n].
inverseSquareRoot[n, digits, debug=false] :=
{
   alter = false
   if n<0
   {
      alter=true
      factor = i
      n = -n
   }

   origPrec = getPrecision[]
   try
   {
      prec = getPrecision[]
      err = 10.^-(ceil[digits]+1)

      // Initial guess
      setPrecision[5]
      if (n<1e+308 and n>1e-323)
      {
         x = 1/(n^(1/2))           // Use floating-point if we can for guess
         err = err * x / 10.
         scale = approxLog2[x] / 3.219   // Approx. log 10
      } else
      {
         x = 1/(2.^(approxLog2[n]/2))// Dumb guess; could use introot function below
         err = err * x
         scale = approxLog2[x] / 3.219  // Approx. log 10
      }

      if scale == 0
         scale = 1
      
      if (debug)
      {
         println["First guess: $x"]
         println["Err is: $err"]
         println["Scale is: $scale"]
      }

      newWorkingPrecision = ceil[min[30,digits+3]]
      if newWorkingPrecision < 30
         newWorkingPrecision = 30
      workingPrecision = 15
      diff = abs[x]

      scaledErr = abs[err*scale]
      
      while (workingPrecision < digits+2) || (diff > scaledErr)
      {
         workingPrecision = newWorkingPrecision
         if debug
            println["precision is $workingPrecision"]
         setPrecision[workingPrecision]
         oldx = x

         diff = 1 - n * x^2
         x = x + diff * x/2

         if debug
         {
            println["x is $x"]
            println["diff is $diff"]
         }
         // This slightly more than doubles  working digits.
         setPrecision[10]
         if diff == 0
            goodDigits = workingPrecision * 2
         else
            goodDigits = -approxLog2[abs[diff]]/3.219+scale
         if debug
            println["Good digits: $goodDigits"]
         if (goodDigits < 30)
            goodDigits = 30
         newWorkingPrecision = min[ceil[digits+3], ceil[goodDigits*1.8]] 
      }

      if (debug)
         println["Final diff is $diff"]

      if alter
      {
         setPrecision[digits+2+1]
         retval = factor * x
      } else
      retval = x
      
      return retval
   }
   finally
      setPrecision[origPrec]
}


// Integer square root--returns the greatest integer less than or equal to the 
// to the square root of n.
// This is Exercise 5.7 in Bressoud with my own modifications for better
// initial guess.
introot[n] :=
{
   a = 2^((bitLength[n]+1) div 2)
   //a = 2^((ceil[approxLog2[n]+1]) div 2)
   b = a - 1

   while b<a
   {
//      println["$a $b"]
      a = b
      b = (a*a + n) div (2*a)
   }

   return a
}


View or download root.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 17957 days, 10 hours, 39 minutes ago.