任意模数NTT(二)
这一版是mx发明的MTT
速度极快,精度基本有保证,在奇技淫巧无效时可以考虑这个东西...
(但是无论如何我都不想用真正的任意模数NTT,那种东西简直毒瘤而且常常数巨大...)
原理:并不关心
#include <cstdio> #include <cmath> #include <cstring> #include <cstdlib> #include <iostream> #include <algorithm> #include <queue> #include <stack> #define ld long double #define ll long long using namespace std; const ld pi=acos(-1.0); const int Siz=(1<<21)+5; const int siz=(1<<15)-1; struct cp { ld x,y; cp (){} cp (ld a,ld b):x(a),y(b){} 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); } }; int to[Siz]; ll p; int n,m; ll a[Siz],b[Siz]; ll ret[Siz] ; void FFT(cp *a,int len,int k) { for(int i=0;i<len;i++)if(i<to[i])swap(a[i],a[to[i]]); for(int i=1;i<len;i<<=1) { cp w0=cp(cos(pi/i),k*sin(pi/i)); for(int j=0;j<len;j+=(i<<1)) { cp w=cp(1,0); for(int o=0;o<i;o++,w=w*w0) { cp w1=a[j+o],w2=a[j+o+i]*w; a[j+o]=w1+w2,a[j+o+i]=w1-w2; } } } } cp A[Siz],B[Siz],C[Siz],D[Siz],E[Siz],F[Siz],Q[Siz],G[Siz],H[Siz]; inline void MTT(ll *f,ll *g,ll *re,int len) { int lim=1,l=0; while(lim<=len)lim<<=1,l++; for(int i=0;i<lim;i++)to[i]=((to[i>>1]>>1)|((i&1)<<(l-1))); for(int i=0;i<=len;i++)A[i]=cp((ld)(f[i]&siz),(ld)(f[i]>>15)); for(int i=0;i<=len;i++)B[i]=cp((ld)(g[i]&siz),(ld)(g[i]>>15)); FFT(A,lim,1),FFT(B,lim,1); for(int i=0;i<lim;i++) { int j=(lim-i)&(lim-1); C[j]=cp(0.5*(A[i].x+A[j].x),0.5*(A[i].y-A[j].y))*B[i]; D[j]=cp(0.5*(A[i].y+A[j].y),0.5*(A[j].x-A[i].x))*B[i]; } FFT(C,lim,1),FFT(D,lim,1); for(int i=0;i<lim;i++) { ll tempA=(ll)(C[i].x/lim+0.5)%p; ll tempa=(ll)(C[i].y/lim+0.5)%p; ll tempB=(ll)(D[i].x/lim+0.5)%p; ll tempb=(ll)(D[i].y/lim+0.5)%p; re[i]=(ll)(tempA+p+(ll)((ll)(tempa+tempB)<<15)%p+p+(ll)((ll)tempb<<30)%p+p)%p; } } int main() { scanf("%d%d%lld",&n,&m,&p); for(int i=0;i<=n;i++)scanf("%lld",&a[i]); for(int i=0;i<=m;i++)scanf("%lld",&b[i]); MTT(a,b,ret,n+m-1); for(int i=0;i<=n+m;i++)printf("%lld ",ret[i]); printf("\n"); return 0; }