【Java】べき乗計算を高速化する(繰り返し二乗法)

ヽ｜∵｜ゝ(Fantom) の開発blog？ホーム » Java »【Java】べき乗計算を高速化する(繰り返し二乗法)

【Java】べき乗計算を高速化する(繰り返し二乗法)

2016/02/04 Thu [edit]

　どちらかというとプロコン用ライブラリという使い方になるかもしれないが、大量のべき乗計算が必要になったとき、高速に計算する方法を知っておけば、イザというとき役に立つかもしれない。

　今回紹介する方法は「繰り返し二乗法」という算法である。数学的に解説してもいいが、ググればたくさん出てくるので、理論はさておき、大雑把な方法論だけ覚えてしまおう(笑)。

●繰り返し二乗法によるべき乗計算

/**<h1>べき乗計算(繰り返し二乗法)</h1>
 * <p>Math.pow() の約5倍速い。</p>
 * @param x ： 底
 * @param n ： べき指数
 * @return<b>long</b> ： x<sup>n</sup>
 */
public static final long pow(long x, long n) {
    long sum = 1;
    while (n > 0) {
        if ((n & 1) == 1) {  //ビットが立っていたら
            sum *= x;        //x^2^k1 * x^2^k2 * x^2^k3 * ・・・
        }
        x *= x;              //((x^2)^2)^2・・・
        n >>= 1;
    }
    return sum;
}

　簡単に解説すると、使うのは主に指数法則である。例えば、べき乗の計算をするとき、for ループなどで x を n 回かける場合、n の値が大きいほど計算時間がかかるので、指数法則を用いて、((x²)²)² ･･･のように２乗を繰り返して、計算回数を減らそうというアルゴリズムである。この場合オーダーは O(log n) となるので、大量の計算でも高速にできる。

　ではその２乗する n をどうやって分解したらいいかというと、２進法を使うのである。というのも２進数というのは元々それぞれのビットが、2⁰=1、2¹=2、2²=4、2³=8、2⁴=16、･･･を表しているので、指数 n のビットを見れば、どの n = 2^k をかければ良いのかわかるのである。言葉で説明すると難しく感じるので、具体的な数値を指数法則で分解したり、２進数と照らしあわせて見て欲しい。

　つまり指数 n のビットが立っているかを下位から判別していき、立っていたら x を 2^k 乗した値をかけていくのである。これが x¹ × x² × x⁴ × x⁸ × x¹⁶ ･･･となっていく。なんとなく数学的によくわからないと思ったら、以下の資料などを見てみると良いだろう。わかってしまえば、そう難しいものでもない。

（参考）
・繰り返し二乗法
・冪乗
・繰り返し自乗法

　ちなみに実行速度は Math.pow() の約５倍ほどの速度で計算できる。

　また、値が long で計算されているが、指数の n 以外を double にしても良い。その場合は、17桁以降が丸められる感じになる（Math.pow() とは少し違う丸め値になる）。double にした場合、Math.pow() の約３倍ほどの速さになるようだ。あと、負の値は考慮してないので、必要あれば例外処理などを加えても良いだろう。

　また、プロコンのような問題では大抵「1000000 の剰余で答えなさい」のようになっていると思うが（普通に計算したらオーバフローするため）、その方法も書いておく。剰余は「余りをとる」「mod をとる」「法とする」などと表記されることもある。

●繰り返し二乗法によるべき乗計算で剰余を返す版

/**<h1>べき乗計算＆剰余(繰り返し二乗法)</h1>
 * <p>Math.pow() の約5倍速い。</p>
 * @param x ： 底
 * @param n ： べき指数
 * @param mod ： 剰余
 * @return<b>long</b> ： x<sup>n</sup> % mod
 */
public static final long modPow(long x, long n, final long mod) {
    long sum = 1;
    while (n > 0) {
        if ((n & 1) == 1) {      //ビットが立っていたら
            sum = sum * x % mod; //x^2^k1 * x^2^k2 * x^2^k3 * ・・・
        }
        x = x * x % mod;         //((x^2)^2)^2 ・・・
        n >>= 1;
    }
    return sum;
}

　long は 19桁なので、9桁×9桁=18桁なら、1000000000 の剰余で計算すれば、下9桁でオーバフローすることはない。検算したければ、BigInteger.modPow() と照合してみれば良いだろう。上の関数は同じ値になる。ちなみに実行速度は BigInteger の約25倍ほど速い。