Download or view arbitrarySqrt.frink in plain text format or (experimental) run in browser
/** This provides arbitrary-precision square root routines that uses
Newton's algorithm (which approximately doubles the number of correct
digits with each iteration).
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, units of measure, 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.
*/
/** Function to take sqrt to current working precision. */
arbitrarySqrt[n, digits=getPrecision[]] := sqrt[n, digits]
/** This function returns the square root of n to the specified number of
digits using Newton's method. */
sqrt[n, digits] :=
{
if (n conforms dimensionless)
units = 1
else
{
[n, units] = factorUnits[n]
units = units^(1/2)
}
// Handle complex numbers
if isComplex[n]
{
r = arbitrarySqrtComplex[n, digits]
origPrec = getPrecision[]
try
{
setPrecision[digits]
return r units
}
finally
setPrecision[origPrec]
}
// Handle negative numbers
if n < 0
{
origPrec = getPrecision[]
try
{
setPrecision[digits]
return sqrt[-n, digits] * i * units
}
finally
setPrecision[origPrec]
}
origPrecision = getPrecision[]
try
{
setPrecision[14]
// Intial estimate. Find smallest power of 2 >= sqrt[n]
if n > 10^308 or n < 10^-308
initial = 2.^ceil[approxLog2[n]/2]
else
{
// The built-in routines will give adequate results for 14 digits or
// less.
if digits <= 14
{
result = sqrt[n]
setPrecision[digits]
return 1. * result * units
} 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) * 0.5
precision = precision * 2
}
setPrecision[digits]
return 1. t units
}
finally
setPrecision[origPrecision]
}
/** Handle square root of a complex number without dimensions. This can be
done in a few ways. The simplest requires an arbitrary-precision exp and
ln operator:
sqrt[z] = exp[ln[z] / 2]
which matches the Common Lisp definition for branch cuts. It may be done
this way in the future. The less simple method that requires three square
roots is:
mag = sqrt[Re[z]^2 + Im[z]^2]
sqrt[z] = sqrt[(mag + Re[z])/2] + i signum[b] sqrt[(mag-Re[z])/2]
*/
arbitrarySqrtComplex[z, digits] :=
{
origPrec = getPrecision[]
workPrec = digits + 3
try
{
setPrecision[workPrec]
a = Re[z]
b = Im[z]
mag = sqrt[a^2 + b^2, workPrec]
real = sqrt[(mag + a) * 0.5, workPrec]
imag = sqrt[(mag - a) * 0.5, workPrec] i signum[b]
setPrecision[digits]
return 1. (real + imag)
}
finally
setPrecision[origPrec]
}
/** This is like the arbitrary-precision square root above but it gives exact
integer (or rational number?) results when possible.
*/
sqrtExact[n, digits] :=
{
if (n conforms dimensionless)
units = 1
else
{
[n, units] = factorUnits[n]
units = units^(1/2)
}
result = sqrt[n, digits]
if isInteger[n]
{
multiplier = 1
if n < 0 // Handle negative numbers
{
result = -i result
n = -n
multiplier = i
}
rounded = round[result]
if rounded^2 == n
return rounded multiplier units
}
return result units
}
Download or view arbitrarySqrt.frink in plain text format or (experimental) run in browser
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