拆系数FFT学习笔记

拆系数FFT

当题目中取模的数不是NTT模数的时候，我们无法利用原根来进行快速数论变换，这个时候就要用到毛啸论文里提到的拆系数FFT。

大致思路

拆系数FFT实际上是将多项式卷积之后的值具体算出来，普通的FFT由于精度误差较大，当然无法胜任。

于是可以考虑将一个较大的数拆成两个较小的数，即将\(x\)拆成\(a \times t + b\)，这样相当于将\(x\)给\(t\)进制划分，在常数\(t\)取根号级别的时候，\(a,b\)的大小也是根号级别的，这样我们再用精度较高的浮点数就可以直接FFT来运算了。

例如：

\[ \begin{aligned} x&=a_1t+b_1\\ y&=a_2t+b_2\\ xy&=(a_1t+b_1)(a_2t + b_2)\\ &=a_1a_2t^2 + (a_1b_2+a_2b_1)t+a_2b_2 \end{aligned} \]

将\(a_1b_1,(a_1b_2+a_2b_1),a_2b_2\)对应算出来再各自乘上相应的系数即可。

这样需要做8次dft，发现中间有两个部分系数相同，对于这两个系数相同的部分我们可以直接点值操作之后再做一次dft即可，这样dft的次数就降到了7次了。

优化

DFT

上面的方法好像常数很大，很多时候需要用到两次DFT合并为一次的方法。

考虑一个多项式A的DFT本质上求的是\(A(\omega_n^k)=\sum_{j=0}^{n-1}A_{j}\times \omega_n^{kj}\)，现在我们需要同时求出A和B，也就是算出\(A(\omega_n^k),B(\omega_n^k)\)。

假设我们知道了\(F_p[k]=A(\omega_n^k)+iB(\omega_n^k),F_q[k]=A(\omega_n^k)-iB(\omega_n^k)\)，我们就可以通过\(F_p[k]\)和\(F_q[k]\)来直接求出\(A(\omega_n^k),B(\omega_n^k)\)了。

于是为了合并两次DFT，我们希望\(F_p\)和\(F_q\)之间能够存在某一种关系，使得求出了某一个之后可以快速求出另外一个。

于是我们来推导\(F_p\)和\(F_q\)的关系，为了方便下面的公式的书写，我们约定\(X=\frac{2\pi kj}{n}\)。

\[ \begin{aligned} F_q[k]&=A(\omega_n^k)-iB(\omega_n^k)\\ &=\sum_{j=0}^{n-1}A_j\omega_n^{kj}-iB_j\omega_n^{kj}\\ &=\sum_{j=0}^{n-1}(A_j-iB_j)(\cos X + i\sin X)\\ &=\sum_{j=0}^{n-1}(A_j\cos X + B_j\sin X)+i(A_j\sin X-B_j\cos X)\\ F_p[n-k]&=A(\omega_n^{-k})+iB(\omega_n^{-k})\\ &=\sum_{j=0}^{n-1}A_j\omega_n^{-kj}+iB_j\omega_n^{-kj}\\ &=\sum_{j=0}^{n-1}(A_j+iB_j)(\cos X -i\sin X)\\ &=\sum_{j=0}^{n-1}(A_j\cos X + B_j\sin X)-i(A_j\sin X-B_j\cos X)\\ \end{aligned} \]

我们发现\(F_p[n-k]\)和\(F_q[k]\)在本质上就是共轭的，于是通过\(F_p\)来计算\(F_q\)就非常方便了。

简单分析一下原因，由于\(\omega_{n}^{n-k}=\omega_n^{-k}\)，而\(\omega_{n}^{k}\)和\(\omega_n^{-k}\)实部和虚部的元素都是相同的，只是符号相反而已，\(F_p[k]\)和\(F_q[k]\)的系数刚好也是符号相反，也就导致\(F_p[k],F_p[n-k],F_q[k],F_q[n-k]\)四个项的它们本质的结果只是一些元素的线性组合，这样从\(F_p\)推出\(F_q\)也就不是那么难了。

IDFT

求出了点值之后我们可以快速将点值乘起来，但是要怎么IDFT回去？

对于两个点值序列\(A,B\)，我们仍然可以用刚才提到的方法，将每一项合并成\(A+iB\)，将新的多项式IDFT之后得到的系数的实部和虚部就对应了真正\(A\)和\(B\)了。

