Line 1: Error: Invalid Blog('by Esehara' )

または私は如何にして心配するのを止めてバグを愛するようになったか

>> Zanmemo

あと何かあれは 「esehara あっと じーめーる」 か @esehara まで

Racket(Scheme)で素因数分解をするようなプログラムを書いてみる

今日の料理

f:id:nisemono_san:20160412170055j:plain

すきやき

概要

ある数が与えられた際に、素因数分解を行なうプログラムを練習がてらに書いてみた。これは数論の本に示唆されていた解法に基づくものである。

今回は、どのようにそのプログラムを実装したかについて、メモをしておく。

はじめに

『はじめての数論』と呼ばれる本が、ピアソン・エデュケーションから出版されている。

この本は数学の一分野である数論の入門書でありながら、ピアソンの本らしく、ところどころにプログラミングとしての課題が乗っている。例えば、今回の「素因数分解する」というプログラムに関してもそうである。この素因数分解するコードを書くという課題について、ちょっと面白そうだったので、Racketで解いてみることにしたのが、今回のエントリである。

どこからどこまで割ってみるか

まず素因数分解といえば、ご存知の通り、{ 18 = 2 * 3^{2} }といったように、なんらかの数を素数の組みあわせに分解することである。このようなプログラムを作るとき、一番単純な方法は、1から求めたい数(以下、これをnとする)まで、順番に割っていけばよい。

とはいえ、直感的には、「n/2以上の数」で割るのは無駄であるように思える。なぜならば、2が最小の素数だからで、この答えにならなければ、以降自身の数になるまで割れないだろうことは、なんとなくわかる。

しかし、もっと範囲を絞ることは可能である。その範囲とは、結論から言ってしまえば、ルートn以下である。

nが素数の場合、そのn自体が答えとなってしまうので、nが何らかの素数の積であるような数であると仮定しよう。このとき、nを割ることができる最低の数をpとし、nをpで割った数をmとする。すると、このとき{ n = pm }という式で表現できるのだが、先程の仮定として、mはpより大きいはずだ(でなければ、最小の数pという仮定に反することになる)。とすると、mとpをかけた数は、pの二乗より大きい筈である。式に直せば{pm \ge p^{2}}ということになる。

ところで、先程{n = pm}とした。つまり、nはpの二乗よりは大きい。そして、式を少し変形させれば、pはnのルートより小さい。

従って、約数として検証する必要のある範囲がわかった。これを関数にする。

(define (range-root x) (range 2 (+ (floor (sqrt x)) 1)))

(check-equal? (range-root 9) '(2 3))

約数のリスト化

さて、範囲がわかったところで、今度は約数をリスト化する必要がある。そこで、約数の候補を与えると、その中から約数だったものだけをリストアップするような関数を考えてみる必要がある。先にテストを提示すると、次のようになる。

(check-equal? (prime-factorization-p 4 '(2) '()) '(2 2))
(check-equal? (prime-factorization-p 10 '(2 3) '()) '(5 2))

ここで考える必要があるのは、「ある数で割りきれたとして、その数でもう一度割りきれる可能性がある」ということと、「リストには候補としてはないものの、割って残った数もまた約数である」ということである。また、「約数が存在しない場合には、その数自体が素数である」ということだ。

恐らく、こういうのは、口で説明するよりも、実際にコードを見せたほうが早いだろう。

(define (prime-factorization-p x y z)
  (if (null? y) (list* x z)
      (let ([y2 (car y)])
        (cond 
          [(= x 1) z]
          [(= 0 (modulo x y2)) (prime-factorization-p (/ x y2) y (list* y2 z))]
          [else (prime-factorization-p x (cdr y) z)]))))

ポイントは、引き数xを「割ったあとの数」として更新していくことと、割り切れないときにのみ、約数のリストから、約数の候補として取り除くこと、そして約数の候補が無くなったときは、約数として、まだ1になっていないxをリストに追加することだろう。

約数をカウントする

とはいえ、この段階では、例えば100などは5 5 2 2といったようなリストになってしまう。普通、素因数分解 を表わす場合、2^2 * 5^2といったように、約数をまとめたかたちに表現したいのではないだろうか。Racketでそのような関数があるか探してみたが、力不足で見つからなかったので、自分でそういった関数を書いてみることにした。

(check-equal? (dubricates-count-p 2 0 '(2 2 5 5) '()) '((2 2) (5 2)))

ポイントは、上記の性質上、同じ約数は連続しているということだ。例えば、5 2 2 5 みたいな変則的な表われ方はしないということである。従って、最初の引き数には数えている数を、そして次には現在の個数をストックするような関数として定義する。

;; dubricate-count-p x c yz
;; 複数あるものをカウントする

(define (dubricates-count-p x c ys zs)
  (if (null? ys) (append zs (list (list x c)))
      (let ([y (car ys)]
            [next-ys (cdr ys)])
        (if (= x y)
            (dubricates-count-p x (+ c 1) next-ys zs)
            (dubricates-count-p (car ys) 1 next-ys (append zs (list (list x c))))))))

;; dubricates-count x
;; 上記関数のラップ版
(define (dubricates-count xs)
  (dubricates-count-p (car xs) 1 (cdr xs) '()))
(check-equal? (dubricates-count '(3)) '((3 1)))

素数のリストを作成する

さて、素因数分解をするにあたって、先ほどの約数候補リストは冗長であるという欠点が存在している。例えば、100の場合、約数候補として4が入っているのだが、これは既に2の段階で分解されているので、わざわざこれを約数候補とする必要はない。従って、あらかじめ、素数のリストを生成し、それを使うほうが望ましい。

;; filter-prime-p xs ys
;; 素数だけのリストを生成する
(define (filter-prime-p xs ys)
  (if (null? xs) (reverse ys)
      (let ([x (car xs)]
            [next-xs (cdr xs)])
        (filter-prime-p (filter 
                         (lambda (y) (not (= (modulo y x) 0))) next-xs)
                        (list* x ys)))))

(check-equal? (filter-prime-p (range 2 10) '()) '(2 3 5 7))

;; filter-prime xs
(define (filter-prime xs)
  (filter-prime-p xs '()))

素数のリスト生成に関しては、よく言及されるエラトステネスのふるいのようなアプローチを利用している。

あるリストが渡されたとき、最初の数は素数だとして、その倍数をリストから除去し、約数のリストに追加する。そして、残ったリストの最初の数も、それ以前には割りきれなかった数なので、素数であるわけだ。それの倍数もまたリストから除去するということをくりかえす。これをわたされたリストが無くなるまでくりかえす。

おわりに

f:id:nisemono_san:20160413083350p:plain

完成品は、こちらのgistに掲載している

素因数分解するプログラムは、考えてしまえば一発なのだが、再帰などの練習をする上において楽しかった。また、素因数分解は、オイラー関数(ある任意の数と素であるような数を求める関数)などにも利用できる。また、なにかの言語を取得するさいにも、練習をするにおいて十分な大きさのプログラムであることがわかった。

はじめての数論 原著第3版 発見と証明の大航海‐ピタゴラスの定理から楕円曲線まで

はじめての数論 原著第3版 発見と証明の大航海‐ピタゴラスの定理から楕円曲線まで