Download or view 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
}
Download or view 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 19966 days, 18 hours, 36 minutes ago.