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 17th 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 🙂
View Jesper Udby's profile
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));
}