Common Lispで遅延評価を作って遊ぶ(3-2) - 素数列で遊ぶ (数列で遊ぶ - その2)

今回の目的

ここまでに作ったcl-lazyライブラリを使って、こまけーことはいいんだよ、の精神で色々な(数)列を作って遊んでみたいと思います、第二回。数といえば自然数自然数といえば素数、ということで素数列を作ってみます。

前回と同じく、冗長なのでREPLによる評価値の出力は(基本的に)表示例から省略します。

前置き:数列の作り方

念のため作成したcl-lazyライブラリの使い方。githubの以下のリポジトリから、quicklispから見える所定の位置にcloneした上で下記のようにすると数列を作成できます。各種例は前回を参照。

eshamster/cl-lazy · GitHub

CL-USER> (ql:quickload :cl-lazy)
CL-USER> (use-package :cl-lazy)
CL-USER> (enable-series-processing-syntax)

; フィボナッチ数列の例
CL-USER> #<a[n] = 0, 1, (+ a[n-1] a[n-2])>

素数列を作る

解1:割と普通の作り方

ある自然数nが素数かを判定プログラムといえば、まず自然数2との大小比較をし、2より大きい場合は、2以上かつnの平方根以下の自然数で割り切れないかをループで調べる、というものになると思います。まずはこの判定式を割りと素直に実装してみます。

ループ部分が若干ひねってありますが、基本の考え方は同じです。aは数列でnはインデックスです。do-seriesは見ての通りだと思いますが、(do-series (x y to z) &body body)とすると、数列yの値をxへ順に捕捉してbodyを実行する、ということをz+1回行います。前回のように(dotimes (i n) (print (lnth i a)))とアクセスすると毎回lnthでリストをたどるため計算量はO(\sum^{n}) = O(n^{2})になりますが、do-seriesはO(n)で計算できるように作っています*1

さて、前回時点である数列から、条件にあった要素をフィルタリングして新しい数列を作るfilter-series関数を作成しました。これを使って、上記のcheck-series-modをフィルタリング関数として自然数列から素数列を取り出してみます。

