フィボナッチ数の話題

前回, この話題では F_n+1F_n-1-F_n² =(-1)ⁿ を証明したが, 実はもっと簡単の単であった.

n=1の時は

1·0 -1·1=-1
つまり(-1)¹ だ.

n+1の時の式は

F_n+2F_n - F_n+1²

このF_n+2とF_n+1の一つを定義により前の2項の和に書き換える.

(F_n+1+F_n)F_n - F_n+1(F_n+F_n-1) = F_n² - F_n+1F_n-1 = -(-1)ⁿ = (-1)ⁿ⁺¹

地図のミウラ折り

ミウラ折りという言葉も知っていたし, 三浦さんに会ったこともあるが, 実際にミウラ折りを折ってみたのは, 最近のことだ.

この度, 国土地理院の地図をミウラ折りしてみたくて, 関連するページを探したら, こういうページ があり, これは正に国土地理院の地形図の大きさ(縦460ミリ, 横580ミリ)に基づいていた. (この種の地図の大きさは昔から殆どかわらず, 私が昭和25年頃に買った地理調査所の地図も同じである. ただこれは紙の大きさで, 地図自体の大きさは, 昔は縦が緯度10分, 横が経度15分であった.)

もとの図は見難いので, 私が書き直した図は次のようだ. (単位はミリ)

下側の図が折り線で, 赤が山折り, 青が谷折りである. 左上の青四角は, 地図の図幅名の場所である.

この図の折り方では, 私の計算では, 疊んだときの寸法が, 縦143.6ミリ, 横105.6ミリになる.

ところでミウラ折りの説明を見ると, 紙を横長に折り, それを傾けながら折り, 山折り谷折りを調整して作るとあるが, 紙を重ねて折ると正確には折れないから, 工夫が必要である.

以下私が今回使った方法を紹介しよう.

まず地図の上に, 後で消せるように鉛筆で上の図の折り線を書き込む.

つぎにこの線に沿って折るわけだが, 山折りは紙の下に木片をいれ, その角を折り線に合せて折り目を付ける. 谷折りは定規を線に合せて, 紙を上に曲げて折り目を付ける.

こうして出来たミウラ折りの地図がこれだ. 下にあるのは, 折り目をつけるのに使った木片と, その上に乗った定規である.

今回, この方法でやってみたところ, 木片の稜を使って山折りの折り目をつけるのは, 下が見えないので, 結構時間がかかった. 一方, 定規による谷折りの線を付けるのは簡単で, 地図を裏返して, 表の折り線が透過するような環境で, 裏から谷折りにするのがいいのではないかと思っている.

また最初の折り線を記入するのも, 正確にするため, ドラフタでも使いたいところである.

フィボナッチ数の話題

TAOCPのフィボナッチ数の話題(1.2.8)には F_n+1F_n-1-F_n² =(-1)ⁿ というこれも気になる式がある.

1·0-1·1=-1
2·1-1·1=1
3·1-2·2=-1
5·2-3·3=1
8·3-5·5=-1
たしかに.

証明したくなった. 式が沢山でてくるので, Texで書いたものが下である.



式0が証明したいもの. F_nはφとそのハットを使って式1の ように書ける.

式2がそれらの定義だ.

φとφハットの積や比は, 式3と4のようだ.

さて, 式0の左辺に式1を代入したのが, 式5である.

とりあえず 1/√を忘れて, 式6のように第1項を掛ける. 第2項の積は 式7のようになる.

6から7を引くと, 式8が得られるが, 先程保留していた分母の√5が2個ずつ あったので, それで5が消え, (-1)ⁿとなる.

フィボナッチ数の話題

ご存じフィボナッチ数の話だ.

F₀=0, F₁=1, F_n+2=F_n+1+F_n
と定義され, 初めの方の値は

0	1	2	3	4	5	6	7	8	9
0	1	1	2	3	5	8	13	21	34

ところでTAOCPの第1巻に,
F_3nは2(=F₃)の倍数
F_4nは3(=F₄)の倍数
F_5nは5(=F₅)の倍数
F_knはF_kの倍数
と書いてあった.

たしかにF₃, F₆, F₉は偶数だし,
F₅と上の表にはないが, F₁₀=55も5の倍数である.

その辺を読んでみると

F_n+m=F_mF_n+1+ F_m-1F_n

という式があり, これが味噌である.

この証明は簡単だ.

F_n+3=F_n+2+F_n+1. F_n+2にF_n+1+F_nを 代入すると F_n+3=2F_n+1+F_n.
この係数の 2はF₃だし, 1はF₂なので, この式は F_n+3=F₃F_n+1+F₂F_nになる.

同様にして

