sitemap →

van NovaLoka

mostly maths

Java BigDecimal sqrt() method - Square Root by Coupled Newton Iteration

```import java.math.*;
import java.util.ArrayList;

/**
* @author Fava
* @version Java 5, 28 September 2007
*/
```
```  /**
* Returns the correctly rounded square root of a positive `BigDecimal`.

* The algorithm for taking the square root of a BigDecimal is most
* critical for the speed of your application. This method performs the fast
* Square Root by Coupled Newton Iteration algorithm by Timm Ahrendt,
* from the book "Pi, unleashed" by Jörg Arndt in a neat loop.
*
* @param squarD   number to get the root from (called "d" in the book)
* @param rootMC   precision and rounding mode (for the last root "x")
*
* @return the root of the argument number
*
* @throws ArithmeticException       if the argument number is negative
* @throws IllegalArgumentException  if rootMC has precision 0
*/

public static BigDecimal bigSqrt(BigDecimal squarD, MathContext rootMC)
{
// Static constants - perhaps initialize in class Vladimir!
BigDecimal TWO = new BigDecimal(2);
double SQRT_10 = 3.162277660168379332;

// General number and precision checking
int sign = squarD.signum();
if(sign == -1)
throw new ArithmeticException("\nSquare root of a negative number: " + squarD);
else if(sign == 0)
return squarD.round(rootMC);

int prec = rootMC.getPrecision();           // the requested precision
if(prec == 0)
throw new IllegalArgumentException("\nMost roots won't have infinite precision = 0");

// 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)                          // add half scales of the root to odds..
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
ArrayList<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))

// The loop of "Square Root by Coupled Newton Iteration" for simpletons
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)
x = x.add(e.multiply(v, nMC));                          // x += e*v     ~ sqrt()
else
{
break;
}

// Next v                                                 // g = 1 - 2*x*v
g = BigDecimal.ONE.subtract(TWO.multiply(x).multiply(v, nMC));

v = v.add(g.multiply(v, nMC));                            // v += g*v     ~ 1/2/sqrt()
}

return x;                        // return sqrt(squarD) with precision of rootMC
}

```

For example:

```----------------------
squarD = 45925227810970180969313328047107279360892708083547808019105709473030650
45569481255426637945266025079673783858555949851093788475843452221325703725424666
61274156519956581249120002755116592960179485715707960253944979126461854585553455
78926810416186541588361998623264959523577851842249324659814465688788133072537667
56359378340224095179392316562356037391060017160108668265033220880325289451251670
45071533268826395981585672508087119433616168639412733006163557735608367682219831
13435868815664032782886652233102152170233490406211339514567772821803566110313553
0500701742115203291371083271861815048739404596117537

n=18     x0=2.14301721437253464E+301
n=22     x=2.143017214372534641897E+301
n=40     x=2.143017214372534641896850098120003621123E+301
n=77     x=2.1430172143725346418968500981200036211228096234110672148875007767407
021022499E+301
n=152    x=2.1430172143725346418968500981200036211228096234110672148875007767407
02102249872244986396757631391716255189345835106293650374290571384628087196915514
9397E+301
n=302     x=21430172143725346418968500981200036211228096234110672148875007767407
02102249872244986396757631391716255189345835106293650374290571384628087196915514
93971496078691355496484619708421492101247422837559083643060929499671638825347975
35118331087892154125829142392955373084335320859663305248773674411336138752
n=602     x=21430172143725346418968500981200036211228096234110672148875007767407
02102249872244986396757631391716255189345835106293650374290571384628087196915514
93971496078691355496484619708421492101247422837559083643060929499671638825347975
35118331087892154125829142392955373084335320859663305248773674411336138752.00000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000001

rootFrac = 1E-300
```

While a more precise calculation affirms this result:

```  rootFrac = 7.69942485265155575166823874397149164934100443256659120312914330331
67270155501489842491326985575190944666586669251750016497777720511707672574425467
4193242873893388799564118248976045937688369E-301
```

posted on Saturday, September 15, 2007 5:20 PM

Feedback

#re: Java BigDecimal sqrt() method - Square Root by Coupled Newton Iteration 9/22/2007 1:22 PM tess

root method heb je ook in de endodontologie ;-)

#re: Java BigDecimal sqrt() method - Square Root by Coupled Newton Iteration 9/23/2007 10:23 AM Frans L.

Worteltrekken bedoel je?

#re: Java BigDecimal sqrt() method - Square Root by Coupled Newton Iteration 9/29/2007 1:10 PM tess

Ha Frans, ja worteltrekken. Hahaha ... Hoe gaat het ermee? Ben je ondergedoken op Java of beter gezegd in Java? Of tussen Java en het schip in? ;-)