CL-USER> (defparameter *primes* (filter-series-using-little 
                                 #'check-series-mod
                                 #<a[n] = (1+ n)>))
CL-USER> (do-series (val *primes* to 19) (format t "~D " val))
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71

この書き方だと、「2以上かつnの平方根以下の自然数」で割るよりも効率が良いです。「nの平方根以下の『素数』」のみで割り切れるかをチェックするためです。これは、check-series-mod関数に渡す数列を自然数列ではなく自分自身、すなわち素数列とすることで実現しています*2

なお、「自然数列から取り出す」という意図に忠実に#<a[n] = (1+ n)>と書きましたが、当然#<a[n] = 2, (1+ (* n 2))>として2, 3, 5, 7, 9,...という数列を対象にした方が効率が良いです。逆に短く書くのであれば#<a[n] = n>で十分です。

解2:エラトステネスの篩の素直な拡張

素数の求め方といえば、エラトステネスの篩ですね。まず2以上の自然数をnまで並べる、次に2で割り切れた(2より大きい)数を消す(=ふるう)、残った中で次に小さい3でさらに残りをふるう、さらに残った中で次に小さい5でふるう、、、と割る側の数がnの平方根以上になるまで繰り返すとn以下の素数がすべて求まる、という手法でした。

上記の解1もよくよく考えるとエラトステネスの篩と同等の計算量でできている、つまり割り切れるかのチェックを素数のみで済ましている(これがエラトステネスの篩の計算量に相当するかは未計算)のですが、もっと篩っている感じを出したいと思いました。そんなわけでエラトステネスの篩のn→∞における自然な拡張として素数列を求めてみたいと思います。先に言いますが、コードがエレガント?で短くなる(あと面白い)代わりに実行効率はだいぶ落ちます。

二次元数列*sieve*では1行目に2から始まる自然数の列を設置し、以降の行ではひとつ上の行から先頭の数で割り切れない値を残して数列を作ります。こうすることで、下記のように2行目に2でふるった数列、3行目にさらに3でふるった数列、4行目にさらに5で振るった数列ができ上がります。あとは先頭の数を取り出していくと素数列ができ上がるというわけです。

CL-USER> (do-series (line *sieve* to 10)
             (do-series (val line to 10)
                   (format t "~4D " val))
               (fresh-line))
   2    3    4    5    6    7    8    9   10   11   12
   3    5    7    9   11   13   15   17   19   21   23
   5    7   11   13   17   19   23   25   29   31   35
   7   11   13   17   19   23   29   31   37   41   43
  11   13   17   19   23   29   31   37   41   43   47
  13   17   19   23   29   31   37   41   43   47   53
  17   19   23   29   31   37   41   43   47   53   59
  19   23   29   31   37   41   43   47   53   59   61
  23   29   31   37   41   43   47   53   59   61   67
  29   31   37   41   43   47   53   59   61   67   71
  31   37   41   43   47   53   59   61   67   71   73

CL-USER> (defparameter *primes* #<a[n] = *sieve*[n][0]>)
CL-USER> (do-series (val *primes* to 19) (format t "~D " val))
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71

だんだんと篩にかける感じが出たと思うのですがどうでしょうか。ただ、先頭の数の2乗まで素数が確定してしまうため、「お、消えてる消えてる」と感じるのはせいぜい最初の2,3行ぐらいで少しがっかりですが…。なお、見ての通り同じ値を大量複製するのでメモリ効率は恐ろしく悪いです。

気持ち改善したバージョンも。初期値数列を2と3以上の奇数の数列にしました。また、毎回割る数を求めるためにa[n-1][0]をたどるコスト(lcdrをn-2回とlcdrを2回)もバカにならないので、2回目以降クロージャに捕捉してO(1)でアクセスできるようにしました。元が効率悪すぎるのでこれでもだいぶ良くはなりますが、たかは知れていますね…。とはいえ、これ以上やるとなんの実装だったか分からなくなりそうなのでこの辺にしておきます。

素数列で遊ぶ

素数といえば素因数分解ということで、素数列を使った素因数分解の実装と、分解ついでに最小公倍数や最大公約数を簡単に求めてみたいと思います。書き始めたら結構面白かったのでライブラリ化しました。cl-lazyと同じくquicklispによるloadに対応しているので、所定の位置にcloneして(ql:quickload :cl-prime-number)で使えます(当然ながらcl-lazyへの依存あり)。なお、この中で使っている素数列*prime-series*は上記の解1改良版(フィルタ対象数列を#<b[k] = 2, (1+ (* k 2))>としたもの)で定義しています。

eshamster/cl-prime-number · GitHub

素因数分解

素因数分解の逆変換

素因数分解を行うfactorize-in-primeとその逆変換を行うinverse-factorize-in-primeです*3。car, cdrとlcar, lcdrが混じっているのが少し気持ち悪いですが、素数列のおかげでリストの再帰操作として統一的に書けている…と思います。factorize-in-primeは下記のように動作します。変換後は主題の遅延評価はなんの関係もなくなっていますが、せっかくなのでもう少し進みます。なお、inverse-…の方は文字通り逆変換なので動作例は省略します。

; ※いちいちprintするのが面倒なのでREPLの出力を表示しています。
CL-USER> (factorize-in-prime 120)
(3 1 1)
CL-USER> (factorize-in-prime 70)
(1 0 1 1)
CL-USER> (factorize-in-prime 50)
(1 0 2)

素数の次数をリスト化した形に分解しています。素数指数表現という呼び名があるようです。これができると、例えば下記のように、最大公約数や最小公倍数を求める問題が単なるリストの操作で解決できるようになります*4

[GCD = 最大公約数を求める]

120                (3 1 1)
 70 --factorize--> (1 0 1 1) --min--> (1 0 1) --inverse--> 10
 50                (1 0 2)

イメージ:(calc-GCD 120 70 50) -> 10
[LCM = 最小公倍数を求める]

120                (3 1 1)
 70 --factorize--> (1 0 1 1) --max--> (3 1 2 1) --inverse--> 4200
 50                (1 0 2)

イメージ:(calc-LCM 120 70 50) -> 4200

実際に、足りない部分であるmin, maxをするためのbundle-lists関数(単なるリスト操作なので実装略)を作って以下のようにcalc-GCD, calc-LCMを実装してみました。上記のイメージ通りに動作します。

整数の世界で計算しているとここまで両計算を共通化するのは難しいと思います。他の活用例では、整数の掛け算や割り算(それぞれ素数指数表現で足し算、引き算)や互いに素の判定というあたりがすぐに思いつくところです。

ひとしきりcl-prime-numberを書いてから「今回逐次評価なリストの世界に持ってきてしまったけど、遅延リストの世界で完結させた方が良かったのでは…」と思い始めました。そうすると、例えば最大公約数を求めるときに全部の数を最後まで素因数分解をする必要がなくなるわけです。実際その方向で書き直し始めたのですが、cl-lazy自体の設計を詰めないときれいに書ききれないのではと思うところがあったので、いったん手を止めて非遅延評価なリスト版で記事にしました。

ちなみに、今回学んだこと:何も考えずに引数を&restを多用すると、&keyと食い合せが悪かったりapplyだらけになったりで嬉しくない

番外編:無限次元列

n次元配列が作れるなら、次元が無限の配列も作れないかと考えてみます。

まず、n次元配列についてあらためて考えてみます。n次元配列aに対して、#{a[x]}のように一次元分取り出すとn-1次元配列が得られます。同様に#{a[x][y]}と2回繰り返すとn-2次元が、n回繰り返すと0次元=値が取り出せます。この延長で考えると、無限次元配列では#{a[x]}とすると無限次元配列が取り出され、#{a[x][y]...}と繰り返すとやっぱり無限次元配列が取り出されることになりそうです。どこまで行っても値を取り出せないので、始める前から役に立たなさそうな雰囲気が漂っていますが実装を考えてみます。

さて、無限長数列が実現できたのは、lcdrを何度実行してもリスト(lcons "値" "再起呼び出し")が返されていたためでした。一方でlcarを実行することで次元が一つ減って値を取り出せていました。ということは、lconsで「値」を入れていた一つ目の要素も「再起呼び出し」に置き換えると何次元目までいっても次の次元を取り出せそうです。

ということで、最小の無限次元列作成関数は(defun f () (lcons (f) (f)))と書けそうです。なお、#<>リーダマクロはlambdaで定義した関数自身を参照する能力を持たないため利用できません*5。このfにもう少し装飾をつけたmake-inf-dimmensionを作って実験してみます。

; ※ここではREPLの評価結果を省略しない
; 無限次元列を作る。routeはそこにたどり着くための経路をリスト化したもの(単なる表示用)
CL-USER> (defun make-inf-dimension (route)
           (format t "~A " (reverse route))
           (lcons (make-inf-dimension (cons 0 route))
                  (make-inf-dimension (cons (1+ (car route)) (cdr route)))))
MAKE-INF-DIMENSION

CL-USER> (defparameter x (make-inf-dimension '(0)))
(0)
X
CL-USER> #{x[2][1][3]}
(0 0) (1) (1 0) (2) (2 0) (3) (2 0 0) (2 1) (2 1 0) (2 2) (2 1 0 0) (2 1 1) (2 1 1 0) (2 1 2) (2 1 2 0) (2 1 3) (2 1 3 0) (2 1 4)
#<COMPILED-LEXICAL-CLOSURE (:INTERNAL MAKE-INF-DIMENSION) #x30200179756F>

; lazy evaluationになっていること(二度目(以降)は評価されずに同じ値が返ること)を確認
CL-USER> #{x[2][1][3]}
#<COMPILED-LEXICAL-CLOSURE (:INTERNAL MAKE-INF-DIMENSION) #x30200179756F>

評価ついでにformatでそこに至る経路を表示しています。例えば、(2 0)であれば、#{x[2][0]}で取り出せる位置にあるということです。なんだか思った以上に出力が多い気がするのでxを再定義して一つ一つ見てみます(下)。こう見ると、次元方向と要素方向にそれぞれ1つずつ余分に評価されていることが分かります。lconsの作りの問題です。

; 思いの外に見づらかったのでREPLの出力はまた省略
CL-USER> (defparameter x (make-inf-dimension '(0)))
(0)
CL-USER> #{x[0]}
(0 0) (1)
CL-USER> #{x[1]}
(1 0) (2)
CL-USER> #{x[2]}
(2 0) (3)
CL-USER> #{x[1][0]}
(1 0 0) (1 1)

無限次元列のまともな定義が分からないので、これでできたと言えるのかもよく分からないですね…。一応、下図のように二分木にマッピング*6して一つ目の「ずれたフィボナッチ数列」(2次元数列の例)と比べてみると、対称性がとれているので大丈夫な気がします。つまり、左に行くと同じ次元の次の要素、右に行くと一つ次元を下りるといった具合です。

↓「一つずれたフィボナッチ数列」(二次元列)を二分木にマッピングf:id:eshamster:20150809200942p:plain:h300

↓無限次元列を二分木にマッピングf:id:eshamster:20150809200941p:plain:h300

ちなみに、lazyマクロを多値返却(multiple-value-bind)に対応させた上で、make-inf-dimension関数の返り値をvaluesを使って多値返却にすることで、一応任意の場所から値を取り出せる無限次元配列が作れそうです。

次回へ

おかしい、なんでこんなに記事が長く…。

次回から実装した数列定義用のリードマクロの話です。実装中に考えていたことのメモ、という記事の趣旨を考えると、こうやって遊んでいる間にそろそろ忘れていそうで怖いです…。

シリーズリンク

*1:実は前回時点で作ってあったのですが、リードマクロの#{x[i][j]}な書き方が割りと気に入っていたのでそちらの書き方で通しました。

*2:自分自身の数列のより若い要素を参照する、という意味で-using-littleという名前にしました。どうもいまいち感のあるネーミングですが…。

*3:両者とも見ていて気になったので記事執筆時点のgithub版とは少し実装が違います(動作確認済み)。factorizeの方はloopマクロの集約機能を忘れていたので使うように直しました。こちらは純粋に改善です。inverseの方は末尾再帰形式で長かったので非末尾再帰再帰で書き直しました。びっくりするぐらい短くなりますね…。当然こちらはトレードオフありです。

*4:素因数分解と逆変換の計算コストを考慮すると総合的には速くないようです。ただ、素数指数表現の世界で色々なことを効率よく実現できるのであれば、変換を前後に1回だけにして極力その世界で計算することで総合的な効率も良くなります。ただ、その仮定の部分が正しいかは数学的素養がないので分かりませんが…。

*5:alambdaを使えば書けそうですが、これだけのために使うのも微妙かなと。仮に対応していれば#<a[n] = (self)>と書けるはずです。もはや何を定義しているのかよく分かりませんが。

*6:無限次元列を二分木にマッピングできることに気づいてすごく頭がスッキリしましたが、一方でもっと摩訶不思議な世界を想像していたので、どうにもガッカリ感を拭えなかったりします。consセルで書いている以上当たり前なんですけどね。