「SOL」TN's Kingdom III

这个标题就出奇的长

---

# 题面

(n) 次多项式 (alpha) 和 (eta) 做循环卷积得到 (gamma)。现在已知 (eta,gamma)，求 (alpha)。

换句话说：

[gamma=sum_{i=0}^{n-1}sum_{j=0}^{n-1}alpha_ieta_jcdot x^{(i+j)mod n} ]

数据规模：(n< 2^{17})。

---

# 解析

直到现在我才知道 FFT 本来就是做的循环卷积……只不过以前用 FFT 做多项式线性卷积都是把多项式次数 (N) 设得特别大，然后体现在系数上就没有循环。

但是我们实现的 FFT 都是分治做的，要求多项式次数是二的幂。就没有办法做任意长度的循环卷积。

Hint.

普通的 FFT 也可以解决一部分循环卷积问题，只需要把多项式次数设得很大，先计算出 $alpha,eta$ 的线性卷积，再手动把 $x^t$ 的系数加到 $x^{tmod n}$ 上去。

只不过这样没有办法倒过来，知道 $alpha imeseta=gamma$ 的 $eta,gamma$ 求 $alpha$。

那就先不管 FFT，来推一下对多项式 (f) 求 DFT 的式子。记 (X_k) 表示 (f(e^{-frac{2pi k}Ni}))（也就是求 DFT 后的第 (k) 项），(f_k) 是 (f) 的 (x^k) 项的系数：

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

Hint.

$e^{frac{2pi} ni}$ 就是单位复根，也可以写作 $omega_n=cosfrac{2pi}n+sinfrac{2pi}ni$（是不是看起来熟悉一点了）。

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

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

于是

[X_k=e^{-frac{k^2pi}{N}}sum_{n=0}^{N-1}e^{frac{(k-n)^2pi}{N}}cdotig(e^{-frac{n^2pi}{N}}f_nig) ]

可以看出 (X_k) 后面的求和式是一个线性卷积，可以用 FFT。

值得注意的是这个卷积的下标 (k-n) 的范围是 ([1-N,N-1])，是可能有负数的。所以先给它平移 (N) 位。

逆变换也基本一样：

[X_k'=e^{frac{k^2pi}{N}}sum_{n=0}^{N-1}e^{-frac{(k-n)^2pi}{N}}cdotig(e^{frac{n^2pi}{N}}f_nig) ]

所以构造出这样两个数列：

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

卷积得到多项式 (C)，别忘了 (C) 是左移了 (N) 位的。

[X_k=e^{-frac{k^2pi}{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
",aryA[i].r);
	return 0;
}

---

THE END

Thanks for reading!

[egin{split} “ &你是沙鸥 你是滴漏\ &你是白昼 你是不朽\ &你是从来都没有 缄默的口\ &你是跌落深谷里 最后的柔\ &你是行走于人间的风\ &你是我遥不可及的梦 ”\ ——& ext{《你是我遥不可及的梦》 By 苍穹} end{split} ]

> Link 你是我遥不可及的梦-Bilibili