SICP § 1.2.6 素数性のテスト

いや〜、わからねぇwww
まだこのセクション理解しきれてないので、以下、途中経過だけ書いておきます。

■除数の探索
覚えていますか素因数分解。「素数^n」の表記を使って自然数を分解していくアレです。
割ることができる素数が見つからない場合、その数は素数って言えるとか言えないとか。（どっちだよ。）
さて、この素因数分解してるときどうやって「これ以上割りきれる数はない」と判断していたか覚えてます？
例えば 131 という数字。

 2・・・ 2 × 65 +  1　ダメ。

 3・・・ 3 × 43 +  2　ダメ。

 5・・・ 5 × 26 +  1　ダメ。

 7・・・ 7 × 18 +  5　ダメ。

11・・・11 × 11 + 10　ダメ。

13・・・13 × 10 +  1　ダメ。

13でそれ以上計算するのやめますよね。
なんでここでやめることができるかというと、これ以上大きな値で探査してもそのペアになる値はそれまで探査してきた小さい値、2 〜 11 の範囲の数の再評価になってしまい、意味がないから。
だから、探査する値の２乗が、素数かどうか調べている数よりも大きくなった時点で打ちきることができるわけです。

11^2 = 121 < 131

13^2 = 169 > 131

それを踏まえてこの条件を終了条件とすると、素数を求める処理は次のように記述できます。

;呼び出すだけ
(define (smallest-divisor n)
  (find-divisor n 2))

;除数を求める処理
(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n) ;これが上の終了条件。
        ((divides? test-divisor n) test-divisor)
        (else (find-divisor n (+ test-divisor 1)))))

;除数かどうかの判別処理
(define (divides? a b)
  (= (remainder b a) 0))

;素数かどうか判別する処理
(define (prime? n)
  (= n (smallest-divisor n)))

これはものすご〜く素直な素数チェックのロジックなので得に問題なく理解できました。
　
　
　
■Fermatテスト
このセクションは読んだだけではさっぱり理解できないので、写経しながら自分の理解度を記述しながら進めてみます。

『Fermatの小定理』
n を素数、a を n より小さい任意の正の整数とすると、a の n 乗は n を法として a と合同である。

いや〜、読むだけだとさっぱりわからん。一つ一つ順番に見ていきます。
まず登場人物の関係を表してみます。

実際の定理の文言を整理すると、こんな感じ？

n: 素数

a: a < n な自然数

a^n ≡ a (※nを法とする）

で、「法」ってなによ？という部分が、その次のカッコ内の記述にある感じ？

（二つの数は、両者を n で割った時、同じ剰余を持てば、n を法として合同という。

　数 a を n で割った剰余を、「n を法とした a の剰余」あるいは単に、

　「n を法とした a」（a modulo n)と言う。）

この部分はFermatの小定理とは無関係で、あくまで「○○を法とした△△」という表現についての説明だと思います。
だから、文中で出てくる n は、「Fermatの小定理の n」とは全然無関係な数だと考えるべきでしょう。
「○○を法とした△△」という表現を自分の脳内で道具として利用できるようにするには、この表現を見たときにすぐにイメージができるようにしておく必要があります。
なので、実はこの部分を理解するのが最優先事項だと思います。俺的には。
というわけで噛み砕いてみますね。

二つの数をそれぞれ a, b、この２つの数を割る値を n、それぞれの剰余を m_a, m_b、そして、（＋とか−とか×とか÷と同じように）剰余を求める数学的演算子として、modulo とすれば、次の関係式が成り立ちます。

a modulo n = m_a

b modulo n = m_b

で、もし、

m_a = m_b

ならば、

a と b は、n を法として合同

と表現すると。ようするに、

「a modulo n = b modulo n」ならば、「a と b は n を法として合同」と言える。

ということですね。 ここまで来るともう一つの重要な表現である、『数 a を n で割った剰余を、「n を法とした a」』 っていうのは、もはやなんてことなくて、単にその数の剰余を求めてますよ、ということになりますね。
今後、「n を法とした a」という文言を発見したら、脳内で「a modulo n の結果のこと」に即変換してやればいいわけです。
その上で、改めてFermatの小定理をもう一度みてみます。

n: 素数

a: a < n な自然数

