平方根の連分数展開

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へ来る. 結果を調べ, 必要なら再起動するようになっていた.

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



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

素因数探し

大きい数が素数かどうか知りたいことがある. 素数と分かればそれでよし. 素数でなければ, 素因数が知りたくなるのが人情だ. 素数であるかは素数性のテストがあるので, それによればよい. 素因数探しはそれに較べ困難である.

Knuth先生のTAOCP第2巻の4.5.4項は素因数に分解する話題である. 最初のアルゴリズム4.5.4Aは2,3,5,...と順に割ってみる方法である. 素数を全部覚えているわけにもいかぬから, 2,3,5のあとは4,2,4,2,...と足して疑似素数を発生する. この辺はもう少し凝れるが, とりあえずはここまで.

その後にあるアルゴリズム4.5.4Cはちょっと面白いので, 今回はその話にしよう. TAOCPによると, この方法は1643年にFermatが使ったものらしい.

分解したい奇数nが与えられた時, このアルゴリズムは下の図の左上のように, 正の整数aとbについて, a²の正方形からb²の正方形を引いた面積をnにするのである. 灰色の部分がnになる.

そういうaとbは, 右上の図のように, aとbの差が1の時, からなず存在する. その隙間の狭い面積をnにするのである. nは奇数だから出来るわけだ. つまりa=(n+1)/2, b=(n-1)/2とすると, a²-b²=(n²+2n+1)/2-(n²-2n+1)/2=4n/4=nである.

このようなaとbが分かれば, n=a²-b²=(a+b)(a-b)だから, a+bとa-bが素因数である. 2つの素因数の差は2b.



素因数がいくつもあると, aとbの組はいくつもあったりする. n=5005なら下の2つの図のように, 71²-6²も73²-18²も5005になる.

求め方はこうだ. 最初a=floor(√n), b=0とする. つまりa²がnに等しいか, ぎりぎりに近いが少し小さ目にする. そしてr=a²-b²を計算し, r<nならaを1増やす. r>nならbを1増やす. r=nならそのaとbが求めるものだ. 素因数を求めるのに割り算をしていない.

たとえばn=9なら, a=3, b=0で決まる. 素因数は3.

n=15ならa=3, b=0, r=9から始め, r<nだからaを4にする. するとr=16になり, r>nだから今度はbを1にする. するとr=15iになり, a=4, b=1に決まる. 素因数は5と3.

n=21ならa=4, b=0, r=16から始める. そろそろ面倒になってきたから, プログラムを書こう.

(define (fermat n)
 (define (try a b)
  (let ((r (- (* a a) (* b b))))
   (display (list a b r)) (newline) ;a,b,rを出力
   (cond ((< r n) (try (+ a 1) b))
         ((= r n) (cons (+ a b) (- a b))) ;素因数が決まる
         ((> r n) (try a (+ b 1))))))
 (try (inexact->exact (floor (sqrt n))) 0))

実行してみると

(fermat 21)
(4 0 16)
(5 0 25)
(5 1 24)
(5 2 21)
=>  (7 . 3)

上の図の下の左は

(fermat 5005)
(70 0 4900)
(71 0 5041)
(71 1 5040)
(71 2 5037)
(71 3 5032)
(71 4 5025)
(71 5 5016)
(71 6 5005)
=> (77 . 65)

これで分かるように, このアルゴリズムはaとbの差の大きい方の解を得る.

TAOCPのアルゴリズム4.5.4Cでは, rの計算を加減算だけで出来るように, うえのaとbの代りにx=2a+1, y=2b+1を使い, rもnと比較するのでなく,　r-nにして0と比較する.

こういうアルゴリズムだ.

C1 x←2(floor(√n)), y←1, r←floor(√n)²-n.
C2 if r=0,終了 n=((x-y)/2)((x+y-2)/2).
C3 r←r+x, x←x+2.
C4 r←r-y, y←y+2.
C5 if r>0 →C4, else →C2.

(TAOCP風の記述では, C2, C3のような各ステップの終わりに, 行き先(→)の指定がなければ, 次のステップへ進むことが了解されている.)

