van NovaLoka

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))
      nPrecs.add(m);


    // 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
      {
        x = x.add(e.multiply(v, rootMC), rootMC);               // root x is ready!
        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