ミラーラビン素数判定法を大雑把に解説してみる

素数判定法

素数判定法には様々なものが知られている。

大雑把に分けると決定的素数判定法と確率的素数判定法がある。

決定的素数判定法は “素数である” と判定した場合には100%素数であるが、確率的素数判定法は “素数である” と判定した場合にも合成数である可能性が存在する。

確率的素数判定法

これだけ見ると確率的素数判定をわざわざ使う必要がないようにも思えるが、実際にはそうではない。

むしろ、現実で使われている素数判定のほとんど全ては確率的素数判定法である。

確率的素数判定法が使われる理由

何故100%判定できる判定法があるのに、そちらを使わないのかというとざっくり言ってしまえば「厳密にチェックするのは時間が掛かるが、大雑把でいいなら早くできるから」ということに尽きる。

実際の暗号で使われる素数は非常に巨大(200桁以上ある!)なので、それが素数かどうかをいちいち厳密にチェックするのはやりとりするメッセージの重要性に大してコストが掛かりすぎるのだ。

「なうめしなう」なんていうメッセージが流出したところで全く困らないでしょ?つまりそういうことです!

確率的素数判定法の信頼度

詳しく話すと長くなるので割愛するが、今回紹介する確率的素数判定法は “素数である” という判定はほんの僅かに誤りである可能性があるが、 “素数ではない” と判定した場合にはそれは100%素数ではないことを保証してくれている。

つまり、確率的ではあるものの判定のうちの一方は決定的であるというわけだ。

さて、確率的素数判定法のおおよその概念がわかったところでこの “ほぼ確実” というのがどのくらい確実なのかが気になるところである。

日常的な尺度でいえば99%あれば十分信頼できそうな感じがするが、確率的素数判定法は一回試して “素数だと判定したが実は素数ではなかった” と間違える確率は約25%である。

つまり、正確性は75%ということになる。

これだと全然ダメじゃないかという気がするが、これは一回試した場合のことである。

例えば、確率的素数判定法を使う機械をたくさん用意して(例えば百台としよう)、その全てが “ほぼ確実に素数である” と判定すれば、一台のときと比べて信頼度が高いのがわかるだろうか?

何故なら、どれか一台でも “素数ではない” と判定すれば絶対に素数でないにもか関わらず、百台全てが “「絶対に素数ではない」とは言い切れない” と言っているからである。

このとき全ての機械が偶然にも誤る確率は非常に低く、この素数と判定された数が実際に素数である確率は、

と見積もることができる。

これは99.999(9が60個つづく)%という、まさに “ほぼ確実” にふさわしい精度を持っているのである。

素数判定アルゴリズム

素数判定アルゴリズムは数多く知られている。

  • 決定的素数判定法
    • 試し割り法
    • エラトステネスの篩
    • 最大公約数
    • AKS素数判定法
  • 確率的素数判定法
    • フェルマーテスト
    • ミラーラビン素数判定

今回はこれらの中から試し割り法とミラーラビン素数判定法を紹介する。

決定的素数判定法

試し割り法

一般的だが、非常に時間のかかるアルゴリズムである。

ただし、チェックしたい数が小さい場合には極めて有効である。

ある数Nを試し割りするアルゴリズムが何回の除算を必要とするかを大雑把に評価すると次のようになる。

全通り試せばN回割り算が必要なので、このときの計算回数は以下のように表せる。

Nまでも平方で割ればいいことに気付けばこれがかなり減らすことができ、

となる。

更に素数リストを使えば計算量を減らすことができる。

ある数Nまでにどのくらい素数があるかというのは素数定理に依存するのがだが、これによるとあるNまでの素数の数π(N)は以下のように近似できる。

今回はこのxにNの平方根が入るので、計算量は以下のようになる。

big Oは定数倍は無視できるので、これは以下のように書き直すことができる。

N通常平方根素数リスト
100001000010011
10000010000031627
10000001000000100072
10000000100000003162196
10000000010000000010000543

