任意模数NTT

Orz myy！

特殊模数

如果模数很特殊，$(2^k+1)|(p-1)$，$(p-1)>$DFT的长度，那么可以求出p的原根g，用$g$代替单位复数根$w_n$。
不会原根？请点这里

一般模数

那就不能用原根代替单位复数根这种方法了。
只能用实数。
但是如果模数很大，比如$1000000007$，用double的话很可能爆精度。
目标：解决精度问题。
1.用中国剩余定理：
对每个模数分别做3次DFT。一般998244353 , 1004535809 ，9985661441
然后将三次的计算结果乘起来，对要求的模数取模即可。
可行条件：三个模数的乘积大于$np^2$
优点：降低错误率。
缺点：做9次DFT，太慢了。常数巨大。

2.7遍DFT
将两个多项式拆一下：
$A_i=a_i*\sqrt{P}+b_i$，$B_i=c_i*\sqrt{P}+d_i$
通过乘法分配律，则
$A_i*B_i=$$P*(a_i*c_i)$$+$$\sqrt{P}(a_i*d_i+b_i*c_i)$$+$$b_i*d_i$
看得出来，做4次点值，分别求出$a_i*c_i,b_i*d_i,a_i*d_i,b_i*c_i$
再做3次插值，分别对三种颜色部分求插值即可。
然而还不优。

重头戏

myy的算法中(DFT+IDFT)的次数只有4次。
看看这究竟是什么奇技淫巧吧。
一个多项式的系数表示是这样的：$A(x)=\Sigma_{i=0}^{n-1}a_ix^i$
在做DFT的时候，$a_i$的虚部一开始都为0。很多人以为这没用。实际上并不然。
这个算法巧妙地利用了虚部的空间。
如果$a_i=k_i+i*b_i$，又会怎么样？
考虑1次DFT搞定原来2次DFT得出的多项式A,B.

凑式子

问题来了，现在既要考虑充分利用虚部空间，又要考虑凑式子。
设$C(x)=\Sigma_{i=0}^{n-1}a_ix^i=\Sigma_{i=0}^{n-1}(k_i+i*b_i)x^i$
由于$x$在FFT里全部都是n次单位复数根，并且在合并$A_0$和$A_1$两个多项式的时候，求$y_k$和$y_{k+n/2}$的式子差不多。
$y_k=A_0(w_n^{2k})+w_n^k*A_1(w_n^{2k})$
$y_{k+n/2}=A_0(w_n^{2k+n})-w_n^k*A_1(w_n^{2k+n})$
某人就想到了共轭复数？这个地方有点难理解。
共轭复数：实部相同，虚部互为相反数的两个复数。
令$a'_i=k_i+i*b_i$
则$D(x)=\Sigma_{i=0}^{n-1}a'_ix^i=\Sigma_{i=0}^{n-1}(k_i-i*b_i)x^i$
则$A(x)$和$B(x)$可以通过简单的加减乘除凑出来。
如果用1次DFT就能求出A,B,C,D，那么
C的某一项和D的某一项是否有直接关系呢？
答案是：有。
求证：$C(w_n^k)$与$D(w_n^{n-k})$互为共轭复数。
证明：
$C(w_n^k)=\Sigma_{i=0}^{n-1}(k_i+i*b_i)w_n^k=\Sigma_{i=0}^{n-1}(k_i+i*b_i)(cos(\theta)+i*sin(\theta))=\Sigma_{i=0}^{n-1}k_i*cos(\theta)-b_i*sin(\theta)+i*(k_i*sin(\theta)+b_i*cos(\theta))$
$D(w_n^k)=\Sigma_{i=0}^{n-1}(k_i-i*b_i)w_n^{-k}=\Sigma_{i=0}^{n-1}(k_i-i*b_i)(cos(\theta)-i*sin(\theta))=\Sigma_{i=0}^{n-1}k_i*cos(\theta)-b_i*sin(\theta)-i*(k_i*sin(\theta)+b_i*cos(\theta))$
很容易看得出结论了。
证毕。
那这样的话，对C求DFT，就相当于对A和B求DFT了。

程序实现

