ZigZagK的博客
[拆系数FFT]洛谷4245【任意模数NTT】题解
2018年5月20日 19:01
洛谷
查看标签

题目概述

求 $A(x)B(x)$ ,系数对任意模数 $p$ 取模。

解题报告

任意模数NTT好像需要三模数NTT搞。然而我不会NTT……所以我来愉快的讲一波任意模数FFT(雾)。

其实洛谷里panda_2134的题解已经说的很详细了……就是把一个数拆成 $ax+b$ ,其中 $x$ 给 $\sqrt{Max\{A,B\}}$ 左右,然后两个拆过的多项式相乘,这样FFT运算过程就不会爆 $int$ 了,最后再模 $p$ 即可。

示例程序

我FFT常数好大啊QAQ。

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
#define fr first
#define sc second
#define mp make_pair
typedef long long LL;typedef long double DB;typedef pair<DB,DB> CP;
const int maxn=100000,maxm=262144;const DB pi=acos(-1);

int n,m,K,REV[maxm+5],a[maxn+5],b[maxn+5],MOD,ans[maxm+5];DB Sin[maxm+5],Cos[maxm+5];
CP A[maxm+5],B[maxm+5],C[maxm+5],D[maxm+5],E[maxm+5],F[maxm+5],G[maxm+5],H[maxm+5];
inline CP operator + (const CP &a,const CP &b) {return mp(a.fr+b.fr,a.sc+b.sc);}
inline CP operator - (const CP &a,const CP &b) {return mp(a.fr-b.fr,a.sc-b.sc);}
inline CP operator * (const CP &a,const CP &b) {return mp(a.fr*b.fr-a.sc*b.sc,a.fr*b.sc+a.sc*b.fr);}

inline void FFT(CP *a,int n,int f){
    for (int i=0;i<n;i++) if (i<REV[i]) swap(a[i],a[REV[i]]);
    for (int k=1;k<n;k<<=1){
        CP wn=mp(Cos[k],f>0?Sin[k]:-Sin[k]),w=mp(1,0),x,y;
        for (int i=0;i<n;i+=(k<<1),w=mp(1,0))
            for (int j=0;j<k;j++,w=w*wn)
                x=a[i+j],y=a[i+j+k]*w,a[i+j]=x+y,a[i+j+k]=x-y;
    }
    if (f<0) for (int i=0;i<n;i++) a[i].fr/=n;
}
inline void Make(int n){
    for (int k=1;k<n;k<<=1) Cos[k]=cos(pi/k),Sin[k]=sin(pi/k);
    for (int i=1;i<n;i++) REV[i]=(REV[i>>1]>>1)|((i&1)*(n>>1));
}
inline void AMOD(int &x,int tem) {if ((x+=tem)>=MOD) x-=MOD;}
inline void MTT(int *a,int *b){
    for (int i=1;i<=n;i++) A[i-1]=mp(a[i]>>15,0),B[i-1]=mp(a[i]&((1<<15)-1),0);
    for (int i=1;i<=m;i++) C[i-1]=mp(b[i]>>15,0),D[i-1]=mp(b[i]&((1<<15)-1),0);
        for (K=1;K<=n-1+m-1;K<<=1);Make(K);FFT(A,K,1);FFT(B,K,1);FFT(C,K,1);FFT(D,K,1);
    for (int i=0;i<K;i++) E[i]=A[i]*C[i],F[i]=A[i]*D[i],G[i]=B[i]*C[i],H[i]=B[i]*D[i];
    FFT(E,K,-1);FFT(F,K,-1);FFT(G,K,-1);FFT(H,K,-1);
    for (int i=0;i<K;i++){
        AMOD(ans[i],(((LL)(E[i].fr+0.5)%MOD)<<30)%MOD);AMOD(ans[i],(((LL)(F[i].fr+0.5)%MOD)<<15)%MOD);
        AMOD(ans[i],(((LL)(G[i].fr+0.5)%MOD)<<15)%MOD);AMOD(ans[i],(LL)(H[i].fr+0.5)%MOD);
    }
}
int main(){
    freopen("program.in","r",stdin);
    freopen("program.out","w",stdout);
    scanf("%d%d%d",&n,&m,&MOD);n++;m++;
    for (int i=1;i<=n;i++) scanf("%d",&a[i]);
    for (int i=1;i<=m;i++) scanf("%d",&b[i]);
    MTT(a,b);printf("%d",ans[0]);
    for (int i=1;i<n+m-1;i++) printf(" %d",ans[i]);
    return puts(""),0;
}
版权声明:本博客所有文章除特别声明外,均采用 CC BY 4.0 CN协议 许可协议。转载请注明出处!