きしだのはてな このページをアンテナに追加 RSSフィード

2008-11-24(月)

[][]BigDecimal版sqrtであらかじめMath.sqrtを使う効果 14:31 BigDecimal版sqrtであらかじめMath.sqrtを使う効果 - きしだのはてな を含むブックマーク

ところで、前回のGauss-Legendre法BigDecimal版sqrtを使いました。で、この中の処理である程度の精度までをdoubleで求めてるのですが、これの効果がどれだけあるか確認しました。


まず、doubleを使ったsqrtで100桁までの円周率を100回求めると、処理時間は0.54秒でした。

これをdoubleを使わず次のようなsqrtを使うと、3.22秒になりました。6倍くらい違います。

    public static BigDecimal sqrt(BigDecimal a, int scale){
        BigDecimal b2 = new BigDecimal(2);
        BigDecimal x = a;
        for(int tempScale = 1; tempScale < scale * 2; tempScale *= 2){
            //x = x - (x * x - a) / (2 * x);
            x = x.subtract(
                    x.multiply(x).subtract(a).divide(
                    x.multiply(b2), 
                    Math.min(scale, tempScale * 2), BigDecimal.ROUND_HALF_EVEN));
        }
        return x;
    }

というわけで、doubleを使って途中までの値を求めておくことには効果があることがわかりました。というか、想像以上に差があったんですが、なんでだろう?


ついでに、マチンの公式による方法で円周率100桁を100回求めると2.87秒でした。Gauss-Legendre法より速い。

これについては、桁数を増やすと挽回するんではないかと思います。

2008-11-23(日)

[][]ガウス・ルジャンドルの方法での円周率 12:35 ガウス・ルジャンドルの方法での円周率 - きしだのはてな を含むブックマーク

もう一個だけ、円周率を求めておく

算術幾何平均というのを使う。

これは、a, bについて

a_1=¥frac{a+b}{2}, b_1=¥sqrt{ab}

とすると

a_{n+1}=¥frac{a_n+b_n}{2}, b_{n+1}=¥sqrt{a_nb_n}

がnが無限のときanとbnが等しくなるというもの。


で、なんかいろいろ調べてたらどこからひっぱってきたか忘れたのだけど、次のようなプログラムで円周率が出せる。

public class PaiGauss {
    public static void main(String[] args){
        double a = 1;
        double b = Math.sqrt(2) / 2;
        double x = 1;
        double t = 1./4;

        for(int i = 0; i < 3; ++i){
            double y = a;
            a = (a + b) / 2;
            b = Math.sqrt(b * y);
            t = t - x * (y - a) * (y - a);
            x = x * 2;
        }
        System.out.println((a + b) * (a + b) / 4 / t);
        System.out.println(Math.PI);


    }
}

結果はこう。

3.141592653589794

3.141592653589793


なんと、3回の繰り返しでdoubleの限界まで求めれた。ただ、これにはちょっと落とし穴があって、sqrtを使っているので、この中にも繰り返しがある。


まあ、とりあえず、BigDecimalを使って100桁まで求めてみる。

これで、次のように100桁もとめれた。ちょっと改行した。

3.1415926535 8979323846 2643383279 5028841971 939937510

5820974944 5923078164 0628620899 8628034825 342117068

ソースはこれ

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

public class PaiGaussBig {
    public static void main(String[] args){
        MathContext mc = new MathContext(100);
        BigDecimal b2 = new BigDecimal(2);
        BigDecimal b4 = new BigDecimal(4);
        BigDecimal a = BigDecimal.ONE;
        BigDecimal b = sqrt(b2, mc.getPrecision()).divide(b2, mc); //Math.sqrt(2) / 2;
        BigDecimal x = BigDecimal.ONE;
        BigDecimal t = BigDecimal.ONE.divide(b4, mc);
        for (int i = 0; i < 6; ++i) {
            BigDecimal y = a;
            a = (a.add(b)).divide(b2, mc);
            b = sqrt(b.multiply(y), mc.getPrecision());
            t = t.subtract(x.multiply((y.subtract(a)).pow(2)));
            x = x.multiply(b2);
        }
        BigDecimal pi = a.add(b).pow(2).divide(b4.multiply(t), mc);
        System.out.println(pi);
    }

