続き

CodeIQというところで出題された結城さんの問題について、話の続きです。問題およびフツーの解答については前編をご覧ください。

高速化を夢見て

概要

前編より、以下の2通りの計算を元にして実装できることが分かっています。

$A_{L}=q_{\lfloor\frac{L+1}{2}\rfloor}-1,\ p_{0}=2,\ q_{0}=1,\ q_{n}=p_{n-1}+q_{n-1},\ p_{n}=p_{n-1}+q_{n}$
$A_{L}= \begin{pmatrix}1\ \ 0\end{pmatrix} \begin{pmatrix}2 \ \ 1 \\ 1 \ \ 1 \end{pmatrix}^{\lfloor\frac{L+1}{2}\rfloor} \begin{pmatrix}1 \\ 1 \end{pmatrix}-1$

それぞれでどれくらい高速に処理ができるのか、CodeIQの実行君 ( ideoneの企業版、実行時間制限 1 秒 ) を使って試してみました。

naiveな実装

こちらは、単純な足し算のループで済む方です。

$A_{L}=q_{\lfloor\frac{L+1}{2}\rfloor}-1,\ p_{0}=2,\ q_{0}=1,\ q_{n}=p_{n-1}+q_{n-1},\ p_{n}=p_{n-1}+q_{n}$

ループをLの半分繰り返すため、見た目の計算量は $O(L)$ ですが、計算が進むほど扱う数字の桁数が増加していくことを考えると、多倍長処理も含めて計算量は $O(L^2)$ です。

前編ではBrainf**k版の移植ということで、Cでchar型配列に10進数の桁を保持する実装を載せましたが、64bitになるべく情報を詰め込むようにして、10¹⁸進数として改めて実装します。

試したところ、L=150,000 を0.88秒で処理できました。

高速化実装

naiveな実装は分かり易いのですが、いかんせん計算量 $O(L^2)$ です。そこで行列計算の方なのですが…。
※ちょっと違う形に変えます。

$A_{L}= \begin{pmatrix}1\ \ 0\end{pmatrix} \begin{pmatrix}2 \ \ 1 \\ 1 \ \ 1 \end{pmatrix}^{\lfloor\frac{L+1}{2}\rfloor} \begin{pmatrix}1 \\ 1 \end{pmatrix}-1= \begin{pmatrix}0\ \ 1\end{pmatrix} \begin{pmatrix}2 \ \ 1 \\ 1 \ \ 1 \end{pmatrix}^{\lfloor\frac{L+3}{2}\rfloor} \begin{pmatrix}1 \\ 0 \end{pmatrix}-1$

ループ回数は $O(\log{L})$ で済むものの、メインとなる演算は多倍長の乗算です ( naiveな方は加算だけで済んでいました )。桁数Nの乗算を単純に実装すると $O(N^2)$ になってしまいますから、全体として $O(L^{2}\log{L})$ 、これでは意味がありません。
そこでいい機会なので、高速版の多倍長演算を組む練習にしてみました。

下準備

行列演算をフルに組んでも良いのですが、今回扱う行列の性質から計算を簡略化することができます。

$\begin{pmatrix}2 \ \ 1 \\ 1 \ \ 1 \end{pmatrix}= \begin{pmatrix}a+b \ \ a \\ \ a \ \ \ \ b \end{pmatrix}\ \ ( 但し a=1,b=1 )$

実は、パラメータ2つで表せる形の行列であり、その積についても同様です。法則は次のようになっています。

$(a,b)\times (c,d)\rightarrow (ac+bc+ad,ac+bd)$
$(a,b)^2 \rightarrow (a^2+2ab,a^2+b^2)$

ということで、行列計算の代わりに、この積を使って実装していきます。

高速な多倍長演算

さて、多倍長演算を高速化するにはどうするか? それは離散フーリエ変換(DFT)を使う方法 ( 参考: 離散フーリエ変換を用いた多倍長乗算が知られています。その高速版アルゴリズム FFTは桁数 ( 要素数 ) nに対して $O(n\log{n})$ ですから、単純な実装の $O(n^2)$ よりも高速です。

ただ、この方法の場合、桁数が増える毎に浮動小数点演算の丸め誤差が蓄積されていくため、正しい値を計算するためには誤差の影響を見積もる必要があります。~~それは面倒くさいため、~~もっと別の方法がないか、ということで、今回はFFTによく似たFMT ( あるいは整数環FFT、正式名は良く分かりません ) を使う方法を試してみました。

要素数Nの離散フーリエ変換が次のような複素数 ( 浮動小数点 ) 演算ですが…

順変換: $X_{k}=\sum_{n=0}^{N-1}\omega^{-kn}x_{n}$
逆変換: $x_{n}=\frac{1}{N}\sum_{k=0}^{N-1}\omega^{nk}X_{k}$
ただし、 $\omega$ は 1 の N乗根 $\omega=\cos{\frac{2\pi}{N}}+i \sin{\frac{2\pi}{N}}$

FMTは、ある整数 p に対する剰余環 Z/pZ ( 要するに mod p の世界 ) に対する整数 ( 固定小数点 ) 演算です。

順変換: $X_{k}\equiv \sum_{n=0}^{N-1}\omega^{kn}x_{n}\ mod\,p$
逆変換: $x_{n}\equiv N^{-1}\sum_{k=0}^{N-1}(\omega^{-1})^{nk}X_{k}\ mod\,p$
ただし、 $\omega$ は 1 の N乗根 $\omega^{N}\equiv 1\ mod\,p$

ほとんど同じ形ですが、こちらは精度を気にすることなく組むことができます。その代わりに、mod p の世界なので、各要素の計算結果が一時的に p を超えるようだと値に狂いが生じます。そこに気を付ける必要があります。

ということで、FMTの実装例 ( 参考:高速剰余変換 ) を元に、Stockham 型 self-sorting FFT ( 参考: 円周率.jp ) ぽく組んでみました。p等のパラメタは整数環FFTにある表から持ってきています。

ちなみにフーリエ変換ですので線形性があり、一旦多倍長整数を変換してしまえば、加算・定数倍も変換したデータの方で済ませてしまうことができます。変換・逆変換をいちいち繰り返すと手間なのでこれは助かる性質です。

実行君で試したところ、L=700,000 を0.99秒ぎりぎりで処理することができました。…意外に遅いですね。まあ、剰余演算を多用するのでしようがないのかもしれません。それに練習なので特にチューニングしたわけでもないですし ( 今回の問題なら分割統治法の方が速いかも )。
ただそれでも計算量はしっかり落ちていて、 $O(L(\log L)^2)$ になっています。

おまけ ( 短いコード )

前編で挙げた次のコード、もっと短く書けると言いましたので一応。

p,q=2,1;1.step(gets.to_i,2){q+=p;p+=q};p q-1

初期値、2,1 ではなくて 1,1 でもいいんですね。そうすれば q の代わりに p に答えが来ます。

p=q=1;1.step(gets.to_i,2){p+=q+=p};p~-p

ああ、いえ。どっちでもいいんですね。どちらでも39文字にはなります。

q=1;1.step(gets.to_i,p=2){p+=q+=p};p~-q

おわりに

まとめるとこんな感じに。

naiveには計算量 $O(L^{2})$
多倍長計算にFMTを使った場合計算量 $O(L(\log{L})^2)$ が可能

結城さんの解説にあったフィボナッチ数列でのアプローチについては取り上げていないですが、1つの問題で色々なネタが楽しめるのは面白いですね。

ange1のブログ

結城浩の「マヨイドーロ問題」問題解答 ( CodeIQ ) 後編

続き