三次变两次 FFT

我们发现：

\[(F(x)+iG(x))^2=F(x)^2-G(x)^2+2iF(x)G(x)
\]

也就是说，我们把 \(F(x)\) 作为实部，\(G(x)\) 作为虚部，那么它的平方的虚部的 \(1/2\) 就是 \(F(x)G(x)\)

可惜精度比较低。

四次 FFT 求任意模数多项式乘法

假设我们要求 \(M(x)\times N(x)\pmod{p}\)，因为如果我们直接 FFT 就会爆 double ，所以我们可以把 \(M(x)\) 拆成 \(kA(x)+B(x)\)，\(N(x)\) 拆成 \(kC(x)+D(x)\)，其中 \(k\approx \sqrt{p}\)，那么，你的值域就大约变为 \(p\times n\) 的了，但是你就需要 \(7\) 次 FFT 了。

我们假设有 \(Q(x)=A(x)+iB(x),E(x)=C(x)+iD(x)\)，设 \(q(x_1)\) 表示 \(Q(x)\) 在 \(x_1\) 的取值，\(A_j\) 表示 \(A(x)\) 第 \(j\) 项系数，那么我们就有：

\[q(\omega_n^x)=\sum_{j=0} A_j\omega_n^{xj}+iB_j\omega_n^{xj}
\]

\[=\sum_{j=0} A_j(\cos{\frac{2\pi xj}{n}}+i\sin{\frac{2\pi xj}{m}})+B_j(i\cos{\frac{2\pi xj}{n}}-\sin{\frac{2\pi xj}{n}})
\]

\[q(\omega_n^{-x})=\sum_{j=0} A_j(\cos{\frac{2\pi xj}{n}}-i\sin{\frac{2\pi xj}{n}})+B_j(i\cos{\frac{2\pi xj}{n}}+\sin{\frac{2\pi xj}{m}})
\]

\[a(\omega_n^{x})=\sum_{j=0} A_j(\cos{\frac{2\pi xj}{n}}+i\sin{\frac{2\pi xj}{n}})
\]

\[b(\omega_n^{x})=\sum_{j=0} B_j(\cos{\frac{2\pi xj}{n}}+i\sin{\frac{2\pi xj}{n}})
\]

我们假设 \(q(x).r\) 表示它的实部，\(q(x).f\) 表示它的虚部，那么我们就可以得到：

\[a(\omega_n^{x}).r=\frac{q(\omega_n^x).r+q(\omega_n^{-x}).r}{2},a(\omega_n^{x}).f=\frac{q(\omega_n^x).f-q(\omega_n^{-x}).f}{2}
\]

\[b(\omega_n^{x}).r=\frac{q(\omega_n^{-x}).f+q(\omega_n^x).f}{2},a(\omega_n^{x}).f=\frac{q(\omega_n^{-x}).r-q(\omega_n^x).r}{2}
\]

然后，我们求出 \(A(x)E(x)\) 和 \(B(x)E(x)\)，我们就可以得到 \(A(x)C(x),A(x)D(x),B(x)C(x),B(x)D(x)\)，然后算就好了。

暂时没有懂 \(3.5\) 次 FFT 的做法，所以不写了。

Code

#include <bits/stdc++.h>

using namespace std;

#define double long double

#define Int register int

#define MAXN 270005

template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c == '-') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}

template <typename T,typename ... Args> inline void read (T &t,Args&... args){read (t);read (args...);}

template <typename T> inline void write (T x){if (x < 0){x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}

struct Complex{

	double x,y;

	Complex(){}

	Complex (double _x,double _y){x = _x,y = _y;}

	Complex operator / (const double &p)const{return Complex{x / p,y / p};}

	Complex operator + (const Complex &p)const{return Complex{x + p.x,y + p.y};}

	Complex operator - (const Complex &p)const{return Complex{x - p.x,y - p.y};}

	Complex operator * (const Complex &p)const{return Complex{x * p.x - y * p.y,x * p.y + p.x * y};}

};

#define pi (double)acos(-1)

int l,lim,rev[MAXN];

void fft (Complex *a,int type){

	for (Int i = 0;i < lim;++ i) if (i < rev[i]) swap (a[i],a[rev[i]]);

	for (Int i = 1;i < lim;i <<= 1){

		Complex Wn(cos(pi / i),type * sin(pi / i));

		for (Int j = 0,r = i << 1;j < lim;j += r){

			Complex w(1,0);

			for (Int k = 0;k < i;++ k,w = w * Wn){

				Complex x = a[j + k],y = w * a[i + j + k];

				a[j + k] = x + y,a[i + j + k] = x - y;

			}

		}

	}

	if (type == -1) for (Int i = 0;i < lim;++ i) a[i] = a[i] / lim;

}

int n,m,mod,a[MAXN],b[MAXN],ans[MAXN];

Complex Q[MAXN],E[MAXN],C[MAXN],D[MAXN];

#define ll long long

signed main(){

	read (n,m,mod),lim = 1;

	while (lim < n + m) lim <<= 1,++ l;

	for (Int i = 0;i < lim;++ i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << l - 1);int up = (1 << 15) - 1;

	for (Int i = 0;i <= n;++ i) read (a[i]),Q[i] = Complex (a[i] >> 15,a[i] & up);

	for (Int i = 0;i <= m;++ i) read (b[i]),E[i] = Complex (b[i] >> 15,b[i] & up);

	fft (Q,1),fft (E,1);

	for (Int i = 0;i < lim;++ i){

		int re = (lim - 1) & (lim - i);

		C[i] = Complex((Q[i].x + Q[re].x) / 2,(Q[i].y - Q[re].y) / 2) * E[i];

		D[i] = Complex((Q[re].y + Q[i].y) / 2,(Q[re].x - Q[i].x) / 2) * E[i];

	}

	fft (C,-1),fft (D,-1);

	for (Int i = 0;i < lim;++ i){

		ll v1 = (ll)(C[i].x + 0.5) % mod,v2 = (ll)(C[i].y + D[i].x + 0.5) % mod,v3 = (ll)(D[i].y + 0.5) % mod;

		ans[i] = ((v1 << 30) + (v2 << 15) + v3) % mod;

	}

	for (Int i = 0;i <= n + m;++ i) write ((ans[i] % mod + mod) % mod),putchar (' ');

	putchar ('\n');

	return 0;

}

巴特西

FFT 的一些技巧

三次变两次 FFT

四次 FFT 求任意模数多项式乘法

Code

最新文章

热门文章