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] :=
{
// 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"]
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 18261 days, 16 hours, 0 minutes ago.