2010年1月24日日曜日

両替問題

研究所の山本君がSICP(Structure and Interpretation of Computer Programs)の両替問題がなんとかといっていたので, 両替問題をおさらいしよう. 英語ではchangeというので, Graham, Knuth, Patashnik著, 有澤他訳の「コンピュータの数学(Concrete Mathematics)」では, 「釣銭の問題」(299ページ)となっている. 正しくは, 空港のCHANGEの窓口と同じで, 外国のお金から, または高額のお金からの両替である. 数学の問題としては, 高額のお金の崩し方の話題である.

この問題はいろいろなところに書いてある. 前述のSICPや, コンピュータの数学の他, G. PolyaのHow to Solve It(如何にして問題を解くか)にもあり, 私がたまたま持っている, PolyaとTarjanの1978年のStanfordでの講義ノート, D.R.Woods, Notes on Introductory Combinatoricd(STAN-CS-79-732)にも出ている.

SICPの記述はこうだ

50セント, 25セント, 10セント, 5 セント, 1セントがあるとして1ドルの両替にはどのくらいの場合があるか. つまり, 金額に対して両替の場合の数を計算する手続きは書けるだろうか.

この問題は再帰的手続きとして単純な解がある. 使える硬貨をある順に並べたとしよう. すると次の関係が成り立つ:

n種類の硬貨を使う, 金額aの両替の場合の数は:
• (a)最初の種類の硬貨以外を使う, 金額aの両替の場合の数, 足す
• (b)dを最初の硬貨の額面金額[denomination]として, n種類の硬貨を使う, 金額a-dの両替の場合の数

つまり1セントと5セントの2種類の硬貨で10セントにするには, 最初の種類の硬貨は5セントであり, 5セント以外の硬貨, つまり1セントで10セントにする場合の数(下の図で赤枠の場合)と, 5セントと1セントで10セントから5セント引いた5セントにする場合の数(青枠の場合)の和である.

また再帰の終わり

• (c)aがちょうど0なら, 両替の場合の数は1
• (d)aが0より少なければ, 両替の場合の数は0
• (e)nが0なら, 両替の場合の数は0
も考えなければならず, 結局プログラムは以下のようになる.


(define (count-change amount)
(cc amount 5))

(define (cc amount kinds-of-coins)
(cond ((= amount 0) 1) ;(c)
((or (< amount 0) (= kinds-of-coins 0)) 0);(d)と(e)
(else (+ (cc amount ;(a)の場合
(- kinds-of-coins 1))
(cc (- amount ;(b)の場合
(first-denomination kinds-of-coins))
kinds-of-coins)))))

(define (first-denomination kinds-of-coins)
(cond ((= kinds-of-coins 1) 1)
((= kinds-of-coins 2) 5)
((= kinds-of-coins 3) 10)
((= kinds-of-coins 4) 25)
((= kinds-of-coins 5) 50)))

10セントの両替で, ccの呼ばれ方を図にすると下のようになる. 頂点の(10 2)はamountが10, kinds-of-coinsが2でccが呼ばれるの意. 左下(10 1)への枝は(a)の下請け, 右下(5 2)への枝は(b)の下請けであり, さらに下請けが次々と呼ばれることを示す. さらに右下方向へ延びる枝が, 図には全体で3本見えるが, この最後が(c)の1通りを返してきて, 結局35回呼ばれ, 3通りになる.


このプログラムの, 100セントの両替では, 関数ccは 15499回呼ばれた. その解析がこの表である.


表の横はamount, 縦はkinds-of-coinsで, 横の見出しの「--」はその間の数, つまり0と5の間なら, 1,2,3,4のことである. amountは硬貨3種類のときに-5で3回呼ばれる. これは(cc 5 3)の時で, 5セントから10セントを引こうとして生じる.

この表の数値を全部足すと, 15499になる. また数値の入っている場所は250個所ある. つまり250通りで呼ばれるのである. 赤枠内のamount=0の欄を足すと, 解の292が得られる.

メモ化

これとかFibonacci数の計算のプログラムのような再帰プログラムでは, 同じ引数で関数を何回も呼ぶことによくなる. 上の計算でも引数の異なる呼出しは, 250回である. こういう時には, メモ化という手法が有効である.

小学生の頃, 長い割り算(商の桁数が多い)をするなら, まず除数の2倍, 3倍,..., 9倍の表を作ってから割り算の取り掛かったが, メモ化は表を新しい引数の組で関数を計算するたびに作るのである.

