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;
	}
};