COGS 2294. [HZOI 2015] 释迦
额，其实就是裸的三模数NTT，上一篇已经说过了
哦，还有一个就是对乘起来炸long long的数取模，用long double之类的搞一下就好，精度什么的，，（看出题人心情？？）
1 #include<cstdio> 2 #include<cstdlib> 3 #include<algorithm> 4 #include<cstring> 5 #define LL long long 6 #define N 300005 7 using namespace std; 8 inline int ra() 9 { 10 int x=0; char ch=getchar(); 11 while (ch<'0' || ch>'9') ch=getchar(); 12 while (ch>='0' && ch<='9') {x=x*10+ch-'0'; ch=getchar();} 13 return x; 14 } 15 16 const int mod=23333333; 17 const int M[]={998244353,1004535809,469762049}; 18 const int G[]={3,3,3}; 19 const LL _M=(LL)M[0]*M[1]; 20 21 LL ksm(LL a, int b, int P) 22 { 23 LL sum=1; 24 for (;b;b>>=1,a=a*a%P) 25 if (b&1) sum=sum*a%P; 26 return sum; 27 } 28 /*LL mul(LL a, LL b, LL p) 29 { 30 a%=p; b%=p; 31 return ((a*b-(LL)((LL)((long double)a/p*b+1e-3)*p))%p+p)%p; 32 } 33 const int m1=M[0],m2=M[1],m3=M[2]; 34 const int inv1=ksm(m1%m2,m2-2,m2),inv2=ksm(m2%m1,m1-2,m1),inv12=ksm(_M%m3,m3-2,m3); 35 int CRT(int a1, int a2, int a3) 36 { 37 LL A=(mul((LL)a1*m2%_M,inv2,_M)+mul((LL)a2*m1%_M,inv1,_M))%_M; 38 LL k=((LL)a3+m3-A%m3)*inv12%m3; 39 return (k*(_M%mod)+A)%mod; 40 }*/ 41 42 inline LL mul(LL a,LL b,LL p){ 43 a%=p; b%=p; 44 return ((a*b-(LL)((LL)((long double)a/p*b+1e-3)*p))%p+p)%p; 45 } 46 47 const int m1=M[0],m2=M[1],m3=M[2]; 48 const int inv1=ksm(m1%m2,m2-2,m2),inv2=ksm(m2%m1,m1-2,m1),inv12=ksm(_M%m3,m3-2,m3); 49 inline int CRT(int a1,int a2,int a3){ 50 LL A=(mul((LL)a1*m2%_M,inv2,_M)+mul((LL)a2*m1%_M,inv1,_M))%_M; 51 LL k=((LL)a3+m3-A%m3)*inv12%m3; 52 return (k*(_M%mod)+A)%mod; 53 } 54 int rev[N]; 55 struct NTT{ 56 int num,P,G; 57 int w[2][N]; 58 void pre(int _P, int _G, int n) 59 { 60 num=n; P=_P; G=_G; 61 int g=ksm(G,(P-1)/num,P); 62 w[1][0]=1; for (int i=1; i<n; i++) w[1][i]=(LL)w[1][i-1]*g%P; 63 w[0][0]=1; for (int i=1; i<n; i++) w[0][i]=w[1][n-i]; 64 } 65 void FFT(int *a, int n, int f) 66 { 67 for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]); 68 for (int i=1; i<n; i<<=1) 69 for (int j=0; j<n; j+=(i<<1)) 70 for (int k=0; k<i; k++) 71 { 72 int x=a[j+k],y=(LL)w[f][num/(i<<1)*k]*a[i+j+k]%P; 73 a[j+k]=(x+y)%P; a[i+j+k]=(x-y+P)%P; 74 } 75 if (!f) for (int i=0,inv=ksm(n,P-2,P); i<n; i++) a[i]=(LL)a[i]*inv%P; 76 } 77 }ntt[3]; 78 79 int n,m; 80 int ans[3][N]; 81 int a[N],b[N],c[N],d[N]; 82 83 int main() 84 { 85 freopen("annona_squamosa.in","r",stdin); freopen("annona_squamosa.out","w",stdout); 86 n=ra(); 87 for (int i=0; i<n; i++) a[i]=ra(); 88 for (int i=0; i<n; i++) b[i]=ra(); 89 for (m=1; m<n<<1; m<<=1); 90 int L=0,x=m>>1; while (!(x&1)) x>>=1,L++; 91 for (int i=0; i<m; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L); 92 // for (m=1;m<=2*(n-1);m<<=1); 93 // int L=0,x=m; while (x>>=1) L++; 94 // for (int i=1;i<=m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)); 95 // for (int i=0; i<m; i++) printf("%d ",rev[i]); cout<<endl; 96 for (int i=0; i<3; i++) ntt[i].pre(M[i],G[i],m); 97 for (int i=0; i<3; i++) 98 { 99 memcpy(c,a,sizeof(int)*(m+5)); memcpy(d,b,sizeof(int)*(m+5)); 100 ntt[i].FFT(c,m,1); ntt[i].FFT(d,m,1); 101 for (int j=0; j<m; j++) c[j]=(LL)c[j]*d[j]%ntt[i].P; 102 ntt[i].FFT(c,m,0); 103 for (int j=0; j<m; j++) ans[i][j]=c[j]; 104 } 105 for (int i=0; i<n; i++) printf("%d ",CRT(ans[0][i],ans[1][i],ans[2][i])); 106 return 0; 107 }
