/** 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] }