中央値関数の排他的論理和

TAOCPには面白い話が満載である. 今回は多数決関数関連の話だ. この本では, 多数決関数のことを中央値(median)関数というので, そう書くことにしよう.

中央値関数は, 奇数個の引数を大きさの順に並べ, その中央の値を返す(平均値でないことに注意).

    x≤y≤zの時 <xyz>=y.

あるグループで流行っている最中限ゲームみたいなものである.

論理関数としては, x,yを偽か真とする時, 常に反対するか常に賛成する偽か真のサクラの入力を加え, <偽xy>がxとyの論理積, <真xy>が論理和である.

実数値については, <-∞xy>はxとyの最小値, <+∞xy>は最大値になる.

こういう話題がいちおう終わってから,

m>0について,
x₀⊕x₁⊕...⊕x_2m=<¬x₀s₁s₁...s_2m>
ここで, s_j=<x₀x_jx_j+1...x_j+m-1¬x_j+m¬x_j+m+1...¬x_j+2m-1>,
ただし, k≥1について, x_2m+k=x_kとする;

と書いてあった.

つまり, 何変数であっても, その排他的論理和がサクラ(定数)なしの2段の中央値関数で得られるのである. 美しい.

そもそも排他的論理和は単調関数ではないから, 中央値関数1段だけでは得られない. 3変数の排他的論理和, つまり全加算器もパラメトロンでは2段必要であった.

上の見事な式は, TAOCPによると, 量子化学が専門の佐々木不可止(ふかし)君が考えたという.

佐々木君は東大物理の学生時代に, 高橋研によく出入りしていて, 多数決論理に関心があったのは知っていたが, こういうことも考えたのかと驚いた. 北海道大学の化学科の教授を定年退職し, 苫小牧駒沢大学に移られたが, そこも退職の後, 惜しいことに2007年に他界された.

あの長々とした式は眼がくらくらするので, TAOCPにあったようにm=1とm=2について書き直すと

x₀⊕x₁⊕x₂=<¬x₀<x₀x₁¬x₂><x₀x₂¬x₁>>

x₀⊕...⊕x₄=<¬x₀<x₀x₁x₂¬x₃¬x₄><x₀x₂x₃¬x₄¬x₁><x₀x₃x₄¬x₁¬x₂><x₀x₄x₁¬x₂¬x₃>>

となる.

この式をテストしてみよう. ます多数決関数mを定義する. nは否定をとる関数.

(define (m . xs)
 (define (n x) (- 1 x))
 (let ((l (length xs)))
  (if (> (apply + xs) (/ (- l 1) 2)) 1 0)))
