Square root of BigDecimal in Java

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?

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.

public static BigDecimal sqrt(BigDecimal A, final int SCALE) {
	BigDecimal x0 = new BigDecimal("0");
    BigDecimal x1 = new BigDecimal(Math.sqrt(A.doubleValue()));
    while (!x0.equals(x1)) {
        x0 = x1;
        x1 = A.divide(x0, SCALE, ROUND_HALF_UP);
        x1 = x1.add(x0);
        x1 = x1.divide(TWO, SCALE, ROUND_HALF_UP);
       
    }
    return x1;
}

This work perfect! Very fast for more than 65536 digits!

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

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)));
}

If you need to find only integer square roots - these are two methods that can be used.

Newton's method - very fast even for 1000 digits BigInteger:

public static BigInteger sqrtN(BigInteger in) {
    final BigInteger TWO = BigInteger.valueOf(2);
    int c;

    // Significantly speed-up algorithm by proper select of initial approximation
    // As square root has 2 times less digits as original value
    // we can start with 2^(length of N1 / 2)
    BigInteger n0 = TWO.pow(in.bitLength() / 2);
    // Value of approximate value on previous step
    BigInteger np = in;

    do {
        // next approximation step: n0 = (n0 + in/n0) / 2
        n0 = n0.add(in.divide(n0)).divide(TWO);

        // compare current approximation with previous step
        c = np.compareTo(n0);

        // save value as previous approximation
        np = n0;

        // finish when previous step is equal to current
    }  while (c != 0);

    return n0;
}

Bisection method - is up to 50x times slower than Newton's - use only in educational purpose:

 public static BigInteger sqrtD(final BigInteger in) {
    final BigInteger TWO = BigInteger.valueOf(2);
    BigInteger n0, n1, m, m2, l;
    int c;

    // Init segment
    n0 = BigInteger.ZERO;
    n1 = in;

    do {
        // length of segment
        l = n1.subtract(n0);

        // middle of segment
        m = l.divide(TWO).add(n0);

        // compare m^2 with in
        c = m.pow(2).compareTo(in);

        if (c == 0) {
            // exact value is found
            break;
        }  else if (c > 0) {
            // m^2 is bigger than in - choose left half of segment
            n1 = m;
        } else {
            // m^2 is smaller than in - choose right half of segment
            n0 = m;
        }

        // finish if length of segment is 1, i.e. approximate value is found
    }  while (l.compareTo(BigInteger.ONE) > 0);

    return m;
}

If you want to calculate square roots for numbers with more digits than fit in a double (a BigDecimal with a large scale) :

Wikipedia has an article for computing square roots: <http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method>

This is my implementation of it:

public static BigDecimal sqrt(BigDecimal in, int scale){
	BigDecimal sqrt = new BigDecimal(1);
	sqrt.setScale(scale + 3, RoundingMode.FLOOR);
	BigDecimal store = new BigDecimal(in.toString());
	boolean first = true;
	do{
		if (!first){
			store = new BigDecimal(sqrt.toString());
		}
		else first = false;
		store.setScale(scale + 3, RoundingMode.FLOOR);
		sqrt = in.divide(store, scale + 3, RoundingMode.FLOOR).add(store).divide(
				BigDecimal.valueOf(2), scale + 3, RoundingMode.FLOOR);
	}while (!store.equals(sqrt));
	return sqrt.setScale(scale, RoundingMode.FLOOR);
}

setScale(scale + 3, RoundingMode.Floor) because over calculating gives more accuracy. RoundingMode.Floor truncates the number, RoundingMode.HALF_UP does normal rounding.

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.)

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;
  }
}

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

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.

Here is very accurate and fast solution, it's based on my BigIntSqRoot solution and the next observation: The square root of A^2B - Is A multiply the root of B. Using this method I can easily calculate the first 1000 digits of square root of 2.

1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079896872533965463318088296406206152583523950547457502877599617298355752203375318570113543746034084988471603868999706990048150305440277903164542478230684929369186215805784631115966687130130156185689872372352885092648612494977154218334204285686060146824720771435854874155657069677653720226485447015858801620758474922657226002085584466521458398893944370926591800311388246468157082630100594858704003186480342194897278290641045072636881313739855256117322040245091227700226941127573627280495738108967504018369868368450725799364729060762996941380475654823728997180326802474420629269124859052181004459842150591120249441341728531478105803603371077309182869314710171111683916581726889419758716582152128229518488472

So here is the source code

public class BigIntSqRoot {
    private static final int PRECISION = 10000;
    private static BigInteger multiplier = BigInteger.valueOf(10).pow(PRECISION * 2);
    private static BigDecimal root = BigDecimal.valueOf(10).pow(PRECISION);
    private static BigInteger two = BigInteger.valueOf(2L);

    public static BigDecimal bigDecimalSqRootFloor(BigInteger x)
            throws IllegalArgumentException {
        BigInteger result = bigIntSqRootFloor(x.multiply(multiplier));
        //noinspection BigDecimalMethodWithoutRoundingCalled
        return new BigDecimal(result).divide(root);
    }

    public static BigInteger bigIntSqRootFloor(BigInteger x)
            throws IllegalArgumentException {
        if (checkTrivial(x)) {
            return x;
        }
        if (x.bitLength() < 64) { // Can be cast to long
            double sqrt = Math.sqrt(x.longValue());
            return BigInteger.valueOf(Math.round(sqrt));
        }
        // starting with y = x / 2 avoids magnitude issues with x squared
        BigInteger y = x.divide(two);
        BigInteger value = x.divide(y);
        while (y.compareTo(value) > 0) {
            y = value.add(y).divide(two);
            value = x.divide(y);
        }
        return y;
    }

    public static BigInteger bigIntSqRootCeil(BigInteger x)
            throws IllegalArgumentException {
        BigInteger y = bigIntSqRootFloor(x);
        if (x.compareTo(y.multiply(y)) == 0) {
            return y;
        }
        return y.add(BigInteger.ONE);
    }

    private static boolean checkTrivial(BigInteger x) {
        if (x == null) {
            throw new NullPointerException("x can't be null");
        }
        if (x.compareTo(BigInteger.ZERO) < 0) {
            throw new IllegalArgumentException("Negative argument.");
        }
        
        return x.equals(BigInteger.ZERO) || x.equals(BigInteger.ONE);
    }
}

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.

public static BigDecimal sqrt( final BigDecimal value )
{
    BigDecimal guess = value.multiply( DECIMAL_HALF ); 
    BigDecimal previousGuess;
    
    do
    {
        previousGuess = guess;
        guess = sqrtGuess( guess, value );
   } while ( guess.subtract( previousGuess ).abs().compareTo( EPSILON ) == 1 );
    
    return guess;
}

private static BigDecimal sqrtGuess( final BigDecimal guess,
                                     final BigDecimal value )
{
    return guess.subtract( guess.multiply( guess ).subtract( value ).divide( DECIMAL_TWO.multiply( guess ), SCALE, RoundingMode.HALF_UP ) );
}

private static BigDecimal epsilon()
{
    final StringBuilder builder = new StringBuilder( "0." );
    
    for ( int i = 0; i < SCALE - 1; ++i )
    {
        builder.append( "0" );
    }
    
    builder.append( "1" );
    
    return new BigDecimal( builder.toString() );
}

private static final int SCALE = 1024;
private static final BigDecimal EPSILON = epsilon();
public static final BigDecimal DECIMAL_HALF = new BigDecimal( "0.5" );
public static final BigDecimal DECIMAL_TWO = new BigDecimal( "2" );

BigDecimal.valueOf(Math.sqrt(myBigDecimal.doubleValue()));