簡単な素数判定法

この記事は2008年に発表されました。

オリジナル記事を読み、私のニュースレターを購読するには、[**ここ**](https://quaily.com/lyric_kiji/p/simple-number-detection-method) でご覧ください。

---

SICPのセクション1.3には素数を判定する方法が記載されています。学生たちに広く知られている単純な素数判定法は、nの最小の約数が1以外で自分自身であれば、nは素数であるというものです。

```lisp
(define (remainder num divisor)
    (- num (* divisor (round (/ num divisor)))))

; nの最小の約数を返す
(define (smallest-divisor n) 
    ; nの約数(2からnまで)を探す
    (define (find-divisor n test-divisor)
        (cond ((< n (sqrt test-divisor)) n)
            ((= (remainder n test-divisor) 0) test-divisor)
            (else (find-divisor n (+ test-divisor 1)))))

    (find-divisor n 2))

; nが素数かどうかをチェックする
(define (prime? n)
    (= n (smallest-divisor n )))
```

しかし、この方法は遅く、複雑度はO(n)です。

SICPでは次にフェルマーの小定理が紹介されています：

> もしnが素数であり、a<nの任意のaについて、a^n≡a(mod n)である。

残念ながら、逆定理は成り立ちませんが、それでも有用です。nより小さい数aを選び、上記のテストを行います。もしa^n≡a(mod n)が成り立つなら、nは非常に高い確率で素数です。より多くのaをnに対してテストするほど、nが素数である可能性が高くなります。したがって、これは確率的なアルゴリズムであり、十分なテスト回数を行えば、nが素数でない確率は非常に小さくなります。

これにより、フェルマーのテストプログラムが導かれます。a^nの計算効率を心配する必要はありません：

考えてみましょう：

```lisp
a^n = (a^(n/2))^2; nが偶数の場合
a^n = a * a^(n-1); nが奇数の場合
```

これは常に成り立つので、ほぼすべての乗算でnが2倍になります。数直線を見てみると、これは基本的に二分探索です。したがって、時間計算量が減少します。これにより、プログラミング言語が提供する`pow`関数が自分で`a*a*a…`と書くよりも常に速い理由が説明されます。

しかし、直接組み込みの`pow`関数を使用すべきではありません。理由は、我々の目標はaのn乗ではなく、aのn乗をmで割った余りがaの余りに等しいかどうかを確認することだからです。簡単に言えば、この式が成り立つかどうかを確認するだけです：`a^n mod m = a`。もし`pow(a, n) % m == a`と直接表現すると、`pow(a, n)`の値が非常に大きくなる可能性があります。したがって、同様の技法を使用する必要があります：

```lisp
a^n mod m = ((a^(n / 2))^2) mod m; nが偶数の場合
a^n mod m = (a * a^(n-1)) mod m; nが奇数の場合
```

この方法では、計算された値はmを超えることはなく、非常に大きな数値での計算を避けることができます。したがって、次のプログラムが簡単に得られます：

```lisp
(define (remainder num divisor)
    (- num (* divisor (round (/ num divisor)))))

(define (expmod base exps m)
    (define (square n) (* n n))
    (cond ((= exps 0) 1)
        ((even? exps)
            (remainder (square (expmod base (/ exps 2) m))
                m))
        (else
            (remainder (* base (expmod base (- exps 1) m))
                m))))

; フェルマーのテスト、nをチェックする。もし(a^n) mod n == aなら、ほぼ間違いなくnは素数
(define (fermat-test n)
    (define (try-it a)
        (= a (expmod a n n)))
        ;(= a (remainder (expt a n) n)))
    (try-it (+ 1 (random (- n 1)))))

; フェルマーのテストで素数を取得する。nが素数かどうかをチェックする。試行回数timesを試す
(define (fermat-prime? n times)
    (cond ((= 0 times) #t)
        ((fermat-test n) (fermat-prime? n (- times 1)))
        (else #f)))
```

しかし、完全ではありません。

フェルマーのテストを欺くことができるカーマイケル数という種類の数が存在します。2から10000の間に、いくつかのこのような数があります：

> 561, 1105, 1729, 2465, 2821, 6601, 8911

これらの頑固な数は合成数ですが、どのaを使ってもフェルマーのテストにおいて素数のように振る舞うことができます。彼らは数が少ないわけではなく、広く分布しており、100,000,000の範囲内に255のカーマイケル数があります。したがって、カーマイケル数のリストを維持したくない場合、素数判定方法を強化することができます。

幸運なことに、エクササイズ1.28でミラーとラビンがミラー・ラビンテスト方法をもたらしました。この方法はフェルマーの小定理の変形に基づいています：

> nが素数なら、a < nの任意のaについて、a^(n-1) ≡ 1(mod n)

MRテストを使用するには、フェルマーのテストの平方ステップで「非自明な平方根1 mod n」を確認する必要があります。もしs = a^2; とs <> 1およびs <> n-1が存在し、s^2 mod m = 1なら、nは素数ではありません。

> [Matrix67](http://www.matrix67.com/blog/archives/234)：ミラー・ラビン素数判定法も確率的アルゴリズムであり、基数aに対するミラー・ラビンテストを通過する合成数を基数aに対する強擬素数と呼びます。基数2に対する最初の強擬素数は2047です。基数2および3に対する最初の強擬素数は1,373,653のように大きいです。

したがって、十分なテスト回数があれば、この方法は工業用途において安全性を保証できます。以下は私の質問1.28に対する答えです：

```lisp
(define (remainder num divisor)
    (- num (* divisor (round (/ num divisor)))))
    
(define (mr-expmod base exps m)
    (define (square x) (* x x))
    (define (check-square-root x n)
        (define (check r)
            (if (and (not (= x 1)) (not (= x (- n 1))) (= r 1)) 0 r))
        (check (remainder (square x) n)))

    (cond ((= exps 0) 1)
        ((even? exps)
            (check-square-root(mr-expmod base (/ exps 2) m) m))
        (else
            (remainder (* base (mr-expmod base (- exps 1) m))
                m))))

; ミラー・ラビンテスト、nをチェックする。もし(a^(n-1)) mod n == 1なら、ほぼ間違いなくnは素数
(define (mr-test n)
    (define (try-it a)
        (= 1 (mr-expmod a (- n 1) n)))
    (try-it (+ 1 (random (- n 1)))))

; ミラー・ラビンテストで素数を取得する。nが素数かどうかをチェックする。試行回数timesを試す
(define (mr-prime? n times)
    (cond ((= 0 times) #t)
        ((mr-test n) (mr-prime? n (- times 1)))


        (else #f)))
```

いくつかのデータをまとめ、[1000, 1100]、[10000, 10100]、[100000, 100100]、および[1000000, 1000100]の範囲でこれらの方法がどのように機能するかを見てみました：

```
    従来の方法        従来の方法の約数2の最適化     ミラー・ラビンテスト方法
    [1000, 1100]                 1 ms                       1 ms                            1 ms
    [10000, 10100]               13 ms                      8 ms                            1 ms
    [100000, 100100]             142 ms                     78 ms                           1 ms
    [1000000, 1000100]           1337 ms                    830 ms                          1 ms
```

これは私のIntel Core 1 1.83GHzマシン（Ubuntu 8.04 + DrScheme）で実行されました。

また、Cバージョンを書き、以下の結果が得られました：

範囲 [2, 10000]：

```
    $ ./prime
    Normal:  return 0 cost 10ms
    Fermat: return 0 cost 24ms
    Miller-Rabin: return 0 cost 28ms
```

範囲 [200000, 1000000]：

```
    $ ./prime
    Normal: return 0 cost 3088ms
    Fermat: return 0 cost 824ms
    Miller-Rabin: return 0 cost 959ms
```

gprof分析結果は概ね予想通りであり、レポートの一部を以下に示します：

```
    各サンプルは0.01秒としてカウントされます。
    %   累積時間   自己時間          自己呼び出し   合計
    時間   秒        秒       呼び出し  us/呼び出し  us/呼び出し  名前
    63.55      2.58     2.58   800000     3.23     3.23  is_prime(int)


    15.76      3.22     0.64   800005     0.80     0.84  fermat_expmod(long, long, long)
    10.10      3.63     0.41   800012     0.51     0.88  mr_expmod(long, long, long)
    7.27      3.92     0.29 14813786     0.02     0.02  chk_unusual_square_root(long, long)
    1.23      3.98     0.05                             main
    0.86      4.01     0.04   800000     0.04     0.93  mr_test(long, long)
    0.74      4.04     0.03 14813658     0.00     0.00  square(long)
    0.25      4.05     0.01   800000     0.01     0.94  mr_is_prime(int)
    0.25      4.06     0.01   800000     0.01     0.85  fermat_is_prime(int)
    0.00      4.06     0.00   800000     0.00     0.84  fermat_test(long, long)
```

小さなn値でのパフォーマンスが低い理由は、Cでの通常の再帰処理のオーバーヘッドによるものであると思います。しかし、データサイズが増えると、このオーバーヘッドは比較的無視できるものになります。誰かが上記のアルゴリズムを反復形式（私にはできません :P）や末尾再帰形式（gccは末尾再帰を反復形式に最適化できます）に変換できれば、パフォーマンスはさらに顕著になると予想されます。