(m 1 0 0) => 0
(m 1 0 1) => 1
(m (n 1）(n 1) 1) => 0

すると5変数(m=2)の排他的論理和はこう書ける.

(define (ex x0 x1 x2 x3 x4)
 (let ((s1 (m x0 x1 x2 (n x3) (n x4)))
       (s2 (m x0 x2 x3 (n x4) (n x1)))
       (s3 (m x0 x3 x4 (n x1) (n x2)))
       (s4 (m x0 x4 x1 (n x2) (n x3))))
  (m (n x0) s1 s2 s3 s4)))
(map (lambda (x) (apply ex (n2bs x 5))) (a2b 0 32)) =>
(0 1 1 0 1 0 0 1 1 0 0 1 0 1 1 0 1 0 0 1 0 1 1 0 0 1 1 0 1 0 0 1)

でどうやらよさそうだ.

パラメトロンは通常3入力止りだが, 5入力が許されるとしてこの回路はこう描ける.



これをさらに納得すべく, s₁,...s₄や, その計算の元になっているx₁+x₂-x₃-x₄
などを, x₀=0である前半の16個について出力してみよう. 自己相対関数なので, 半分やればいいことになっている.

xの否定は1-xと計算し, (l-1)/2の2と比較したから, 閾値と較べる値は, 本来は
    x₁+x₂+(1-x₃)+(1-x₄)=2+x₁+x₂-x₃-x₄.

従って, 先頭の2を外して,
    x₁+x₂-x₃-x₄
の和sumについてs₁はsum≤0なら偽(x₀=0を除外したから, 等号が要る), sum>0なら真である.

s₁と添字がmだけ違うs₃は, 足す項と引く項がちょうど反対だから, 和も正負逆で同じ値になる. こういう対を相棒といおう. ところで, 真が奇数個だと, 和が0になることはないので, 相棒同士の真理値は反対になり, x₀がキャスティングボートをとって, exは真になる. すなわち排他的論理和である.

ところが, 真が偶数個だと, どこかに相棒同士で和が0になるものがあり, x₀は0だから, sは2つとも偽になる. 先頭の¬x₀が1を主張しても, 和が0でない他の相棒同士は打ち消し合い, この2個(あるいは4個)の0連合には勝てず, exは0になる.

   + s1 + s2 + s3 + s4 ex
0( 0 0  0 0  0 0  0 0  0)
1(-1 0 -1 0  1 1  1 1  1)
2(-1 0  1 1  1 1 -1 0  1)
3(-2 0  0 0  2 1  0 0  0)
4( 1 1  1 1 -1 0 -1 0  1)
5( 0 0  0 0  0 0  0 0  0)
6( 0 0  2 1  0 0 -2 0  0)
7(-1 0  1 1  1 1 -1 0  1)
8( 1 1 -1 0 -1 0  1 1  1)
9( 0 0 -2 0  0 0  2 1  0)
a( 0 0  0 0  0 0  0 0  0)
b(-1 0 -1 0  1 1  1 1  1)
c( 2 1  0 0 -2 0  0 0  0)
d( 1 1 -1 0 -1 0  1 1  1)
e( 1 1  1 1 -1 0 -1 0  1)
f( 0 0  0 0  0 0  0 0  0)

左端は行数である. 各行のかっこ内の9個の値は, 左からx₁+x₂-x₃-x₄, s₁, x₂+x₃-x₄-x₁, s₂, x₃+x₄-x₁-x₂, s₃, x₄+x₁-x₂-x₃, s₄, exである. 偶数個の場合, どの対で和が0になるかを調べたのが下だ. 赤の-は, x₁とx₂の和がx₃とx₄の和に等しい(つまりs₁s₃が0になる)場合, 青の-は, x₂とx₃の和がx₄とx₁の和に等しい(つまりs₂s₄が0になる)場合である.


0(0 0-0-0-0-)
1(0 0 0 0 1 )
2(0 0 0 1 0 )
3(0 0 0-1 1-)
4(0 0 1 0 0 )
5(0 0-1-0-1-)
6(0 0-1 1-0 )
7(0 0 1 1 1 )
8(0 1 0 0 0 )
9(0 1-0 0-1 )
a(0 1-0-1-0-)
b(0 1 0 1 1 )
c(0 1 1-0 0-)
d(0 1 1 0 1 )
e(0 1 1 1 0 )
f(0 1-1-1-1-)


和が0になる対は本当にあるか. こういう図を描いてみた.



いま真というか1の数は偶数である. s₁とs_m+1の範囲で1の数が同じであれば, 一発で和が同じところは見つかったわけでハッピーだ.

一方, s₁がオール0でs_m+1がオール1なら, 左の密度は0, 右の密度は1. そこで図の下へs₂, s₃のように加算の範囲をずらしていけば, どこかで2つの範囲の1の数が一致する, 同じ密度になる場所があり, そこで赤や青の-が引ける対が見つかるはずである(平均値の定理を習ったのは高校生の時だっけ).

こういうわけで, 偶数の場合はexが偽になるのである.

「なるほど」であった.

実はサクラを使うと, どんな関数でも2段の中央値関数で作れるが, そのことはいずれ.

平方根の連分数展開

Norman Lehmerの論文に, factor stencilはNの連分数展開でも使えると書いてあった.

その始めの方. まずx₁=√N=q₁+1/x₂とおく. 但し, q₁はNの平方根の最大整数. 同様にx₂=(√N+P₂)/Q₂=q₂+1/x₃のようにする. (ははぁ. TAOCPの連分数展開のUとVがPとQになっているわけだ.)

順に進んで, x_n=(√N+P_n)/Q_n=q_n+1/x_n+1になるというまでは明快だが, その後に
「従って
P_n=q_n-1Q_n-1-P_n-1,
Q_n=q_n-1(P_n-1-P_n)+Q_n-2
がexpediouslyに得られる」
と書いてある.

上の方の式は11月30日のブログの(4)と(5)の間にあるU_n=A_nV_n-U_n-1と同じだが, 下の方の式はどう導くか.

今回のブログはその計算法である.

(0)はx_nの標準形である. (1)はx_n-1で, そのx_nに(0)を代入し, 前回のように計算を続けると(1)の最後になる. (1)の最初と分母を等しいとしたのが(2)で, ただちに(3)と書ける. 分子の√Nの後に続く項を等しいとしたのが(4)で, 途中で(3)を利用する. P_nとP_n-1を交換したのが上の式(5)だ.

一方, (3)の添字を1減らしたのが(6)で, (3)-(6)が(7),(8)である. 1/Q_n-1は, (5)を移項した(9)から(10)のように得られ, これを(8)に代入すると(11)を経て, 下の方の式(12)になるのが分かる.

漸化式はP₁=0, Q₁=1から始まるが, 問題は, Q_n-2があるので, Q₀が必要なことだ.

実はQ₀=Nなのだが, この計算は次の通り.

これで, 連分数のP,Q,qが次々と求まることになる. なんだかまた久しぶりにこのような計算をした.

素因数探し

11月6日の本ブログにある自転車チェーンの篩を作ったのはHenry Lehmerだ. その父親のNorman LehmerもCalifornia大学Berkeley校の数学の教授であった. 父のLehmerが素因数探しの道具として, Factor Stencilを考案したという話がある.

それがどういうものかが分かってきた. 今回はその話をしよう.

aが素数pの平方剰余であるとは, x²=a mod pとなる何かのxがあることなのは周知のとおり. こういうa(≠0)の時, Legendreの記法で(a/p)=+1とする. (htmlなので, 分子分母が横並びだが, 通常は上下に書く.)

xがないときは, aを平方非剰余といい, (a/p)=-1とする.

p=7の時, 剰余は(0,1,2,4), 非剰余は(3,5,6)である. ここで, 0を除外し, 非剰余同士を掛けると, 3*3=2 mod 7, 3*5=1 mod 7, 5*5=4 mod 7のように剰余になり, 剰余と非剰余を掛けると2*3=6 mod 7のように非剰余のまま, 剰余同士を掛けると4*4=2 mod 7で剰余であり, +1と-1で掛算になる.

ところで, ある整数Nを素因数分解する場合, 任意の整数aについて, a²-N=Rは, Nがpで割れるなら, mod pで考えると, Nの項はないのと同じだから, Rはこのpの平方剰余になる.

aを√Nの程度にとると, Rは小さくなり, 平方剰余の乗算に規則でもっと分解すると, さらに小さい剰余か非剰余が得られ, それから素因数の形が決るという.

たとえば, 平方剰余に-1があるなら, すなわちp-1なら, 素因数は4n+1, 2があるなら, 素因数は8n±1の形, というふうに, 素因数はlinear formになっている. それを組み合せて, 素因数の探索の範囲を狭めたい. ところが, 平方剰余が段々大きくなってくると, linear formの形も複雑になり, 式の形を決めにくくなる. Normanの論文によると, Rが113と199の場合, formの形は11088にもなるそうだ.

そこで, 各平方剰余について, 候補になる素数の表をあらかじめ作っておこうというのが, この発想である.

小さい範囲で, テストしてみよう. 下がその剰余(縦方向)と素数(横方向)の表である. 10月26日のブログの図の素数のところを抜き出したものと思えばよい.

    3 5 7111317192329313741434753596167717379838997
-1    X     X X     X   X X     X   X     X     X X
 2      X     X   X   X   X   X         X X X   X X
 3        X X     X     X     X   X X   X X   X   X
 5        X     X   X X   X       X X   X   X   X
 7  X           X   X X X     X X X           X
11    X X       X       X   X   X           X X X X
13  X         X   X X       X   X   X       X
17          X   X           X X X X   X       X X
19  X X       X       X           X X X X X X

Normanの素数表には1があるというので, この左に1と2の列があるかも知れないが, 除数としては不要なので, とりあえず3からの表だ. 7の下を見ると, 2と11にXがあり, 2と11=4 mod 7が平方剰余であることを示す. 他の素数についても同じだ. たとえば,

(quadres 97)=>(0 1 2 3 4 6 8 9 11 12 16 18 22 24 25 27 31 32
 33 35 36 43 44 47 48 49 50 53 54 61 62 64 65 66 70 72 73 75
 79 81 85 86 88 89 91 93 94 95 96)

のうち, 19までの素数, 2,3,11と, 96(= -1 mod 97)があるから-1とにXがある.

また, 横向きに眺めると, -1の行は, 5,13,17,29,37,41,53,61,73,89,97で, すべて4n+1. 2の行の7,23,31,47,71,79は8n-1, 17,41,73,89,97は8n+1である.

さて, これを使い, N=2279(=43x53) の素因数を探してみる. (sqrt N) => 47.738873049120045 だから, aを48から順に増やしながら, a²-Nを計算し, その素因数分解もしておく.

 (map (lambda (a) (factorize (- (* a a) n))) (a2b 48 80))
=> ((5 5) (61 2) (17 13) (23 7 2) (17 5 5) (53 5 2) (13 7 7)
 (373 2) (857) (97 5 2) (31 7 5) (601 2) (1321) (103 7 2)
 (313 5) (13 13 5 2) (79 23) (139 7 2) (67 31) (17 13 5 2)
 (67 7 5) (73 17 2) (2621) (1381 2) (83 7 5) (61 5 5 2)
 (139 23) (239 7 2) (269 13) (73 5 5 2) (761 5) (283 7 2))

ここに並んでいるのは, 掛けると+1になるもの同士である. (5 5)からは, 5が+1かも, -1かも知れないということが分かるだけだ. (61 2)では, 61が上の表の剰余にないから, 無視. (17 13)は, とりあえず13と17は+1と+1か, -1と-1が分かっただけ. その次の(17 5 5)はラッキーだ. 5と5は打ち消すから, 17が+1と判明した. したがって, (17 13)から, 13も+1と決る.

(13 7 7)からも13が剰余と確信出来る. さらに進むと, (13 13 5 2)や(17 13 5 2)から, 2と5は+1同士か-1同士か, いずれにしろ同じ仲間ということが分かるが, 役にたつかどうかは不明.

しかし, 13と17が得られたので, 13と17の行にXのある素数を探すと, ブラボー! 43と53が見つかった. ついでだが, 43と53の列では, 2と5は空白で, 同じ組なことが示せる. いまは1組だけが見つかったが, 通常は候補がたくさん得られ, 実際に割ってみる必要がある.

Norman Lehmerが1914年に出版したFactor Stencilは, 剰余が±238まで, 素数が48593までの素数表であって, 各剰余ごとに1枚で, Xの代りに素数の位置に孔をあけてある. 今の例では13と17の紙を取り出してきて重ねると, 条件に合う素数のところから光が洩れる仕掛けであった. 当時としては大変な労作であった. いまでは, MacBookで簡単に実験が出来, ありがたい時代だ.

平方根の連分数展開

TAOCPの第4.5.3項の演習問題12は, quadratic irrationalityという話題である. よく分からぬが, どうやら平方数でない整数の平方根の連分数展開をしているらしい.

そのアルゴリズムを, さっそくSchemeに書き直した.

(define (cf d n)
 (let* ((x (sqrt d)) (v 1)
   (a (inexact->exact (floor x))) (u a))
  (define (loop)
   (display (list 'v v 'a a 'u u)) (newline)
   (set! v (/ (- d (* u u)) v))
   (set! a (inexact->exact (floor (/ (+ x u) v))))
   (set! u (- (* a v) u))
   (set! n (- n 1))
   (if (> n 0) (loop) 'ok))
  (loop)))

関数cfは, dの平方根の連分数の部分商をn回計算するものである. 31の平方根は

(cf 31 13)
(v 1 a 5 u 5)
(v 6 a 1 u 1)
(v 5 a 1 u 4)
(v 3 a 3 u 5)
(v 2 a 5 u 5)
(v 3 a 3 u 4)
(v 5 a 1 u 1)
(v 6 a 1 u 5)
(v 1 a 10 u 5)
(v 6 a 1 u 1)
(v 5 a 1 u 4)
(v 3 a 3 u 5)
(v 2 a 5 u 5)

となる. このaの値の5,1,1,3,5,3,1,1,10,1,1,3,5が, √31の連分数の部分商である. √31=5+1/(1+1/(1+1/(3+1/(5+....)))).

アルゴリズムに戻って, 変数はv,a,uの3個(nは繰返し数). let*で設定するvの初期値は1, aは√dの整数部分, つまり最初の部分商, uもaであって,
v_n+1=(d-u_n²)/v_n,
a_n+1=floor((√d+u_n)/v_n+1),
u_n+1=a_n+1v_n+1-u_n
を繰り返す. aに√dがあるが, 整数部分をとるので下の方はいらない. vもuも整数のままで, 誤差が入らないから, いつまでも計算が続く.

この演習問題の次の設問は, ある回数から先は, 0<u_n<√d, 0<v_n<2√dを証明せよというのだ. 私は証明しようとも思わないが, そういうことはこの計算が周期的になるわけで, 平方根は無理数なのに, 連分数にすると繰り返すというのは面白い.

さて, このアルゴリズムは何をやっていたか調べてみる. √dをx₀, 最初の部分商をa₀とすると, (1)のように書ける. (TEXのプログラムでは, D,A,U,Vは大文字になっている.) ただし, u_-1, v₀は(0)のようだ. uの添字がvのより1だけ小さいのに注意. x=(√d+u)/vを各xの標準形とする.



(1)の添字をnだけ増やすと, 一般形x_nが(2)のように得られる. a_nを左辺に移すと, 1/x_n+1が(3)である. 上下を引っくり返せば(4)のx_n+1になるが, 分母に平方根があるので, 有理化し, 左上に√dが来るように, v_nで割る. これを(4)の最後のような標準形だと思うと, 対応する項から(5)の漸化式を得る. u_nを書いたら, 添字を1増やすとu_n+1が出来る. v_n+1も作れるが, u_nの式を入れると, vの漸化式になり, たしかにschemeで書いたプログラムが得られる.

というわけで, 今回もアルゴリズムが解明出来た.

素因数探し

TAOCPの素因数探しの最初のアルゴリズム4.5.4Aは, 2,3,...と素数で次々を割ってみるやつだ.

A1. t←0, k←0, n←N.
A2. if n=1, 終了.
A3. q←floor(n/dk), r←n mod dk.
A4. if r≠0 →A6.
A5. t←t+1, pt←dk, n←q, →A2.
A6. if q>dk, (k←k+1, →A3).
A7. t←t+1, pt←n, 終了.

ここでd₀=2,d₁=3, ....は割ってみるべき, 次々の素数で, d₂=5の次は, 2,4,2,4,...と足していくという方法があるとこのアルゴリズムには書いてある. つまり, d=2,3,5,7,11,13,17,19,23,25,29,31,35,...とする. これをみると, 25,35以外は素数で, 一見うまく行きそうだが, 要するに2と3で割れない数を仮に素数としている. 多くの素数を記憶するわけにはいかぬので, 疑似素数列を生成するわけだ.

上のプログラムをSchemeでコーディングしてみる. p_tは素因数の列の配列なので, Schemeではリストpにする. A5とA7でtを1増やすのは, 配列pの添字を進めるので, リストではいらない. t←t+1, p_t←nは(set! p (cons n p))となる. kは約数の配列の添字だ. A6でkを増やすのは次の約数にするのだから, ここは(nextd)として, 後で考える.

 (define (a2) (if (= n 1) p (a3)))
 (define (a3) (set! q (quotient n d))
  (set! r (modulo n d)) (a4))
 (define (a4) (if (not (= r 0)) (a6) (a5)))
 (define (a5) (set! p (cons d p)) (set! n q) (a2))
 (define (a6) (if (> q d) (begin (nextd) (a3)) (a7)))
 (define (a7) (cons n p))

a4の(not (= r 0))は(> r 0)でよいが, その辺にはこだわらぬ. dは初期値を2とし, 1,2,2,4と順に足せば3,5,7,11が得られる. そこでこのリスト(1 2 2 4)の最後を(...2 _^ 2...)の間の ^ のところに繋げれば無限リスト(1 2 2 4 2 4 2 4 ...)が作れて, その後も2,4,2,4と足せるはずだ. このリストをddといおう. 次のようにして作る.

 (define dd '(1 2 2 4))
 (set-cdr! (list-tail dd 3) (cddr dd))

そうすれば, (nextd)は

 (define (nextd)
  (set! d (+ d (car dd))) (set! dd (cdr dd)))

でよい.
その上と下に

(define (algorithm454a n)
 (let ((p '()) (d 2) (q 0) (r 0))
と
 (a2)))

をつければ完成で, (algorithm454a 3628800)でよべば, (7 5 5 3 3 3 3 2 2 2 2 2 2 2 2)が得られる.

このプログラムでは, 1をfactorizeしようとすると, a2でいきなり止り, 空リストが返る. TAOCPを読み直すと, 「every positive integer n」は p_kを素数, p₁ ≤ p₂ ≤...≤ p_tとして,

n=p1p2...pt

のように一意に表わせる. n=1の時はt=0で成り立つ とあるので, 空リストは当然だ. n=0ならどうか. 0をpositive integerとするかどうかだが, 上のプログラムは止らない.

ところで, 50年も前, パラメトロン計算機でこういう計算をしたころ, 疑似素数列に現れる非素数を, 我々は臨時素数と呼んでいた.

2,4,2,4と足すということは, 6の幅の間に2つの数をテストするから, 疑似素数の全整数に対する頻度は1/3=0.333...である.

2と3と5で割れない疑似素数列なら, その頻度はどうなるか. TAOCPは, 計算時間は20%節約になると説明する.

0から2,3,5のLCM, 30の前までで, 29までで, 2でも3でも5でも割れない数は, 1,7,11,13,17,19,23,29の8個なので, 8/30=0.266... . これが1/3の何%かを計算すると, (8/30)/(1/3)=8/10だから, 計算時間は80%になり, 20%の節約になる.

ついでに, TAOCPに7までにすると, さらに14%節約とあるのを確認しよう.

2*3*5*7=210個のうち, 2でも3でも5でも7でも割れないものは,

(filter (lambda (n) (> (* (modulo n 2) (modulo n 3)
(modulo n 5) (modulo n 7)) 0))  (a2b 0 (* 2 3 5 7))) =>
(1 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79
83 89 97 101 103 107 109 113 121 127 131 137 139 143 149
151 157 163 167 169 173 179 181 187 191 193 197 199 209)

このリストのlengthは48なので, 頻度は48/210. (48/210)/(8/30)=6/7=0.857... これを0.86と思えばTAOCPのいうように節約は14%である.

ところで, 2でも3でも5でも割れない数は, 0から29までに8個あると上に述べた. 計算方法はこうだ. 0から29までの数を書き, 2で割れるもの, 3で割れるもの, 5でわれるものに印をつけると図のようになる.



2で割れる15個の数には赤線を引いた. 3の10個は緑で, 5の6個は青である. これらを30から引くと, 30 - 15 - 10 - 6=-1. 6や10のように2重線のものは, 2回引いたから, 実は引きすぎである. これは1回ずつ戻さなければならない. それは6で割れる0,6,12,18,24の(30/6=)5個, 10で割れる0,10,20の(30/10=)3個, 15で割れる0と15の(30/15=)2個で, -1+5+3+2=9になる. しかし, 0は今度は3回戻されてしまった. ここは本来は1回引く場所であるが, -3+3をやってしまった. そこで(30/30=)1回を引いて, 結局四角で囲んだ8個が残る.

こういう計算法をprinciple of inclusion and exclusionという. 日本語では包除原理とか和積の原理とかいうらしい.

この原理により, 2,3,5,7,11までの頻度を計算してみよう.

(define t (* 2 3 5 7 11))
t=>2310

(- t
(+ (/ t 2) (/ t 3) (/ t 5) (/ t 7) (/ t 11))
(- (+ (/ t 2 3) (/ t 2 5) (/ t 2 7) (/ t 2 11) (/ t 3 5)
      (/ t 3 7) (/ t 3 11) (/ t 5 7) (/ t 5 11) (/ t 7 11)))
(+ (/ t 2 3 5) (/ t 2 3 7) (/ t 2 3 11) (/ t 2 5 7)
   (/ t 2 5 11) (/ t 2 7 11) (/ t 3 5 7) (/ t 3 5 11)
   (/ t 3 7 11) (/ t 5 7 11))
(- (+ (/ t 2 3 5 7) (/ t 2 3 5 11) (/ t 2 3 7 11)
      (/ t 2 5 7 11) (/ t 3 5 7 11)))
(+ (/ t 2 3 5 7 11)))
=> 480

(/ 480 2310)=>16/17=.20779...

大体1/5である.

上に書いたSchemeのプログラムは, dやpの列をリストにした以外はTAOCPのプログラムの焼き直しである. もう少しSchemeらしくしたのが, 次のプログラムだ.

(define (factor n)
 (define dd '(1 2 2 4))
 (set-cdr! (list-tail dd 3) (cddr dd))
 (define (loop n d p)
  (define (nextd)
   (set! d (+ d (car dd))) (set! dd (cdr dd)))
  (define (next)
   (let ((q (quotient n d)))
    (cond ((= (modulo n d) 0) (loop q d (cons d p)))
          ((> q d) (nextd) (next))
          (else (cons n p)))))
  (if (= n 1) p (next)))
 (loop n 2 '()))

(map factor (a2b 1 30))
=> (() (2) (3) (2 2) (5) (3 2) (7) (2 2 2) (3 3) (5 2) (11)
(3 2 2) (13) (7 2) (5 3) (2 2 2 2) (17) (3 3 2) (19) (5 2 2)
(7 3) (11 2) (23) (3 2 2 2) (5 5) (13 2) (3 3 3) (7 2 2) (29))

プログラムは本質的にループだから, loopというプログラムを末尾再帰で呼ぶことにする. ある約数dで割れなければ, 次の約数を(nextd)で準備し, nextを呼ぶ. 素因数分解のたびにddを作り直すのも癪だが, そう長くはないから我慢する. さらに多くの素数の倍数をスキップするddの作り方については, またいつか書こう.

素因数探し

昨日のブログ(篩を使った素因数探し)は, TAOCPのアルゴリズムを理解するのが目的であったので, 同書のアルゴリズムの特徴であるgoto文をそのまま反映していた.

しかし, もっとSchemeらしくするにはどうするか.

前のアルゴリズムでは, 絶えずmoduloを取っているのが問題であった. あのアルゴリズムでは, 篩が2重の配列になっているので, 添字を配列の範囲に収めるためにmoduloを取るのである.

しかし, Scheme風にすると, 配列はリストになり, リストとなれば, 自転車のチェーンのように無限リストが作れる.

すなわち, (define foo '(0 1 2 3)) と設定し(図の上), (set-cdr! (list-tail foo 3) foo)
とすると, (list-tail foo 3)でfooのcdrを3回とり, 3のセルに達する. そのcdrのnilをfooに書き換えると(図の下), 3の次が0になり, 無限ループが出来る.



第2の改良点は, 篩の要素を0と1にせず, #fと#tにする. そうするとandが一発でとれる. TAOCPでは[述語]というIverson blacketを多用していて, これは述語が真のとき1, 偽のとき0になるものなので, 前回のプログラムもそうなっていた. これを(1 1 1... 1)とequal?で真偽を判定した.



lispのandは先頭からみて, 偽をみつけると途端に終了するから, この方が早いのである.

第3は, 無限リストを次々とcdr downするのに, いちいち代入するのではなく, 引数として末尾再帰で渡すことである.

このようにして書直したのが, 次のプログラムである.

(define (algorithm454d n)
 (define (makesieve m n)
  (let* ((x (a2b 0 m))
    (x2 (map (lambda (a) (modulo (* a a) m)) x))
    (s (map (lambda (b)
     (if (member (modulo (- (* b b) n) m) x2) #t #f)) x)))
   (set-cdr! (list-tail s (- (length s) 1)) s)
  s))
 (define s3 (makesieve 3 n))
 (define s5 (makesieve 5 n))
 (define s7 (makesieve 7 n))
 (define s11 (makesieve 11 n))
 (define s13 (makesieve 13 n))
 (define s17 (makesieve 17 n))
 (define s19 (makesieve 19 n))
 (call-with-current-continuation
  (lambda (throw)
   (define (next x s3 s5 s7 s11 s13 s17 s19)
    (if (and (car s3) (car s5) (car s7) (car s11)
         (car s13) (car s17) (car s19))
     (let ((y (sqrt (- (* x x) n))))
      (if (integer? y) (throw (cons (+ x y) (- x y))))))
    (next (+ x 1) (cdr s3) (cdr s5) (cdr s7) (cdr s11)
      (cdr s13) (cdr s17) (cdr s19)))
   (let ((x (inexact->exact (ceiling (sqrt n)))))
    (next x
     (list-tail s3 (modulo x 3))
     (list-tail s5 (modulo x 5))
     (list-tail s7 (modulo x 7))
     (list-tail s11 (modulo x 11))
     (list-tail s13 (modulo x 13))
     (list-tail s17 (modulo x 17))
     (list-tail s19 (modulo x 19)))))))

解が見つかったとき, 脱出するのにcall-with-current-continuationを使っているが, これが結局一番簡単なようである.

引数をぞろぞろ引き摺っていくのは, 素数の個数を変更するのに困るわけだが, とりあえずはこれでさくさく動く.

篩全体をリストにするプログラムも書いてはみたが, mapをとったりするので, 上のプログラムより遅かった. プログラムを動的に生成するという考えもあるが, 分かり難くもなり, 今はためらっている.

素因数探し

前回のブログはTAOCPのAlgorithm4.5.4Cが話題であった.

TAOCPのその次はAlgorithm4.5.4Dで, 篩を使うものである. 今回はその説明をしたい.

以下の例で, 素因数を探す数Nは, 23番のMersenne数M₂₃=2²³-1
=8388607=178481×47である. また, この方法では複数の素数を利用する. それに3,5,7,11を使おう.

TAOCPの説明の通りに進めると, まず下のような表を作る. 一番左が法にする素数mである. 次にその法に現れる数x, つまり0からm-1を書く. 更にその右は, xの自乗のmod mである. つまりmを法として, 自乗の数にはこれしか現れないことを確認する.



最後は0〜m-1のxについて, (x²-N)mod mを書く. この中で, 隣りの自乗の表にあるものだけが, 考慮に値するので, それを赤で示す. 例えばm=3だと, 2,0,0のうち, 自乗の表には0と1しかないので, 2は黒のまま, 0は赤にする. するとその2つの赤なので, xの1と2の場合に対応する.

m=5だと, 赤の位置, x=1か4ならよい.

上の表の最後の数列は,

(define n (- (expt 2 23) 1))
(map (lambda (m)
(let* ((x (a2b 0 m))
 (x2 (map (lambda (x) (modulo (* x x) m)) x))
 (s (map (lambda (x) (modulo (- (* x x) n) m)) x)))
(display (list x x2 s)))) '(3 5 7 11))

のように計算した.

このようにして, あるxについて, xがすべての素数の法で赤の位置に対応したとき, つまり, この表で篩われた時に, x²-Nがあるyの自乗かどうかを調べるのである.

TAOCPのアルゴリズムはごちゃごちゃしているが, わたし流にSchemeで書き直すとこうなる.

(define (algorithm454d n)
 (define (makesieve m n)
  (let* ((x (a2b 0 m)) 
   (x2 (map (lambda (a) (modulo (* a a) m)) x))
   (s (map (lambda (b) 
    (if (member (modulo (- (* b b) n) m) x2) 1 0)) x)))
   s))
 (define ms '(3 5 7 11 13 17 19))
 (define m1 (map (lambda (x) 1) ms)) ;ms length 1's
 (define ss (map (lambda (m) (makesieve m n)) ms))
 (display ss)
 (define x 0) (define ks '())

 (define (d1)
  (set! x (inexact->exact (ceiling (sqrt n)))) (display x)
  (set! ks (map (lambda (m) (modulo x m)) ms))
  (display ks) (d2))

 (define (d2)
  (if (equal? (map (lambda (s k) (list-ref s k)) ss ks) m1)
   (d4) (d3)))

 (define (d3)
  (set! x (+ x 1))
  (set! ks (map (lambda (k m) (modulo (+ k 1) m)) ks ms))
  (d2))

 (define (d4)
  (let* ((d (- (* x x) n)))
    (if (integer? (sqrt d)) (let ((y (sqrt d)))
        (list (+ x y) (- x y)))
        (d3))))
 (d1))

もう少しScheme風にも改良できそうだが, 今はこの辺で. 途中でss, xとksの初期値を出力している. それらは次の通りである.

ss=((0 1 1) (0 1 0 0 1) (1 0 1 0 0 1 0)
(1 1 0 0 1 0 0 1 0 0 1) (0 0 0 1 1 0 1 1 0 1 1 0 0)
(1 0 1 1 1 1 0 0 0 0 0 0 1 1 1 1 0)
(1 0 1 1 1 0 1 0 0 0 0 0 0 1 0 1 1 1 0))
x=2897
kx=(2 2 6 4 11 7 9)

TAOCPのアルゴリズムでは, kを(-x)mod mと計算しているが, 上の表から分かるように, 篩の値は, 0を除いて対称的なので, マイナスにする理由がなく, 私のアルゴリズムでは, x mod mで計算する.

素数をたくさん使うと, d4に来る回数が少なくなるのは当然である. ms=(3 5 7 11 13 17 19)だと523回だが, (3 5 7 11 13 17 19 23 29 31 37 41 43 47)では11回であった.

ところで, 1926年に, Henry Lehmerが自転車のチェーンをつかって篩を作った話は有名である. Mountain ViewのComputer History Museumにはその複製が展示されている.



小さい素数はチェーンが短いので, 何回か繰り返て実装されている. 篩われる数がチェーンの輪の上端に来ると, 設定してあるピンで回路が切れ, 回転が止る. つまり上のアルゴリズムのstep D4へ来る. 結果を調べ, 必要なら再起動するようになっていた.

上の例で, ピンがどうなっているかを示したのが次の図である.



篩のピンが対称的なのがよく分かるではないか.