MTT 学习笔记

很久以前就听说了这东西，一直没空学。最近重学多项式，就重新搞了一下。

MTT 主要解决的是任意模数（或者说是没有模数）的多项式乘法，可以用于应对专门恶心人的毒瘤题。

首先，假设多项式次数 $10^5$,值域 $10^9$，那么朴素 FFT 光答案就高达 $10^{23}$，就算你不损失精度也转不成 long long，况且中间值比这还要高几个量级。

MTT 采用拆系数的方法，即对于多项式 $A(x)$，拆成两个多项式 $M{A}_0(x)+A_1(x)$ 来减小精度压力。其中 $M$ 为一较大整数，一般取 $\sqrt{MOD}$。要注意模数 $2e9$ 级别的时候不要偷懒用位运算，$2^{15}$ 或 $2^{16}$ 都容易炸，并且复杂度瓶颈不在这，不用卡这点常。

然后你只需要求 $(M{A}_0(x)+A_1(x))(M{B}_0(x)+B_1(x))=M^2{A}_0(x)B_0(x)+M(A_0(x)B_1(x)+A_1(x)B_0(x))+A_1(x)B_1(x)$ 就可以了，这样答案是 $10^{14}$ 级别，开 long double 完全不虚。

这样一共需要 $7$ 次 DFT，一般问题不大。接下来是利用高超的玄学技巧压到 $4$ 次 DFT 的方法。

前置知识：普通 FFT 三次变两次

对于一般的实系数多项式乘法 $A(x)\cdot B(x)$，有一种优化方法是构造一个多项式 $A(x)+i\cdot B(x)$，因为是实系数，所以直接把 $B(x)$ 的系数挂到虚部上就可以了。

对它做一次 DFT 然后自己乘自己，得到的是 $A^2(x)-B^2(x)+2i\cdot A(x)B(x)$

这样算出来的虚部除以 $2$ 就是答案。

该算法将大量用到该思想。

DFT

考虑我们对两个实系数多项式 $A(x),B(x)$ 做 DFT，因为是实系数，我们把虚部在那里空着显得很浪费，尝试把这个优化成一次 DFT。

构造多项式 $A(x)+iB(x)$，花费一次 DFT 得到 $\mathrm{DFT}(A(x)+i(B(x)))$。因为 DFT 的本质是点值，所以这个也等于 $\mathrm{DFT}(A(x))+i\mathrm{DFT}(B(x))$。

这个显然不是实部和虚部的关系了。现在我们需要不用 DFT 把这两个分离开。

考虑共轭虚多项式 $A(x)-iB(x)$，如果它的 DFT 可以快速得到，我们就可以简单地求出 $A(x)$ 和 $B(x)$ 的 DFT。

然后又是愉快的推式子时间。

设两个多项式的系数为 $a_i,b_i$，代入一个单位根 ${\omega }_n^k$

$$\mathrm{DFT}(A(x)+iB(x))=A({\omega }_n^k)+iB({\omega }_n^k)$$

$$=\sum _{j=0}^{n-1}{a}_j{\omega }_n^{kj}+i{b}_j{\omega }_n^{kj}$$

推不动了，考虑硬算。令 ${\omega }_n^{kj}=x_j+i{y}_j$

$$\mathrm{DFT}(A(x)+iB(x))=\sum _{j=0}^{n-1}{a}_j(x_j+i{y}_j)+i{b}_j(x_j+i{y}_j)$$

$$=\sum _{j=0}^{n-1}(a_j{x}_j-b_j{y}_j)+i(a_j{y}_j+b_j{x}_j)$$

然后这不是三角函数公式吗？

所以 DFT 的某一点值的本质就是两个多项式对应项组成的向量 $(a_j,b_j)$ 旋转同一角度后的向量和。

对于另一边 $\mathrm{DFT}(A(x)-iB(x))$，其对应的向量是 $(a_j,-b_j)$ ，也就是关于 $x$ 轴对称。那么我们把上面那个角度倒着旋转，得到的向量和也是对称的。

翻译成人话， DFT 后翻转再共轭就得到了共轭虚多项式的 DFT。

注意是长度为 $n$ 的翻转，所以 $0$ 的位置需要特判。

这样，我们用两次 DFT 得到了 $A_0(x),A_1(x),B_0(x),B_1(x)$ 的 DFT。

