絶滅

どうでもいい

バイナリ法

むかしむかし,あるところにおじいさんとおばあさんがいました.おじいさんは山へ芝刈りに,そしておばあさんは 3 の 1000000000 乗 を 1000000007 で割った余りを求めることにしました.




方法 1: 単純な for ループを回す

一番単純な方法は 3 を 1000000000 回かける方法で,O(n) なので計算機でもまあまあ時間がかかり,おばあさんは死にます.

typedef long long lint;

//n^p % m
lint bigModExpNaive(lint n, lint p, lint m)
{
    lint res = 1;
    for (lint i = 1; i <= p; i++) { // ループ
        res *= n; // n をかける
        res %= m; // m で割った余りを求める
    }
    return res;
}
int main()
{
    clock_t start, end;
    start = clock();

    cout << "ans = " << bigModExpNaive(3, 1000000000, 1000000007) << endl;

    end = clock();
    printf("%d msec.\n", end - start);

    return 0;
}
ans = 235939645
23206 msec.






方法 2: バイナリ法を使う


乗算を工夫することで時間計算量 O(log n) で済むバイナリ法という方法がある.

指数 p が 2 のべき乗,つまり p = 2k と表される場合は k 回の乗算で済む.

(k = 0) 320 = 3
(k = 1) 321 = (320)2 = 32 = 9
(k = 2) 322 = (321)2 = 92 = 81
(k = 3) 323 = (322)2 = 812 = 6561
(k = 4) 324 = (323)2 = 65612 = 43046721

青字が再利用部分.3 を 2 乗して 9,9 を 2 乗して 81,81 を 2 乗して 6561,と繰り返すことで 32k の値を次々に求めることができる.


p = 1000000000 の場合について,

p = 100000000010
 = 0011 1011 1001 1010 1100 1010 0000 00002
 = 229 + 228 + 227 + 225 + 224 + 223 + 220 + 219 + 217 + 215 + 214 + 211 + 29

なので,f(k) = 32k とすると,

3p = 3229 + 228 + 227 + 225 + 224 + 223 + 220 + 219 + 217 + 215 + 214 + 211 + 29
  = f(29) * f(28) * f(27) * f(25) * f(24) * f(23) * f(20) * f(19) * f(17) * f(15) * f(14) * f(11) * f(9)

が成り立つ.

このとき乗算回数は f(29) までを求める 29 回と f(k) 同士の積をとる 12 回の合計たったの 41 回.よかったですね.





以上を一般化し, n を p 乗して m で割った余りを求める bigModExp のコードは以下.

// n を p 乗して m で割った余りを求める
lint bigModExp(lint n, lint p, lint m) {
    
    //返り値
    lint res = 1;

    // i = 1 = 2^0 (二進表記で ...0000 0000 0001) から始めて
    // ビットを一つずつ左シフトしていく.
    for (lint i = 1; i <= p; i <<= 1) {

        // k 回目のループにおいて
        // p の二進表記の 2^(k - 1) の位にビットが立っているか?
        if (p & i) {
            //f(k - 1) を返り値に乗算する
            res = res * n % m;
        }
        //f(k) を求める.
        n = n * n % m;
    }
    return res;
}


で,これをちょいひねると再帰バージョンもできる.(乗算および剰余演算の回数はまったく同じ)

// n を p 乗して m で割った余りを求める (再帰)
lint bigModExpRec(lint n, lint p, lint m)
{
    //n^0 は 1
    if (p == 0)  return 1;

    // n を (p / 2) 乗して m で割った余りを求める
    lint ret = bigModExpRec(n, p >> 1, m);

    // n を (p / 2) * 2 乗して m で割った余り
    ret = ret * ret % m;

    // p の最下位ビットが 1,つまり奇数の時は
    // (p / 2) * 2 == p - 1 なので
    // n を乗算してから返す
    if (p & 1)
        return ret * n % m;
    else
        return ret;
}


3p % 1000000007 の計算を 0 ≦ p ≦ 10000000 で回して実行時間を計測してみる.

#define MOD 1000000007

int main()
{
    lint N = 10000000;
    clock_t start, end;

    start = clock();
    for (lint p = 0; p <= N; p++)
        bigModExp(3, p, MOD);
    end = clock();
    printf("非再帰 : %d msec.\n", end - start);

    start = clock();
    for (lint p = 0; p <= N; p++)
        bigModExpRec(3, p, MOD);
    end = clock();
    printf("再帰 : %d msec.\n", end - start);

    return 0;
}
非再帰 : 6522 msec.
再帰 : 14274 msec.

関数呼び出しによるオーバーヘッドの影響か再帰は二倍ぐらい遅いっぽい.





おまけ:

f:id:my316g:20170831081652j:plain:w500