F_n+4=F₄F_n+1+F₃F_nになる.

この4をmにすると上の式になるわけだ.

証明はこんな具合である.

F_n+m+1={Fibonacci数の定義}F_n+m+F_n+m-1. このそれぞれに上の式を代入する.

F_mF_n+1+F_m-1F_n
F_m-1F_n+1+F_m-2F_n
F_m+F_m-1=F_m+1だし F_m-1+F_m-2=F_mなので, 上の式が得られた.

ところでF_nkがF_kの倍数である証明も 簡単であった.

F_k=F_kだからF_kの倍数.
F_k+k=F_kF_k+1+ F_k-1F_kだからそれぞれの項にF_kがあり, F_kの倍数.

F_nkがF_kの倍数とすると
F_k+nk=F_nkF_k+1+ F_nk-1F_kとなり, 前の項にはF_nk, 後の項にはF_kがあり, それぞれF_kの 倍数なので, F_(n+1)kもF_kの倍数である.

面白い.

予定表

徐々に年末が近付きつつある. 年末独特の心配事は他にもあるが, まずは来年の予定の記帳だ. 通常の手帳は12月(かせいぜい1月かそれも第1週だけとか)で終りであり, 次の会合が翌年の始めだったりすると, 書き込む場所がない.

また例年に送られてきて使っている企業の手帳も年末にならないと届かない.

そういうわけで, この数年, もしかしたら十数年は10月ころになると, 次のような簡単な予定表をA4の紙に印刷し, 手帳に挟んで利用していた.

下に示すのは2017年版である. 左上の01S, 02M, ...は1月1日(日曜), 1月2日(月曜), ... を表す. 右端の列の31Fは3月31日(金曜)に対応する.

1日分が非常に小さいが, この予定表を作るころは, 新春の予定も少いのであまり困ることはなかった. 絵から分るように日曜だけは太字にしてある.

もちろん新しい手帳が届くまでの使命で, 手帳が来るとそちらへ転記し, お役御免になる.



ところが, 世の中では皆さんがスマートフォンで予定を管理するようになり, 学生のころから毎年使っていた手帳は, 企業でも今年(2016)の分から遂に作るのを止めた.

今年はとりあえず市販の手帳に予定を記入してきたが, やはり使い易いとは思わなかった. そこでこの暮は, 来年の予定表を自作することにした. その方が楽しくもある.



上の写真がこうして出来た個人用予定表の表紙である. 右側が表紙, 左側が裏表紙である. 大きさはB7版になる. (こうして開いた状態はB6版である.)

裏表紙はすぐ分るように, 2017年の七曜表である. しかし, 最上段が左から1月2月3月, 次の段が4月5月6月, ... ときて, よく見ると5段あり, 最下段は実は2018年の1月2月3月なのである.

この予定表の特徴は, 2017年1月から2018年3月までになっていることだ. つまり3ヶ月のバッファーを準備している.

従って, 表紙の次は



のように2017年1月になり, 裏表紙の前は



のように2018年3月である.

1日分が狭いのではないかといわれそうだが, 最近は予定表に記入する事項も少くなり, たとえ沢山書き込む日があっても, 欄外に書いて線で結ぶリンク構造にすれば楽々収まる.

そういう次第で来年はこの予定表を活用しようと考えている.

いい遅れたが, 表紙の西暦年号の下は曜日(上段)と月齢(下段)を計算する定数(枠内は左から1月2月3月 | 4月5月6月 | 7月8月9月 | 10月11月12月)である.

どちらも私のブログの話題になったもので, m月d日の曜日は左欄外の4にm月の定数とdを足し, 7で除した法をとる. 月齢は左欄外の1に足し, 30で除した法をとる. 右欄外は2018年の定数である.

曜日の定数が1月と2月で赤になっているのは, 閏年には1を引くことの注意である. 閏年には欄外の定数も赤で表示するようにプログラムしてある.

こういう予定表を作るには紙の中央で綴じる必要がある. 中綴じホチキスを探したが, 「マックス多機能スケール ナカトジール」というものを見付けそれを購入して使った.

ただ, 上の写真で分るように, 丁度折り目のところで留めるのは結構むずかしい.

予定表を吐き出すPostScriptのプログラムを書くのは, これまでカレンダーのプログラムを何度も書いたことがあるので, 簡単である. 来年, 再来年の祝日はインターネットで探すとすぐに見付かる.

ループ検出

前回のこのタイトルのブログで, Gosper流のループ検出法の説明をしたが, これでループが発見できる理由があれだけではいまいちである.

私がどう納得しているかも書いておく必要があろう. x_iが次々と置かれていく様子を添字だけで再掲する.