Nが1億になっても素数リストはおよそ543回の除算で済むことがわかり、単純な平方根と比べても1/20程度であることがわかる。

確率的素数判定法

ミラーラビン素数判定法

高速に素数かどうかを確率的に判定できる方法であるが、まずはこのミラーラビン素数判定法とはなんぞやらというのを明らかにしていこう!

ミラーラビン素数判定法

ミラーラビン素数判定法はフェルマーの小定理の対偶を応用した確率的素数判定アルゴリズムである。

フェルマーの小定理を利用するということはまずはこっちを理解せねばなるまい!

フェルマーの小定理

pを素数とし、aをpと互いに素な数とするとき以下が成立する。

これをフェルマーの小定理という。

今回使うのはこの対偶なので、aと互いに素なある数nが、

を満たせばnは素数ではないということが言いたいわけである。

このフェルマーの小定理の対偶を利用した素数判定アルゴリズムをフェルマーテストなどと言ったりする。

これは必要条件は与えているが、十分条件は与えていないので計算結果が1でなければ “間違いなく素数でない” ということは保証してくれるが、計算結果が1だったものが素数であることは保証してくれない。

もちろん互いに素な数aをいろいろ試しても “素数ではない” という結果が返ってこなければその数nは素数である可能性が高くなるが、如何なる互いに素な数aを用いても誤って素数と判定される数がある。

それらの数はカーマイケル数と呼ばれ、最小の数が561で割合は少ないものの無限に存在することがわかっている。

カーマイケル数は10^21までの範囲でみるとおよそ5*10^13個に1つの割合らしい。割合で見るとめちゃくちゃ少ないのだが20,138,200個もあると考えるとあんまり安心できない気もする。

出典: フリー百科事典『ウィキペディア(Wikipedia)』

自明な平方根

まず、有限体 Zpの単位元の平方根についての補題を考える。ここで、pは奇素数である。pを法とした剰余(mod p)において、1と-1の二乗は常に1となる。これらを1の「自明な」平方根と呼ぶ。pは素数なので1の「自明でない」平方根は存在しない。これを示すため、1 mod p の非自明な平方根を x とする。

さて、1の自明でない平方根とはなんだろうか?

まず、自明な平方根を考える。

平方根を考える前にもっと簡単な以下の式を考えよう。

上式を満たすxは何が考えられるだろうか?

少し考えればそれが1とn+1であることがわかるはずだ。

xを有限体Zpから選べば1しかないことがわかる。

わかったところで自明な平方根を考えよう。

これは式変形して因数分解すると次のように変形できる。

剰余が0ということはx+1かx-1がnで割り切れるということを意味する。

しかし、nが素数であればそのようなx+1やx-1が存在するはずがないので、1と-1であることがわかる。

ここでいう-1とはn-1のことを意味している。

これはnがどんな数(素数や合成数)であっても必ず平方根となりうるので “自明” というわけである。

例えば541の自明な1の平方根は1と540である。

nが素数の場合は非自明な平方根が存在しない、というのがミラーラビン素数判定法のキモである。

このあたりはまた後日詳しく解説したい。

フェルマーの小定理の利用

ではここで、調べたい奇素数をnとすると以下のように式変形できる。

このとき、次のどちらかが成立する。

これは証明が長くなるので割愛するが、要するに奇素数であればこれが成立するのである。

ただし、aは以下の式を満たす数である。

素数判定のアルゴリズム

これまでの流れはこうである。

  • 素数の自明な平方根は1と-1(その数自身-1)しかない。
  • フェルマーの小定理を利用すると、奇素数は二つの式のうちどちらかを満たす(証明は省略)

つまり、これの対偶をとれば「二つの式のどちらも満たさなければ与えられた数は素数ではない」と言えるのである。

確かめるべき二つの式

ある数nに対していろいろなaを試しても “素数ではない” という結果が返ってこなければその数が素数である可能性が高くなるという仕組みである。

