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 ðŸ™‚

1st comment: for f(x)=sqrt(x) and an approximation x_n > sqrt(x_n), the next approximaton x_(n+1) is > sqrt(x) as well, there is no oscillation. Only when you start with a value below the corret one, you will -once only- end up with a value larger than the correct one, and it will stay larger.

2nd Comment : Newton-Ralphson is a method with quadratic convergence as long as the start value is close to the correct solution. I.e. doing n approximations for a precision of n digits is far far far too much. If you start with say 4 valid digits, your next step will bring you 8 digits, the next 16, asf.

3rd comment: Your code works for me, at least as long as I am not too far from an argument of order 1, since x/2 as a start for sqrt(x) is not exactly close…Havent tried MathContext.DECIMAL64 though, what I do is

private MathContext mc = new MathContext(50);

and

public void changeAccuracy(int accuracy) {

mc = new MathContext(accuracy);

}

on request within my calculator class.

What you should probably do is

– don’t test on equals

– scale your argument down to be sure that x=fact *y where fact is easy get the root (because it is something like 10**2n or 2**2m or both, and y is of the order of 1. Then use the double sqrt as a start value with 15 valid digits instead of x/2. Doing so the iteration, one step brings you to ~30 digits, the next to ~ 60 asf.) Of course only, if you want to have this precision!

Cheers, Henrik

Apropos, in case you wonder about the code for the sqrt as sketched before…

here it is:

public BigDecimal bigSqrt(BigDecimal d, MathContext mc) {

// 1. Make sure argument is non-negative and treat Argument 0

int sign = d.signum();

if(sign == -1)

throw new ArithmeticException(“Invalid (negative) argument of sqrt: “+d);

else if(sign == 0)

return BigDecimal.ZERO;

// 2. Scaling:

// factorize d = scaledD * scaleFactorD = scaledD * (sqrtApproxD * sqrtApproxD)

// such that scalefactorD is easy to take the sqare root ( based on the

BigInteger bigI=d.unscaledValue();

int bigS=d.scale();

int bigL = bigI.bitLength();

BigInteger scaleFactorI;

BigInteger sqrtApproxI;

if ((bigL%2==0)){

scaleFactorI=BigInteger.ONE.shiftLeft(bigL);

sqrtApproxI=BigInteger.ONE.shiftLeft(bigL/2);

}else{

scaleFactorI=BigInteger.ONE.shiftLeft(bigL-1);

sqrtApproxI=BigInteger.ONE.shiftLeft((bigL-1)/2 );

}

BigDecimal scaleFactorD;

BigDecimal sqrtApproxD;

if ((bigS%2==0)){

scaleFactorD=new BigDecimal(scaleFactorI,bigS);

sqrtApproxD=new BigDecimal(sqrtApproxI,bigS/2);

}else{

scaleFactorD=new BigDecimal(scaleFactorI,bigS+1);

sqrtApproxD=new BigDecimal(sqrtApproxI,(bigS+1)/2);

}

BigDecimal scaledD=d.divide(scaleFactorD);

// 3. This is the core algorithm:

// Newton-Ralpson for scaledD : In case of f(x)=sqrt(x),

// Heron’s Method or Babylonian Method are other names for the same thing.

// Since this is scaled we can be sure that scaledD.doubleValue() works

// for the start value of the iteration without overflow or underflow

System.out.println(“ScaledD=”+scaledD);

double dbl = scaledD.doubleValue();

double sqrtDbl = Math.sqrt(dbl);

BigDecimal a = new BigDecimal(sqrtDbl, mc);

BigDecimal HALF=BigDecimal.ONE.divide(BigDecimal.ONE.add(BigDecimal.ONE));

BigDecimal h = new BigDecimal(“0”, mc);

// when to stop iterating? You start with ~15 digits of precision, and Newton-Ralphson is quadratic

// in approximation speed, so in roundabout doubles the number of valid digits with each step.

// This fmay be safer than testing a BigDecifmal against zero.

int prec = mc.getPrecision();

int start = 15;

do {

h = scaledD.divide(a, mc);

a = a.add(h).multiply(HALF);

start *= 2;

} while (start <= prec);

// 3. Return rescaled answer. sqrt(d)= sqrt(scaledD)*sqrtApproxD :

return (a.multiply(sqrtApproxD));

}