『列島強靭化論』を読んだ

ミッションクリティカルなシステムでは、普通は何らかの冗長構成を採って、サーバーが一台落ちたりハードディスクが一個壊れたりしても何ともないようになっている。落ちたら困る度合いに応じてお金を掛けて、一台ですむところを冗長化したりその他の安全のための工夫をしている。まあ、何てことない普通の話だと思う。

プログラマとして藤井聡教授の『列島強靭化論』を読むと、こうした冗長化などの考えを日本の国土に対して適用し、「日本」というシステムの Availability（可用性）を高めようというアイデアが見えてくる。ついこの間、自民党の国土強靭化基本法案が提出されたところなので、元ネタと思われるこれを読んでみた。

システム開発者として意訳すると、曰く、日本というシステムは言うまでもなく日本人にとってもの凄くクリティカルなのだけど、にもかかわらずアーキテクチャとして余りにも脆弱である。ここで、もちろん Single Point of Failure は東京圏なわけだけど、30年以内に直下型地震が 70%以上の確率で起きる事が既に分かっている。

異音を立て始めたハードディスクを何の冗長化もなくミッションクリティカルなシステムに使ってるようなもので、これは恐ろしい。

つうわけで、東京に集まっている経済機能、都市機能を、太平洋ベルトに入らない日本海側や九州、北海道の諸都市にも分散しようという話だけど、システム屋ならすんなり納得できるんじゃないだろうか。クリティカルなシステムに冗長性を与えようって話だから、むしろ言われる前にやっておけと。

しかも藤井教授によると、そのためのインフラ整備の公共事業がデフレ・ギャップ解消／経済成長にも役立つから一石二鳥。それに先立つ東北の復興も加えると一石三鳥という事になる。10年で200兆って事だけど、通読するとそれほど法外なお金じゃない事も分かり易く説かれている。

（藤井教授はデフレ解消をもの凄く重視しているけど、この辺りは学者によって賛否両論あるらしい。自分は藤井教授の考えに一理あると思う。まあ、新聞なんかでは単に「公共事業は悪」ってイメージだけの批判もあったようだけど、そういうのは論外。）

震災直後の1ヶ月で書き上げたとの事で、妙に臨場感があって、読むと当時の気分が蘇ってくる。

終章が面白い。もしも日本人が本書で提示された「強靭化」を日本の国土に対して施さなかった場合、この「不作為の罪」によって、どんな風に日本が終わってしまうかというストーリーが提示というか予言されている。これがもの凄くリアル。

発行から１年以上たった今、デフレ下では絶対ダメと本書で厳禁されていたはずの消費増税も肯定されてしまいそうな勢いだし、TPPも交渉参加だし、これは悪い方の予言が当たって、「30年以内に70%の確率」で日本終了のお知らせって事になっちゃうかもしれん・・・

余りが１になる割り算から元の数を求める方法

1000003 × B ÷ 1000000 = q 余り 1 となる最小のBは何か？

これを素朴に求めると、B＝1から開始して初めて余りが 1になるまでBをインクリメントするような事になるが、ずっと効率的なやり方があったのでメモっておく。

以下、euler 134 のネタバレあり

==== ==== ==== ====

project euler に、与えられた素数 p1, p2 について、末尾が p1 で、なおかつ p2 で割りきれる数 S をもの凄く大量に計算する問題があった。

合同式で書くとこうなる

S ≡ 0 (mod p2)
S ≡ p1 (mod m) (m はp1を越える最小の10の正数乗)

実は合同式自体あまり詳しくないので、"連立 合同"でググってみると、連立合同式の解法を紹介しているサイトがあった。

これをよく読んで、上の問題に当てはめてみると、

  S ≡ p1*p2*B (mod p2*m)
  ただし、B は p2*B ≡ 1 (mod m) を満たす整数)

を解けば良い事がわかる。

試しにやってみると、

p1 = 19, p2 = 23
23*B ≡ 1 (mod 100) より、B = 87
(19 × 23 × 87) ÷ (23 × 100) = 16 余り 1219
答え： 1219

となり、一応計算できてることが分かる。

ただし困ったのが、B を求めるやり方で、素朴にやると以下のようになるけど（Haskell）

ghci> find (\a ->1==mod (23*a) 100) [1..]
Just 87

これだと冒頭であげた1000000近辺の数だとかなり時間がかかる。

そこで更に調べてみると、「逆元（inverse）」なるものを求める方法があって、どうやらこれが使えそう。

ググりまわって Wikipedia でこの擬似コードを見つけた

remainder[1] := f(x)
remainder[2] := a(x)
auxiliary[1] := 0
auxiliary[2] := 1
i := 2
while remainder[i] > 1
    i := i + 1
    remainder[i] := remainder(remainder[i-2] / remainder[i-1])
    quotient[i] := quotient(remainder[i-2] / remainder[i-1])
    auxiliary[i] := -quotient[i] * auxiliary[i-1] + auxiliary[i-2]
inverse := auxiliary[i]

これを、以下のように解釈して、表計算ソフトに打ち込んでみる。

n	r	a	q
1	m	0	-
2	p2	1	-
・・・
i	MOD(ri-2, ri-1)	= - qi * ai-1 + ai-2	= DIV( ri-2, ri-1)

で、r が1 になったときの a が求める数字らしい。

試しに 1000003 でやってみると以下のようになる。

1000000	0
1000003	1
1000000	0	0
3	1	1
1	-333333	333333

-333333 が負数になっているので、mod (-333333,1000000)としてみると 666667 が得られる。 1000003 × 666667 ＝ 666669000001で、確かに 1000000 で割ると余りが 1 になる事が分かる。

もの凄く計算量が減ってる。ははは、これは面白い。

