代码

#include<bits/stdc++.h>
using namespace std;
int const N=(1<<22)+5;
namespace polynomials
{
    int const p=998244353,g=3;
    int w[N],iv[N],r[N],last;
    typedef vector<int> vec;
    inline int mod(int x){return x>=p?x-p:x;}
    int pw(int x,int y)
    {
        int res=1;
        while(y)
        {
            if(y&1)res=1ll*res*x%p;
            x=1ll*x*x%p;
            y>>=1;
        }
        return res;
    }
    void init(int n)
    {
        int lim=1;
        while((lim<<1)<n)lim<<=1;
        iv[1]=1;
        for(int i=2;i<=lim;i++)iv[i]=p-1ll*p/i*iv[p%i]%p;
        for(int i=1;i<=lim;i<<=1)
        {
            int wn=pw(g,(p-1)/(i<<1)),j=0,ww=1;
            for(j=0,ww=1;j+7<i;j+=8)
            {
                w[i+j]=ww;
                ww=1ll*ww*wn%p;
                w[i+j+1]=ww;
                ww=1ll*ww*wn%p;
                w[i+j+2]=ww;
                ww=1ll*ww*wn%p;
                w[i+j+3]=ww;
                ww=1ll*ww*wn%p;
                w[i+j+4]=ww;
                ww=1ll*ww*wn%p;
                w[i+j+5]=ww;
                ww=1ll*ww*wn%p;
                w[i+j+6]=ww;
                ww=1ll*ww*wn%p;
                w[i+j+7]=ww;
                ww=1ll*ww*wn%p;
            }
            for(;j<i;j++)w[i+j]=ww,ww=1ll*ww*wn%p;
        }
    }
    void ntt(vec &f,int n,int op)
    {
        if(last!=n)
        {
            for(int i=1;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)?(n>>1):0);
            last=n;
        }
        for(int i=1;i<n;i++)if(i<r[i])swap(f[i],f[r[i]]);
        for(int i=1;i<n;i<<=1)
            for(int j=0,k=0,x,y;j<n;j+=i<<1)
            {
                for(k=0;k+7<i;k+=8)
                {
                    x=f[j+k],y=1ll*f[i+j+k]*w[i+k]%p;
                    f[j+k]=mod(x+y);f[i+j+k]=mod(x-y+p);
                    x=f[j+k+1],y=1ll*f[i+j+k+1]*w[i+k+1]%p;
                    f[j+k+1]=mod(x+y);f[i+j+k+1]=mod(x-y+p);
                    x=f[j+k+2],y=1ll*f[i+j+k+2]*w[i+k+2]%p;
                    f[j+k+2]=mod(x+y);f[i+j+k+2]=mod(x-y+p);
                    x=f[j+k+3],y=1ll*f[i+j+k+3]*w[i+k+3]%p;
                    f[j+k+3]=mod(x+y);f[i+j+k+3]=mod(x-y+p);
                    x=f[j+k+4],y=1ll*f[i+j+k+4]*w[i+k+4]%p;
                    f[j+k+4]=mod(x+y);f[i+j+k+4]=mod(x-y+p);
                    x=f[j+k+5],y=1ll*f[i+j+k+5]*w[i+k+5]%p;
                    f[j+k+5]=mod(x+y);f[i+j+k+5]=mod(x-y+p);
                    x=f[j+k+6],y=1ll*f[i+j+k+6]*w[i+k+6]%p;
                    f[j+k+6]=mod(x+y);f[i+j+k+6]=mod(x-y+p);
                    x=f[j+k+7],y=1ll*f[i+j+k+7]*w[i+k+7]%p;
                    f[j+k+7]=mod(x+y);f[i+j+k+7]=mod(x-y+p);
                }
                for(;k<i;k++)
                {
                    x=f[j+k],y=1ll*f[i+j+k]*w[i+k]%p;
                    f[j+k]=mod(x+y);f[i+j+k]=mod(x-y+p);
                }
            }
        if(op==-1)
        {
            reverse(&f[1],&f[n]);
            int i=0;
            for(i=0;i+7<n;i+=8)
            {
                f[i]=1ll*f[i]*iv[n]%p;
                f[i+1]=1ll*f[i+1]*iv[n]%p;
                f[i+2]=1ll*f[i+2]*iv[n]%p;
                f[i+3]=1ll*f[i+3]*iv[n]%p;
                f[i+4]=1ll*f[i+4]*iv[n]%p;
                f[i+5]=1ll*f[i+5]*iv[n]%p;
                f[i+6]=1ll*f[i+6]*iv[n]%p;
                f[i+7]=1ll*f[i+7]*iv[n]%p;
            }
            for(;i<n;i++)
                f[i]=1ll*f[i]*iv[n]%p;
        }
    }
}
using namespace polynomials;
vec A,B;
void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c)
{
    n++,m++;
    int lim=1;
    while(lim<(n+m-1))lim<<=1;
    init(lim+5);
    A.resize(lim),B.resize(lim);
    for(int i=0;i<n;i++)A[i]=a[i];
    for(int i=0;i<m;i++)B[i]=b[i];
    ntt(A,lim,1),ntt(B,lim,1);
    int i=0;
    for(i=0;i+7<lim;i+=8)
    {
        A[i]=1ll*A[i]*B[i]%p;
        A[i+1]=1ll*A[i+1]*B[i+1]%p;
        A[i+2]=1ll*A[i+2]*B[i+2]%p;
        A[i+3]=1ll*A[i+3]*B[i+3]%p;
        A[i+4]=1ll*A[i+4]*B[i+4]%p;
        A[i+5]=1ll*A[i+5]*B[i+5]%p;
        A[i+6]=1ll*A[i+6]*B[i+6]%p;
        A[i+7]=1ll*A[i+7]*B[i+7]%p;
    }
    for(;i<lim;i++)A[i]=1ll*A[i]*B[i]%p;
    ntt(A,lim,-1);
    for(int i=0;i<n+m-1;i++)c[i]=A[i];
}

评测结果