    private 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), Math.min(scale, tempScale * 2),
                    BigDecimal.ROUND_HALF_EVEN));
        }
        return x;
    }
}

2008-11-22(土) BigDecimalで平方根を求めてみる

[][]BigDecimalで平方根を求めてみる 13:31 BigDecimalで平方根を求めてみる - きしだのはてな を含むブックマーク

唐突に平方根を求めてみます。

平方根を求めるときは¥sqrt{a}を求めるとして

x^2-a=0

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


方程式の解はニュートン法で求めます。

404 Not Found

そうすると次のような計算を繰り返すと平方根が求まります。

a_2=a_1-(¥frac{a_1^2-a}{2a_1})


ということで、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略

nowokaynowokay 2008/11/22 17:07 divideの第二引数のscaleはtempScale*2にしたほうが効率がいいかも

nowokaynowokay 2008/11/23 11:50 Math.min(scale, tempScale * 2)のほうがいいかな

2008-11-16(日)

エビス・ザ・ホップ

[][]さらに円周率を求めてみる 21:07 さらに円周率を求めてみる - きしだのはてな を含むブックマーク

こないだ積分の練習に円周率を求めてみたんだけど、もうちょっとちゃんと円周率を求めてみる

とりあえず、モンテカルロでってことなので、やってみる。

円周率をモンテカルロで求めるときには、ランダムに選んだ点が円に入るかどうかで判定します。

つまりランダムにx、yを選んで

¥sqrt{x^2+y^2}<1

を満たす割合を求めます。

public class PaiMontecarlo {
    public static void main(String[] args){
        Random r = new MTRandom();
        int in = 0;
        int total = 100000;
        for(int i = 0; i < total; ++i){
            double x = r.nextDouble();
            double y = r.nextDouble();
            if(x * x + y * y < 1) ++in;
        }
        System.out.println((double)in / total * 4);
    }
}

そしたらこう。

3.14888

3桁精度出すのがやっと。計算結果として桁数を増やすにも実行回数の桁数を増やす必要がある。あと、乱数がだいたいちゃんとばらついていると円周率出るので、あんまり乱数の性能は関係ない気がする。


ということで、まじめに級数展開したもので計算。

0+¥frac{0+¥pi}{4}=1-¥frac{1}{3}+¥frac{1}{5}-¥frac{1}{7}+

(はてなtexモジュールの動きが変なので、0足してます。数学的な意味はないです。)

public class PaiNormal {
    public static void main(String[] args){
        double d = 0;
        double sig = 1;
        for(int i = 1; i < 1000000; i += 2){
            d += sig / i;
            sig *= -1;
        }
        System.out.println(d * 4);
        System.out.println(Math.PI);
    }
}

そしたらこう。

3.141590653589692

3.141592653589793

50万回計算してやっと小数点以下5桁がもとまった。


これはこれでちょっと収束遅すぎなので、オイラーの加速を使って収束を速くしてみる。

π/2=1+2/3×1/2+2・4/3・5×1/2^2+2・4・6/3・5・7×1/2^3+・・・

(もうどうやってもtexモジュールが表示されない)

¥frac{¥pi}{2}=1+¥frac{2}{2 ¥cdot 3}+¥frac{2¥cdot 4}{2^2 ¥cdot 3¥cdot 5}+¥frac{2 ¥cdot 4 ¥cdot 6}{2^3 ¥cdot 3 ¥cdot 5 ¥cdot 7}+¥cdots

public class PaiEularAcc {
    public static void main(String[] args){
        double pi = 1;
        double a = 1;
        double b = 1. / 2;
        for(int i = 2; i < 99; i += 2){
            a *= i / (double)(i + 1);
            pi += a * b;
            b /= 2;
        }

        System.out.println(pi * 2);
        System.out.println(Math.PI);
    }
}