上の問題に戻って、p1 = 999983, p2 = 1000003 としてみると、mod (999983 * 1000003 * 666667) (1000003000000) = 666662999983 となり、1000003 で割りきれて末尾が 999983 になる数が確かに得られた。

上手く行きそうなので、表計算ソフトに書いた事を以下のような Haskell コードに置き換えてみる。

inverse a m = mod (f m 0 0 a 1 0) m where
    f 1  a  _  _  _  _  = a
    f r1 a1 q1 r2 a2 q2 = f r a q r1 a1 q1
      where r = mod r2 r1
            q = div r2 r1
            a = -q*a1 + a2

euler 134 に適用してみると、1秒で正答が得られた。（ちなみに「余り１になるB」をインクリメントしながら計算すると、4時間くらいかかった。）

ピタゴラス数の小ネタ

今までプログラマやってきて特に気にすることもなかったけど、project euler を始めてみるとピタゴラス数に関連する問題に出会うことがある。で、いろいろ調べてみると、ちょっと面白いネタというか Tips が見つかった。

■ 原始ピタゴラス数の無限リスト

このサイトで原始ピタゴラス数（primitive Pythagorean triple (PPT)）を次々に生成する行列が紹介されている。こんな行列。

（ -1　 -2　 2）
（ -2 　-1　 2）
（ -2　 -2　 3）

この行列に符号を3通りに変えたピタゴラス数を掛けると、子のピタゴラス数が3つ生成される。例えば(3,4,5)から始めると以下のようになる。

(-1 -2 2) (-3)   ( 5)
(-2 -1 2) ( 4) = (12) 
(-2 -2 3) ( 5)   (13)

(-1 -2 2) ( 3)   (15)
(-2 -1 2) (-4) = ( 8) 
(-2 -2 3) ( 5)   (17)

(-1 -2 2) (-3)   (21)
(-2 -1 2) (-4) = (20) 
(-2 -2 3) ( 5)   (29)

算出された子のピタゴラス数から、さらに３つずつの孫を作ることができ、これを再帰的に繰り返すとピタゴラス数の無限リストを作ることができる。Haskell で書くとこんな感じになる。

pyt1 (a,b,c) = ((-1)*a+(-2)*b+2*c,(-2)*a+(-1)*b+2*c,(-2)*a+(-2)*b+3*c)
pyt2 (a,b,c) = [pyt1(-a,b,c), pyt1(-a,-b,c), pyt1(a,-b,c)]
pyt3 ts = ts ++ pyt3 (concatMap pyt2 ts)

以下のような結果になる。

ghci< pyt3 [(3,4,5)]
[(3,4,5),(5,12,13),(21,20,29),(15,8,17),(7,24,25),(55,48,73),(45,28,53),(39,80,89),(119,120,169),(77,36,85),(33,56,65),(65,72,97),(35,12,37),(9,40,41),(105,88,137),(91,60,109),(105,208,233),(297,304,425),(187,84,205),(95,168,193),(207,224,305),(117,44,125),(57,176,185),(377,336,505),(299,180,349),(217,456,505),(697,696,985),(459,220,509),(175,288,337),(319,360,481),(165,52,173),(51,140,149),(275,252,373),(209,120,241),(115,252,277),(403,396,565),(273,136,305),(85,132,157),(133,156,205),(63,16,65),(11,60,61),(171,140,221),(153,104,185),(203,396,445),(555,572,797),・・・

■ 高さ＝底辺の二等辺三角形に漸近する漸化式

このサイトでは、ピタゴラス数を直角二等辺三角形に近づけていく漸化式が紹介されていた。直角を構成する辺 a と b の差を1に こんな式。

an+2 = 6*an+1 - an + 2  a1 = 3, a2 = 20
cn+2 = 6*bn+1 - bn - 2  b1 = 4, b2 = 21
cn+2 = 6*bn+1 - cn      c1 = 5, c2 = 29

で、ここから類推して高さ＝底辺の二等辺三角形（の片側の三角形）に漸近する漸化式もきっとあるだろうなと思って探してみた（きっかけは euler 138）。

まず直角を構成する2辺のうち長い辺と短い方を二倍したものとの差が1であるものを何個か探す。これは上述の無限リストに対して、条件に合うものをフィルタリングした。以下のようなものが見つかった（これ以上は数が大きくなりすぎて難しい）。

(15,8,17)
(273,136,305)
(4895,2448,5473)
(87841,43920,98209)
(1576239,788120,1762289)
(28284465,14142232,31622993)
(507544127,253772064,567451585)
(9107509825,4553754912,10182505537)

で見つかったピタゴラス数から連立方程式を立ててこれを解いた。手で計算するのは面倒なので、ネット上の連立方程式を解くプログラムを使った。例えば、斜辺ならこんな式を解けば良い。

305*x + 17y + z = 5473
5473*x + 305*y + z= 98209
98209*x + 5473*y + z =1762289

解）x=18, y=-1, z=0

で、c₁ = 17, c₂ = 305, c_n+2 = 18*c_n+1 - c_n という漸化式が得られる。Haskell で書くとこうなる。

l = 17:305:zipWith (\a b->18*b-a) l (tail l)

同様に、直角を構成する2辺についても漸化式が得られる。 a1 = 15, a2 = 273, a_n+3 = 17*a_n+2 + 17*a_n+1 - a_n (bも同様)

試しに20番目の要素を見てみると、a= 10151021471800938910964641, b = 5075510735900469455482320, c = 11349187026003431978487841)で、a² + b² = c², 1 = |2b-a| が確かに成立する。

==== ==== ==== ====

さっき problem 138 を解いて、euler project を150問解いた事になる。簡単な問題を選んでやってきたが、そろそろ難しくなってきた。かなり頑張ったけど、まだ半分にも至ってないのか…