このプログラムの意外なのは, C2でr=0でなければxとyを増やしてしまう点だ.　xを増やした結果はr=0にならなず, r>0になるから, yも同時に増やしている.　n=19(素数)でトレースしてみる. 赤字のrはnより小さく, aを増やす時を示す.

(fermat 19)
(4 0 16)
(5 0 25)
(5 1 24)
(5 2 21)
(5 3 16)
(6 3 27)
(6 4 20)
(6 5 11)
(7 5 24)
(7 6 13)
(8 6 28)
(8 7 15)
(9 7 32)
(9 8 17)
(10 8 36)
(10 9 19)
=> (19 . 1)

nが平方数なら, n=9の例のように一発で決る. そうでないなら, aは√nより小さいからr<nになり, aを増やす. その後bを増やすとrは減って, nに等しくなるか(つまり停止するか), r<nになりaを増やす. 先ほどはr>nだったrからb²を引いてr<nになったところへ, bより大きいa²を足すのだから, b²を引くまえのrより大きくなり, r=nとなるはずはないのである. (上の結果の赤字の上下のrの値を見較べると, 下の方が大きいのが分かる.)

前のプログラムも

(define (fermat n)
 (define (try a b)
  (let ((r (- (* a a) (* b b))))
   (display (list a b r)) (newline) ;a,b,rを出力
   (cond ((< r n) (try (+ a 1) (+ b 1)))
         ((= r n) (cons (+ a b) (- a b))) ;素因数が決まる
         ((> r n) (try a (+ b 1))))))
 (try (inexact->exact (floor (sqrt n))) 0))

と改良できて,

(fermat 19)
(4 0 16)
(5 1 24)
(5 2 21)
(5 3 16)
(6 4 20)
(6 5 11)
(7 6 13)
(8 7 15)
(9 8 17)
(10 9 19)
=> (19 . 1)

たしかにこの方がスマートだ.

平方剰余

MathworldのQuadratic Residueを見ると, 10月1日のブログUlam Spiralのような, なんやらフラクタル風の図がある. これもやはり自分でも描くべしと, 調べてみた.

Quadratic Residue, つまり平方剰余に興味を持つ人には時々出会う.

まずaがpの平方剰余であるとは, 0<x<pのxについて, x²=a (mod p)となるxがあることだ. 従って, この範囲のについて計算すると, たとえばp=10として

(define p 10)
(map (lambda (x) (modulo (* x x) p)) (a2b 1 p))
=>(1 4 9 6 5 6 9 4 1)

従って10に対しては, 1,4,5,6,9が平方剰余であり, ここに現れない2,3,7,8が平方非剰余(Quadratic Nonresidue)である. 上の計算で見る通り, 1からp-1のxに対して, 現れる平方剰余は対称なので, 真ん中まで計算すれば良い.

従って

(define (quadratic-residue p)
 (map (lambda (x) (modulo (* x x) p))
  (a2b 1 (+ (quotient p 2) 1))))

(quadratic-residue 10) => (1 4 9 6 5)
(quadratic-residue 11) => (1 4 9 5 3)
(quadratic-residue 12) => (1 4 9 4 1 0)
(quadratic-residue 13) => (1 4 9 3 12 10)

0は平方剰余と言わないらしいから, 0を除き, nubを使って重複を省き, 大きさの順にするには,

(define (quadratic-residue p)
 (sort (nub (filter (lambda (x) (> x 0))
  (map (lambda (x) (modulo (* x x) p))
   (a2b 1 (+ (quotient p 2) 1)))
   )) <))

(map (lambda (p) (quadratic-residue p)) (a2b 10 14))
=> ((1 4 5 6 9) (1 3 4 5 9) (1 4 9) (1 3 4 9 10 12))

これを見るとどれにも1,4,9があるが, xが1,2,3の時のもので当然だ.

そこで, p=2から40までの平方剰余を計算してみると

(for-each (lambda (p) (display p) 
 (display (quadratic-residue p)) (newline))
 (a2b 2 41))
