BigDecimalで平方根を求めてみる

唐突に平方根を求めてみます。
平方根を求めるときはを求めるとして

という方程式の解を求めればいいことになります。

方程式の解はニュートン法で求めます。
http://www.akita-pu.ac.jp/system/elect/comp1/kusakari/japanese/teaching/Programming/2005/note/6/Slide03.html
そうすると次のような計算を繰り返すと平方根が求まります。

ということで、doubleでやってみる。

public class Sqrt {
    public static void main(String[] args){
        System.out.println(sqrt(2)+"^2="+sqrt(2)*sqrt(2));
        System.out.println(
            Math.sqrt(2)+"^2="+Math.sqrt(2)*Math.sqrt(2));
    }

    public static double sqrt(double a){
        double x = a;
        for(int i = 0; i < 5; ++i){
            x = x - (x * x - a) / (2 * x);
        }
        return x;
    }
}

Math.sqrtと同じ値が求まりました。

1.4142135623730951^2=2.0000000000000004
1.4142135623730951^2=2.0000000000000004

で、doubleでMath.sqrtと同じ値が求めれてもあまりうれしくないので、BigDecimal版を作ります。ここで、ニュートン法は繰り返しごとに精度が倍になるので、16桁まではdoubleで求めておくと最初の数回の繰り返しが減らせます。

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

public class SqrtBig {
    public static void main(String[] args){
        BigDecimal s2_50 = sqrt(new BigDecimal(2), 50);
        System.out.println(s2_50);
        System.out.println(s2_50.multiply(s2_50));
    }

    public static BigDecimal sqrt(BigDecimal a, int scale){
        //とりあえずdoubleのsqrtを求める
        BigDecimal x = new BigDecimal(
                Math.sqrt(a.doubleValue()), MathContext.DECIMAL64);
        if(scale < 17) return x;

        BigDecimal b2 = new BigDecimal(2);
        for(int tempScale = 16; tempScale < scale; tempScale *= 2){
            //x = x - (x * x - a) / (2 * x);
            x = x.subtract(
                    x.multiply(x).subtract(a).divide(
                    x.multiply(b2), scale, BigDecimal.ROUND_HALF_EVEN));
        }
        return x;
    }
}

計算結果はこう。50桁まで求まっています。

1.41421356237309504880168872420969807856967187537695
2.000000000000000000000000000000000000000000000000005449略