模板题：洛谷4245。
注意一点：传统写法：$w=w*w_n^k$，乘法运算的次数过多，会使得精度有大问题，建议将n个n次单位复数根用数组存储。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define N 262210
#define LL long long
#define DB double
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
const DB pi=acos(-1);
struct Cpx{
    DB r,i;
    Cpx(){}
    Cpx(DB R,DB I){r=R,i=I;}
    Cpx operator +(const Cpx&x)const{return Cpx(r+x.r,i+x.i);}
    Cpx operator -(const Cpx &x)const{return Cpx(r-x.r,i-x.i);}
    Cpx operator *(const Cpx &x)const{return Cpx(r*x.r-i*x.i,r*x.i+i*x.r);}
}f1[N],f2[N],w[N];
int i,j,k,L,n,m,ans;
int Rev[N];
Cpx df1[N],df3[N],df4[N];
LL x2,mo;
void FFT(Cpx *h,int delta){
    int i,k,m,Half;
    fo(i,0,L-1)if(i<Rev[i])swap(h[i],h[Rev[i]]);
    for(m=2;m<=L;m<<=1){
        Half=m>>1;
        for(i=0;i<L;i+=m){
            int v=0;
            fo(k,i,i+Half-1){
                Cpx u=h[k],t=w[v]*h[k+Half];
                h[k]=u+t;
                h[k+Half]=u-t;
                v=(v+delta*(L/m)+L)&(L-1);
                //w=w*wn;//精度差很多 乘法运算次数太多 
            }
        }
    }
    if(delta==-1)fo(i,0,L-1)h[i].r/=L;
}

void Mulpoly(Cpx *x,Cpx *y,Cpx *z){
    int i,j;
    static Cpx a[N],b[N],c[N];
    fo(i,0,L-1)a[i]=Cpx((LL)x[i].r&32767,(LL)x[i].r>>15);
    fo(i,0,L-1)b[i]=Cpx((LL)y[i].r&32767,(LL)y[i].r>>15);
    FFT(a,1);
    FFT(b,1);
    fo(i,0,L-1){
        static Cpx d1,d2,d3,d4;
        j=(L-i)&(L-1);
        d1=(a[i]+Cpx(a[j].r,-a[j].i))*Cpx(0.5,0);
        d2=(a[i]-Cpx(a[j].r,-a[j].i))*Cpx(0,-0.5);
        d3=(b[i]+Cpx(b[j].r,-b[j].i))*Cpx(0.5,0);
        d4=(b[i]-Cpx(b[j].r,-b[j].i))*Cpx(0,-0.5);
        df1[i]=d1*d3;
        df3[i]=d1*d4+d2*d3;
        df4[i]=d2*d4;
    }
    fo(i,0,L-1)a[i]=df1[i]+df4[i]*Cpx(0,1);
    //巧妙利用虚部空间，一次IDFT完成红色和蓝色的插值。
    fo(i,0,L-1)b[i]=df3[i];
    FFT(a,-1);FFT(b,-1);
    fo(i,0,L-1){
        static LL d1,d3,d4;
        d1=(LL)(a[i].r+0.5)%mo;
        d3=(LL)(b[i].r+0.5)%mo;
        d4=(LL)(a[i].i/L+0.5)%mo;
        z[i].r=(((LL)d4<<30)%mo+((LL)d3<<15)%mo+d1)%mo;
    }
}
int main(){
    scanf("%d%d%d",&n,&m,&mo);
    for(L=1;L<n+m+1;L<<=1);
    fo(i,0,L-1)Rev[i]=Rev[i>>1]>>1|(i&1)*(L>>1);
    fo(i,0,n){
        scanf("%lld",&x2);
        f1[i]=Cpx(x2,0);
    }
    fo(i,0,m){
        scanf("%lld",&x2);
        f2[i]=Cpx(x2,0);
    }
    fo(i,0,L-1)w[i]=Cpx(cos(2*pi*i/L),sin(2*pi*i/L));
    Mulpoly(f1,f2,f1);
    fo(i,0,n+m)printf("%lld ",(LL)f1[i].r);
    return 0;
}