i	T0	T1	T2	T3
1	1
2	1	2
3	3	2
4	3	2	4
5	5	2	4
6	5	6	4
7	7	6	4
8	7	6	4	8
9	9	6	4	8
10	9	10	4	8

これらの表中のx_iに, 次のxがあるかどうかを調べるのだが, そこで添字の差だけの表にすると

i	T0	T1	T2	T3
1	-1
2	-2	-1
3	-1	-2
4	-2	-3	-1
5	-1	-4	-2
6	-2	-1	-3
7	-1	-2	-4
8	-2	-3	-5	-1
9	-1	-4	-6	-2
10	-2	-1	-7	-3

つまり1つ前2つ前に同じものがあると, それらはT₀で見付かり, 3つ前4つ前にあると, それらはT₁で見付かり,...というわけだ. Tの添字が大きいほど, 繰り返し開始時のμの値の曖昧さが大きいことがこの表から判明する.

TAOCPでは演習問題3.1-7にこの話があり, その解答には

x_nが表T_ρnに格納されると, それはその後x_n+1, x_n+2, ..., x_n+2^ρn+1と比較される.

と記載してあるが, 上の説明の方が分り易いと私は思っている.

Gosperの元の記事はHAKMEM 132 (http://home.pipeline.com/~hbaker1/hakmem/flows.html#item132)にある.

ループ検出

数列x_iが例えば32ビットの語の範囲の値しかとれなくて, x_i+1=f(x_i)を繰り返していくと, いつかは昔のある値が現れ, 以後同じ数列を辿ることになって, ループを回り始める.

このループの開始位置とループ長を求めるのがループ検出である.

そのひとつの方法にBill Gosperの考案したものがあり, あちこちに説明があるが, あまり真面目に読んだことはなかった. しかし気にはなっていたので, この数日プログラムを書いてみたりしながら, その仕掛けを調べてみた.

いまx₁, x₂, x₃, ... が1,2,3,4,5,6,7,8,9.10,6,7,8,9,10,6,7,8,...であったとしよう.(x₀から始める流儀もあり, また私はそれが好きだが, 今回の添字は1から始めることにした.) The Art of Computer Programming(TAOCP)の演習問題3.1,6の記述を借用すれば, 「x₁,x₂,...,x_μ,..., x_μ+λ-1は互いに相異なるが, n ≥μの時はx_n=x_n+λであるμとλがある.」

このμが開始位置, λがループ長で, 上の例ではμ=6, λ=5だ.

記憶場所も時間もふんだんにある時は, xの値をすべて順に記憶しておき, 新しいxについて, それを過去のすべてのxの値のそれぞれと比較する. そこに同じものがあれば, 前の値のあったところがループの開始場所であり, そこから今の値までの長さがループ長である.

しかし, こんなことはとてもやっていられないから, あの手この手を考えることになる. その一つにGosperの提案したものがある. 以下は添字などを私流に多少修正したものである.

これは簡単にいえば, 繰り返しが「現れるや否や」発見するのは諦める, 最新のxの1つ前, 2つ前, 4つ前, 8つ前, ... などを覚えておき, ループに入り同じ値が現れ始めたらなるべく早目にそれを発見しようという方針である.

「1つ, 2つ, 4つ, 8つ,...」から推察されるように, n個のxを調べた時は, log₂ n個程度の配列Tを使う.

ところで, n>0を二進法で表した時, もっとも右の1の位置, もっとも左の1の位置をTAOCPに倣い, それぞれρ(n), &lambda(n)とする. ビットの位置は右端を0とする. (上のループ長の定義にもλがあって申し訳ない.)

ρ(1)=0, λ(1)=0; 
ρ(2)=1, λ(2)=1; 
ρ(3)=0, λ(3)=1; 
ρ(4)=2, λ(4)=2; 
ρ(5)=0, λ(5)=2; 
ρ(6)=1, λ(6)=2; 
ρ(7)=0, λ(7)=2; 
ρ(8)=3, λ(8)=3; 

ρの方は別名ルーラー関数(物差関数)という.



十進法の物差しとは様子が違い, ひと目盛りが半分ずつに分割されていて, インチの物差しが大体こういう感じである.

ρとλを計算する関数 rhoとlamは次のようだ.

ρは

(define (rho x)
 (define (r x rh)
   (if (= (modulo x 2) 0) (r (quotient x 2) (+ rh 1)) rh))
 (if (= x 0) 'error (r x 0)))

(map rho (a2b 1 16)) => (0 1 0 2 0 1 0 3 0 1 0 2 0 1 0)

またλは関数logがあれば

(define (lam n)
 (inexact->exact (floor (/ (log n) (log 2))))

(map lam (a2b 1 16))=>
(0 1 1 2 2 2 2 3 3 3 3 3 3 3 3)

ρのように自力で計算するには

(define (lam x)
 (cond ((= x 0) 'error)
       ((= x 1) 0)
       (else (+ (lam (quotient x 2)) 1))))

さて配列をT₀, T₁, T₂, ... とする. (こちらは0から始める) そしてx_iをT_ρ(i)に順に置くのである.

x_iが次々と置かれていく様子は次のようだ.

i	T0	T1	T2	T3
1	x1
2	x1	x2
3	x3	x2
4	x3	x2	x4
5	x5	x2	x4
6	x5	x6	x4
7	x7	x6	x4
8	x7	x6	x4	x8
9	x9	x6	x4	x8
10	x7	x10	x4	x8

注意すべきは, この表のうち見えているのは最後の行だけということである. つまりi=8の時には, x₁, x₂, x₃, x₅はすでに消えている.

さて, 上の例のようなループがあった時の同様な配列のx₁₀=10までの変化の様子は,

i	T0	T1	T2	T3
1	1
2	1	2
3	3	2
4	3	2	4
5	5	2	4
6	5	6	4
7	7	6	4
8	7	6	4	8
9	9	6	4	8
10	9	10	4	8

のようだ. ここでx₁₁=6になった時, この配列は9, 10, 4, 8であり, この中に6はないから, 6を(ρ(11)=0だから)T₀のところに置く. したがって配列は6, 10, 4, 8になる. 次, x₁₂=7も配列にないから, 7をT₂に置き, 配列は6, 10, 7, 8になる.

次はx₁₃=8だが, i=8の時に作ったT₃に8があるので, 循環が始まっていたのがわかる.

そこで問題はiの行のT_kを置いたのは誰かということだ. 最初の表からxを削除した表は次のようだが

i\k	0	1	2	3
1	1
2	1	2
3	3	2
4	3	2	4
5	5	2	4
6	5	6	4
7	7	6	4
8	7	6	4	8
9	9	6	4	8
10	7	10	4	8

これで例えばi=7, k=2ならこの4を置いたのはi=4であったということが知りたい.

いろいろやってみると, iからiの下からk+1ビットをとったものの, 先頭の0,1を逆にしたものを引くといいらしいことが分かる.

iの下からk+1ビットをとるには, 2^k+2-1で1のビット列を作ってマスクする.

i &(2^k+2-1)

これを 2^kでxorすると先頭のビットが反転する.

それをiから引くのである. Schemeで書くと

(define (m i k)
   (- i (fix:xor (fix:and i (- (fix:lsh 1 (+ k 1)) 1)) 
    (fix:lsh 1 k))))

であって, 実際に計算すると

(do ((i 1 (+ i 1))) ((= i 11))
(display (cons i 
 (map (lambda (k) (m i k)) (a2b 0 (+ (lam i) 1)))))
(newline))

(1 1)
(2 1 2)
(3 3 2)
(4 3 2 4)
(5 5 2 4)
(6 5 6 4)
(7 7 6 4)
(8 7 6 4 8)
(9 9 6 4 8)
(10 9 10 4 8)

となって上の表と一致する.

従ってi=13の時のx_13>=8がi=12, k=3にあり, (m 12 3)=8が分ったから, λ=13-8=5であった.

プログラムを書いてやってみよう. 最初の数列x_iは, μとλから

(define (f i mu lm)
 (if (< i mu) i (+ (modulo (- i mu) lm) mu)))

(map (lambda (i) (f i 6 5)) (a2b 1 19)) => 
(1 2 3 4 5 6 7 8 9 10 6 7 8 9 10 6 7 8)

と作る. プログラム全体と実行結果を下に示す.

(define (iloop i x)
 (define (kloop k)
  (if (> k (lam (- i 1))) (iloop i x)
  (begin 
   (if (= x (vector-ref tab k)) (list i k (m (- i 1) k))
     (kloop (+ k 1))))))
(vector-set! tab (rho i) x)
(display i) (display (take (+ (lam i) 1) (vector->list tab))) 
(newline)
(set! i (+ i 1))
(set! x (f i mu lm))
(kloop 0))

(define mu 6) (define lm 5)
(iloop 1 (f 1 mu lm))

1(1)
2(1 2)
3(3 2)
4(3 2 4)
5(5 2 4)
6(5 6 4)
7(7 6 4)
8(7 6 4 8)
9(9 6 4 8)
10(9 10 4 8)
11(6 10 4 8)
12(6 10 7 8)

=> (13 3 8)

つまりi=13のとき, k=3の所に同じものがあり, その時のiは8であった. 従ってλは13-8=5であった.