満たせば素数ではない
100以下のn, r, dのリスト

しかも、これはnがある値よりも小さければそれほど多くのaを試す必要がないことがわかっている。

na
n < 20472
n < 1,373,6532, 3
n < 9,080,19131, 73
n < 25,326,0012, 3, 5
n < 3,215,031,7512, 3, 5, 7
n < 4,759,123,1412, 7, 61
n < 1,122,004,669,6332, 13, 23, 1662803

なんとaとして2, 7, 61の3つを試せば47億以下の全ての数が素数かどうかを正しく判定できるのだ。

47億というのは暗号強度的には全くダメだが、現実問題として素数かどうか調べるプログラムとしては十分すぎる大きさである。

そもそもint型がせいぜい4.3億くらいまでしか扱えないというのもあるが…

ちなみにもっと大きいアドレス空間(例えばint_64など)であってもaとして3, 5, 7, 11, 13, 17, 19, 23を用いれば3,825,123,056,546,413,051までの数が判定できるので安心である。

From Wikipedia, the free encyclopedia

まとめ

今回はnの範囲を制限することでミラーラビン素数判定法を決定的アルゴリズムに変えるというのを目標にしている。

n < 4,759,123,141においては以下が成立する。

これの対偶をとると、

notの論理演算を行なえば以下の式が成立するので、

これを用いると、

のように変形できる。

コードを書いてみる

ミラーラビン素数判定で最も計算時間がかかるのは累乗の剰余を計算する箇所である。

そのまま実装するとあっという間にオーバーフローするのでバカ正直に書いてはいけない。

int nmodpow(int a, int d, int n)
{
    int k = a % n;
    for (int i = 1; i < d; i++)
    {
        k = k * a % n;
    }
    return k;
}

上記のコードのように毎回剰余をとる方法もあるが、ここはバイナリ法を使うと高速に実装できる。

バイナリ法のアルゴリズム

  1. a^dを計算したいとする。
  2. 一時計算用の値p=aと返り値v=1を設定する
  3. dの最下位桁が1ならv=v*pとする
  4. n/=2とし、p=p*pとして2.に戻る

要するに指数を二進数展開し各ビットをみて、1であれば返り値をどんどん大きくしていくという仕組みである。

これは指数をnとすると、二進数展開したときの桁数くらいしか乗算を行わないので効率的である。

int bpow(int a, int d, int n)
{
    int p = a;
    int v = 1;

    // 指数が0でない間くり返す
    while (d)
    {
        // 最下位ビットを見る
        if (d &amp; 1)
        {
            v = v * p;
        }
        // 半分にする
        d >>= 1;
        p *= p;
    }
    return v;
}

ここでいろいろ高速化のコードが活かされている。

例えば、二進数展開したときに最下位ビットが1かどうか確認するというのは “奇数かどうか調べる” ということに等しいので、

d % 2 == 1

という書き方もできるのだが、これは剰余計算をしているために遅い。

このコードのように1と論理積をとればdの最下位ビットが1ならば1、違えば0が返るという仕組みである。

コンピュータは根本的には二進数で処理をしているのでビット演算をするのが最も手っ取り早い。

また、2で割って新たなdを設定するところも、

d /= 2;

d *= 0.5

のような書き方もできるが、やはりこれは乗算や除算命令を伴うために遅い。

2で割るというのは右に1ビットシフトするのと同じなので、あのように書いたほうが断然速い。

賢いコンパイラなら上手く解釈してくれるかもしれないが…

これを剰余をとるbmodpow()に変更する。

__int64 bmodpow(int a, int d, int n)
{
    __int64 p = a;
    __int64 v = 1;

    while (d)
    {
        if (d &amp; 1)
        {
            v = (v * p) % n;
        }
        p = (p * p) % n;
        d >>= 1;
    }
    return v;
}

では、このバイナリ法がどのくらい速いか考えてみよう。

