Download or view sqrtWayne.frink in plain text format
/** This provides arbitrary-precision square root routines that use both
Newton's algorithm (which approximately doubles the number of correct
digits with each iteration) and Halley's algorithm (which approximately
triples the number of correct digits with each iteration, but has more
overhead.)
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, 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] :=
{
// Handle negative numbers
if n < 0
return sqrt[-n, digits] * i
origPrecision = getPrecision[]
try
{
setPrecision[20]
// Intial estimate. Find smallest power of 2 >= sqrt[n]
// TODO: Make this return exact integers or rational numbers when possible.
if n > 10^308 or n < 10^-308
initial = 2.^ceil[approxLog2[n]/2]
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
}
finally
setPrecision[origPrecision]
}
/** This uses Halley's algorithm to find the square root of a number to the
specified number of digits. Halley's algorithm approximately triples the
number of correct digits with each iteration. It switches to Newton's
algorithm for the final iterations, as this is found to be more efficient.
*/
sqrtHalley[n, digits] :=
{
// Handle negative numbers
if n < 0
return sqrtHalley[-n, digits] * i
origPrecision = getPrecision[]
try
{
setPrecision[15]
// Intial estimate. Find smallest power of 2 >= sqrt[n]
// TODO: Make this work with arbitrary precision? Or faster?
if n > 10^308 or n < 10^-308
initial = 2.^ceil[approxLog2[n]/2]
else
initial = 2.^ceil[log[n, 2] / 2]
// println["initial is $initial"]
x = initial
precision = 1
minError = 10.^-digits
maxPrecision = ceil[digits * 1.05]
do
{
if precision * 4 < maxPrecision
{
newPrecision = 1 + 3 * precision
if newPrecision > maxPrecision
newPrecision = maxPrecision
println["Using Halley for $newPrecision"]
setPrecision[newPrecision]
x = (x^3 + 3 n x) / (3x^2 + n)
precision = precision * 3
} else
{
newPrecision = 1 + 2 * precision
if newPrecision > maxPrecision
newPrecision = maxPrecision
precision = newPrecision
println["Using Newton for $newPrecision"]
setPrecision[newPrecision]
x = (n/x + x) / 2.
precision = precision * 2
}
calcError = undef
if precision > digits
{
println["Calculating error..."]
calcError = abs[x^2 - n]
println["calcError is " + formatSci[calcError,1,3]]
}
// println["Error: " + formatSci[other - x, 1, 3] + "\t" + (calcError != undef ? formatSci[calcError, 1, 3] : "") + "\tPrecision = " + newPrecision]
} while (calcError == undef) or (calcError > minError)
// println["Done. calcError is " + formatSci[calcError,1,3]]
setPrecision[digits]
return 1. x
}
finally
setPrecision[origPrecision]
}
Download or view sqrtWayne.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, eliasen@mindspring.com