Squareroot (sqrt) with BigDecimal

Recently I was refactoring some code from using double to using BigDecimal and suddenly needed a square root method.

I remembered years ago I was taught a simple “successive approximation” method; believe it was grammar school.

I searched the net and found the “Babylonian method” (which is actually a “Newton-Raphson method” implementation). This is the method I was taught and I present a BigDecimal-implementation here.

First a little theory
We start out with a guess g, that in this case could be just about any number (this is not always the truth when using the “Newton-Raphson method”, see the references).
From that we calculate a result r (where x is the number we want to calculate the square root of):

r = (x/g + g)/2

Assume g is a good approximation, we can say:

g ~ sqrt(x) =>
 
r ~ (x/sqrt(x) + sqrt(x))/2 <=>
 
r ~ (sqrt(x) + sqrt(x))/2 <=>
 
r ~ sqrt(x)

In other words, if g is a good approximation of sqrt(x) then r will be too. If g is higher or lower than sqrt(x), the approximation gets better in the “opposite” direction; below is shown if g < sqrt(x)

g < sqrt(x) =>
 
x/g > x/sqrt(x) =>
 
r = (x/g + g)/2 > sqrt(x)

In other words, if the guess is above the right answer, the next approximation is lower and vice-versa. Whatever the guess is, the next guess will be better 🙂

This is the basic idea behind Newton-Raphson approximations.

So we implement the above approximation, assuming we get closer to the right result for each iteration, and stop when the current approximation equals the previous.

One problem here is oscillations: situations where the result “oscillates” between 2 different values. In those cases the approximation never gets a definite answer. In the solution given, we assume each iteration calculates a new digit; i.e. if the MathContext dictates 16 decimals precision, we break out of the loop after the 17^th iteration.

The [Java] code

private static final BigDecimal TWO = BigDecimal.valueOf(2L);
 
public static BigDecimal sqrt(BigDecimal x, MathContext mc) {
	BigDecimal g = x.divide(TWO, mc);
	boolean done = false;
	final int maxIterations = mc.getPrecision() + 1;		
	for (int i = 0; !done && i < maxIterations; i++) {
		// r = (x/g + g) / 2
		BigDecimal r = x.divide(g, mc);
		r = r.add(g);
		r = r.divide(TWO, mc);
		done = r.equals(g);
		g = r;
	}
	return g;
}

This works out nicely. If MathContext.DECIMAL64 is used for math-context, we would expect to get same behavior as using a double. This is unfortunately not the case. If the above method is used to calculate sqrt(x) for x = 0.01..10000.00 there are 352261 situations where the BigDecimal method is wrong; or just about 1 out of 3.

Using the above method with MathContext.DECIMAL64 is not better than using e.g. StrictMath.sqrt() for double.

In other words, if calculating the square root of a BigDecimal with precision comparable with double, one should use the method below:

public static BigDecimal sqrt(BigDecimal x) {
	return BigDecimal.valueOf(StrictMath.sqrt(x.doubleValue()));
}

Doing a little performance testing, we see that calculating sqrt using double as intermediate type is about 7 times faster and is more precise than the above mentioned sqrt method; see output of my test program below:

started
timing (1):75.753043015
6.666661664588234E7
timing (2):11.431122455
6.666661664588234E7
timing (1):74.843013516
6.666661664588234E7
timing (2):11.661304554
6.666661664588234E7
For REAL
timing (1):73.35738404700001
6.666661664588234E7
timing (2):11.169711972
6.666661664588234E7

In other words – the BigDecimal.sqrt() implementation described above has academic interest only, when doing calculations with double precision 🙂

References

• BigDecimal sqrt computations source
• wiki: Newtons method
• wiki: Babylonian method