/** 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. */ /** 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 negative numbers if n < 0 return sqrt[-n, digits] * i * units 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) / 2. precision = precision * 2 } setPrecision[digits] return 1. t units } finally setPrecision[origPrecision] } /** 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 }