萌新五次 FFT WA 50pts 求助

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

Moeebius @ 2022-08-19 11:08:49

RT,\mathtt{long\ double+\_\_int128} 无果,求 debug……

#include<bits/stdc++.h>
using namespace std;

#define il inline
#define mkp make_pair
#define pii pair<int,int>
#define lll __int128
#define ll __int128_t // 注意这一行!!
#define For(i,j,k) for(int i=(j); i<=(k); ++i)
#define ForDown(i,j,k) for(int i=(j); i>=(k); --i)
#define pb push_back
#define init(filename) freopen(filename ".in" ,"r",stdin);freopen(filename ".out" ,"w",stdout)
template<typename T> 
il void read(T &x){ x=0;int f=1;char c=getchar();while(!isdigit(c)){if(c=='-')f=-1;c=getchar();}while(isdigit(c)){x=x*10+c-'0';c=getchar();}x*=f;}
template<typename T, typename ... Args>
il void read(T &x, Args &... y){ read(x);read(y...); }

ll n,m,mod;

namespace MTT_algorithm
{
    struct Complex
    {
        long double real,imag;
        Complex(double __r=0, double __i=0) : real(__r),imag(__i) {}
    };
    il Complex operator * (const Complex &a, const Complex &b)
    {
        return Complex(a.real*b.real-a.imag*b.imag,a.real*b.imag+a.imag*b.real);
    }
    il void operator *= (Complex &a, const Complex &b)
    {
        a=Complex(a.real*b.real-a.imag*b.imag,a.real*b.imag+a.imag*b.real);
    }
    il Complex operator + (const Complex &a, const Complex &b)
    {
        return Complex(a.real+b.real,a.imag+b.imag);
    }
    il Complex operator - (const Complex &a, const Complex &b)
    {
        return Complex(a.real-b.real,a.imag-b.imag);
    }

    static const int MAXN=5e6+6;
    static int rev[MAXN];
    static const long double PI=acos(-1);
    il void change(Complex *A, int N)
    {
        For(i,0,N-1)
        {
            rev[i]=rev[i>>1]>>1;
            if(i&1) rev[i]|=N>>1;
        }
        For(i,0,N-1)
        {
            if(i<rev[i])
            {
                swap(A[i],A[rev[i]]);
            }
        }
    }
    void FFT(Complex *A, int N, int tp)
    {
        change(A,N);
        for(int len=2; len<=N; len<<=1)
        {
            Complex W(cos(2.0*PI/len),sin(2.0*tp*PI/len));
            for(int i=0; i<N; i+=len)
            {
                Complex w(1.0,0.0),x,y;
                for(int j=i; j<i+(len>>1); j++)
                {
                    x=A[j],y=A[j+(len>>1)]*w;
                    A[j]=x+y;
                    A[j+(len>>1)]=x-y;
                    w*=W;
                }
            }
        }
        if(tp==-1)
        {
            For(i,0,N-1)
            {
                A[i].real/=N;
                A[i].imag/=N;
            }
        }
    }
    Complex F[MAXN],G[MAXN],H[MAXN];
    ll get(long double x) { return (ll)(x+0.5); }
    void mul(ll *A, ll *B, int N, int M, ll *R)
    {
        For(i,0,N) F[i]=A[i]&32767,G[i]=A[i]>>15;
        For(i,0,M) H[i]=Complex(B[i]&32767,B[i]>>15);
        int len=1; while(len<=N+M) len<<=1;
        FFT(F,len,1); FFT(G,len,1); FFT(H,len,1);
        For(i,0,len-1) F[i]*=H[i],G[i]*=H[i];
        FFT(F,len,-1); FFT(G,len,-1);
        For(i,0,N+M) (R[i]=get(F[i].real)+((get(F[i].imag)+get(G[i].real))<<15)%mod+(get(G[i].imag)<<30)%mod)%=mod; 
    }
};

const int MAXN=5e6+6;
ll f[MAXN],g[MAXN],r[MAXN];

signed main()
{
    read(n,m,mod);
    For(i,0,n) read(f[i]);
    For(i,0,m) read(g[i]);
    MTT_algorithm::mul(f,g,n,m,r);
    For(i,0,n+m) printf("%lld ", (long long)r[i]);
    return 0;
}

by thomaswmy @ 2022-08-19 11:18:44

%%%


by Ruiqun2009 @ 2022-11-18 14:51:45

预处理单位根试试?这样既不用开 long double,也不用开 __int128


by Ruiqun2009 @ 2022-11-18 14:51:58

@Xiaohuba


|