root.frink

// 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
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
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
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
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
}

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 18321 days, 9 hours, 11 minutes ago.