View or download sqrtWayne.frink in plain text format
/** This is an attempt to convert 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.
https://introcs.cs.princeton.edu/java/92symbolic/ArbitraryPrecisionSqrt.java.html
Thanks to Jeremy Roach for the research.
*/
sqrt[n, digits] :=
{
// Handle negative numbers
if n < 0
return sqrt[-n, digits] * i
// 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
origPrecision = getPrecision[]
try
{
while (precision <= 2 * digits)
{
setPrecision[1 + 2 * precision]
t = (n/t + t) / 2.
precision = precision * 2
}
setPrecision[digits]
return 1. t
}
finally
setPrecision[origPrecision]
}
sqrtHalley[n, digits] :=
{
// Intial estimate. Find smallest power of 2 >= sqrt[n]
// TODO: Make this work with arbitrary precision? Or faster?
initial = 2.^ceil[ln[n] / ln[2] / 2]
// println["initial is $initial"]
x = initial
precision = 1
// other = sqrt[n, digits]
origPrecision = getPrecision[]
try
{
setPrecision[15]
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]
}
View or download 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 was born 18370 days, 13 hours, 58 minutes ago.