多项式乘法卡常技巧与任意模数NTT

一.普通多项式乘法.

作为最基础的多项式运算，这里不再介绍，具体参见快速傅里叶变换FFT与快速数论变换NTT入门.

这里的多项式乘法卡常技巧不会过于深入，且只适用于FFT.


二.多项式乘法优化1.

这是一种非常简单的多项式乘法常数优化，对于两个系数为实数的多项式$A(x),B(x)$，我们写成一个多项式$P(x)$：
$$P(x)=A(x)+iB(x)$$

然后将$P(x)$平方得到：
$$\begin{gathered} P^2(x)=(A(x)+iB(x){)}^2 \\ =A^2(x)-B^2(x)+2iA(x)B(x) \end{gathered}$$

我们发现只需要取出$P(x)$的虚部除$2$就是两个多项式卷积的值了，这样我们需要计算DFT的次数从$3$次变成了$2$次.

代码如下：

comp ta[N+9];

void Poly_mul(int *a,int *b,int n){
  Get_len(n);
  for (int i=0;i<len;++i) ta[i]=comp(a[i],b[i]),b[i]=0;
  FFT(ta,len,0);
  for (int i=0;i<len;++i) ta[i]=ta[i]*ta[i]*0.5;
  FFT(ta,len,1);
  for (int i=0;i<len;++i) a[i]=(int)(ta[i].y+0.5),ta[i]=comp();
}

但是这种算法只能直接得到乘积，并不能得到每一个多项式对应的点值，这导致它很不实用.


三.多项式乘法优化2.

对于多项式$A(x),B(x)$，我们设：
$$\begin{gathered} P(x)=A(x)+iB(x) \\ Q(x)=A(x)-iB(x) \end{gathered}$$

设$F_P(k)$表示对多项式$P(x)$代入${\omega }_n^k$后的点值，那么有：
$$\begin{gathered} F_A(k)=\frac{F_P(k)+F_Q(k)}{2} \\ {F}_B(k)=-i\frac{F_P(k)-F_Q(k)}{2} \end{gathered}$$

现在我们把$A(x),B(x)$两个毫不相关的多项式的点值求解转化为了$P(x),Q(x)$这两个看起来非常有关的多项式点值求解，这意味着有可能只要求出$P(x)$的点值就能直接得到$Q(x)$的点值.

来试试看，设$X=\frac{2\pi jk}{n}$，用$conj(z)$表示$z$的共轭复数，那么有：
$$\begin{gathered} F_Q(k)=\sum _{j=0}^{n-1}{\omega }_n^{jk}(a_j-i{b}_j) \\ =\sum _{j=0}^{n-1}(\cos X+i\sin X)(a_j-i{b}_j) \\ =\sum _{j=0}^{n-1}((a_j\cos X+b_j\sin X)+i(a_j\sin X-b_j\cos X)) \\ =conj\left(\sum _{j=0}^{n-1}((a_j\cos X+b_j\sin X)-i(a_j\sin X-b_j\cos X))\right) \\ =conj\left(\sum _{j=0}^{n-1}((a_j\cos (-X)-b_j\sin (-X))+i(a_j\sin (-X)+b_j\cos (-X)))\right) \\ =conj\left(\sum _{j=0}^{n-1}(a_j+i{b}_j)(\cos (-X)+i\sin (-X))\right) \\ =conj\left(\sum _{j=0}^{n-1}{\omega }_n^{-jk}(a_j-i{b}_j)\right) \\ =conj(F_P(n-k)) \end{gathered}$$

突然发现真的有关系，那么我们只需要做$1$次FFT就可以得到$2$个多项式的DFT了，仍然是从$3$次DFT变为$2$次DFT.

代码如下：

comp ta[N+9],tb[N+9];

void Poly_mul(int *a,int *b,int n){
  Get_len(n);
  for (int i=0;i<len;++i) ta[i]=comp(a[i],b[i]),b[i]=0;
  FFT(ta,len,0);
  tb[0]=comp(ta[0].x,-ta[0].y);
  for (int i=1;i<len;++i) tb[len-i]=comp(ta[i].x,-ta[i].y);
  for (int i=0;i<len;++i) ta[i]=comp(0,0.25)*(ta[i]+tb[i])*(tb[i]-ta[i]),tb[i]=comp();
  FFT(ta,len,1);
  for (int i=0;i<len;++i) a[i]=(int)(ta[i].x+0.5),ta[i]=comp();
}

四.三模数NTT求解任意模数NTT.

对于一个模数$p$，若要直接用一般NTT求在模$p$意义下的多项式卷积，由于原根与分治过程需要的特性，模数必须满足可以被分解为$p=2^{k}c+1$，其中$k$必须足够大，到达$2^k$超过多项式项数的程度.

