2011年12月27日火曜日

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

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

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

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

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

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

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

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

m>0について,
x0⊕x1⊕...⊕x2m=<¬x0s1s1...s2m>
ここで, sj=<x0xjxj+1...xj+m-1¬xj+m¬xj+m+1...¬xj+2m-1>,
ただし, k≥1について, x2m+k=xkとする;

と書いてあった.

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

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

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

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

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

x0⊕x1⊕x2=<¬x0<x0x1¬x2><x0x2¬x1>>

x0⊕...⊕x4=<¬x0<x0x1x2¬x3¬x4><x0x2x3¬x4¬x1><x0x3x4¬x1¬x2><x0x4x1¬x2¬x3>>

となる.

この式をテストしてみよう. ます多数決関数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入力が許されるとしてこの回路はこう描ける.



これをさらに納得すべく, s1,...s4や, その計算の元になっているx1+x2-x3-x4
などを, x0=0である前半の16個について出力してみよう. 自己相対関数なので, 半分やればいいことになっている.

xの否定は1-xと計算し, (l-1)/2の2と比較したから, 閾値と較べる値は, 本来は
    x1+x2+(1-x3)+(1-x4)=2+x1+x2-x3-x4.

従って, 先頭の2を外して,
    x1+x2-x3-x4
の和sumについてs1はsum≤0なら偽(x0=0を除外したから, 等号が要る), sum>0なら真である.

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

ところが, 真が偶数個だと, どこかに相棒同士で和が0になるものがあり, x0は0だから, sは2つとも偽になる. 先頭の¬x0が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個の値は, 左からx1+x2-x3-x4, s1, x2+x3-x4-x1, s2, x3+x4-x1-x2, s3, x4+x1-x2-x3, s4, exである. 偶数個の場合, どの対で和が0になるかを調べたのが下だ. 赤の-は, x1とx2の和がx3とx4の和に等しい(つまりs1s3が0になる)場合, 青の-は, x2とx3の和がx4とx1の和に等しい(つまりs2s4が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の数は偶数である. s1とsm+1の範囲で1の数が同じであれば, 一発で和が同じところは見つかったわけでハッピーだ.

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

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

「なるほど」であった.

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

2011年12月4日日曜日

平方根の連分数展開

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

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

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

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

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

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

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

漸化式はP1=0, Q1=1から始まるが, 問題は, Qn-2があるので, Q0が必要なことだ.

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

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

2011年12月3日土曜日

素因数探し

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

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

aが素数pの平方剰余であるとは, x2=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について, a2-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から順に増やしながら, a2-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で簡単に実験が出来, ありがたい時代だ.

2011年11月30日水曜日

平方根の連分数展開

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であって,
vn+1=(d-un2)/vn,
an+1=floor((√d+un)/vn+1),
un+1=an+1vn+1-un
を繰り返す. aに√dがあるが, 整数部分をとるので下の方はいらない. vもuも整数のままで, 誤差が入らないから, いつまでも計算が続く.

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

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



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

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

2011年11月19日土曜日

素因数探し

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, 終了.

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

上のプログラムをSchemeでコーディングしてみる. ptは素因数の列の配列なので, Schemeではリストpにする. A5とA7でtを1増やすのは, 配列pの添字を進めるので, リストではいらない. t←t+1, pt←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」は pkを素数, p1 ≤ p2 ≤...≤ ptとして,

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の作り方については, またいつか書こう.

2011年11月7日月曜日

素因数探し

昨日のブログ(篩を使った素因数探し)は, 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をとったりするので, 上のプログラムより遅かった. プログラムを動的に生成するという考えもあるが, 分かり難くもなり, 今はためらっている.

2011年11月6日日曜日

素因数探し

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

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

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

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



最後は0〜m-1のxについて, (x2-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がすべての素数の法で赤の位置に対応したとき, つまり, この表で篩われた時に, x2-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へ来る. 結果を調べ, 必要なら再起動するようになっていた.

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



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

2011年10月28日金曜日

素因数探し

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

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について, a2の正方形からb2の正方形を引いた面積をnにするのである. 灰色の部分がnになる.

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

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



素因数がいくつもあると, aとbの組はいくつもあったりする. n=5005なら下の2つの図のように, 712-62も732-182も5005になる.

求め方はこうだ. 最初a=floor(√n), b=0とする. つまりa2がnに等しいか, ぎりぎりに近いが少し小さ目にする. そしてr=a2-b2を計算し, 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)2-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からb2を引いてr<nになったところへ, bより大きいa2を足すのだから, b2を引くまえの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)

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

2011年10月26日水曜日

平方剰余

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

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

まずaがpの平方剰余であるとは, 0<x<pのxについて, x2=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未満の数で割った剰余の線であった.

2011年10月1日土曜日

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に点を打つ.

2011年9月25日日曜日

再帰曲線

このブログにdragon曲線のことを書いたのは, 2年以上も前であった.

dragon曲線は, 紙を同じ向きに半分半分と折り, 折り目を90度になるように広げたものであった.

最近, この折る向きを毎回逆にしたらどうなるかと思った. 紙を折ってやってみるのは, やはり3,4回が限度であり, 手元の計算機の威力を借りたくなった.

下の図を見て欲しい.



左端のAは折る前である. 下の黒丸は出発点を示す. その上の方に右向きの矢印が示すように, 最初は右に折る. するとBのようになる. 紙の長さが倍にになっているが, 今は折り方が問題なので, 分かりやすく描いてある.

出発点から途中右折したので, その角には+が付けてある.

Bを今度は左へ折る. 最初の折り目の+は左の下へ移動する. Cの出発点から辿ると, 最初は左折, 次がBで出来た右折, さらに右折する.

Cを次は右へ折るとDになる. 上の方の+や-は, 今回出来たのもで, 下の方のは左がCの下にあった折り目, 右がCの上にあったものだ.

出発点からの右折左折を抜き出すと

B +
C - + +
D + - - + + + -
この伝でつづけると次は
E - + + - - - + + - + + + - - +
となる.

一段上の列の要素を間に+と-を交互に入れたものになっているが, dragon曲線の場合と同じで, このパターンに着目する.

EをX3, DをX3とすると, X3の左半分はX2の+-を反対にしたものだ. 右半分は中央だけが違って右半分を殆んど同じである.

そこで, 'で+-の反転を表わすとすると,
X3=X2'+Y2'.

結局
X0=+, Y0=-,
Xn=Xn-1'+Yn-1',
Yn=Xn-1'-Yn-1'
となることが分かる.

Schemeでプログラムしてみる. この引数のnは上の漸化式の2n-1になっている.

