Orz说好的long double能过呢???

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

小菜鸟 @ 2018-08-08 22:53:48

50分QAQ

#include<cstdio>
#include<math.h>
#include<algorithm>
using std::sin;
using std::cos;
using std::swap;
using std::max;

const int maxl(1<<20);
const long double PI=3.1415926535897932384626433832795028841971693993751;
int P;
int rev[maxl];

char buf[1<<22],*p1=buf,*p2=buf;
namespace quickio
{
    static char getchar(){return ((p1==p2)&&(p2=(p1=buf)+fread(buf,1,1<<22,stdin)),p1==p2?EOF:*p1++);}
    template<class T>
    static T read(T& x)
    {
        using quickio::getchar;
        T n=0;
        bool sym=0;
        char c=getchar();
        while(c<48||c>57)sym|=(c==45),c=getchar();
        while(c>=48&c<=57)n=(n<<1)+(n<<3)+(c^48),c=getchar();
        return x=sym?~n+1:n;
    }
    template<class T>
    static void write(T x)
    {
        if(x<0)putchar(45),x=~x+1;
        if(x>=10)write(x/10);
        putchar((x%10)^48);
    }
};
using namespace quickio;

template<class T>
class complex
{
    T r,i;
    public:
    complex(const T& x=0,const T& y=0){r=x;i=y;}
    ~complex(){}
    T real()const{return r;}
    T imag()const{return i;}
    static T abs(const complex& x){return sqrt(x.r*x.r+x.i*x.i);}
    static T arg(const complex& x){return atan2(x.i,x.r);}
    static complex polar(const T& _r,const T& Theta){return complex(_r*cos(Theta),_r*sin(Theta));}
    complex conj(){return complex(r,-i);}
    complex operator+=(const T& b){r+=b;return *this;}
    complex operator+=(const complex& b){r+=b.r;i+=b.i;return *this;}
    template<class U>
    complex operator+(const U& b)const
    {
        complex c=*this;
        c+=b;
        return c;
    }
    friend complex operator+(const T& a,const complex& b){return b+a;}
    complex operator-=(const T& b){r-=b;return *this;}
    complex operator-=(const complex& b){r-=b.r;i-=b.i;return *this;}
    template<class U>
    complex operator-(const U& b)const
    {
        complex c=*this;
        c-=b;
        return c;
    }
    friend complex operator-(const T& a,const complex& b){return complex(a)-b;}
    complex operator*=(const T& b){r*=b,i*=b;return *this;}
    complex operator*=(const complex& b){complex c=*this;r=b.r*c.r-b.i*c.i;i=b.r*c.i+b.i*c.r;return *this;}
    template<class U>
    complex operator*(const U& b)const
    {
        complex c=*this;
        c*=b;
        return c;
    }
    friend complex operator*(const T& a,const complex& b){return b*a;}
    complex operator/=(const T& b){r/=b,i/=b;return *this;}
    complex operator/=(const complex& b){complex c=b.conj();*this/=(b*c).r;*this*=c;return *this;}
    template<class U>
    complex operator/(const U& b)const
    {
        complex c=*this;
        c/=b;
        return c;
    }
    friend complex operator/(const T& a,const complex& b){return complex(a)/b;}
    friend void fft(complex*,int,int,bool);
};
typedef complex<long double> cd;

void getrev(int bit)
{
    for(int i=0;i<(1<<bit);++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
void fft(cd *a,int n,int dft=1)
{
    for(int i=0;i<n;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int step=1;step<n;step<<=1)
    {
        cd wn(cos(PI/step),sin(dft*PI/step));
        for(int j=0;j<n;j+=step<<1)
        {
            cd wnk(1);
            for(int k=j;k<j+step;++k)
            {
                cd x=a[k];
                cd y=a[k+step]*wnk;
                a[k]=x+y;
                a[k+step]=x-y;
                wnk*=wn;
            }
        }
    }
    if(dft==-1)for(int i=0;i<n;++i)a[i]/=n;
}

void mtt(int deg,int *l,int *r,int *res)
{
    cd A[maxl],B[maxl],C[maxl],D[maxl],E[maxl],F[maxl],G[maxl],H[maxl];
    int len=1,bit=0;
    while(len<deg<<1)len<<=1,++bit;
    getrev(bit);
    for(int i=0;i<len;++i)
    {
        l[i]%=P;r[i]%=P;
        A[i]=l[i]>>15,B[i]=l[i]&0xffff;//get hbits and lbits of F(x)
        C[i]=r[i]>>15,D[i]=r[i]&0xffff;//get hbits and lbits of G(x)
    }
    fft(A,len),fft(B,len);
    fft(C,len),fft(D,len);
    for(int i=0;i<len;++i)
    {
        E[i]=A[i]*C[i];//hbits*hbits
        F[i]=B[i]*C[i];//lbits*hbits
        G[i]=A[i]*D[i];//hbits*lbits
        H[i]=B[i]*D[i];//lbits*lbits
    }
    fft(E,len,-1);
    fft(F,len,-1);
    fft(G,len,-1);
    fft(H,len,-1);
    for(int i=0;i<len;++i)res[i]=(
        (((long long)(E[i].real()+0.5)%P)<<30)%P//hbits
        +(((long long)(F[i].real()+0.5)%P)<<15)%P//mbits
        +(((long long)(G[i].real()+0.5)%P)<<15)%P//mbits
        +(long long)(H[i].real()+0.5)%P//lbits
        )%P;
}

int n,m,deg,F[maxl],G[maxl],H[maxl];

int main()
{
    read(n),read(m),read(P);
    for(int i=0;i<=n;++i)read(F[i]);
    for(int i=0;i<=m;++i)read(G[i]);
    mtt(max(n,m),F,G,H);
    for(int i=0;i<=n+m;++i)write(H[i]),putchar(' ');
}

@panda_2134 求出题巨佬帮我看看啊

@jiyutian
@Stilwell
@xumingyang
@a3_de_liefde_zz
所有我认识的julao都来帮帮我啊。。。


by かねき けん @ 2018-08-09 00:02:49

这谁能看的懂啊


by かねき けん @ 2018-08-09 00:03:02

逗我呢


by 姬小路秋子 @ 2018-08-09 00:44:42

天天打模板有意思吗?

不如做点别的


by 小菜鸟 @ 2018-08-10 19:36:02

@a3_de_liefde_zz julao帮我看看啊
智商不够只好练模板啊。。。
毕竟靠着FFT我过了唯一一道我能过的ZJOI,,,


|