TopCoder SRM 436 Div1 Hard CircularShifts
問題
二つの長さnの整数列A,Bが与えられる。
Bに対して次のような、シフト操作を何度でも行うことが出来る。
- B[n-1]を新たなB[0]とし、B[0],B[1],...,B[n-2]を新たなB[1],B[2],..,B[n-1]とする
このとき、A[0]*B[0] + A[1]*B[1] + …… + A[n-1]*B[n-1]の最大値を求めよ。
制約条件
n≦60000
A[i],B[i]≦100
方針
ディジタル法またはFFT(高速フーリエ変換)を用いた高速な畳み込み乗算。
前者はO(n^1.585)、後者はO(n log n)の時間計算量になる。
後者を実装した。考え方は次の通り。
数列の各項を桁としてみた時、A,Bは多倍長の巨大整数とみなせる。
(ただし、Bは桁を逆転させておく)
これら二つの数の積を考える。
A[0],A[1],A[2]…,A[n-1] * B[n-1],B[n-2],...,B[0] --------------------------
すると掛け算の結果は、
下から1番目の桁はA[n-1]*B[0],
下から2番目の桁はA[n-1]*B[1]+A[n-2]*B[0]
下から3番目の桁はA[n-1]*B[2]+A[n-2]*B[1]+A[n-3]*B[0]
……
となる。
下から1番目の桁と1+n番目の桁を足すと、これは
A[n-1]*B[0]+A[0]*B[1]+A[1]*B[2]+……+A[n-2]*B[n-1]
下から2番目の桁と2+n番目の桁を足すと、これは
A[n-2]*B[0]+A[n-1]*B[1]+A[0]*B[2]+……+A[n-3]*B[n-1]
のようになり、掛け算の答えが、BをシフトしたときのΣA[i]B[i]になっていることがわかる。
この多倍長の掛け算を高速に計算する方法がわかればいい。
フーリエ変換によってA[i]を列a[i]に変換し、B[i]をb[i]に変換する。
各項をa[i]*b[i]とした数列c[i]というのは、
この多倍長の掛け算を各項とする数列C[i]をフーリエ変換したものになっている。
なので、フーリエ変換によりa[i],b[i]を求めた後、c[i]を掛け算により求め、
c[i]をフーリエ逆変換によりC[i]に戻すことで多倍長の掛け算が全体でO(nlog n)の計算時間でできることがわかる。
ソースコード
FFTのコードはspaghetti sourceのコードなので省略。
Complex a[131072],b[131072]; double res[60000]; class CircularShifts { public: int maxScore(int N, int Z0, int A, int B, int M) { memset(a,0,sizeof(a)); memset(b,0,sizeof(b)); memset(res,0,sizeof(res)); a[0].real()=Z0%M%100; rep(i,N-1){ Z0=(Z0*(ll)A+B)%M; a[i+1].real()=Z0%100; } rep(i,N){ Z0=(Z0*(ll)A+B)%M; b[N-i-1].real()=Z0%100; } int n=1; while(n<2*N)n*=2; fft(n,2*PI/n,a); fft(n,2*PI/n,b); rep(i,n)a[i]=a[i]*b[i]; fft(n,-2*PI/n,a); rep(i,n)res[i%N]+=a[i].real(); int ans=0; rep(i,N)ans=max(ans,(int)(res[i]/n+0.5)); return ans; } };