}
// Initial precision is that of double numbers 2^63/2 ~ 4E18
int BITS = 62; // 63-1 an even number of number bits
int nInit = 16; // precision seems 16 to 18 digits
MathContext nMC = new MathContext(18, RoundingMode.HALF_DOWN);
// Iteration variables, for the square root x and the reciprocal v
BigDecimal x = null, e = null; // initial x: x0 ~ sqrt()
BigDecimal v = null, g = null; // initial v: v0 = 1/(2*x)
// Estimate the square root with the foremost 62 bits of squarD
BigInteger bi = squarD.unscaledValue(); // bi and scale are a tandem
int biLen = bi.bitLength();
int shift = Math.max(0, biLen - BITS + (biLen%2 == 0 ? 0 : 1)); // even shift..
bi = bi.shiftRight(shift); // ..floors to 62 or 63 bit BigInteger
double root = Math.sqrt(bi.doubleValue());
BigDecimal halfBack = new BigDecimal(BigInteger.ONE.shiftLeft(shift/2));
int scale = squarD.scale();
if (scale % 2 == 1) {
root *= SQRT_10; // 5 -> 2, -5 -> -3 need half a scale more..
}
scale = (int) Math.floor(scale/2.); // ..where 100 -> 10 shifts the scale
// Initial x - use double root - multiply by halfBack to unshift - set new scale
x = new BigDecimal(root, nMC);
x = x.multiply(halfBack, nMC); // x0 ~ sqrt()
if (scale != 0) {
x = x.movePointLeft(scale);
}
if (prec < nInit) { // for prec 15 root x0 must surely be OK
return x.round(rootMC); // return small prec roots without iterations
}
// Initial v - the reciprocal
v = BigDecimal.ONE.divide(TWO.multiply(x), nMC); // v0 = 1/(2*x)
// Collect iteration precisions beforehand
List<Integer> nPrecs = new ArrayList<Integer>();
assert nInit > 3 : "Never ending loop!"; // assume nInit = 16 <= prec
// Let m be the exact digits precision in an earlier! loop
for (int m = prec + 1; m > nInit; m = m/2 + (m > 100 ? 1 : 2)) {
nPrecs.add(m);
}
// The loop of "Square Root by Coupled Newton Iteration"
for (int i = nPrecs.size() - 1; i > -1; i--) {
// Increase precision - next iteration supplies n exact digits
nMC = new MathContext(nPrecs.get(i), (i%2 == 1) ? RoundingMode.HALF_UP :
RoundingMode.HALF_DOWN);
// Next x // e = d - x^2
e = squarD.subtract(x.multiply(x, nMC), nMC);
if (i != 0) {