/*=======================================
 * Author : ylsoi
 * Time : 2019.3.14
 * Problem : luogu4245
 * E-mail : ylsoi@foxmail.com
 * ====================================*/
#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
#define double long double
typedef long long ll;

using namespace std;

void File(){
    freopen("luogu4245.in","r",stdin);
    freopen("luogu4245.out","w",stdout);
}

template<typename T>void read(T &_){
    _=0; T fl=1; char _c=getchar();
    for(;!isdigit(_c);_c=getchar())if(_c=='-')fl=-1;
    for(;isdigit(_c);_c=getchar())_=(_<<1)+(_<<3)+(_c^'0');
    _*=fl;
}

const int maxn=1e5+10;
const double pi=acos(-1);
int n,m,mod,a[maxn<<2],b[maxn<<2],c[maxn<<2];

namespace Poly{
    struct cp{
        double x,y;
        cp(double xx=0,double yy=0){x=xx,y=yy;}
        cp conj(){return cp(this->x,-this->y);}
        friend cp operator + (cp a,cp b){return (cp){a.x+b.x,a.y+b.y};}
        friend cp operator - (cp a,cp b){return (cp){a.x-b.x,a.y-b.y};}
        friend cp operator * (cp a,cp b){return (cp){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
        friend cp operator * (cp a,double b){return cp(a.x*b,a.y*b);}
        friend cp operator / (cp a,double b){return cp(a.x/b,a.y/b);}
    }om[maxn<<2],aa[maxn<<2],bb[maxn<<2],cc[maxn<<2],dd[maxn<<2];
    int lim,cnt,dn[maxn<<2];
    inline void dft(cp *A,int ty){
        if(ty==-1)reverse(A+1,A+lim);
        REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
        for(int len=1;len<lim;len<<=1){
            cp w=om[len<<1];
            for(int L=0;L<lim;L+=len<<1){
                cp wk=cp(1,0);
                REP(i,L,L+len-1){
                    cp u=A[i],v=A[i+len]*wk;
                    A[i]=u+v;
                    A[i+len]=u-v;
                    wk=wk*w;
                }
            }
        }
        if(ty==-1)REP(i,0,lim-1)A[i]=A[i]/lim;
    }
    void mtt(int *A,int *B,int *C,int la,int lb){
        REP(i,0,lim-1){
            dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
            aa[i]= i<=la ? cp(A[i]>>15,A[i]&32767) : cp(0,0);
            bb[i]= i<=lb ? cp(B[i]>>15,B[i]&32767) : cp(0,0);
        }
        dft(aa,1),dft(bb,1);
        REP(i,0,lim-1){
            int j=(lim-i)&(lim-1);
            cp a1=(aa[i]+aa[j].conj())*cp(0.5,0);
            cp a0=(aa[i]-aa[j].conj())*cp(0,-0.5);
            cp b1=(bb[i]+bb[j].conj())*cp(0.5,0);
            cp b0=(bb[i]-bb[j].conj())*cp(0,-0.5);
            cc[i]=a1*b1+a1*b0*cp(0,1);
            dd[i]=a0*b1+a0*b0*cp(0,1);
        }
        dft(cc,-1),dft(dd,-1);
        REP(i,0,lim-1){
            int a1=(ll)(cc[i].x+0.5)%mod;
            int b1=(ll)(cc[i].y+0.5)%mod;
            int c1=(ll)(dd[i].x+0.5)%mod;
            int d1=(ll)(dd[i].y+0.5)%mod;
            C[i]=(((ll)a1<<30)+((ll)(b1+c1)<<15)+(ll)d1)%mod;
        }
    }
}

int main(){
    File();

    read(n),read(m),read(mod);
    REP(i,0,n)read(a[i]);
    REP(i,0,m)read(b[i]);

    using namespace Poly;

    lim=1,cnt=0;
    while(lim<=n+m)lim<<=1,++cnt;
    if(!cnt)cnt=1;
    REP(i,0,lim-1)dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
    for(int i=lim;i;i>>=1)
        om[i]=cp(cos(2.0*pi/i),sin(2.0*pi/i));

    Poly::mtt(a,b,c,n,m);
    REP(i,0,n+m)printf("%d ",(c[i]+mod)%mod);

    return 0;
}