そしたらこう。

3.1415926535897922

3.141592653589793

50回くらいでdoubleの限界まで来た。速い。


まあ、でもここまでならSICPにも書いてあるので、あまり面白くない。

ということでマチンの公式を使う。

¥frac{¥pi}{4}=4 tan^{-1} ¥frac{1}{5}-tan^{-1}¥frac{1}{239}

public class PaiMachin {
    public static void main(String[] args){
        double sig = 1;
        double a = 0;
        double b = 0;
        double s1 = 1. / 5;
        double s2 = 1. / 239;
        for(int i = 1; i < 22; i += 2){
            a += sig * s1 / i;
            b += sig * s2 / i;
            s1 /= 5 * 5;
            s2 /= 239 * 239;
            sig *= -1;
        }
        System.out.println(a * 16 - b * 4);
        System.out.println(Math.PI);
    }
}

そしたらこう。

3.141592653589794

3.141592653589793

10回程度の繰り返しでdoubleの限界まできた。


じゃあ、ということで、もちょっと桁数を増やす。doulbeの限界を超えるので、BigDecimalを使うことにする。

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

public class PaiMachinBig {
    public static void main(String[] args){
        MathContext mc = new MathContext(100);
        BigDecimal sig = BigDecimal.ONE; // 1
        BigDecimal a = BigDecimal.ZERO;  // 0
        BigDecimal b = BigDecimal.ZERO;  // 0
        BigDecimal d5 = new BigDecimal(5, mc); // 5
        BigDecimal d239 = new BigDecimal(239, mc); // 239
        BigDecimal s1 = BigDecimal.ONE.divide(d5, mc); // 1 / 5
        BigDecimal s2 = BigDecimal.ONE.divide(d239, mc); // 1 / 239
        BigDecimal d5_2 = d5.multiply(d5); // 5*5
        BigDecimal d239_2 = d239.multiply(d239); // 239 * 239
        System.out.println("開始");
        long start = System.nanoTime();
        for(int i = 1; i < 140; i += 2){
            BigDecimal bi = new BigDecimal(i, mc);
            a = a.add(s1.divide(bi, mc).multiply(sig)); // a + s1 / i * sig
            b = b.add(s2.divide(bi, mc).multiply(sig)); // b + s2 / i * sig
            s1 = s1.divide(d5_2, mc);   // s1 / 5^2
            s2 = s2.divide(d239_2, mc); // s2 / 239^2
            sig = sig.negate(); // -sig
        }
        System.out.println((System.nanoTime() - start) / 1000 / 1000 / 1000.);
        a = a.multiply(new BigDecimal(16, mc));
        b = b.multiply(new BigDecimal(4, mc));

        System.out.println(a.subtract(b, mc));
    }
}

こんな感じで出た(10桁で区切ってます)

3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 342117068

75回程度でここまで来た。マチンすごい!

女性ホルモンの分泌を増やしてバストアップに載ってるのと同じ数字が出た。

ループの条件をi < 300くらいにして150回繰り返すようにすると、200桁までちゃんと求まる。

けど、BigDecimalのプログラムめんどすぎ泣ける。Scala使ったほうがいいね。

ところで、円周率100万桁のページを表示するとFireFox3のCPU使用率が100%になる。どうにかならんもんか

2008-11-13(木)

エルディンガー オクトーバーフェスト

[][]フェルマーテストは乱数の性能にあんまり影響されない 17:37 フェルマーテストは乱数の性能にあんまり影響されない - きしだのはてな を含むブックマーク

試しにフェルマーテストが乱数の質で性能に違いがあるかやってみました。

1千万までの合成数について、何回の判定で合成数と認識されるかというのを、100回やってみました。

1回 - 433428013

2回 - 97607

3回 - 7565

4回 - 2207

5回 - 1160

中略

213回 - 1

223回 - 2

