拆系数FFT样例都没过求助/dk

P4245 【模板】任意模数多项式乘法

YamadaRyou @ 2022-08-11 21:51:02

全输出 0/fad

// 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

草,第 71 行把取模写成除,但是还是不对


|