「SOL」TN's Kingdom III - Assassination(POJ)

這個標題就出奇的長

---

# 題面

\(n\) 次多項式 \(\alpha\) 和 \(\beta\) 作循環卷積獲得 \(\gamma\)。如今已知 \(\beta,\gamma\)，求 \(\alpha\)。

換句話說：

\[\gamma=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\alpha_i\beta_j\cdot x^{(i+j)\bmod n} \]

數據規模：\(n< 2^{17}\)。

---

# 解析

直到如今我才知道 FFT 原本就是作的循環卷積……只不過之前用 FFT 作多項式線性卷積都是把多項式次數 \(N\) 設得特別大，而後體如今係數上就沒有循環。

可是咱們實現的 FFT 都是分治作的，要求多項式次數是二的冪。就沒有辦法作任意長度的循環卷積。

Hint.

普通的 FFT 也能夠解決一部分循環卷積問題，只須要把多項式次數設得很大，先計算出 $\alpha,\beta$ 的線性卷積，再手動把 $x^t$ 的係數加到 $x^{t\bmod n}$ 上去。

只不過這樣沒有辦法倒過來，知道 $\alpha\times\beta=\gamma$ 的 $\beta,\gamma$ 求 $\alpha$。

那就先無論 FFT，來推一下對多項式 \(f\) 求 DFT 的式子。記 \(X_k\) 表示 \(f(e^{-\frac{2\pi k}Ni})\)（也就是求 DFT 後的第 \(k\) 項），\(f_k\) 是 \(f\) 的 \(x^k\) 項的係數：

\[\begin{aligned} X_k&=\sum_{n=0}^{N-1}f_ne^{-\frac{2\pi kn}{N}i} \end{aligned} \]

Hint.

$e^{\frac{2\pi} ni}$ 就是單位復根，也能夠寫做 $\omega_n=\cos\frac{2\pi}n+\sin\frac{2\pi}ni$（是否是看起來熟悉一點了）。

接下來就是最重要的部分：

\[\begin{aligned} &nk=\frac{(n-k)^2-n^2-k^2}2\\ &e^{-\frac{2\pi nk}Ni}=e^{\frac{-k^2\pi}N}\cdot e^{\frac{(k-n)^2\pi}{N}}\cdot e^{\frac{-n^2\pi}{N}} \end{aligned} \]

因而

\[X_k=e^{-\frac{k^2\pi}{N}}\sum_{n=0}^{N-1}e^{\frac{(k-n)^2\pi}{N}}\cdot\big(e^{-\frac{n^2\pi}{N}}f_n\big) \]

能夠看出 \(X_k\) 後面的求和式是一個線性卷積，能夠用 FFT。

值得注意的是這個卷積的下標 \(k-n\) 的範圍是 \([1-N,N-1]\)，是可能有負數的。因此先給它平移 \(N\) 位。

逆變換也基本同樣：

\[X_k'=e^{\frac{k^2\pi}{N}}\sum_{n=0}^{N-1}e^{-\frac{(k-n)^2\pi}{N}}\cdot\big(e^{\frac{n^2\pi}{N}}f_n\big) \]

因此構造出這樣兩個數列：

\[A=\sum_{i=0}^{2N-1}e^{\frac{(i-N)^2\pi}{N}}x^i\\ B=\sum_{i=0}^{N-1}f_ie^{-\frac{i^2\pi}{N}}x^i \]

卷積獲得多項式 \(C\)，別忘了 \(C\) 是左移了 \(N\) 位的。

\[X_k=e^{-\frac{k^2\pi}{N}}C_{k+N} \]

這樣就實現了循環卷積。整個算法思路歸爲 bluestein，是一種把循環卷積轉化爲線性卷積的方法。

---

# 源代碼

/*Lucky_Glass*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;

struct Complex{
	double r,i;
	Complex(const double &_r=0,const double &_i=0):r(_r),i(_i){}
	Complex operator *(const Complex &v)const{
		return Complex(r*v.r-i*v.i,r*v.i+i*v.r);
	}
	Complex operator +(const Complex &v)const{
		return Complex(r+v.r,i+v.i);
	}
	Complex operator -(const Complex &v)const{
		return Complex(r-v.r,i-v.i);
	}
	Complex operator /(const double &k)const{
		return Complex(r/k,i/k);
	}
	Complex operator /(const Complex &v)const{
		return Complex((i*v.i+r*v.r)/(v.r*v.r+v.i*v.i),(i*v.r-r*v.i)/(v.r*v.r+v.i*v.i));
	}
};
Complex omega(const double &k){return Complex(cos(k),sin(k));}

const int N=(2<<17)+10;
const double PI=acos(-1.0);

int rev[N<<2];

void FFT(Complex *ary,int n,int typ){
	for(int i=1;i<n;i++){
		rev[i]=rev[i>>1]>>1|((i&1)*(n>>1));
		if(rev[i]<i) swap(ary[rev[i]],ary[i]);
	}
	for(int i=1,ii=2;i<n;i<<=1,ii<<=1){
		Complex wn=omega(PI/i*typ);
		for(int j=0;j<n;j+=ii){
			Complex wi(1,0);
			for(int k=j;k<j+i;k++,wi=wi*wn){
				Complex tmp=wi*ary[k+i];
				ary[k+i]=ary[k]-tmp;
				ary[k]=ary[k]+tmp;
			}
		}
	}
	if(typ==-1) for(int i=0;i<n;i++) ary[i]=ary[i]/n;
}
void bluestein(Complex *ary,int n,int typ){
	static Complex aryA[N<<2],aryB[N<<2];
	memset(aryA,0,sizeof aryA);
	memset(aryB,0,sizeof aryB);
	for(int i=0;i<n;i++)
		aryA[i]=omega(-typ*PI*i*i/n)*ary[i];
	for(int i=0;i<(n<<1);i++)
		aryB[i]=omega(typ*PI*(i-n)*(i-n)/n);
	int len=1;
	while(len<(n<<2)) len<<=1;
	FFT(aryA,len,1),FFT(aryB,len,1);
	for(int i=0;i<len;i++) aryA[i]=aryA[i]*aryB[i];
	FFT(aryA,len,-1);
	for(int i=0;i<n;i++){
		ary[i]=aryA[i+n]*omega(-typ*PI*i*i/n);
		if(typ==-1) ary[i]=ary[i]/n;
	}
}

int n;
Complex aryA[N],aryB[N];

int main(){
	// freopen("input.in","r",stdin);
	scanf("%d",&n);
	for(int i=0;i<n;i++) scanf("%lf",&aryA[i].r);
	for(int i=0;i<n;i++) scanf("%lf",&aryB[i].r);
	bluestein(aryA,n,1),bluestein(aryB,n,1);
	for(int i=0;i<n;i++) aryA[i]=aryB[i]/aryA[i];
	bluestein(aryA,n,-1);
	for(int i=0;i<n;i++) printf("%.4f\n",aryA[i].r);
	return 0;
}

---

THE END

Thanks for reading!

\[\begin{split} 「\ &你是沙鷗\ 你是滴漏\\ &你是白晝\ 你是不朽\\ &你是歷來都沒有\ 緘默的口\\ &你是跌落深谷裏\ 最後的柔\\ &你是行走於人間的風\\ &你是我高不可攀的夢\ 」\\ ——&\text{《你是我高不可攀的夢》 By 蒼穹} \end{split} \]

> Link 你是我高不可攀的夢-Bilibili