19. 多項式クラス(12)
例えば、
1 + x + x16
を因数分解するとき、
1次式、2次式、で割れるかどうかチェックしていき、
最後に8次式で割れるかをチェックします。
このとき8つの異なる整数 k について
f(k)
を計算しなければなりません。
オーバーフローを避けるためにはなるべく絶対値が小さな
k について計算したいところですが、4か-4を使うことは避けられません。
どちらについても、
416
を計算しなければなりませんが、
これはすでに32ビットの int ではオーバーフローです。
しょうがないので、64ビットの int を導入したのが次のプログラムです。
factorize1g.cpp
これで、
1 + 3x + x15
を因数分解したところ、2.6秒とやや遅くなりました。
仕方のないところですね。
肝心の16次の多項式についても計算してみましょう。
2 + 2x + x16
で、24.2秒かかりました。
これに若干の変更を加えたのが、
factorize1h.cpp
です。
あとからさらに表面化してきますが、
f(k) の素因数分解に時間がかかることがあります。
因数定理で因数分解する部分と素因数分解の部分の順序を逆にして、
f(k) = 0 になったら素因数分解しなくてもいいようにしました。
あと、前回最後に分数計算を一部整数計算にしましたが、それを全面展開しました。
これで、
2 + 2x + x16
に17.6秒かかりました。
あと、これから次元が高くなってくると、
爆発的にループの回数が多くなっていくため、
これを減らす機構をループ内に入れます。
そのためループ内の処理が複雑になり、
1回のループにかかる時間が長くなります。
これを短くするためにループ内ではなるべく
STL を避けて静的な配列を使うようにしました。
上で計算したように、
32ビットの整数を使うと16次以上の因数分解ができないことがわかりましたが、
64ビットだと、
624 < 263 < 625
となって、25次までの因数分解ができないことがわかります。
24次までしか扱わないことを前提に配列の大きさを前もって決めることができます。
ただ、約数の配列はそういうわけにはいかないので新たにクラスを作りました。
factorize1i.cpp
同じ多項式の因数分解に、11.9秒かかりました。
問題の逆行列の最後の列の要素は全て整数になります。
b0 = a00c0 + ... + a0ncn
...
bn = an0c0 + ... + anncn
と書くと、a
in は全て整数です。
よって、b
in の整数チェックに
c
n は関係ないことになります。
今まで c
n も回していましたが、
整数チェックを通ったときのみ、c
n を回すことにしましょう
（ここの説明とソースとでは表記が違う、b → a、c → b となってます）。
factorize1j.cpp
同じ多項式の因数分解に、6.2秒かかりました。
当たり前ですが、半分になっています。
前のページにも少し書いたのですが、
n-1行目が最初に整数になるかどうかチェックしています。
bn-1 = an-1 0 c0 + ... + an-1 n cn
ここで、
an-1i = (-1)n-1-in-1Ci/(n-1)! (i < n)
だから、係数を整数計算にするために、
b'n-1 = a'n-1 0 c0 + ... + a'n-1 n cn
a'n-1i = n!an-1i
として計算しているのでした。
だから単純に考えると、ここの整数チェックをクリアする確率は 1/(n-1)! で、
上の例だと、
1/(8-1)! = 1/5040
になるはずです。実際に調べてみると、16777216回のループに対し、
チェックをクリアしたのは5628回と、約1/3000になっています。
要するにこのチェックを速くやることが重要で、
あとはほとんど時間はかからないということです。
だからこの先はこのチェックをいかに速くするかのみを考えます。
bn-1 = a'n-1 0 c0 + ... + a'n-1 n cn
の計算をするとき、c
i を回しますが、
c
0 を回しているときは、次は不変で、
a'n-1 1 c1 + ... + a'n-1 n cn
c
i を回しているときは、
a'n-1i+1ci+1 + ... + a'n-1ncn
が不変です。
だから、この計算結果は配列を作ってそこに格納しておけばよいです。
こうして作ったのが、
factorize1k.cpp
これで、今までの因数分解が、2.9秒となりました。