指数をNとするとnmodpow()は単純にN回の乗算と剰余計算を行うので、計算量は以下のように表せる。

それに対してbmodpow()は指数Nの二進数展開の桁数程度しか乗算と剰余計算を行わないので、

で済む計算になる。

バイナリ法はどのくらい速いか

int a = 9999;
    int sum = 0;
    clock_t st = clock();
    for (int d = 3; d <= 50000; d++)
    {
        sum += npowmod(a, d, d - 1);
    }
    cout << clock() - st << " sum: " << sum << endl;

    sum = 0;
    st = clock();
    for (int d = 3; d <= 50000; d++)
    {
        sum += bpowmod(a, d, d - 1);
    }
    cout << clock() - st << " sum: " << sum << endl;
    return 0;

9999の3乗から10000乗までの剰余の和をとるプログラムである。

nmodpow()bmodpow()
処理時間13.37s0.01s
621500601621500601

というわけでバイナリ法のほうが1300倍も速い結果になりました。

コード実装

実際に動くコードがこちら。

#include <iostream>
#include <vector>
#include <time.h>

#define N 1000000

using namespace std;

__int64 bmodpow(int a, int d, int n)
{
    __int64 p = a;
    __int64 v = 1;

    __int64 x = d &amp; 1;

    while (d)
    {
        if (d &amp; 1)
        {
            v = (v * p) % n;
        }
        p = (p * p) % n;
        d >>= 1;
    }
    return v;
}

bool mr(int a, int d, int s, int n)
{
    size_t x = bmodpow(a, d, n);

    if (x == 1)
        return true;
    while (s)
    {
        if (x == n - 1)
            return true;
        x = (x * x) % n;
        --s;
    }
    return false;
}

int main()
{
    double st = clock();
    int p = 2;
    int n;

    for (n = 3; p <= N; n += 2)
    {
        if (n == 7 || n == 61)
        {
            ++p;
            continue;
        }
        __int64 d = n >> 1;
        __int64 s = 1;
        while (!(d &amp; 1))
        {
            d >>= 1;
            ++s;
        }
        if (mr(2, d, s, n) &amp;&amp; mr(7, d, s, n) &amp;&amp; mr(61, d, s, n))
            ++p;
    }
    cout << n - 2 << " is " << p - 1 << " th prime number." << endl;
    cout << (clock() - st) / CLOCKS_PER_SEC << endl;
    return 0;
}

mr()の内部でbmodpow()が動作しており、mr()はs回のループをする。

sはどんなに大きくてもlogNで抑えることができるので、ミラーラビン判定法は大雑把に、

程度の計算量と見積もることができる。

これはnが47億以下という条件の下での計算量であり、実際にはaは2, 7, 61以外にももっと多くの数を試さなければいけない。

拡張リーマン予想が真であるなら(数学者の多くはリーマン予想は恐らく正しいと考えている)、素数かどうかを確実にチェックできるaの集合が(lnN)^2より小さい元から必ずつくることができる。

拡張リーマン予想が真であるなら、ミラーラビン素数判定法は決定的に(確率的ではない!)素数判定が可能ということになる。

このときの計算量は以下の計算式で大雑把に(logN)^4で抑えられることになり、これはwikiの値と一致する。

これは、数が十分小さいうちは素数リストのほうが速いことを意味するが、あるNでミラーラビン判定法のほうが高速になる。

そのNの値は以下の不等式を解くことで得られ、

これを解くとおよそ、

という値が得られる。

つまり、調べたい数が2400万を超えた辺りからミラーラビン判定法の方が高速になるということである。

実際には1000万番目の時点でわずかにミラーラビンが素数リストを上回った。

まあ、オーダーは定数倍は無視するのでこの辺りは誤差の範囲だろう。

最後に

ぼくは数学の専門家ではないので、論理に間違いがあったりした場合は優しく(←大事)間違いを教えてください!!