WA50分求助

Lacrymabre @ 2020-03-26 19:48:23

//借鉴了MTT模板
#include<iostream>
#include<cstdio>
#include<iomanip>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<queue>
#include<stack>
#include<vector>
#include<map>
#include<set>

#define ll long long
#define db double
#define MAX 0x7fffffff
#define init inline int
#define INF 0X3fffffff
#define MAXX 1e-6

using namespace std;

inline long long read()
{
    ll f=1,s=0;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')s=(s<<1)+(s<<3)+ch-'0',ch=getchar();
    return s*f;
}

const int N=4e5+5;
const double pi=acos(-1.0);
const int M=32768;

struct cp
{
    double x,y;
    cp(double x=0,double y=0):x(x),y(y){}
};

cp operator+(cp a,cp b) {return cp(a.x + b.x , a.y + b.y);}
cp operator-(cp a,cp b) {return cp(a.x - b.x , a.y - b.y);}
cp operator*(cp a,cp b) {return cp(a.x*b.x - a.y*b.y , a.x*b.y + a.y*b.x);}

ll lim,l;
ll r[N];
ll n,m,mod;
cp a1[N],b1[N],a2[N],b2[N];
ll ans[N]; 
cp x[N];

ll mypow(ll n,ll m,ll mod) {
    ll ret=1;
    while(m){
        if (m&1) ret=(ret*n)%mod;
        n=(n*n)%mod;
        m>>=1;
    }
    return ret%mod;
}

void FFT(cp*f,int op)//fft
{
    for(int i=0;i<lim;++i)
    if(i<r[i]) swap(f[i],f[r[i]]);
    for(int j=1;j<lim;j<<=1)
    {
        cp F(cos(pi/j),op*sin(pi/j)) ;
        for(int k=0;k<lim;k+=(j<<1)){
            cp t(1,0);
            for(int l=0;l<j;l++,t=t*F){
                cp x=f[k+l],y=t*f[k+j+l] ;
                f[k+l]=x+y;
                f[k+j+l]=x-y;
            }
        }
    }
}

void back(cp *a,cp *b,ll res,ll mod){//乘回去 
    for(int i=0;i<lim;++i) x[i]=a[i]*b[i];
    FFT(x,-1);
    for(int i=0;i<=n+m;++i)(ans[i]+=(ll)(x[i].x/lim+0.5)%mod*res%mod)%=mod;
}

void MTT(cp *a,cp *b,cp *c,cp *d,ll mod){//MTT
    FFT(a,1);FFT(b,1);FFT(c,1);FFT(d,1);
    back(a,c,M*M%mod,mod);back(a,d,M%mod,mod);
    back(c,b,M%mod,mod);back(b,d,1,mod);
} 

int main()
{
    n=read();m=read();mod=read();
    for(int i=0,x;i<=n;++i) x=read(),a1[i].x=x/M,b1[i].x=x%M;//分开来 
    for(int i=0,x;i<=m;++i) x=read(),a2[i].x=x/M,b2[i].x=x%M;
    for(lim=1,l=0;lim<=n+m;lim<<=1,l++);
    for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    MTT(a1,b1,a2,b2,mod);
    for(int i=0;i<=n+m;++i) cout<<ans[i]<<" ";
}

by Lacrymabre @ 2020-03-26 19:51:03

我谔谔

by Smile_Cindy @ 2020-03-26 19:53:02

@爆零_自动机 试下long double，不行的话可以看下这份代码

虽然常数大，但是可以保证精度，还好理解。

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const long double PI=acos(-1);
const int MAX_N=1<<18;
const int base=1<<15;
int n,m,Mod; 
struct Complex
{
    long double x,y;
    Complex(){}
    Complex(long double _x,long double _y)
    {
        x=_x;
        y=_y;
    }
};
Complex operator + (Complex w1,Complex w2)
{
    return Complex(w1.x+w2.x,w1.y+w2.y);
}
Complex operator - (Complex w1,Complex w2)
{
    return Complex(w1.x-w2.x,w1.y-w2.y);
}
Complex operator * (Complex w1,Complex w2)
{
    return Complex(w1.x*w2.x-w1.y*w2.y,w1.x*w2.y+w1.y*w2.x);
}
int rev[MAX_N];
void FFT(Complex A[],int n,int flag)
{
    for(int i=0;i<n;i++)if(rev[i]>i)swap(A[i],A[rev[i]]);
    for(int i=1;i<n;i<<=1)
    {
        Complex w1=Complex(cos(2*PI*flag/(i<<1)),sin(2*PI*flag/(i<<1)));
        for(int j=0;j<n;j+=(i<<1))
        {
            Complex w2=Complex(1,0);
            for(int k=0;k<i;k++)
            {
                Complex t1=A[j+k],t2=A[i+j+k]*w2;
                A[j+k]=t1+t2;
                A[i+j+k]=t1-t2;
                w2=w2*w1;
            }
        }
    }
    if(flag==-1)for(int i=0;i<n;i++)A[i].x/=n;
}
int len=1,L=0;
void init()
{
    while(len<=n+m)
    {
        len<<=1;
        L++;
    }
    for(int i=1;i<len;i++)rev[i]=rev[i>>1]>>1|(i&1)<<(L-1);
}
Complex A1[MAX_N],A2[MAX_N],B1[MAX_N],B2[MAX_N];
Complex A[MAX_N],B[MAX_N],C[MAX_N];
int F[MAX_N],G[MAX_N];
int ans[MAX_N];
signed main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
    cin>>n>>m>>Mod;
    for(int i=0;i<=n;i++)cin>>F[i];
    for(int i=0;i<=m;i++)cin>>G[i];
    init();
    for(int i=0;i<=n;i++)A1[i].x=F[i]/base,B1[i].x=F[i]%base;
    for(int i=0;i<=m;i++)A2[i].x=G[i]/base,B2[i].x=G[i]%base;
    FFT(A1,len,1);
    FFT(A2,len,1);
    FFT(B1,len,1);
    FFT(B2,len,1);
    for(int i=0;i<len;i++)A[i]=A1[i]*A2[i];
    for(int i=0;i<len;i++)B[i]=A1[i]*B2[i]+A2[i]*B1[i];
    for(int i=0;i<len;i++)C[i]=B1[i]*B2[i];
    FFT(A,len,-1);
    FFT(B,len,-1);
    FFT(C,len,-1);
    for(int i=0;i<len;i++)
    {
        ans[i]=(LL)(A[i].x+0.5)%Mod;
        ans[i]=(1LL*ans[i]*base+(LL)(B[i].x+0.5))%Mod;
        ans[i]=(1LL*ans[i]*base+(LL)(C[i].x+0.5))%Mod;
    }
    for(int i=0;i<=n+m;i++)cout<<ans[i]<<" \n"[i==n+m];
    return 0;
}

by Lacrymabre @ 2020-03-26 19:53:56

@Alpha 谢谢大佬！我先去试试。

by Lacrymabre @ 2020-03-26 19:55:13

@Alpha 能不能请大佬解释下普通double会爆的原因？？