ルート・パワー問題の高速化

riverplus氏のルート・パワー問題。
(1+sqrt(2)+sqrt(3)+sqrt(5))**nを展開した時の整数項を答える問題。
log(n)ステップ必要で、行列解法だと(3つの素数が登場するので2**3をTとすると)1ステップごとに、O(T**3)が必要である。これをO(T**2)にすることができる。

以下、第n項と書く時、0-indexであるとする。

(1+sqrt(2)+sqrt(3)+sqrt(5))**nを展開した時に、A+Bsqrt(2)+Csqrt(3)+Dsqrt(6)+Esqrt(5)+Fsqrt(10)+Gsqrt(15)+Hsqrt(30)と表すことにする。
ここで、sqrt(5)の項をあえて第4項、すなわち2の冪番目に置いていることに注意。

第n項において、nの各ビットが立っていればその素数を因数に持つことを示す。
sqrt(10)は第5項なので、1+4。10は第1項の2と第4項の5の積ということ。
素因数分解の結果得られた集合を2進数にマッピングしているとも言える。第n項の値の素因数の個数はpopcnt(n)個。

例えばFsqrt(10)とGsqrt(15)の積は、5FGsqrt(6)となる。
sqrt(10)は第5項すなわち集合(2,5)、sqrt(15)は第6項すなわち集合(3,5)である。
5FGsqrt(6)の5は(2,5)と(3,5)の共通集合を表すので第5 and 6項、6は集合(2,3)すなわち和集合から共通集合を除いたものなので第5 xor 6項の値ということができる。

これは一般化できるので、乗算は次のように定義できる。

T=[1,2,3,6,5,10,15,30]
M=10000000

def mult(x,y)
	a=[0]*T.size
	T.size.times{|i|T.size.times{|j|
		a[i^j]=(a[i^j]+T[i&j]*x[i]*y[j])%M
	}}
	a
end

あとは繰り返し二乗法を用いれば良い。

Ruby 初版
Ruby 素数一覧からTを自動計算する版
C++ 素数一覧からTを自動計算する版

Rubyの自動計算版でA=[2,3,5,7];N=11とした結果がWolframAlphaと一致したので大丈夫だと思います。。

追伸

1. 素数以外の平方根が問題に登場する場合にこの乗算が定義できるかどうかはまだ確かめていません。
2. neko_the_shadowさんの答案はmult()を書き下したものです。
3. 平方根の性質から考えてビット演算に行き着くのは自然な流れだと思うのですが、違うのでしょうか…。実は行列を持ち出す考えには至ってないです。行列の第1列は私のプログラムで言うeに一致するのですが、理由がよくわかっていません。
4. (a&b)+(a^b)==(a|b)であることにいまさら気が付きました。