Square root of BigDecimal in Java

Question

Can we compute the square root of a BigDecimal in Java by using only the Java API and not a custom-made 100-line algorithm?

Answer 1

As it was said before: If you don't mind what precision your answer will be, but only want to generate random digits after the 15th still valid one, then why do you use BigDecimal at all?

Here is code for the method that should do the trick with floating point BigDecimals:

    import java.math.BigDecimal;
    import java.math.BigInteger;
    import java.math.MathContext;



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 square root
    // you use scale and bitlength for this, and if odd add or subtract a one
    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));
}

As a test, try to repeatedly square a number a couple of times than taking the repeated square root, and see how close you are from where you started.

Answer 2

I came up with an algorithm that doesn't just take the square root, but does every root below an integer of every BigDecimal. With the big advantage that it doesn't do a search algorithm, so it's quite fast with a 0,1ms - 1ms runtime.

But what you get in speed and versatility, it lacks in accuracy, it averages 5 correct digits with a deviancy of 3 on the fifth digit. (Tested with a million random numbers and roots), although the test ran with very high roots, so you can expect a bit more accuracy if you keep the roots below 10.

The result only holds 64 bits of precision, with the rest of the number being zeroes, so if you need very high levels of precision, don't use this function.

It's made to handle very large numbers, and very large roots, not very small numbers.

public static BigDecimal nrt(BigDecimal bd,int root) {
//if number is smaller then double_max_value it's faster to use the usual math 
//library
    if(bd.compareTo(BigDecimal.valueOf(Double.MAX_VALUE)) < 0) 
        return new BigDecimal( Math.pow(bd.doubleValue(), 1D / (double)root ));

    BigDecimal in = bd;
    int digits = bd.precision() - bd.scale() -1; //take digits to get the numbers power of ten
    in = in.scaleByPowerOfTen (- (digits - digits%root) ); //scale down to the lowest number with it's power of ten mod root is the same as initial number

    if(in.compareTo(BigDecimal.valueOf( Double.MAX_VALUE) ) > 0) { //if down scaled value is bigger then double_max_value, we find the answer by splitting the roots into factors and calculate them seperately and find the final result by multiplying the subresults
        int highestDenominator = highestDenominator(root);
        if(highestDenominator != 1) {
            return nrt( nrt(bd, root / highestDenominator),highestDenominator); // for example turns 1^(1/25) 1^(1/5)^1(1/5)
        }
        //hitting this point makes the runtime about 5-10 times higher,
        //but the alternative is crashing
        else return nrt(bd,root+1) //+1 to make the root even so it can be broken further down into factors
                    .add(nrt(bd,root-1),MathContext.DECIMAL128) //add the -1 root and take the average to deal with the inaccuracy created by this
                    .divide(BigDecimal.valueOf(2),MathContext.DECIMAL128); 
    } 
    double downScaledResult = Math.pow(in.doubleValue(), 1D /root); //do the calculation on the downscaled value
    BigDecimal BDResult =new BigDecimal(downScaledResult) // scale back up by the downscaled value divided by root
            .scaleByPowerOfTen( (digits - digits % root) / root );
    return BDResult;
}
private static int highestDenominator(int n) {
    for(int i = n-1; i>1;i--) {
        if(n % i == 0) {
            return i;
        }
    }
    return 1;
}

It works by using a mathematical property that basicly says when you are doing square roots you can change x^0.5 to (x/100)^0,5 * 10 so divide the base by 100 take the power and multiply by 10.

Generalized it becomes x^(1/n) = (x / 10^n) ^ (1/n) * 10.

So for cube roots you need to divide the base by 10^3, and for quad roots you need to divide with 10^4 and so on.

The algorithm uses that functions to scale the input down to something the math library can handle and then scale it back up again based how much the input was scaled down.

It also handles a few edge cases where the input can't be scaled down enough, and it's those edge cases that adds a lot of the accuracy problems.

Answer 3

As of Java 9 you can! See BigDecimal.sqrt().

Answer 4

By using Karp's Tricks, this can be implemented without a loop in only two lines, giving 32 digits precision:

public static BigDecimal sqrt(BigDecimal value) {
    BigDecimal x = new BigDecimal(Math.sqrt(value.doubleValue()));
    return x.add(new BigDecimal(value.subtract(x.multiply(x)).doubleValue() / (x.doubleValue() * 2.0)));
}

Answer 5

There isn't anything in the java api, so if double is not accurate enough (If not, why use BigDecimal at all?) then you need something like the below code.)

From http://www.java2s.com/Code/Java/Language-Basics/DemonstrationofhighprecisionarithmeticwiththeBigDoubleclass.htm

import java.math.BigDecimal;

public class BigDSqrt {
  public static BigDecimal sqrt(BigDecimal n, int s) {
    BigDecimal TWO = BigDecimal.valueOf(2);

    // Obtain the first approximation
    BigDecimal x = n
        .divide(BigDecimal.valueOf(3), s, BigDecimal.ROUND_DOWN);
    BigDecimal lastX = BigDecimal.valueOf(0);

    // Proceed through 50 iterations
    for (int i = 0; i < 50; i++) {
      x = n.add(x.multiply(x)).divide(x.multiply(TWO), s,
          BigDecimal.ROUND_DOWN);
      if (x.compareTo(lastX) == 0)
        break;
      lastX = x;
    }
    return x;
  }
}

Answer 6

I've used this and it works quite well. Here's an example of how the algorithm works at a high level.

Edit: I was curious to see just how accurate this was as defined below. Here is the sqrt(2) from an official source:

(first 200 digits) 1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147

and here it is using the approach I outline below with SQRT_DIG equal to 150:

(first 200 digits) 1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206086685

The first deviation occurs after 195 digits of precision. Use at your own risk if you need such a high level of precision as this.

Changing SQRT_DIG to 1000 yielded 1570 digits of precision.

private static final BigDecimal SQRT_DIG = new BigDecimal(150);
private static final BigDecimal SQRT_PRE = new BigDecimal(10).pow(SQRT_DIG.intValue());

/**
 * Private utility method used to compute the square root of a BigDecimal.
 * 
 * @author Luciano Culacciatti 
 * @url http://www.codeproject.com/Tips/257031/Implementing-SqrtRoot-in-BigDecimal
 */
private static BigDecimal sqrtNewtonRaphson  (BigDecimal c, BigDecimal xn, BigDecimal precision){
    BigDecimal fx = xn.pow(2).add(c.negate());
    BigDecimal fpx = xn.multiply(new BigDecimal(2));
    BigDecimal xn1 = fx.divide(fpx,2*SQRT_DIG.intValue(),RoundingMode.HALF_DOWN);
    xn1 = xn.add(xn1.negate());
    BigDecimal currentSquare = xn1.pow(2);
    BigDecimal currentPrecision = currentSquare.subtract(c);
    currentPrecision = currentPrecision.abs();
    if (currentPrecision.compareTo(precision) <= -1){
        return xn1;
    }
    return sqrtNewtonRaphson(c, xn1, precision);
}

/**
 * Uses Newton Raphson to compute the square root of a BigDecimal.
 * 
 * @author Luciano Culacciatti 
 * @url http://www.codeproject.com/Tips/257031/Implementing-SqrtRoot-in-BigDecimal
 */
public static BigDecimal bigSqrt(BigDecimal c){
    return sqrtNewtonRaphson(c,new BigDecimal(1),new BigDecimal(1).divide(SQRT_PRE));
}

be sure to check out barwnikk's answer. it's more concise and seemingly offers as good or better precision.