Round 1A-C-Large

GoogleCodeJam

Google Code Jam Round 1AのC-Largeの問題「3 + sqrt(5) の20億乗の整数部を1000で割ったあまりを求めよ」(sqrtは平方根)という問題の解き方を教えてもらったので、自分の理解を確認するのもかねて解説を書いてみる。あと、数学的な記号の使い方に抵抗がある人も多いみたいだからプログラミング言語的な表現にこだわって書いてみる。

        • -

まず背景知識として

式1: a[0] == A0
式2: a[1] == A1
式3: a[i + 2] == s * a[i + 1] + t * a[i]

という条件が成り立つ数列aの一般項を求める方法について解説する。つまり、a[10000]の値は下のような9998回のループで求まるが、それをせずにいきなり求める方法について解説する。

a[0] = A0
a[1] = A1
for i in range(9998):
    a[i + 2] = s * a[i + 1] + t * a[i]

もし式3: a[i + 2] == s * a[i + 1] + t * a[i] が二つの定数p, qを使って
式4: a[i + 2] - p * a[i + 1] == q * (a[i + 1] - p * a[i])
と表すことができると仮定する。

そうすると、
a[2] - p * a[1] == q * (a[1] - p * a[0])
a[3] - p * a[2] == q * (a[2] - p * a[1]) == q * q * (a[1] - p * a[0])
...と掛け算の繰り返しに落とし込むことができるので、累乗**を使って
式5: a[n + 1] - p * a[n] == (q ** n) * (a[1] - p * a[0])
と表現することができる。

また、式4: a[i + 2] - p * a[i + 1] == q * (a[i + 1] - p * a[i])の右辺の括弧を外すと
a[i + 2] - p * a[i + 1] == q * a[i + 1] - q * p * a[i]
となり、左辺の2番目の項を右辺に移動すると
a[i + 2] == (p + q) * a[i + 1] - q * p * a[i]
になる。この式はpとqについて対称(pとqを入れ替えても成立する)なことから想像できるのだけど、式4のpとqを入れ替えた式
式6: a[n + 1] - q * a[n] == (p ** n) * (a[1] - q * a[0])
が成立することも式5を導くために行った式の変形をもう一度することで示すことができる。

ここでp != qだと仮定すると
式5: a[n + 1] - p * a[n] == (q ** n) * (a[1] - p * a[0])
式6: a[n + 1] - q * a[n] == (p ** n) * (a[1] - q * a[0])
この2式の両辺で差をとって
(q - p) * a[n] == (q ** n) * (a[1] - p * a[0]) - (p ** n) * (a[1] - q * a[0])
つまりa[n]の値は式7: 

a[n] == ((q ** n) * (a[1] - p * a[0]) - (p ** n) * (a[1] - q * a[0])) / (q - p)

となる。

いったんまとめ。
与えられた漸化式 a[i + 2] == s * a[i + 1] + t * a[i] が二つの異なる定数p, qを使って a[i + 2] - p * a[i + 1] == q * (a[i + 1] - p * a[i]) と表現できるなら、aの任意の項は式7を使って求めることができる。

        • -

さてようやく背景知識の説明が終わったが全然 (3 + sqrt(5)) ** n を求めろという問題が解ける気配がない。ここからが魔術的(笑)
ある1より小さい数smallを使って 式8: (3 + sqrt(5)) ** n + small ** n という式を考える。このときsmall ** n < 1なので式8がもし整数であるなら、(3 + sqrt(5)) ** nの整数部は式8の結果ひく1。そこで式8が整数になるようなsmallを探そう。

といっても small = 3 - sqrt(5) ぐらいしかそんな値は思いつかないのですけど。それぞれn乗した後での和が整数になってほしいので、sart(5)が1回だけ現れる項を消したくて、「和をとったときに消える」っていうとそりゃ符号を逆にしてみるしか。運のいいことに 3 - sqrt(5) は 1 より小さい(出来レース！)

というわけでめでたく問題は (3 + sqrt(5)) ** n + (3 - sqrt(5)) ** n を求める問題に帰着しました。しかもこの式は整数であることがわかっているので、それを利用した効率化が将来的には出来るはず！

さて式7: 

a[n] == ((q ** n) * (a[1] - p * a[0]) - (p ** n) * (a[1] - q * a[0])) / (q - p)

を思い出してみよう。
下の2つの式が成り立つようなa[0], a[1]があれば a[n] = q ** n + p ** n の形になってうれしい。

q - p == a[1] - p * a[0]
q - p == -(a[1] - q * a[0])

p = 3 - sqrt(5), q = 3 + sqrt(5)を代入して、二つの式を両辺足してみよう。
4 * sqrt(5) == 2 * sqrt(5) * a[0]. よってa[0] == 2。
ついでこのa[0]を代入して2 * sqrt(5) == a[1] - (3 - sqrt(5)) * 2. 整理するとa[1] == 6がでる。
あとは s == p + q, t == - p * qにも代入して、s == 6, t == -4。

ここまでのまとめ:
a[0] == 2, a[1] == 6, a[n + 2] = 6 * a[n + 1] - 4 * a[n]
を満たすようなa[n]を求めれば、答えがわかる！

        • -

あとは「1000で割ったあまり」は1000個しかないので、その値a[n], a[n + 1]個を使って次のペアa[n + 1], a[n + 2]を求める処理はかならず1000 * 1000回以内で元の数に戻ってくる。なるべく短い周期のループであることを期待しながら順番に計算して調べてみると、a[103]以降はa[3]以降の繰り返しであることがわかる。a[2000000000]のようなまじめに計算するとすごく時間のかかる値でも、a[100]と一致することがわかるのでまじめに計算しなくていい。a[102]までの値を計算すれば、あとはどの値でも一瞬で回答できる。

        • -

感想。漸化式の一般項を求める方法で高速に計算するのかと思いきや、逆に使って累乗の計算を漸化式に落とした上でその漸化式がループすることを利用して高速化するわけですね。

余談だけどはてなでは「かっこかっこ」を特殊な記法として扱うので書きにくいな。行中のコードもあまり読みやすくないし。こういう解説を書くならやはりTeXで書いてPDFで公開するのが一番なのかなぁ。