231回 - 1

320回 - 1

平均1.001判定 分散 0.0130 偏差 0.1141


Javaの標準Randomと、メルセンヌ・ツイスターで試してみましたが、だいたいおんなじ感じでした。実行時間も、ちょっとメルセンヌ・ツイスターが遅いくらい。

結局、フェルマーテストでは、ほとんどの合成数が一回の判定で認識できてしまっているので、あまり乱数の質の影響を受けないのではないかと思います。また、200回を超えても判定できてないものがいくつかあるので、フェルマーテストでは判定回数が数百回程度では信頼性があまり高いとはいえないです。


実際には、数回の判定をすりぬけるのはカーマイケル数で、そのカーマイケル数も多くは19までの素数の倍数になっているので、19までの素数の倍数になっていないかチェックすることと、19までの素数を約数に持たないカーマイケル数を配列に持たせて判定することで、intの範囲では実用的になると思います。

まあ、intの範囲なら、√21億≒46500までの素数を求めておいて、その全部で割れるか試してみればいいだけの話ですが。


ということで、ソースはこれ。Core Duo L2500でだいたい15分くらいかかりました。

手軽に試すには、eratosの引数を100000くらい、iのループを10回くらいにするといいと思います。

import java.util.*;

public class FermetEfficient {
    public static void main(String[] args){
        r.nextLong();
        boolean[] comp = eratos(10000000);
        Map<Integer, Integer> counter = new TreeMap<Integer, Integer>();
        for(int i = 0; i < 100; ++i){
            for(int n = 9; n < comp.length; n += 2){
                if(!comp[n]) continue;//素数はとばす
                int count = fermetCheckCount(n);
                if(!counter.containsKey(count)) counter.put(count, 0);
                counter.put(count, counter.get(count) + 1);
            }
        }
        int sum = 0;//判定回数
        int total = 0;//試行回数
        for(int i : counter.keySet()){
            int c = counter.get(i);
            System.out.printf("%d回 - %d%n", i, c);
            total += c;
            sum += i * c;
        }
        double avg = (double)sum / total;//平均
        //分散を求める
        double disp = 0;
        for(int i : counter.keySet()){
            int c = counter.get(i);
            disp += (avg - i) * (avg - i) * c;
        }
        System.out.printf("平均%.3f判定 分散 %.4f  偏差 %.4f%n", 
            avg, disp / total, Math.sqrt(disp / total));
    }
    static Random r = new Random();

    /** エラトステネスの篩により合成数ならtrue/素数ならfalseの配列 */
    private static boolean[] eratos(int n){
        boolean[] ret = new boolean[n + 1];
        ret[0] = true;
        ret[1] = true;
        for(int i = 4; i <= n; i += 2){
            ret[i] = true;
        }
        for(int i = 3; i <= n; i += 2){
            if(ret[i]) continue;
            for(int j = i * 2; j <= n; j += i){
                ret[j] = true;
            }
        }
        return ret;
    }

    /** フェルマーテストで合成数と判定されるまでの回数 */
    private static int fermetCheckCount(int n){
        for(int i = 1; ; ++i){
            if(n % 2 == 0) return i;//偶数なら合成数
            int a = r.nextInt(n - 2) + 2;//2以上n未満の値
            int g = (int) gcd(n, a);
            if(g != 1) return i;//最大公約数がみつかったなら合成数
            if(modPow(a, n - 1, n) != 1) return i;
        }
    }

    /** (a ^ b) % m を計算 */
    static int modPow(int a, int b, int m){
        if(a == 0) return 0;

        int result = 1;
        a %= m;
        while(b != 0){
            if(b % 2 == 1){
                result = (int)(((long)result * a) % m);
            }
            a = (int)(((long)a * a) % m);
            b = b >> 1;
        }

        return result;
    }

    /** 最大公約数を求める(a > b) */
    private static long gcd(long a, long b){
        if(b == 0) return a;
        long m = a % b;

        return gcd(b, m);
    }
}