SICPでは, 3.3.3節の最後にメモ化の話題がある. それに倣って, memo-ccを作ってみた.


(define (make-table) ;和訳SICP 159ページにあるmake-table
(let ((local-table (list '*table*)))
(define (lookup key-1 key-2)
(let ((subtable (assoc key-1 (cdr local-table))))
(if subtable
(let ((record (assoc key-2 (cdr subtable))))
(if record
(cdr record)
false))
false)))
(define (insert! key-1 key-2 value)
(let ((subtable (assoc key-1 (cdr local-table))))
(if subtable
(let ((record (assoc key-2 (cdr subtable))))
(if record
(set-cdr! record value)
(set-cdr! subtable
(cons (cons key-2 value)
(cdr subtable)))))
(set-cdr! local-table
(cons (list key-1
(cons key-2 value))
(cdr local-table)))))
'ok)
(define (dispatch m)
(cond ((eq? m 'lookup-proc) lookup)
((eq? m 'insert-proc!) insert!)
(else (error "Unknown operation -- TABLE" m))))
dispatch))

(define (memoize f) ;和訳SICP 160ページのmemoizeを2次元化
(let ((table (make-table)))
(define get (table 'lookup-proc))
(define put (table 'insert-proc!))
(lambda (x y)
(let ((previously-computed-result (get x y)))
(or previously-computed-result
(let ((result (f x y)))
(put x y result)
result))))))


(define memo-cc ;和訳SICP 22ページのccをメモ化
(memoize (lambda (a k)
(set! cccount (+ cccount 1))
(cond ((= a 0) 1)
((or (< a 0) (= k 0)) 0)
(else (+ (memo-cc a (- k 1))
(memo-cc (- a (first-denomination k)) k)))))))

これでやってみると, memo-ccが呼ばれたのはちょうど250回であった.

2010年1月23日土曜日

微分解析機

遊星歯車の動き方のアニメーションを作りたいと思っていたが, やっとなんとかなったので, 公開することにした. 試作の順に番号をつけたので, planetary2planetary4とがある.

planetary2では


またplanetary4では


のような絵が得られる.

とりあえず動きを見るために作ったのがplanetary2で, 歯車を円で近似した. しかし, 歯車の廻るところも見たくなり, planetary4を書いた. ただ歯車は早く廻転すると, ストロボ効果で, 逆転のように見えるので, 結構遅く廻転させている.

下にある赤, 緑, 橙, 青の四角をクリックすると, いずれかの行動をする. そのパターンが

である. 赤線は太陽歯車の廻転角, 緑線は内歯車の廻転角, 青線は衛星キャリヤの廻転角である.

2010年1月8日金曜日

2010

Knuth先生のホームページ
Did you know that 2009 is (1/2+3)*(4+5*(6*7+8*9))?
と書いてある. (年があけてからはHappy 2010 (i.e., MMIX++) to all!になった.)

今年は2010年. 今度はどうなるかと思ってプログラムを書いた.
結論からいうと

(- (- (+ 1 2) (* (* (* (- (- 3 4) 5) 6) 7) 8)) 9)
(* (/ (* 1 2) 3) (- (* (* (* (+ 4 5) 6) 7) 8) 9))
(+ (- 1 (* (+ 2 (* (* (- (- 3 4) 5) 6) 7)) 8)) 9)
(* (- (- 1 2) 3) (- (/ (+ 4 5) 6) (* (* 7 8) 9)))
(* (* (/ (* 1 2) 3) (+ (- 4 5) (* (* 6 7) 8))) 9)
(* (+ (+ 1 2) (* 3 4)) (+ (- 5 6) (* (+ 7 8) 9)))
(* (+ 1 (* 2 (+ 3 4))) (+ (- 5 6) (* (+ 7 8) 9)))
(* (+ (* 1 2) (* (- 3 (* 4 (- (/ 5 6) 7))) 8)) 9)
(* (* 1 (+ 2 (* (- 3 (* 4 (- (/ 5 6) 7))) 8))) 9)
(+ (+ 1 2) (* (- (* (+ (* 3 (+ 4 5)) 6) 7) 8) 9))
(+ (+ 1 2) (* (+ 3 (* 4 (+ (+ 5 (* 6 7)) 8))) 9))
(+ (+ 1 (* (* 2 (- (* (* 3 4) (+ 5 6)) 7)) 8)) 9)
(* (+ (+ 1 2) 3) (- (* 4 (+ (* 5 6) (* 7 8))) 9))
(* (* (* 1 2) 3) (- (* 4 (+ (* 5 6) (* 7 8))) 9))
(+ (+ 1 2) (* (+ 3 (* 4 (+ (- 5 6) (* 7 8)))) 9))
(* (+ (- (- 1 2) 3) (* 4 (+ (/ 5 6) (* 7 8)))) 9)
(+ (+ 1 2) (* 3 (+ (* (* 4 (+ 5 6)) (+ 7 8)) 9)))

の17通りが2010になる. 11番目と最後の加算と乗算だけの式は美しい.
規則としては, 数字は1から9までこの順に使う. 演算子は+, -, *, /だけである.

実はTAOCPの第4巻, 分冊4の7.2.1.6小項の演習問題122は,

100=1+2*3+4*5-6+7+8*9=(1+2-3-4)*(5-6-7-8-9)
=((1/((2+3)/4-5+6))*7+8)*9

などの例の後で, このような100の作り方はどれだけあるかというのである. 候補の式の数は, 10646016あり, 100になるのは1640通りあるそうだ.

100の美しい式は 1+2+3+4+5+6+7+8*9 である.

このような式で表せない最小の正の整数は, 2284だそうだから, まだ当分楽しめる.

上の17通りの全解探索に要した計算時間は, MacBookAirで50分弱であった.
MIT Schemeを使ったが, 分数が使えるので, こういう計算には便利である.


後記
折角プログラムを書いたので, 100になる式をやってみた. ところが1640個ではなく1641個が得られた. Knuth先生にメイルしたところ, 珍しく早々にメイルで返事が来た. もちろん先生の秘書から来る.
それによると
2004年に書いたプログラムにはミスがあり, 1641が正しい.
先頭に単項の-をつけると2010になる式には他に26通りある.
ということだ.

2010年1月4日月曜日

微分解析機

微分解析機の加算装置に, 車のデフ型でなく, 遊星歯車型の差動歯車を使う話を書いたが, 正月休みにもう少し定量的に理解したいと思った.

遊星歯車の動き方のアニメーションがここここにある. でもこのくらいは自分で描きたい.

まず各歯車の動き方を調べた. 下の図で, 周囲の大きい円が外側の内歯車(歯が内側を向いている), その半分の半径の中央の円が太陽歯車Sである. 右の小型の円が衛星歯車で, 真横にあるのが最初の位置Pである.



太陽と衛星歯車は, AとBで, 衛星歯車と内歯車は, CとDで接している. 太陽歯車がA'まで廻転し(回転角をφ), 衛星キャリヤがP'まで来た時, キャリヤの廻転角をθとする. この時, 遊星歯車と内歯車はC"とD"で接している. その移動中に接した歯車は円弧DD"とC'C"でこの長さは等しい. 一方AとBはすでにA'とB'まで移動していて, 今はA"とB"で接している. この方の移動中に接した歯車は円弧A"A'とB"B'でこれも前の長さに等しい.
内歯車, 太陽歯車, 遊星歯車の半径をそれぞれ, R, r, ρ; B"P'B'の角をαとする.
DD"=Rθ=C'C"=ρα=B"B'=A"A'=r(φ-θ)
従ってキャリアの廻転角は θ=φr/(R+r);
衛星歯車の廻転角は α=θR/ρ=φRr/(R+r)ρとなる.

一方, 太陽歯車を固定し, 内歯車をψだけ廻転すると, キャリヤの廻転角θも同様に計算出来る.

θ=ψ*R/(R+r);
α=-ψ*Rr/(R+r)ρとなる.

この両者を合わせ, 太陽歯車の廻転角 φと内歯車の廻転角 ψから, キャリヤの廻転角θは
θ=φ*r/(R+r)+ψ*R/(R+r);
遊星歯車の廻転角αは
α=(φ-ψ)*Rr/ρ(R+r).

前回の歯車の図や, 上の図ではR=200,r=100,ρ=50なので,
θ=φ/3+2*ψ/3;
α=4*(φ-ψ)/3となる.

例えば, 太陽歯車を90度進めると, φ=90からθ=30度, α=120度(上の図). 次に内歯車を90度進めて追い付くと, ψ=90からθ=60度, α=-120度で, キャリヤも90度まで進み, 衛星歯車も最初の接していた両側の点で再び接するようになる.

通信機用のバリコンのつまみは, 太陽歯車に相当し, バリコンはキャリヤ軸についていたとすると, 1/3の減速になるのであった.