a^n ≡ a (※nを法とする）

１行目、２行目はまぁいいとして、３行目！ これは次のように書き換えられますね。

n: 素数

a: a < n な自然数

a^n modulo n = a modulo n

やっと理解できた。。では続きを写経。
あ、その前に、「俺が」混乱しないように今やろうとしていることを整理しておきます。
今議論している内容は、

「Fermatの小定理を証明すること」ではなく、

「Fermatの小定理を利用して素数を求めるにはどうするか」

です。あぶねー。勘違いして時間を無駄にするところだったぜ。。

n が素数でない時、一般に a < n なる殆どの数について上の関係は成り立たない。

このことから次の素数性テストのアルゴリズムができる。
与えられた数 n に対し、 a < n なるランダムな数の n を法とした a^n の剰余を計算する。

結果が a でなければ、n は確実に素数ではない。

結果が a ならば、n が素数である確率は高い。

別の数 a をランダムにとり、同じようにテストする。

また等式が成り立つなら、n が素数であることに更に確信が持てる。

更にいろいろな a の値で確かめることで結果の確信度が増す。
このアルゴリズムはFermatテストとして知られている。

ここで話題にしてるのは、Fermatの小定理を使った、「Fermatテスト」なるものについて。 「法」の概念がわかりゃ、上の分かりにくい文章もすっきり整理できるっしょ。

主役： n

脇役： a (< n)

脚本： a^n modulo n

結末： m (= a^n modulo n)

       ? m ≠ a　→ n は素数ではない。

       ? m ＝ a　→ n は「素数である確率が高い」。

          この場合、ランダムに a を選んで同じ評価を繰り替えしていけば、

　　　　　より素数である確率が高くなる。

ここまでくれば、次のschemeの実装も理解できるかも？

Fermatテストを実装するには、ある数のべき乗の別の数を法とした剰余を計算する手続きがいる。

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

「ある数のべき乗の別の数を法とした剰余を計算する手続き」という表現に、文字変数も一緒に記述して欲しいなぁ。。
この定義式のアルゴリズムが理解できないが、定義式中の登場人物、base, exp, m の役割を推測してみると、こんな感じでしょうか。

base ： ある数（べき乗される数）？

exp  ： べき乗（する数）？

m    ： 別の数？

さて、この定義式のアルゴリズムについてですが・・・悲しくなる程理解できません。。
普通に考えたら再帰して「べき乗演算」を実施してからその結果に対して剰余算するのが素直なやり方だと思うんですが、この定義式は、自身を呼び出して、remainderまで含めて再帰している。。

・・・根本的に何かを見落としていそうだなぁ。 ちょっと数式で考えて行ってみましょう。

　

  base^exp modulo m

= (base * base * ... * base) modulo m  (※ (base * base * ... * base)は exp 回繰り返し）

modulo演算は四則演算と同じような位置づけだが、各種演算子に対してどういう法則（分配法則とか結合法則とか）が成り立つかどうかがミソっぽいですね。
ようするに、

= base * ((base * base * ... * base) modulo m)      (※ (base * base * ... * base)は exp-1 回繰り返し)

= base * base * ((base * base * ... * base) modulo m)  (※ (base * base * ... * base)は exp-2 回繰り返し)

　　・

　　・

　　・

= base * base * ... (base modulo m)              (※ base * base * ... は exp-1 回繰り返し)

のように変換ができると何となく再帰で処理できそうな・・・っていやいやいや、これは違うなwww
答えの数字がどんどんデカくなってっちゃうもの。

仕方ない、自分の脳みそで予想できないなら、定義式のアルゴリズムから何をしているかを解析してみましょう。
元の定義式の逐次平方を削除して簡略化して、正規順序評価してみます。
簡略化するとこんな感じ。

(define (expmod base exp m)
  (if (= exp 0)
      1
      (remainder (* base (expmod base (- exp 1) m))
                 m)))

で、例としてこんなパラメータでやってみます。

base = 7

exp  = 5

m    = 3

if構文は逐次評価してやると、こんなツリー構造になる。

<pre>(remainder
 (* 7 
    (remainder
     (* 7 
        (remainder
         (* 7 
            (remainder
             (* 7
                (remainder
                 (* 7 
                    1)
                 3))
             3))
         3))
     3))
 3)

こいつを数式にするとこんな風に記述できます。

((7 * ((7 * ((7 * ((7 * ((* 7 1) modulo 3)) modulo 3)) modulo 3)) modulo 3)) modulo 3)

え！？modulo演算てこんな風に展開されるの？？マジで？

シラネーヨw

ギブギブ！！これはもうそういうモンだと解釈して話を進めようw

と思ってたらこんなんありました。

冪剰余（Wikipedia）

この中の「途中で剰余を取る」って項目。
これ、証明できるんだろうなぁ。どうしようか。。