=>
2(1)
3(1)
4(1)
5(1 4)
6(1 3 4)
7(1 2 4)
8(1 4)
9(1 4 7)
10(1 4 5 6 9)
11(1 3 4 5 9)
12(1 4 9)
13(1 3 4 9 10 12)
14(1 2 4 7 8 9 11)
15(1 4 6 9 10)
16(1 4 9)
17(1 2 4 8 9 13 15 16)
18(1 4 7 9 10 13 16)
19(1 4 5 6 7 9 11 16 17)
20(1 4 5 9 16)
21(1 4 7 9 15 16 18)
22(1 3 4 5 9 11 12 14 15 16 20)
23(1 2 3 4 6 8 9 12 13 16 18)
24(1 4 9 12 16)
25(1 4 6 9 11 14 16 19 21 24)
26(1 3 4 9 10 12 13 14 16 17 22 23 25)
27(1 4 7 9 10 13 16 19 22 25)
28(1 4 8 9 16 21 25)
29(1 4 5 6 7 9 13 16 20 22 23 24 25 28)
30(1 4 6 9 10 15 16 19 21 24 25)
31(1 2 4 5 7 8 9 10 14 16 18 19 20 25 28)
32(1 4 9 16 17 25)
33(1 3 4 9 12 15 16 22 25 27 31)
34(1 2 4 8 9 13 15 16 17 18 19 21 25 26 30 32 33)
35(1 4 9 11 14 15 16 21 25 29 30)
36(1 4 9 13 16 25 28)
37(1 3 4 7 9 10 11 12 16 21 25 26 27 28 30 33 34 36)
38(1 4 5 6 7 9 11 16 17 19 20 23 24 25 26 28 30 35 36)
39(1 3 4 9 10 12 13 16 22 25 27 30 36)
40(1 4 9 16 20 24 25 36)

これを使って描いたのが次の図である.


p=200までを描くと


これは, 最初に述べた, Mathworldにあった図と同じである.

PostScriptのプログラムは

40 40 translate
/d 2 def
/dot {2 dict begin /b exch def /a exch def a d mul
b d mul moveto d 0 rlineto 0 d rlineto d neg 0 rlineto
closepath fill end} def

1 1 200{/p exch def
1 1 p 1 sub{/x exch def
/a x x mul p mod def
a 0 gt {p a dot} if} for} for

のようになっている.

高さ1,4,9のような横線の他に, 右下がりの線も目立つ. つまり座標でいうと(5,4)(6,3)(7,2)(8,1)の線; (9,7)(10,6)(11,5)(12,4)(13,3)(14,2)(15,1)の線などだ. 最初の線はxとyの和が9, 次のでは和が16なのに気づく.

最初の(5,4)は, 5を法として, 4は平方剰余であるということだが, 4に法の5を足してみると9になり, 3掛ける3を5で割った余りが4なのである. その次の(6,3)は同じ9は6を法として3であるということで, この線は9を9未満の数で割った剰余の線. 次は16を16未満の数で割った剰余の線であった.

Ulam Spiral

WikipediaにUlam Spiralという項目があった.

1から図のように螺旋状に自然数を配置し, それが素数ならその場所に点を打つというものである.





私は100万くらいまでの素数のビット表を持っているから, こういう図を再現するのは何でもない.

Wikipediaには縦横200ドットの図があるので, 同じものを書いてみた.

まったく同じ図が出来て安心する.




PostScriptのプログラムはこんな具合いだ.

/n 1 def /x 0 def /y 0 def /d 2.5 def /p 1 def
/q {n primep {x y dot}if /n n 1 add def} def
/x+ {q /x x d add def} def
/x- {q /x x d sub def} def
/y+ {q /y y d add def} def
/y- {q /y y d sub def} def
100 {
p {x+} repeat p {y+} repeat /p p 1 add def
p {x-} repeat p {y-} repeat /p p 1 add def} repeat

primepは引数nが素数なら真, 合成数なら偽を返す. dotはx, yに点を打つ.