// Function for taking roots of a number using Newton's method.
//
// NOTE: sqrt can be optimized by iterating x = 1/2 (x + n/x)
root[n, p, digits] :=
{
if p==2
x = sqrt[n, digits]
else
{
err = 10.^-(ceil[digits]+1)
// Initial guess (really really dumb)
x = 2.^(approxLog2[n]/p)
do
{
//println[x]
errterm = (x^p - n)/ (p x^(p-1));
x = x - errterm;
} while (abs[errterm/n] > err)
}
r = round[x]
if (abs[r-x] < err)
if (r^p == n)
return r
return x
}
// Optimization for 2
sqrt[n, digits] :=
{
prec = getPrecision[]
err = 10.^-(ceil[digits]+1)
// Initial guess
setPrecision[5]
if (n<10.^308)
{
x = sqrt[n] // Use floating-point if we can for guess
err = err * x / 10.
scale = log[x]
} else
{
x = 2.^(approxLog2[n]/2) // Dumb guess
err = err * x
scale = approxLog2[x] / 3.219 // Approx. log 10
}
// println["First guess: $x"]
// println["Err is: $err"]
newWorkingPrecision = ceil[min[30,digits+2]]
if newWorkingPrecision < 30
newWorkingPrecision = 30
workingPrecision = 15
diff = abs[x]
while (workingPrecision < digits+2) || (diff > err*scale)
{
workingPrecision = newWorkingPrecision
// println["precision is $workingPrecision"]
setPrecision[workingPrecision]
oldx = x
x = (x + n/x) / 2
diff = abs[oldx-x]
// println["diff is $diff"]
// This slightly more than doubles working digits.
setPrecision[10]
goodDigits = -approxLog2[diff]/3.219+scale
// println["Good digits: $goodDigits"]
if (goodDigits < 30)
goodDigits = 30
newWorkingPrecision = min[digits+2, ceil[goodDigits*2.1]]
}
setPrecision[prec]
return x
}
// 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)
b = a - 1
while b<a
{
// println["$a $b"]
a = b
b = (a*a + n) div (2*a)
}
return a
}
View or download 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 14705 days, 3 hours, 9 minutes ago.