(define (r s) (if (eq? s '+) '- '+))

(define (x n i) (cond ((= i n) '+)
((< i n) (r (x (/ n 2) i)))
((> i n) (r (y (/ n 2) (- i n))))))
(define (y n i) (cond ((= i n) '-)
((< i n) (r (x (/ n 2) i)))
((> i n) (r (y (/ n 2) (- i n))))))

(map (lambda (i) (x 8 i)) (a2b 1 16))
=> (- + + - - - + + - + + + - - +)


これで右折左折の情報が得られたから, いよいよ交互折りの曲線を描くことする.

そして出来たのが次の図である. X1からX8までが, 色を変えながら描いてある. 左中ほど上の, 直角に下の曲がった青の線がB(X1)である.



左から右向きに出発するのが, 出発点である. その下の緑の線がC(X2)に対応する.

折角描いてみたが, 竜にも鳳凰にもならず, 単に三角形が出来ただけであった. 思うにdragon曲線を考えた人も, これもやっては見たが, つまらない結果だったので, このことは書いて置かなかったのかもしれない.

つまらない絵しか描けないことが分かっただけでも, 一応の知見が得られたというべきか.

2011年9月2日金曜日

整数立方根

平方根に較べ, 立方根はおおいに疎外されている感があるが, 立方根を計算したいこともある.

その場合, 実用的には
(expt 9270 (/ 1 3)) => 21.00680051861007
なことで済ませてしまう.

さて, Warrenの「ハッカーの楽しみ」を眺めていたら, 整数立方根という話題があった(訳書の226ページ). こういうプログラムが書いてある(32ビット用だ).



同書の巻頭の「推薦の辞」に私が書いたように, 本書にはそのアルゴリズムでよい理由が殆んど書いてない.

このアルゴリズムは下の説明のように出来ているのだ.

例によってSchemeに書き直すと




最初のsの設定が30でなく18なのは, そんなに大きい数でテストすることもないからだ. また最後にyだけでなく, xも返すのは, 余りも欲しいからである.

途中の経過を見るdisplayが3行ある. 9270の立方根を計算してみる.



すなわち, 立方根は21で, 剰余が9. 整数の範囲で開立をやめ, 剰余が得られるから, 整数立方根というわけである. 同じものを十進法で書くと,



つまりbの始めの方は, xより小さくなるまで218, 215, 212,... と減っていき, 4096がxより初めて小さくなった時点でyを1とする. yは最後は二進で10101だが, その最初の1(=24)が決ったのである.

xの方はそれを引いて, 9270-4096=5174になった.

次はyのその下のビットを1にするか, 0にするかを決めることである. y=16の時, その下のビットは8だから(16+8)3=163+3*162*8+3*16*82+83=13824

しかし, 最初の4096はすでに引いてあるので, 比較するbは, 上の4段目の13824-4096=9728である. こうするよりは, y=2, (y+1)3-y3=3y(y+1)+1を作り, 左シフトする方が楽というのが, プログラムの趣旨であろう.

8のビットは立たなかった. 次は4のビットで, (16+4)3-163=3904である. これはxより小さいからyは1増えて101になった.

このようにしてyのビットを上から次々を立てていき, 1ビットごとにきめていく. 最後に残ったxが開立の剰余である.

このプログラムの場合は, 方針は最初から見え見えであったが.

2011年8月28日日曜日

多面体描画道楽

Archimedes多面体13種のうち, rhomb(斜方とか菱形)という形容詞のつくものが2つある. 斜方立方八面体と斜方二十十二面体である. 何ゆえに斜方といわれるのか.





これらの図に示すように, この2種類の多面体には, それぞれ3種類の面がある.

斜方立方八面体では, 青は元々の立方体の面, 緑は正八面体の面である. 斜方二十十二面体では, 青は元々の正十二面体の面, 緑は正二十面体の面である.

では, それぞれの赤の面はなにか.

赤の面を延長して, 緑の面の上方での交点を決め, 赤の正方形に外接する菱形を考えることが出来る. 赤は元々そういう面の多面体の面である.

この多面体がどちらも菱形になるので, 斜方といわれる所以である. つまり, 斜方(赤)立方(青)八(緑)面体, 斜方(赤)二十(緑)十二(青)面体という命名であったわけだ.

ではその斜方多面体を描いてみよう.


上は斜方立方八面体である. x, y, zの対称軸も描いてある. z軸(青)のまわりに少し回転し, y軸(緑)のまわりにも僅かに回転したようになっている.


そのz軸, y軸まわりの回転をやめ, x軸方向から眺めたのが, この上の図だ. そこで, 最初の赤い正方形を含む面を考えてみると, この図の赤線で示す菱形が得られる. その座標が分かれば, 次の図が描けるわけだ.



これは, 対角線の比が1:√2の菱形が12個で構成されているので, 斜方十二面体(rhombic dodecahedron)という.


斜方二十十二面体の方は, 面が多いせいか, 多少手ごわい.



これが元の斜方二十十二面体で, 上と同様にz軸, y軸について回転してある.



回転角を0とし, x方向から眺めると, このように見えるはずだ. 赤線はそこにかぶせた菱形である. この方は, 菱形30枚で構成され, 斜方三十面体(triacontahedron)という.



菱形の対角線の比は1:φ(黄金比)である.

同一の菱形だけで出来る多面体はこの2つだけらしい. またKeplerはこの2つの多面体のあることを知っていたといわれる.

2011年8月27日土曜日

3シリンダ機関車

ずいぶん前のこのブログに3シリンダ機関車のことを書いた. そこで引用しているURLが, 最近なくなっているのに気づいたので, もう一度3シリンダ機関車の弁装置について書くことにする.

今でも蒸気機関車は人気があり, 各地で復活されているが, その構造はあまり理解されていないに違いない.

簡単にいってしまえば, 機関車の左右の前方につけた2個のシリンダ(気筒)の中にピストンを入れ, ピストンの両側から交互に蒸気を補給してピストンを前後に動かす. ピストンの往復運動を主連棒で繋ぎ, 動輪を回転して機関車を進める. ピストンがどちらかの端に寄っていると, 蒸気の補給が出来ず, 死点といって力が出せないが, 反対側のピストンが最も力の出る中央にいるように, 左右のピストンの位相を90度づらせて配置してある. このことは, 鉄道博物館で実物の機関車を見ると確認出来る.

ピストンの前後に交互に蒸気の補給をしているのが, シリンダと並んだ蒸気室で, その中に滑り弁があり, これが蒸気室の中を前後に滑りながら, 蒸気をピストンの一方に送り込む.

Walschaertの走り装置の図はここにあるので参照されたい.

滑り弁は, 動輪に動力を伝える主連棒と90度づれたエキセン棒から, 逆転装置を経由した弁棒で押し引きされて運動する. この辺の装置は大変に込み入っているから, 説明は省略!

ところで, 表題の3シリンダ機関車は, 2つのシリンダでは力不足の時, 第3のシリンダを両シリンダの中央に搭載し, 合わせて3本の主連棒で動輪を駆動する方式である. 動輪軸の中央はクランク状になっていて, そこに中央の主連棒が接続されている.

問題は, 左右の動輪のエキセン棒のようなものが, 中央のクランクからも出ているかということだが, そういう複雑な走り装置にはなっておらず, 中央のシリンダの弁の制御は両側の弁の制御から構成していたというのが前回のブログであった.

下の図を見て欲しい. Cylinder AとCylinder Bと書いたのが, 両端のシリンダの蒸気室で, 左(機関車の後方)から弁棒が入り, 途中に滑り弁があり, 弁棒の延長が右方に突き出している. 機関車の先頭にAOD, BDCのような梃子があり, 支点AはCylinder Aの弁棒の先に, 支点BはCylinder Bの弁棒の先に固定されている. 梃子だから, 図で左右に揺れると支点の上下移動もあるはずだが, 今は各支点は左右にだけ動くとしている.

支点Oは固定点である. またOD=OA/2, つまりAが動くとDはAの半分の振幅で逆位相に動く. そのDを支点として, DB=DCの梃子があり, CがCylinder Cの弁棒になっている. そうすると, CはAとBが120度の位相差で動くと, さらに120度の位相差で動き, 3番目の滑り弁が望み通りに動くのである.

このアニメーションがここにある.

アニメの下でくるくる回る3芒星は, 120度の位相で回転することを示す.

2011年8月26日金曜日

ベクトルの和

「とっておきの数学パズル」に, 多面体の各面の, 直交外向きで, 面積に比例する大きさのベクトルの和が0になることを証明せよというのがあった. Q5.13

多分そうなりそうだが, 証明となると....

「いかにして問題を解くか」には, 「一様な四面体の重心を求めよ」という問題が解けなければ, 平面の類似の問題「一様な三角形の重心を求めよ」を考えてみるという教訓がある.

その伝でいくなら, 多角形の各辺の, 直交外向きで, 辺の長さに比例する大きさのベクトルの和が0になることを証明せよとなる.

たとえば, 下の図の左は, 最も簡単な多角形の三角形の場合で, 辺の長さをa,b,cとする. (5:3:4になっている.) 従って, 外向きのベクトルは, 黒い矢印のように出来る. ベクトルの和を計算するのだから, aの先端までbの付け根を移動し(青矢印), その先端までcの付け根を移動する(青矢印)と, 3つのベクトルで構成される三角形は, もとの三角形と相似になり, cの先端はaの付け根に一致し, 和が0なことが分かる.



しかし, こういう方法はどうだろうか. (図の右)

元の三角形の左右方向をx軸. 上下方向をy軸とする. 各ベクトルのx, y成分を考える. bはx成分しか, cはy成分しか持たない. aは, 図にax, ayと書いて示したx成分とy成分がある. ベクトルaの大きさが, 斜面aの長さに比例するなら, axはaの斜めの面の, x方向からみた断面(青線で示す)の長さに比例するわけで, この断面の長さはbの辺の長さに等しく, 従ってaxの大きさはbの大きさを同じで方向が反対である.従って, ベクトルのx成分の和は0.

同様にして, y方向の成分の和も0. 従って, ベクトル全体の和も0になる.

これを3次元の問題でも使ってみよう. ある面の, その面積に比例した直交外向きベクトルaのx成分は, ベクトルとx軸のなす角をθとすると, a cosθ, この面のx方向から見た断面積もa cosθのはずである.



したがって, この多面体の直交外向きのベクトルのx+方向の和は, x+方向から見た断面積になり, x-方向の和は同じ断面積で反対向きになり, x方向全体は打ち消す.

y,z方向も同じわけだから, 結局和は0というので証明出来たのではないか.

ところで, ある時気づいた, 多面体を水の中に置くという考えも棄て難い. ある面に立てた, 面積に比例する直交外向きのベクトルは, その面に働く水圧に相当する. 水中の多面体が流れないのは, その水圧の総和が0だからだと思うと, この証明は終りらしいが.

2011年8月12日金曜日

3での整除

TAOCPに, レジスタ上の64ビットの数が3で整除出来るかどうかの判定プログラムをMMIXで書けという問題7.1.3-215がある.

十進の場合, 各桁の和が3で整除出来るなら, 元の数も3で整除出来るのは知られている.

理由は1も10も100も..., つまり10nが3を法として1であるからだ.

二進の場合も同様に考えられる. 二進の1, 10, 100,...は, 交互に3を法として1と-1だから, ビットに右から0,1,2,...と番号をつければ, 偶数番目の1の数と奇数番目の1の数の差が3の倍数の時に3で整除出来る.

これは十進法で11で整除出来る判定と同じだ. 3003や1023や5060は11で整除出来る.

だから, (nu x)をxの横方向の和をとる関数, (band x y)をxとyとのビットごとのandをとる関数とし, 32ビットでお許し頂くと

(= (modulo (- (nu (band x #x55555555)) (nu (band x #xaaaaaaaa))) 3) 0)

ということになる. 最後にまた3での剰余をとるのは違反っぽいが, 十進の場合の各桁の和の3の整除と同じで, 大きい数をいきなり割るのよりは簡単であり, とりあえずはこういう考えでよかろう.

この辺でTAOCPの解答をみると, 一見妙なプログラムである. MMIXのアセンブリコードに慣れていない向きには申し訳ないが, コメントをつけてある.



最後の3で割った剰余をみるのに, 001001...を左シフトし, 符号ビットが1になるのを利用している. これはうまい.

ではその前は何をしたか. レジスタ0は偶数ビットと奇数ビットの横方向の和(これをr0とする). レジスタ1は偶数ビットの横方向の和(r1とする). ところで, 奇数ビットの横方向の和(r2とする)も欲しい. r2=r0-r1のはず. レジスタ0はその後32-r0になり, 2ADDUでレジスタ1の2倍に足される.

つまりレジスタ1はr1+r1+(32-r0)=32+r1-r2になり, #x249...をこれだけ左シフトすることになる. bの定数は32ビット目, 8文字目の次が#x9なので, r1=r2のとき符号ビットが1, 負数になり, 3で整除出来ることが分かる.

巧妙なプログラムであった. 斯くの如きプログラムの解読は楽しい.

2011年8月11日木曜日

多面体描画道楽

立方体は正方形6枚で出来ている. 正方形は1つの内角が90度. 従って正方形全体で内角の和は360度. それがF=6枚あるから, 2160度, つまり12π(パイ)になる.



一方頂点の数Vは8. それに対して内角の総和は(2V-4)π=12πという関係をEulerは知っていたらしい.

別の例. 正十二面体は正五角形がF=12枚. 頂点はHamilton世界巡礼でお馴染のV=20. 正五角形の内角は108度. 従って内角の総和は, 108*5*12=6480. 6480/180=36(π). 2V-4も36だ.

さらに2種類の面のある, 切頭二十面体(サッカーボール)は, 正五角形が12, 正六角形が20, 頂点は60. 内角の総和は, 108*5*12+120*6*20=20880, 20880/180=116. 2*60-4=116となる.

こう考えてみた.

Eulerの式に, 面の数F, 辺の数E, 頂点の数Vとすると, F+V=E+2という有名なのがある.

いま, ある多面体がn1角形, n2角形,...,nF角形のFの面で出来ていたとする.

n角形の内角の和は(n-2)πであったから, 計算したい内角の総和Sは
S=(n1-2)π+(n2-2)π+...+(nF-2)π=(n1+n2+...+nF)π-2Fπ.

一方, この多面体の辺の数を計算すると, 各多角形の辺は2回ずつ寄与しているから, E=(n1+n2+...+nF)/2.

従って, n1+n2+...+nF=2E. これをSの式に入れると,

S=2Eπ-2Fπ. Eulerの式から, E-F=V-2.

ゆえにS=2(V-2)π.

なるほど.

また, Descartesの定理というのもある.



左下は, 正四面体の頂点は, 正三角形が3個集まっていて, その内角の和が, 360度より足りない度合い(deficiency)が180度であることを示す. 正四面体には頂点が4個あるから,deficiencyの総和は720度. 4πである.

これを他の正多面体でやって見ると, やはり, どれも720度になるのである.

サッカーボールではどうか.

この図のように, 120度, 120度, 108度に挟まれた角は12度. それが60個あるから, やはり720度である.




Eulerの方は, サッカーボールの図でいうと, 348度の方を総和している. 60倍すれば先ほどの20880度になる.

V個のradianで表わしたdeficiencyをd1, d2,...,dVとする.

Eulerの方の内角の和は, (2π-d1)+(2π-d2)+...(2π-dV)
=2πV-(d1+d2+...+dV) と足す.

これが2(V-2)πだったから, d1+d2+...+dV=2πV-2(V-2)π=4π.

これもなるほどであった.

正多面体では, 面の数の順に正四面体, 立方体(正六面体), 正八面体, 正十二面体, 正二十面体と並べる. 正四面体の頂点は尖っているから, あとに行くほど丸くなるかと思うのは誤解であって, 正八面体も正二十面体も結構尖っている. それは正八面体は立方体より, 正二十面体は正十二面体よりも, deficiencyが大きく, 頂点の数が少ないためであろう.

サッカーボールがよく転がるのは, deficiencyが12度で, V=60もあるからである.

2011年8月4日木曜日

3つの円

「とっておきの数学パズル」に, 「3つの円」という問題があった. 5.14



平面上で離れている, 半径の違う2つの円の, 円の間で交わらない(外側の)共通接線の交点を, それらの円の焦点という. 3つの円があると, 焦点は3個出来るわけだが, その3点は一直線上にあることを証明せよというのである.

こういう図で考えてみる.


円O0は中心の座標をx0,y0, 半径をr0, 円O1は座標をx1,y1, 半径をr1とする. また共通接線の1本と, 円との接点をT0, T1, 共通接線の交点P01の座標を, x01, y01とする.

三角形O0P01T0と三角形O1P01T1とは相似で, 大きさの比はr0 : r1である. 従って
(x01-x0) : (x01-x1)=r0 : r1,   (y01-y0) : (y01-y1)=r0 : r1.

これから, x01=(r0x1-r1x0)/(r0-r1),   y01=(r0y1-r1y0)/(r0-r1).

x12, y12, x20, y20についても, 同様な式が得られる.

ここで, 3点(x0,y0), (x1,y1), (x2,y2)が一直線上にあることの検査は, このブログの2008年4月24日の「同一直線上の3点」にあるように, 3点で作る三角形の面積が0かどうかで判定出来る.

今回もP01, P12, P20の三角形の面積を求めれば良い.

まず座標は
x01=(r0x1-r1x0)/(r0-r1)
y01=(r0y1-r1y0)/(r0-r1)
x12=(r1x2-r2x1)/(r1-r2)
y12=(r1y2-r2y1)/(r1-r2)
x20=(r2x0-r0x2)/(r2-r0)
y20=(r2y0-r0y2)/(r2-r0)
だから, 三角形の面積の2倍は
x01y12+x12y20+x20y01-x12y01-x20y12-x01y20

通分して分母を払い, 計算を続けると

  (r0x1-r1x0)(r1y2-r2y1)(r2-r0)+(r1x2-r2x1)(r2y0-r0y2)(r0-r1)
+(r2x0-r0x2)(r0y1-r1y0)(r1-r2)-(r1x2-r2x1)(r0y1-r1y0)(r2-r0)
-(r2x0-r0x2)(r1y2-r2y1)(r0-r1)-(r0x1-r1x0)(r2y0-r0y2)(r1-r2)

展開すると,
 (r0r1x1y2-r0r2x1y1-r1r1x0y2+r1r2x0y1)(r2-r0)
+(r1r2x2y0-r1r0x2y2-r2r2x1y0+r2r0x1y2)(r0-r1)
+(r2r0x0y1-r2r1x0y0-r0r0x2y1+r0r1x2y0)(r1-r2)
-(r1r0x2y1-r1r1x2y0-r2r0x1y1+r2r1x1y0)(r2-r0)
-(r2r1x0y2-r2r2x0y1-r0r1x2y2+r0r2x2y1)(r0-r1)
-(r0r2x1y0-r0r0x1y2-r1r2x0y0+r1r0x0y2)(r1-r2)

さらに展開する.
+r0r1r2x1y2-r0r2r2x1y1-r1r1r2x0y2+r1r2r2x0y1
-r0r0r1x1y2+r0r0r2x1y1+r0r1r1x0y2-r0r1r2x0y1
+r0r1r2x2y0-r0r0r1x2y2-r0r2r2x1y0+r0r0r2x1y2
-r1r1r2x2y0+r0r1r1x2y2+r1r2r2x1y0-r0r1r2x1y2
+r0r1r2x0y1-r1r1r2x0y0-r0r0r1x2y1+r0r1r1x2y0
-r0r2r2x0y1+r1r2r2x0y0+r0r0r2x2y1-r0r1r2x2y0
-r0r1r2x2y1+r1r1r2x2y0+r0r2r2x1y1-r1r2r2x1y0
+r0r0r1x2y1-r0r1r1x2y0-r0r0r2x1y1+r0r1r2x1y0
-r0r1r2x0y2+r0r2r2x0y1+r0r0r1x2y2-r0r0r2x2y1
+r1r1r2x0y2-r1r2r2x0y1-r0r1r1x2y2+r0r1r2x2y1
-r0r1r2x1y0+r0r0r1x1y2+r1r1r2x0y0-r0r1r1x0y2
+r0r2r2x1y0-r0r0r2x1y2-r1r2r2x0y0+r0r1r2x0y2

同じrの項をまとめると,
r0r0r1(-x1y2-x2y2-x2y1+x2y1+x2y2+x1y2)
r0r0r2(+x1y1+x1y2+x2y1-x1y1-x2y1-x1y2)
r0r1r1(+x0y2+x2y2+x2y0-x2y0-x2y2-x0y2)
r0r2r2(-x1y1-x1y0-x0y1+x1y1+x0y1+x1y0)
r1r1r2(-x0y2-x2y0-x0y0+x2y0+x0y2+x0y0)
r1r2r2(+x0y1+x1y0+x0y0-x1y0-x0y1-x0y0)
r0r1r2(+x1y2-x0y1+x2y0-x1y2+x0y1-x2y0
  -x2y1+x1y0-x0y2+x2y1-x1y0+x0y2)

となり, かっこ内はプラスマイナス打ち消し, めでたく0になった.

実際に上のように書いて計算したわけではない. 変数の順を覚えておき, 引数だけを使って計算した. 最後に, ブログ用に変数名を追加するには, utilispのプログラムを書いて実行した.

ところで, 数式処理システムrisa/asirを使ってみたら,

X01=(r0*x1-r1*x0)/(r0-r1);
X12=(r1*x2-r2*x1)/(r1-r2);
X20=(r2*x0-r0*x2)/(r2-r0);

Y01=(r0*y1-r1*y0)/(r0-r1);
Y12=(r1*y2-r2*y1)/(r1-r2);
Y20=(r2*y0-r0*y2)/(r2-r0);

X01*Y12+X12*Y20+X20*Y01-X12*Y01-X20*Y12-X01*Y20;

を入力すると, 瞬時に0が返り, 安堵する.

さて, 例によって, 幾何学的に証明できないかと考える. 幾何の証明で大切なのは, 図をなるべく正確に描くことだ. 幸い, このブログの最初の図を描いたPostScriptのプログラムがあるから, それを多少修正すると以下の図が得られる.

今度は3つの円O0, O1, O2の中心をA, B, Cとしている. また, 焦点をE, F, Gとした. BからEFに平行線を引き, AFとの交点をDとする.

円の中心O0, O1と焦点P01の距離の比は, 半径r0, r1の比であった.

EはA,Bの焦点だから, BE/AE=r1/r0
FはA,Cの焦点だから, CF/AF=r2/r0

BD//EFにしたから, DF/AF=r1/r0
従って, CF/DF==r2/r1

また, GがB,Cの焦点だから, CF/DR=r2/r1=CG/GB
従って, FG//BD//EF

QED.

この方が楽だなぁ.

2011年8月2日火曜日

条件付きの倍数

三角関数表, 対数表をはじめとして, 数表の存在感の薄い昨今である. なにしろ手元の計算機で関数値が一発で計算出来る世の中だからだ. だが今でも, 理科年表には, 10.0から100.9までの四桁の対数表があり, なんとなく空しい.

計算機が手元にない, 古き良き時代には, 数表は重要であった. Charles Babbageが階差機関を考案したのも, 数表を自動作成したいからであった.

江戸時代に歩いて日本地図を作成した伊能忠敬に関連した展示会にいってみると, 漢数字で書かれた各種の数表も並べられていて, そんな昔から数表が使われていたことに驚く.

ところで対数表の話だ. 高木先生の解析概論には気になる脚注がある. 「Wolframノ表ニハ1000以下ノ素数ノ自然対数ノ50桁ノ表ガ掲ゲラレテヰル. コノ表ハ既ニ少年がうすガ愛用シタモノデアル. 」だ.

なるほど, ある数の対数は, その素因数の対数の和なので, 素数の対数表だけあれば充分なのであった. しかし50桁も足すのは大変であったろう. また常用対数に変換するには, ln 10で割らなければならず, これも一仕事であるが, 昔の人は計算は得意であったからだとまずは納得する.

常用対数表は, スコットランド人のJohn Napierが1614年に作成, 出版したといわれている. そこで, 2014年に, 出版400年を記念してエディンバラで会合を開こうという話が進んでいる.

Napierの後, 何人もが常用対数表を計算し, 刊行したが, 先頃, Edward Sangの対数表の話をウェブで見つけた. Sangは1875年に, 素数について28桁の常用対数表を完成させた. その計算方法は次のようだ.

小さい素数の方から順に計算していく. さて, 素数1619の常用対数を計算するには, その1619を何倍かし, 下の方の桁を...999999のようにする. 9は何桁でもよいが, 多い方がよい. 例えば

714021*1619=1155999999

を作る. この計算が, このブログのタイトルの, ...9999で終るという「条件付きの倍数」である.

この倍数の作り方は次の通り.

まず1619の1から9までの倍数表を用意する.

(map (lambda (x) (display (list x (* x 1619))) (newline))
(a2b 1 10))
(1 1619)
(2 3238)
(3 4857)
(4 6476)
(5 8095)
(6 9714)
(7 11333)
(8 12952)
(9 14571)

被乗数が2と5以外の素数だから, 最後の桁には1から9が出揃う.

次にこの表を見ながら, 下の桁から積の各桁が9になるように1619の倍数を足していく. まず1倍の1619が9で終るから, 1619と書く. 下から2桁目は, 1だから, 最後が8で終る倍数, 3238をそこに足す. すると33999になる. 4桁目は3だから, 最後が6で終る倍数, 6476を足す. 今度は6509になる. 0のところを9にしなければならにから, 9で終る1619をもう一度足す. 2269になる. 6のところを9にするので, 11333を足す. こうしてここまでで1155999999が得られる.

1619 1
3238 2
----------
33999
6476 40
----
6509
1619 1
-----
2269
11333 7
-----
11559

この数を素因数分解すると, 1155999999=3×7×11×11×281×1619 である. いま常用対数を求めようとしている1619以外の素因数は, すべてこれより小さいから, それらの常用対数は分かっているはずだ.

従って, 1155999999の常用対数が得られれば, 1619の常用対数は,
log 1619=log 1155999999 - log 3 - log 7 - log 11 - log 11 - log 281
で得られるわけである.

ところで, log 1155999999=log(1156×106-1).
1156=17×17×2×2だから, この常用対数は既知である.

一方, x=1/(1156×106)とおいて,

log(1156×106-1)=log(1156×106×(1-x))=log(1156×106)+log(1-x).

自然対数をlnと書くと, ln(1+x)のTaylor展開は

ln(1+x)=x-x2/2+x3/3-x4/4+...    (|x|≤1, x≠-1)

だったから,

log(1-x)=ln(1-x)/ln 10=(-x-x2/2-x3/3-x4/4-...)/ln 10

Schemeで計算してみよう. Schemeのlog関数は自然対数lnのことである.

(define x (/ 1. (* 1156 1000000)))

x=>8.650519031141868e-10

(/ (* x x) 2)=>3.741573975407383e-19

(/ (* x x x) 3)=>2.157770458712447e-28

(log 10)=>2.302585092994046

常用対数のlog 1156×106

(/ (log (* 1156 (expt 10 6))) (log 10))=>9.06295783408451

従って log 11559999として

(- 9.06295783408451
(/ (+ x (/ (* x x) 2) (/ (* x x x) 3)) 2.302585092994046))
=> 9.062957833708824

が得られる. そこでlog 1619は

(- 9.062957833708824
(/ (+ (log 3) (log 7) (log 11) (log 11) (log 281)) (log 10)))
=>3.2092468487533754

log 1619をSchemeで直接計算すると

(/ (log 1619) (log 10)) => 3.2092468487533736

最後の2桁が違うがまぁよかろう.

こんな数値計算をするもの久しぶりだ.

2011年7月30日土曜日

条件付きの倍数

「とっておきの数学パズル」に0,1,2という問題があった. 2.2

「十進表記で0と1のみからなる自然数nの0でない倍数が存在する」というのである.

計算してみると0から9のnについて, 次のようになる. つまり3は37倍すると1だけの十進数111となる.

n 倍 十進表記
0 1 0
1 1 1
2 5 10
3 37 111
4 25 100
5 2 10
6 185 1110
7 143 1001
8 125 1000
9 12345679 111111111

この9の場合はタイガー計算器でよく遊んだ計算である.

十進表記の求め方を考えてみた. %で剰余を表すと, 3の場合は, 1%3=1, 10%3=1, 100%3=1だから, 111ではその和が3, つまり剰余が0, よって111は3の倍数である.

7の場合は, 1%7=1, 10%7=3, 100%7=2, 1000%7=6だから, 和が7の倍数になるように1と1000を選び, 1001が7の倍数となるわけだ.

そこで早速プログラムを考えてみる.

10のべきをnで割った剰余のべき集合の和のnによる剰余が0のものを求めたい.

まず(0)なるリストを用意する. これは空のべき集合である. 次に10の0乗, 1をnで割った剰余1を, これまでのリストの各要素に足し, (1)が出来, これを(0)の前にappendする. (1 0)が出来る.

次は10の1乗, 10をnで割った剰余を, これまでのリストの各要素に足す. nが7なら, 剰余3を足し, (4 3)が出来る. appendすると(4 3 1 0)になる. それぞれ(11 10 1 0)を7で割った剰余である.

同様に100を7で割った剰余2を足すと(6 5 3 2)が出来, appendして(6 5 3 2 4 3 1 0)になる. この段階では, 最後の0以外に7の倍数はないので, まだ続ける.

1000でやると剰余は10の時の剰余3と100の時の剰余2の積で6. よって新しいリストは(12 11 9 8 10 9 7 6)となり, べき集合に7の倍数7が現れた.

7は新しいリストにあるから1000番台は1, リストの右半分にあるから, 100番台は0, そのまた右半分にあるから, 10番台も0, その(7 0)では左にあるから, 1番台は1で, 結局1001が7の0と1による倍数(の最小のもの)であった.

Schemeのプログラムは次のとおり.

(define (zeroone n e rs)
(let* ((r (modulo e n))
(rs1 (map (lambda (x) (modulo (+ x r) n)) rs)))
(if (member 0 rs1)
(- (length (member 0 (append rs1 rs))) 1)
(zeroone n (* e 10) (append rs1 rs)))))

引数のeは10の順次のべきで, 最初の起動では1, rsはリストで最初は(0). let*で用意するrはこの桁の剰余, rs1は新しいリストである. let*の後の本体のifは, 新しいリストに0があるかを見, あればその位置を返す. なければ次の桁へ進む.

起動は次のようにする.

(zeroone 13 1 '(0))

9が結果の値だが, これを二進で1001としたのが求める十進表記になる. n=2から20までの値は

(map (lambda (x) (zeroone x 1 '(0))) (a2b 2 21)) =>
(2 7 4 2 14 9 8 511 2 3 28 9 18 14 16 29 1022 25 4)

二進で示すと

(map (lambda (x) (number->string (zeroone x 1 '(0)) 2))
(a2b 2 21)) =>
("10" "111" "100" "10" "1110" "1001" "1000" "111111111"
"10" "11" "11100" "1001" "10010" "1110" "10000" "11101"
"1111111110" "11001" "100")


ところでこういう倍数が存在する理由はなにか. 10のべきをnで割った剰余を表示してみると:

(define (foo n)
(map (lambda (x) (modulo (expt 10 x) n)) (a2b 0 20)))

(foo 2) =>
(1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0)
(foo 3) =>
(1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1)
(foo 4) =>
(1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0)
(foo 5) =>
(1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0)
(foo 6) =>
(1 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4)
(foo 7) =>
(1 3 2 6 4 5 1 3 2 6 4 5 1 3 2 6 4 5 1 3)
(foo 8) =>
(1 2 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0)
(foo 9) =>
(1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1)

のようで, 0があるものはその桁を1つとれば完成. 0がないものは剰余だから当然同じ数が何回も現れる. その桁をn個足せば当然nの倍数が得られるわけだ.

面白い問題であった. ところでこういう条件付き倍数の生成を対数表の計算で利用した人がいた. Edward Sangのその話は次回に.

2011年7月28日木曜日

板チョコ分割

数日前に書いた私のツィッター

割れ目でm×nのピースになる板チョコ. 割れ目に沿って割りピースにまで分ける時の割る回数は, 割る順に関係なく一定とパズル本にあった. 証明は易しい. トーナメントの試合数を思い出す.

について書いてみたい.

これは先頃頂いた, ピーター・ウィンクラー著, 坂井, 岩沢, 小副川訳「とっておきの数学パズル」にあったものだ. 8.11

やってみよう.



2×3だと, たしかにすべて5回だ.



3×4だと, 11回. ということは, m×n-1かもしれない.

そこで数学的帰納法のお出ましとなる.



m×nに分割する回数をf(m,n)と書き, m×n-1であろうとして証明する. f(1,1)は分割するまでもなく1区画になっているから, 1×1-1=0だ.

n≥1,n'≥1,n+n'>1について, nとn'の幅に分割したとする. 分割に1回かかるから, f(1,n+n')=f(1,n)+f(1,n')+1={仮定により}(n-1)+(n'-1)+1=n+n'-1. 横1列の板チョコの分割はf(1,n)=n-1で良さそうである.




今度は縦がm≥1,m'≥1でm+m'列ある板チョコを, m列とm'列に分割したとする. 分割に今度も1回かかるから,

f(m+m',n)=f(m,n)+f(m',n)+1={仮定により}(m×n-1)+(m'×n-1)+1
=(m+m')×n-1

確かにf(m,n)=m×n-1であった.


こんなに簡単な式だと, もっと別のエレガントな考え方がありそうだ.




板チョコは最初 ひと(1)塊である. 1回分割すると, ふた(2)塊になる. そのいずれかを分割すると, 3塊になる. つまり, t回分割すると, t+1塊が得られる. 斯くしてm×n塊が得られるまでには, m×n-1回の分割が必要なのであった.


問題を解くとき, とりあえず絶対に答がえられそうな, 正面からの道をとる. それで答が得られてから, 次にもっとうまい手がないか探すのが, 私の普段のやり方だ. プログラムを書くときもそうで, きたなくてもよいから, まず動くプログラムを書く. すると見通しがよくなり, さらにきれいで速いプログラムが得られるのである.

2011年7月23日土曜日

ラスターグラフィック

TAOCPV4F1はビット操作が前半の話題で, 当然ラスターグラフィックの話も登場する. そのうち, 細線化アルゴリズムの話は前に書いた. 今回は, 塗り潰しについて.

ビットマップが閉曲線で囲まれている時, その内部を塗り潰したい. PostScriptなら, 閉曲線を描き. strokeの代りにfillとすればよい. ProcessingならFill(色); としておくと, その色で塗り潰してくれるから, 簡単である.

それを自分でやる方法が, filling algorithmである.

TAOCPによると, 正方形ピクセルの2次元配列を, 円錐曲線が正方形の中心のいずれの側を通るかにより, 縦横の線の境界で分離するのは簡単だそうである.

下の図はTAOCPの7.1.3(171)を私が描き直したもので, 黒の曲線で示す, 円錐曲線の代表選手, 長軸20, 短軸10の楕円の長軸の右端の辺りのピクセルの様子である. 右下の黒丸の座標が(20,0)だ.



ピクセルの格子点は, x, yが整数値のところに対応している. 座標の整数値の間にピクセルがあるから, ピクセルの中央の座標には, 1/2が付くわけだが, それはいやらしいから, ピクセルの右上の隅の座標で指定する.

従って, ピクセル(20,0)は, 黒丸の左下のピクセルである.

さて, このピクセルを塗るべきか否かは, ピクセルの中央における, 楕円の式の値の正負で決める. 丁度0のところを楕円が通る. 楕円の式は(x/20)2+(y/10)2-1=0だが, ピクセルの中央の値は, xとyからそれぞれ1/2を引いて計算する.

つまり, ピクセル(x,y)の中央の値は,
((x-1/2)/20)2+((y-1/2)/10)2-1
=(x2-x+1/4)/400+(y2-y+1/4)/100-1

1/4や1/100がついて回って, これもやめたいから整数係数にする

400倍して
x2-x+1/4+4y2-4y+1-400
4倍して
4x2-4x+1+16y2-16y+4-1600
=4x2+16y2-4x-16y+1+4-1600
=4x2+16y2-4x-16y+-1595
これをQ(x,y)とする. このx,yに20,0を入れると, その左下のピクセルの中央の値, Q(20,0)=-75が得られる.

こういう値を, この図のピクセルについて計算した結果が, ピクセルの下の方の数値である. これが分かると, 正と負のピクセルの境界が, 塗る塗らぬの境界になる.

この図には, すべてのピクセルの値が計算してあるが, 曲線が左上へ単調に伸びると分かっていれば, 上隣りのピクセルの値を計算し, 負なら上へ進み, 正なら左へ進めば良い. すなわち, 今黒丸の点にいて, 上のピクセルが-75なので, 赤線のように境界を上に延ばす. その上も-43なので, 上へ延ばす. その上は21だから左へ行く. その上は-131だから上へ. というように進むと, 左上へ単調に進む第1象限の部分の境界が得られる.

ところで, 一般的な円錐曲線の式Q(x,y)=ax2+bxy+cy2+dx+ey+fの順次の計算は, Qx(x,y)=2ax+by+d, Qy(x,y)=bx+2cy+eを保持していると, 簡単に求まるということを見つけた知恵者がいて, これを3レジスタ法という. それを使ったアルゴリズムがTAOCPにあるが, 例によってgoto文だらけなので, Scheme風に書いたプログラムが以下である. 関数定義の後のTと番号のコメントは, TAOCPのアルゴリズムTの番号である.

(define (algorithm713t x x1 y y1 a b c d e f)
(define (q x y) (+ (* a x x) (* b x y)
(* c y y) (* d x) (* e y) f))
(define (qx x y) (+ (* 2 a x) (* b y) d))
(define (qy x y) (+ (* b x) (* 2 c y) e))
(let ((qreg 0) (rreg 0) (sreg 0))

(call-with-current-continuation
(lambda (exit)
(define (left) (display "left\n")
(bit-string-xor!
(list-ref hs (+ y 11)) (- xmax x)))
(define (right) (display "right\n")
(bit-string-xor!
(list-ref hs (+ y 11)) (- xmax (+ x 1))))
(define (up) (display "up\n"))
(define (down) (display "down\n"))
(define (initset x y a c)
(set! qreg (q x y)) (set! rreg (+ (qx x y) a))
(set! sreg (+ (qy x y) c)))
(define (update rs ab bc)
(set! qreg (+ qreg rs)) (set! rreg (+ rreg ab))
(set! sreg (+ sreg bc)))

(define (hfinish) ;T10
(cond ((< x x1) (right) (set! x (+ x 1)) (hfinish))
((> x x1) (left) (set! x (- x 1)) (hfinish))
(else (exit 'ok))))
(define (vfinish) ;T11
(cond ((< y y1) (up) (set! y (+ y 1)) (vfinish))
((> y y1) (down) (set! y (- y 1)) (vfinish))
(else (exit 'ok))))

(define (moveup) (up) (set! y (+ y 1)) ;T6
(if (= y y1) (hfinish) (update sreg b (* 2 c))))
(define (movedown) (down) (set! y (- y 1)) ;T7
(if (= y y1) (hfinish)
(update (- sreg) (- b) (- (* 2 c)))))
(define (moveleft) (left) (set! x (- x 1)) ;T8
(if (= x x1) (vfinish)
(update (- rreg) (- (* 2 a)) (- b))))
(define (moveright) (right) (set! x (+ x 1)) ;T9
(if (= x x1) (vfinish) (update rreg (* 2 a) b)))
(define (rightup) ;T2
(if (< qreg 0) (moveright) (moveup)) (rightup))
(define (downright) ;T3
(if (< qreg 0) (movedown) (moveright)) (downright))
(define (upleft) ;T4
(if (< qreg 0) (moveup) (moveleft)) (upleft))
(define (leftdown) ;T5
(if (< qreg 0) (moveleft) (movedown)) (leftdown))

(if (= x x1) (vfinish))
(if (= y y1) (hfinish))
(if (< x x1)
(if (< y y1)
(begin (initset (+ x 1) (+ y 1) a c) (rightup))
(begin (initset (+ x 1) y a (- c)) (downright)))
(if (< y y1)
(begin (initset x (+ y 1) (- a) c) (upleft))
(begin (initset x y (- a) (- c)) (leftdown))))
))))

こう定義してから, 順次の象限について

(algorithm713t 20 0 0 10 4 0 16 -4 -16 -1595)
(algorithm713t 0 -20 10 0 4 0 16 -4 -16 -1595)
(algorithm713t -20 0 0 -10 4 0 16 -4 -16 -1595)
(algorithm713t 0 20 -10 0 4 0 16 -4 -16 -1595)

のように呼ぶと, 下のような境界が得られる.



しかし, 境界だけでは塗り潰したことにならない. 塗り潰しの実装はこのようにする. 境界の横線の真上のピクセルを1に, それ以外を0のビット列の配列hs[-11..10], bs[-11..10]を用意する. 私流のプログラムでは,

(define hs
(map (lambda (x) (make-bit-string 42 #f)) (a2b -11 11)))
(define bs
(map (lambda (x) (make-bit-string 42 #f)) (a2b -11 11)))

これはMIT Schemeにある長さ42ビットのbit-stringで, 初期値はオール0である. リストの長さは22だ.

上のプログラムの, 下請け関数leftとrightのところは,

(define (left) (display "left\n")
(bit-string-xor! (list-ref hs (+ y 11)) (- xmax x)))
(define (right) (display "right\n")
(bit-string-xor! (list-ref hs (+ y 11)) (- xmax (+ x 1))))

となっていて, 横向きの境界に対応するビットを0から1にしている. その結果, 1周すると, hsとして(最下行がリストの先頭)

#*000000000000000111111111111000000000000000
#*000000000011111000000000000111110000000000
#*000000001100000000000000000000001100000000
#*000000110000000000000000000000000011000000
#*000011000000000000000000000000000000110000
#*000100000000000000000000000000000000001000
#*001000000000000000000000000000000000000100
#*000000000000000000000000000000000000000000
#*010000000000000000000000000000000000000010
#*000000000000000000000000000000000000000000
#*000000000000000000000000000000000000000000
#*000000000000000000000000000000000000000000
#*010000000000000000000000000000000000000010
#*000000000000000000000000000000000000000000
#*001000000000000000000000000000000000000100
#*000100000000000000000000000000000000001000
#*000011000000000000000000000000000000110000
#*000000110000000000000000000000000011000000
#*000000001100000000000000000000001100000000
#*000000000011111000000000000111110000000000
#*000000000000000111111111111000000000000000
#*000000000000000000000000000000000000000000

のようなビット列が得られる. hsが出来たら,

(do ((i 1 (+ i 1))) ((= i 22))
(list-set! bs i (bit-string-xor (list-ref bs (- i 1))
(list-ref hs i))))

のプログラムで, bsが作れる. bsはこういう値である.

#*000000000000000000000000000000000000000000
#*000000000000000111111111111000000000000000
#*000000000011111111111111111111110000000000
#*000000001111111111111111111111111100000000
#*000000111111111111111111111111111111000000
#*000011111111111111111111111111111111110000
#*000111111111111111111111111111111111111000
#*001111111111111111111111111111111111111100
#*001111111111111111111111111111111111111100
#*011111111111111111111111111111111111111110
#*011111111111111111111111111111111111111110
#*011111111111111111111111111111111111111110
#*011111111111111111111111111111111111111110
#*001111111111111111111111111111111111111100
#*001111111111111111111111111111111111111100
#*000111111111111111111111111111111111111000
#*000011111111111111111111111111111111110000
#*000000111111111111111111111111111111000000
#*000000001111111111111111111111111100000000
#*000000000011111111111111111111110000000000
#*000000000000000111111111111000000000000000
#*000000000000000000000000000000000000000000

タイプの文字の縦横比は1:1でなく, これは正しい形ではない. 再びMIT Schemeの機能graphicsを利用すると,

(define mydevice (make-graphics-device 'x))
(graphics-set-coordinate-limits mydevice -24 -24 24 24)
(do ((i 0 (+ i 1))) ((= i 22))
(let ((b (list-ref bs i)))
(do ((j 0 (+ j 1))) ((= j 42))
(if (bit-string-ref b j)
(graphics-operation mydevice 'fill-circle
(- 21 j) (- i 11) 0.5)))))

によって, 次のようなパターンが得られる.

2011年7月9日土曜日

多面体描画道楽

Archimedesに興味をもって, その多面体を描いてきた. 今回はその総括である. この多面体は13種ある. まとめると以下の通り.

ます, 英語名と日本語名を並べる. 次の行は(pi npi)のリストで, 正pi角形の面がnpi枚あることを示す. 多面体屋なら, この情報から稜の数と頂点の数は計算できるので, その情報をここには記載しない. その次は各頂点について, 何角形が集まっているかを示すリストである. これがこの立方体の実在正を示すキューになる.

日付がある項は, その日のこのブログに絵があったことを示す.

a. truncated tetrahedron 切頭四面体
((3 4) (6 4)) (3 6 6)
2008年10月29日

b. cub octahedron 立方 八面体
((3 8) (4 6)) (3 4 3 4)

c. truncated octahedron 切頭八面体
(4 6) (6 8)) (4 6 6)
2008年10月29日

d. truncated cube 切頭立方体
((3 8) (8 6)) (3 8 8)

e. rhomb cub octahedron 斜方立方 八面体
((3 8) (4 18)) (3 4 4 4)

f. truncated cub octahedron 切頭立方 八面体
((4 12) (6 8) (8 6)) (4 6 8)

g. icosi dodecahedron 二十 十二面体
((3 20) (5 12)) (3 5 3 5)

h. truncated icosahedron 切頭二十面体
((5 12) (6 20)) (5 6 6)
2008年10月29日

i. truncated dodecahedron 切頭十二面体
((3 20) (10 12)) (3 10 10)

j. snub cube 変形立方体
((3 32) (4 6)) (3 3 3 3 5)
2011年6月30日

k. rhomb icosi dodecahedron 斜方二十 十二面体
((3 20) (4 30) (5 12)) (3 4 5 4)
2011年6月19日

l. truncated icosi dodecahedron 切頭二十 十二面体
((4 30) (6 20) (10 12)) (4 6 8)
2011年6月29日

m. snub dodecahedron 変形十二面体
((3 80) (5 12)) (3 3 3 3 5)
2011年7月7日

それ以外のを以下に示す. 描いた時が違うので, 描画法がばらばらだが, ご容赦を. また時間のある時に揃えることもあろう.



立方 八面体

切頭立方体


斜方立方八面体


切頭立方八面体


二十 十二面体


切頭十二面体


ところで私はこれらの立体の頂点の座標を計算し, PostScriptでプログラムして描いただけだが, 東京理科大学 近代科学資料館には菱田為吉氏による, 木製の多面体模型がある. 私は資料館に行くたびに, うっとりと眺めてくる.

Platonの5個の正多面体, Archimedesの13個の多面体は, 立体マニアの興味を引き続けてきた. 私もそれに取り込まれた1人である.

2011年7月7日木曜日

多面体描画道楽

Archimedesの多面体の最後は, snubdodecahedron, 日本語では変形正十二面体という.

ものの本によると, この立体は, 三角形が80, 五角形が12で構成される. 従って稜は(3×80+5×12)/2=150, 頂点は150+2-(80+12)=60である.

正十二面体は, 正五角形の面が12, 頂点はHamilton巡礼で良く知られた20, 稜(辺)はEulerの式から12+20-2=30あった. 従って, 変形正十二面体では, 正五角形は正十二面体の面に対応し, 三角形は頂点+2×稜である.

さて, 下の図を見て欲しい. 大きい正五角形が3つ描いてある. 左の上と下は, x軸の方から見て正面になるものだ. 右の寝そべったのは, y軸方向から見えるものだ.



それぞれの正五角形の中に, 適当に縮小した正五角形を描き, それをθだけ傾ける. 縮小しただけでは, 斜方正十二面体になる. 今回は傾けるのが味噌だ. 縮小したことで 元は3つの正五角形が1点に集まっていたものが, 図のE, H, Fのように分かれ, 三角形EHFが出来る.

また元は1つの辺だったものが, AEとGFに分かれる. 傾けたため, AとFが近づき, それを結ぶと, 三角形AFGと三角形FAEが図のように出来る.

つまり1つの頂点が1つの三角形, 1つの辺が2つの三角形になるから, 前述のように20+30*2=80の三角形になる.

従って, この立体の図を描くためには, 赤で示した線分の長さが等しくなるように, 縮小比と回転角を決めなければならない.

正面上の辺AEを持つ正五角形の各頂点の立体座標が分かれば, 下の辺FGを持つ正五角形の各頂点の立体座標は, 上の座標をx軸について, 180度回転すれば得られるということから始まり, すべての面の頂点の座標はなんとか軸対称回転で計算出来, 思い通りの変形正十二面体の図が得られるはずである.

いつものように, 1辺が2の正五角形で出来た正十二面体で計算する. 正五角形の中心までの高さは, tan 54°=1.3763819204711734; 上の頂点までの高さは, tan72°=3.0776835371752527; 半径はこの差だから 1.7013016167040793. その中に半径rに内接する正五角形を作る.

その中心から見た頂点の座標をまず計算する. 72n+36°の正弦と余弦をsn0やcs0のように書くことにする.

rとθが与えられると, 以下のプログラムは, 3次元に変換し, 各点のx,y,z座標を計算する. うっとうしいが全部書くと



これらの点から, AF, FE, EH, HFの距離とAEの距離の差の絶対値の和を最小にするrとθを決めたい. こんな式は, これ以上計算する気にならないので, 今回は山登り法で最適値を探すようにした.

とりあえず, rを0.8, 0.85, ... 1.25, thetaを0, 0.5, ...0.45と変え, この関数の値をプロットすると, この範囲に最小値があるらしい.



大体r=0.95, θ=0.21 radianあたりが解らしい. PostScriptで描画するから, この程度の値で充分である.



今回の描画プログラムは, 見えない面の破線は, 一度しか描かないようにしているから, 破線がきれいに見えている.

2011年7月3日日曜日

八方睨みの猫

ご近所の知友 金子和夫氏の個展があった. 鉛筆画の「八方睨みの猫」が評判である. チラシからその一部をスキャンさせていただいた. この猫ちゃんには, 金子家で会ったことがある. 可愛い猫だ.



ところで八方睨みとは, この猫をどの方向から見ても, 猫は見られている方を向いているように見えるということである. このチラシを, 右左に傾けても, 上下から眺めても, 確かに視線はこちらを向いている.

実はこの猫だけでなく, 多くの肖像画の目付きも, 大体はこうなっている.

歩きながら月をみると, 歩くにつれて, いつまでも月が付けて来るように, 絵の前を通り過ぎても, 肖像の眼の視線はずーっと自分を追ってくる. 私は若い頃は, これを不思議に思っていた. また, それについて, おそらく黒目が眼の中央にあるからだろうと想像してた.

大学院生の頃, 研究室でこれが話題になったことがあり, 高橋秀俊先生は「黒目が真ん中だから」とこともなげにいわれ, 私としては, 決着がついている.

この猫の黒目は, 上下左右の真ん中にある. 従って四方八方を睨むことになる.

私は, ご存じのようにPostScriptで絵を描くのが大好き人間だ. で, さっそくやってみた. それぞれの絵では, 上から黒目が右より, 中央, 左よりに描いてある. やはり中央のが八方睨みに見える.



横方向から見ても, (つまり横方向を縮小しても)



のようになろう. 真昼の猫の眼のように, 細くなっても同じだろうか.