YamadaRyou @ 2022-08-11 21:51:02
全输出
// Problem: P4245 【模板】任意模数多项式乘法
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P4245
// Memory Limit: 500 MB
// Time Limit: 2000 ms
#include<stdio.h>
#include<ctype.h>
#include<cmath>
namespace fasti{
char buf[1<<21],*p1=buf,*p2=buf;
#define getc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*(p1++))
inline void read(int&x){
char c=getc();
for(;!isdigit(c);c=getc());
for(x=0;isdigit(c);c=getc())x=x*10+(c^48);
}
inline void read(double&x){
char c=getc();
for(;!isdigit(c);c=getc());
for(x=0;isdigit(c);c=getc())x=x*10+(c^48);
}
}
using fasti::read;
namespace fasto{
char obuf[1<<21];int o_p=-1;
inline void flush(){fwrite(obuf,1,o_p+1,stdout),o_p=-1;}
inline void putc(const char &c){if(o_p==(1<<21)-1)flush();obuf[++o_p]=c;}
void print(const int&x){if(x>9)print(x/10);putc((x%10)^48);}
}
using fasto::print;
using fasto::putc;
using fasto::flush;
inline int sum(int x,int y,int p){return (p-x)>y?x+y:y-(p-x);}
inline int dif(int x,int y,int p){return x<y?p-(y-x):x-y;}
inline int pro(int x,int y,int p){return (long long)x*y%p;}
constexpr double PI=std::acos(-1);
int rev[2097152];
struct complex{
double a,b;
complex():a(0),b(0){}
complex(double _a,double _b):a(_a),b(_b){}
complex operator +(const complex&o)const{return complex(a+o.a,b+o.b);}
complex operator -(const complex&o)const{return complex(a-o.a,b-o.b);}
complex operator *(const complex&o)const{
double t1=a*o.b,t2=b*o.a,t3=(a+b)*(o.a-o.b);
return complex(t1-t2+t3,t1+t2);
}
complex operator ~()const{return {a,-b};}
}omega[2097152];
void FFT(complex*a,int op,int n){
for(int i=1;i<n;++i)if(i<rev[i]){complex t=a[i];a[i]=a[rev[i]];a[rev[i]]=t;}
for(int i=1;i<n;i<<=1){
for(int j=0;j<n;j+=(i<<1))
for(int k=j;k<i+j;++k){
complex x=a[k],y=omega[k-j?(op?n-n/(i<<1)*(k-j):n/(i<<1)*(k-j)):0]*a[i+k];
a[k]=x+y,a[i+k]=x-y;
}
}
if(op)for(int i=0;i<n;++i)a[i].a=(a[i].a/n),a[i].b=(a[i].b/n);
}
void init(int n,int l){
for(int i=1;i<n;++i)rev[i]=(i&1)<<l-1|rev[i>>1]>>1;
for(int i=0;i<n;++i)omega[1]=complex(std::cos(2.0*i*PI/n),std::sin(2.0*i*PI/n));
}
void MTT(complex*a,complex*b,int n,int l,int p){
init(n,l);
static complex c[2097152],d[2097152];
int m=std::floor(std::sqrt(p));
for(int i=0;i<n;++i)
a[i].a=(int)(a[i].a)%p,b[i].a=(int)(b[i].a)/p,
c[i]={(int)(a[i].a)%m,(int)(a[i].a)/m},d[i]={(int)(b[i].a)%m,(int)(b[i].a)/m};
FFT(c,0,n),FFT(d,0,n);
complex z1,z2,z3,z4;
for(int i=0;i<n;++i)
z1=(c[i]+(~c[i?n-i:0]))*complex(0.5,0),
z2=(c[i]-(~c[i?n-i:0]))*complex(0,-0.5),
z3=(d[i]+(~d[i?n-i:0]))*complex(0.5,0),
z4=(d[i]-(~d[i?n-i:0]))*complex(0,-0.5),
a[i]=z1*z4+z2*z3,c[i]=z1*z3+z2*z4*complex(0,1);
FFT(a,1,n),FFT(c,1,n);
// for(int i=0;i<n;++i)printf("%d %d %d\n",(int)((c[i]-(~c[i?n-i:0])).b+0.5)/2,(int)(a[i].a+0.5),(int)((c[i]+(~c[i?n-i:0])).a+0.5)/2);
for(int i=0;i<n;++i)
a[i]=complex(sum(sum((int)(((c[i]-(~c[i?n-i:0])).b+0.5)/2)%p,pro((int)(((c[i]+(~c[i?n-i:0])).a+0.5)/2)%p,pro(m,m,p),p),p),pro((int)((int)(a[i].a+0.5))%p,m,p),p),0);
}
complex a[2097152],b[2097152];
int main(){
int n,m,p,t,l;
scanf("%d%d%d",&n,&m,&p);
for(int i=0;i<=n;++i)scanf("%lf",&a[i].a);
for(int i=0;i<=m;++i)scanf("%lf",&b[i].a);
t=n+m+1;l=(t&-t)==t?__builtin_ctz(t):32-__builtin_clz(t);t=1<<l;
MTT(a,b,t,l,p);
for(int i=0;i<=n+m;++i)print((int)a[i].a),putc(' ');
flush();
return 0;
}
by esquigybcu @ 2022-08-11 22:06:31
orz
by harvey2019 @ 2022-08-11 22:38:49
orz
by fjy666 @ 2022-08-12 09:17:22
orz
另:您应该知道我想说啥,所以我不说了
by YamadaRyou @ 2022-08-12 09:21:56
@fjy666 求调/dk
by Register_int @ 2022-08-12 09:59:38
stCCCCCCCCCO
by YamadaRyou @ 2022-08-12 10:05:16
草,我 init
那里把 omega[i]
写成了 omega[1]
,但是改了还是没过
by YamadaRyou @ 2022-08-12 10:13:03
草,第