但是若$p$并不满足这个性质该如何解决问题？

一般来说，$p$在$10^9$级别左右，$n$在$10^5$级别左右，所以最终答案序列的每一位会达到$10^{23}$级别大小，显然是不太可以用C++中内置的整数和实数类型操作的（__int128是C++11的）.

一种方法是找到三个$10^9$级别的NTT模数$A,B,C$，分别求出乘积的每一位系数在模$A,B,C$意义下的值，设为$x_1,x_2,x_3$，那么对于真正的系数$x$有：
$$\begin{gathered} x\equiv x_1(modA) \\ x\equiv x_2(modB) \\ x\equiv x_3(modC) \\ x<ABC \end{gathered}$$

先合并$x_1,x_2$，有：
$$\begin{gathered} x_1+k_{1}A\equiv x_2(modB) \\ {k}_1\equiv \frac{x_2-x_1}{A}(modB) \end{gathered}$$

于是我们现在求得$k_1$，并且有形式$x\equiv x_1+k_{1}A(modAB)$，记$x_4=x_1+k_{1}A$，继续合并$x_4,x_3$：
$$\begin{gathered} x_4+k_{4}AB\equiv x_3(modC) \\ {k}_4\equiv \frac{x_3-x_4}{AB}(modC) \end{gathered}$$

那么有$x\equiv x_4+k_{4}AB(modABC)$，由于$x<ABC$，所以$x=x_4+k_{4}AB=x_1+k_{1}A+k_{4}AB$.

直接在模$p$意义下做这个式子即可.

由于这个方法并不是非常好用（至少我这么觉得），所以就不贴代码了.


五.拆位FFT求解任意模数NTT.

一般情况下，double存$10^{18}$级别的数还是没什么问题的，所以我们可以尝试拆位分别做FFT.

我们将多项式$A(x)$拆成$A(x)=2^{15}{A}_0(x)+A_1(x)$的形式，那么多项式$A_0(x),A_1(x)$的系数均在一个相对原来小了不少的范围内.

对$A(x)$和$B(x)$这么处理后，我们写出最后乘积的式子：
$$A(x)B(x)=2^{30}{A}_0(x)B_0(x)+2^{15}(A_0(x)B_1(x)+A_1(x)B_0(x))+A_1(x)B_1(x)$$

接下来的问题就是分别计算这$4$个卷积了，看起来这$4$个卷积依次实现需要$12$次DFT啊…

但是我们不需要依次实现，事实上，我们只要分别跑出$4$个多项式的DFT，计算出每一个卷积后再IDFT回来就行了，其中系数为$2^{15}$的两个卷积可以直接加好再DFT回来，这样只需要$7$次DFT.

然而看起来这和上面的三模数NTT的次数差不了多少啊，而且FFT还有复数运算的常数因子在…

但是别忘了，FFT是可以进行两两合并的！利用多项式乘法优化2中的技巧，DFT的次数降成了$4$次，常数远远小于三模数NTT（实测是一半）！

代码如下：

comp ta[N+9],tb[N+9],tc[N+9],te[N+9];

void Split_fft(comp *a,comp *b){
  b[0]=comp(a[0].x,-a[0].y);
  for (int i=1;i<len;++i) b[i]=comp(a[len-i].x,-a[len-i].y);
  for (int i=0;i<len;++i){
	comp x=a[i],y=b[i];
	a[i]=(x+y)*0.5;b[i]=comp(0,0.5)*(y-x);
  }
}

void Poly_mul(int *a,int *b,int n){
  Get_len(n);
  for (int i=0;i<len;++i){
	ta[i]=comp(a[i]>>15,a[i]&(1<<15)-1);
	tc[i]=comp(b[i]>>15,b[i]&(1<<15)-1);
  }
  FFT(ta,len,0);
  FFT(tc,len,0);
  Split_fft(ta,tb);
  Split_fft(tc,te);
  for (int i=0;i<len;++i){
	ta[i]=ta[i]*tc[i]+comp(0,1)*ta[i]*te[i];
	tb[i]=tb[i]*tc[i]+comp(0,1)*tb[i]*te[i];
  }
  FFT(ta,len,1);
  FFT(tb,len,1);
  for (int i=0;i<len;++i){
	LL x=(LL)(ta[i].x+0.5)%mod,y=(LL)(ta[i].y+tb[i].x+0.5)%mod,z=(LL)(tb[i].y+0.5)%mod;
	a[i]=((x<<30)+(y<<15)+z)%mod;
	ta[i]=tb[i]=comp();
  }
}