然后我们就可以求出 $A_0(x)B_0(x),A_0(x)B_1(x),A_1(x)B_0(x),A_1(x)B_1(x)$ 的点值了。

IDFT

还是用一样的方法抱团变换。

考虑构造多项式 $P(x)=A_0(x)B_0(x)+i{A}_0(x)B_1(x)$，注意是 IDFT 之后的，也就是这里设的是我们要求的答案。

因为我们知道 $A_0(x)B_0(x),A_0(x)B_1(x)$ 在所有单位根上的点值，所以可以直接求出 $P(x)$ 的点值，然后做 IDFT，实部虚部就是 $A_0(x)B_0(x)$ 和 $A_0(x)B_1(x)$。

另外两个同理，需要 $2$ 次 IDFT。

总共做了 $4$ 次 DFT，比三模数快了 $8$ 倍。

注意输出的时候答案已经很大了，要先模了再乘 $M$。

#include <iostream>
#include <cstring>
#include <cctype>
#include <cstdio>
#include <cmath>
#define MAXN 262144+5
#define double long double
using namespace std;
inline int read()
{
	int ans=0;
	char c=getchar();
	while (!isdigit(c)) c=getchar();
	while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
	return ans;
}
typedef long long ll;
const double Pi=acos(-1.0);
int p;
struct complex
{
	double x,y;
	inline complex(const double& x=0,const double& y=0):x(x),y(y){}	
};
inline complex operator +(const complex& a,const complex& b){return complex(a.x+b.x,a.y+b.y);}
inline complex operator -(const complex& a,const complex& b){return complex(a.x-b.x,a.y-b.y);}
inline complex operator *(const complex& a,const complex& b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
inline complex rev(const complex& a){return complex(a.x,-a.y);}
complex rt[2][20];
int l,lim,r[MAXN];
inline void init(){lim=1<<l;for (int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));}
inline void fft(complex* a,int type)
{
	for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);
	for (int L=0;L<l;L++)
	{
		int mid=1<<L,len=mid<<1;
		complex Wn=rt[type][L+1];
		for (int s=0;s<lim;s+=len)
		{
			complex w(1,0);
			for (int k=0;k<mid;k++,w=w*Wn)
			{
				complex x=a[s+k],y=w*a[s+mid+k];
				a[s+k]=x+y,a[s+mid+k]=x-y;
			}
		}
	}
	if (type) for (int i=0;i<lim;i++) a[i].x/=lim,a[i].y/=lim;
}
complex A[MAXN],B[MAXN],C[MAXN],D[MAXN];
int main()
{
	freopen("test.in","r",stdin);
	freopen("test.out","w",stdout);
	for (int i=0;i<20;i++)
	{
		double a=2*Pi/(1<<i);
		rt[0][i]=rev(rt[1][i]=complex(cos(a),sin(a)));
	}
	int n,m;
	n=read(),m=read(),p=read();
	for (l=0;(1<<l)<=n+m;l++);
	init();
	for (int i=0;i<=n;i++)
	{
		int x=read();
		A[i]=complex(x>>15,x&((1<<15)-1));
	}
	for (int i=0;i<=m;i++)
	{
		int x=read();
		C[i]=complex(x>>15,x&((1<<15)-1));
	}
	fft(A,0),fft(C,0);
	B[0]=rev(A[0]),D[0]=rev(C[0]);
	for (int i=1;i<lim;i++) B[i]=rev(A[lim-i]),D[i]=rev(C[lim-i]);
	for (int i=0;i<lim;i++) 
	{
		complex a=A[i],b=B[i],c=C[i],d=D[i];
		A[i]=(a+b)*complex(0.5,0),B[i]=(a-b)*complex(0,-0.5);
		C[i]=(c+d)*complex(0.5,0),D[i]=(c-d)*complex(0,-0.5);
	}
	for (int i=0;i<lim;i++)
	{
		complex a=A[i],b=B[i],c=C[i],d=D[i];
		A[i]=a*c+a*d*complex(0,1);
		B[i]=b*c+b*d*complex(0,1);
	}
	fft(A,1),fft(B,1);
	for (int i=0;i<=n+m;i++)
	{
		ll ans=0;
		ll a=A[i].x+0.5,b=A[i].y+B[i].x+0.5,c=B[i].y+0.5;
		ans=(ans+((a%p)<<30))%p;
		ans=(ans+((b%p)<<15))%p;
		ans=(ans+c)%p;
		printf("%lld ",ans);
	}
	return 0;
}