代码

// Radix-4 DIF forward + DIT inverse (no bitrev), AVX2 Montgomery.
#pragma GCC target("avx2")
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <immintrin.h>
typedef unsigned long long u64; typedef unsigned int u32;
const u32 MOD=998244353u,NINV=998244351u,R2=932051910u,G=3u;
const u32 IM=911660635u, IINV=MOD-911660635u;
const int MAXN=1<<18;
alignas(32) static u32 a[MAXN],b[MAXN];
alignas(32) static u32 tw1[MAXN],tw2[MAXN],tw3[MAXN];
alignas(32) static u32 tw1i[MAXN],tw2i[MAXN],tw3i[MAXN];
static inline u32 modpow(u32 base,u64 e){u64 r=1,bb=base%MOD;for(;e;e>>=1){if(e&1)r=r*bb%MOD;bb=bb*bb%MOD;}return (u32)r;}
static inline u32 mont_s(u32 a,u32 b){u64 t=(u64)a*b;u32 m=(u32)t*NINV;u32 u=(u32)((t+(u64)m*MOD)>>32);if(u>=MOD)u-=MOD;return u;}
static inline __m256i mont8(__m256i x,__m256i y,const __m256i&ninv,const __m256i&modv,const __m256i&modm1,const __m256i&shuf){
    __m256i t_lo=_mm256_mullo_epi32(x,y);
    __m256i m=_mm256_mullo_epi32(t_lo,ninv);
    __m256i pe=_mm256_mul_epu32(x,y);
    __m256i po=_mm256_mul_epu32(_mm256_srli_epi64(x,32),_mm256_srli_epi64(y,32));
    __m256i me=_mm256_mul_epu32(m,modv);
    __m256i mo=_mm256_mul_epu32(_mm256_srli_epi64(m,32),modv);
    __m256i ue=_mm256_srli_epi64(_mm256_add_epi64(pe,me),32);
    __m256i uo=_mm256_srli_epi64(_mm256_add_epi64(po,mo),32);
    __m256i e=_mm256_shuffle_epi8(ue,shuf);
    __m256i o=_mm256_shuffle_epi8(uo,shuf);
    __m256i u=_mm256_unpacklo_epi32(e,o);
    __m256i mask=_mm256_cmpgt_epi32(u,modm1);
    u=_mm256_sub_epi32(u,_mm256_and_si256(mask,modv));
    return u;
}
static inline __m256i addm(__m256i a,__m256i b,const __m256i&modv){__m256i s=_mm256_add_epi32(a,b);return _mm256_min_epu32(s,_mm256_sub_epi32(s,modv));}
static inline __m256i subm(__m256i a,__m256i b,const __m256i&modv){__m256i d=_mm256_sub_epi32(a,b);return _mm256_min_epu32(d,_mm256_add_epi32(d,modv));}

static void dif4(u32*x,int n,const u32*TW1,const u32*TW2,const u32*TW3,u32 Im){
    const __m256i ninv=_mm256_set1_epi32(NINV),modv=_mm256_set1_epi32(MOD),modm1=_mm256_set1_epi32(MOD-1);
    const __m256i shuf=_mm256_setr_epi8(0,1,2,3,8,9,10,11,0,0,0,0,0,0,0,0, 0,1,2,3,8,9,10,11,0,0,0,0,0,0,0,0);
    const __m256i imv=_mm256_set1_epi32(Im);
    for(int len=n;len>1;len>>=2){int s=len>>2;
        for(int i=0;i<n;i+=len){u32*y=x+i;int p=0;
            for(;p+8<=s;p+=8){
                __m256i x0=_mm256_loadu_si256((__m256i*)(y+p));
                __m256i x1=_mm256_loadu_si256((__m256i*)(y+p+s));
                __m256i x2=_mm256_loadu_si256((__m256i*)(y+p+2*s));
                __m256i x3=_mm256_loadu_si256((__m256i*)(y+p+3*s));
                __m256i t0=addm(x0,x2,modv);
                __m256i t1=subm(x0,x2,modv);
                __m256i t2=addm(x1,x3,modv);
                __m256i t3=mont8(subm(x1,x3,modv),imv,ninv,modv,modm1,shuf);
                __m256i z0=addm(t0,t2,modv);
                __m256i z2=subm(t0,t2,modv);
                __m256i z1=addm(t1,t3,modv);
                __m256i z3=subm(t1,t3,modv);
                __m256i w1=_mm256_loadu_si256((__m256i*)(TW1+p));
                __m256i w2=_mm256_loadu_si256((__m256i*)(TW2+p));
                __m256i w3=_mm256_loadu_si256((__m256i*)(TW3+p));
                _mm256_storeu_si256((__m256i*)(y+p),z0);
                _mm256_storeu_si256((__m256i*)(y+p+s),mont8(z1,w1,ninv,modv,modm1,shuf));
                _mm256_storeu_si256((__m256i*)(y+p+2*s),mont8(z2,w2,ninv,modv,modm1,shuf));
                _mm256_storeu_si256((__m256i*)(y+p+3*s),mont8(z3,w3,ninv,modv,modm1,shuf));
            }
            for(;p<s;p++){
                u32 X0=y[p],X1=y[p+s],X2=y[p+2*s],X3=y[p+3*s];
                u32 T0=X0+X2;if(T0>=MOD)T0-=MOD;
                u32 T1v=X0-X2;if(T1v>=MOD)T1v+=MOD;
                u32 T2=X1+X3;if(T2>=MOD)T2-=MOD;
                u32 T3=mont_s((u32)(X1-X3+MOD*(X1<X3)),Im);
                u32 Z0=T0+T2;if(Z0>=MOD)Z0-=MOD;
                u32 Z2=T0-T2;if(Z2>=MOD)Z2+=MOD;
                u32 Z1=T1v+T3;if(Z1>=MOD)Z1-=MOD;
                u32 Z3=T1v-T3;if(Z3>=MOD)Z3+=MOD;
                y[p]=Z0;y[p+s]=mont_s(Z1,TW1[p]);y[p+2*s]=mont_s(Z2,TW2[p]);y[p+3*s]=mont_s(Z3,TW3[p]);
            }
        }
        TW1+=s;TW2+=s;TW3+=s;
    }
}
static void dit4(u32*x,int n,const u32*TW1,const u32*TW2,const u32*TW3,u32 Im){
    const __m256i ninv=_mm256_set1_epi32(NINV),modv=_mm256_set1_epi32(MOD),modm1=_mm256_set1_epi32(MOD-1);
    const __m256i shuf=_mm256_setr_epi8(0,1,2,3,8,9,10,11,0,0,0,0,0,0,0,0, 0,1,2,3,8,9,10,11,0,0,0,0,0,0,0,0);
    const __m256i imv=_mm256_set1_epi32(Im);
    for(int len=4;len<=n;len<<=2){int s=len>>2;
        for(int i=0;i<n;i+=len){u32*y=x+i;int p=0;
            for(;p+8<=s;p+=8){
                __m256i x0=_mm256_loadu_si256((__m256i*)(y+p));
                __m256i x1=mont8(_mm256_loadu_si256((__m256i*)(y+p+s)),_mm256_loadu_si256((__m256i*)(TW1+p)),ninv,modv,modm1,shuf);
                __m256i x2=mont8(_mm256_loadu_si256((__m256i*)(y+p+2*s)),_mm256_loadu_si256((__m256i*)(TW2+p)),ninv,modv,modm1,shuf);
                __m256i x3=mont8(_mm256_loadu_si256((__m256i*)(y+p+3*s)),_mm256_loadu_si256((__m256i*)(TW3+p)),ninv,modv,modm1,shuf);
                __m256i t0=addm(x0,x2,modv);
                __m256i t1=subm(x0,x2,modv);
                __m256i t2=addm(x1,x3,modv);
                __m256i t3=mont8(subm(x1,x3,modv),imv,ninv,modv,modm1,shuf);
                __m256i z0=addm(t0,t2,modv);
                __m256i z2=subm(t0,t2,modv);
                __m256i z1=addm(t1,t3,modv);
                __m256i z3=subm(t1,t3,modv);
                _mm256_storeu_si256((__m256i*)(y+p),z0);
                _mm256_storeu_si256((__m256i*)(y+p+s),z1);
                _mm256_storeu_si256((__m256i*)(y+p+2*s),z2);
                _mm256_storeu_si256((__m256i*)(y+p+3*s),z3);
            }
            for(;p<s;p++){
                u32 X0=y[p];
                u32 X1=mont_s(y[p+s],TW1[p]);
                u32 X2=mont_s(y[p+2*s],TW2[p]);
                u32 X3=mont_s(y[p+3*s],TW3[p]);
                u32 T0=X0+X2;if(T0>=MOD)T0-=MOD;
                u32 T1v=X0-X2;if(T1v>=MOD)T1v+=MOD;
                u32 T2=X1+X3;if(T2>=MOD)T2-=MOD;
                u32 T3=mont_s((u32)(X1-X3+MOD*(X1<X3)),Im);
                u32 Z0=T0+T2;if(Z0>=MOD)Z0-=MOD;
                u32 Z2=T0-T2;if(Z2>=MOD)Z2+=MOD;
                u32 Z1=T1v+T3;if(Z1>=MOD)Z1-=MOD;
                u32 Z3=T1v-T3;if(Z3>=MOD)Z3+=MOD;
                y[p]=Z0;y[p+s]=Z1;y[p+2*s]=Z2;y[p+3*s]=Z3;
            }
        }
        TW1+=s;TW2+=s;TW3+=s;
    }
}

static char inbuf[1<<20]; static size_t inpos=0,inlen=0;
static inline int readbyte(){if(inpos>=inlen){inlen=fread(inbuf,1,1<<20,stdin);inpos=0;if(inlen==0)return -1;}return (unsigned char)inbuf[inpos++];}
static inline int readint(){int c=readbyte();while(c==' '||c=='\n'||c=='\r'||c=='\t')c=readbyte();int x=0;while(c>='0'&&c<='9'){x=x*10+(c-'0');c=readbyte();}return x;}
static char outbuf[1<<21]; static size_t outpos=0;
static inline void putc(char c){outbuf[outpos++]=c;}
static inline void putint(int x){if(x==0){putc('0');return;}char t[12];int k=0;while(x){t[k++]='0'+(x%10);x/=10;}while(k)putc(t[--k]);}

int main(){
    int n=readint(), m=readint();
    int na=n+1, nb=m+1;
    for(int i=0;i<na;i++)a[i]=(u32)readint();
    for(int i=0;i<nb;i++)b[i]=(u32)readint();
    int size=1; while(size<na+nb-1)size<<=1;
    {int s=size,lg=0;while(s>1){s>>=1;lg++;}if(lg&1)size<<=1;}
    u32 wstep=(MOD-1)/(u32)size;
    u32 w=modpow(G,wstep), iw=modpow(w,MOD-2);
    {u32*p1=tw1,*p2=tw2,*p3=tw3;
     for(int len=size;len>1;len>>=2){int s=len>>2,step=size/len;
       u32 ws=modpow(w,(u64)step); u32 wsm=mont_s(ws,R2);
       u32 ws2m=mont_s(wsm,wsm), ws3m=mont_s(ws2m,wsm);
       p1[0]=p2[0]=p3[0]=mont_s(1,R2);
       for(int p=1;p<s;p++){p1[p]=mont_s(p1[p-1],wsm);p2[p]=mont_s(p2[p-1],ws2m);p3[p]=mont_s(p3[p-1],ws3m);}
       p1+=s;p2+=s;p3+=s;}}
    {u32*p1=tw1i,*p2=tw2i,*p3=tw3i;
     for(int len=4;len<=size;len<<=2){int s=len>>2,step=size/len;
       u32 ws=modpow(iw,(u64)step); u32 wsm=mont_s(ws,R2);
       u32 ws2m=mont_s(wsm,wsm), ws3m=mont_s(ws2m,wsm);
       p1[0]=p2[0]=p3[0]=mont_s(1,R2);
       for(int p=1;p<s;p++){p1[p]=mont_s(p1[p-1],wsm);p2[p]=mont_s(p2[p-1],ws2m);p3[p]=mont_s(p3[p-1],ws3m);}
       p1+=s;p2+=s;p3+=s;}}
    u32 Im=mont_s(IM,R2), Iinv=mont_s(IINV,R2);
    for(int i=0;i<size;i++)a[i]=mont_s(a[i],R2);
    for(int i=0;i<size;i++)b[i]=mont_s(b[i],R2);
    dif4(a,size,tw1,tw2,tw3,Im);
    dif4(b,size,tw1,tw2,tw3,Im);
    {const __m256i ninv=_mm256_set1_epi32(NINV),modv=_mm256_set1_epi32(MOD),modm1=_mm256_set1_epi32(MOD-1);
     const __m256i shuf=_mm256_setr_epi8(0,1,2,3,8,9,10,11,0,0,0,0,0,0,0,0, 0,1,2,3,8,9,10,11,0,0,0,0,0,0,0,0);
     for(int i=0;i<size;i+=8){__m256i x=_mm256_loadu_si256((__m256i*)(a+i));__m256i y=_mm256_loadu_si256((__m256i*)(b+i));_mm256_storeu_si256((__m256i*)(a+i),mont8(x,y,ninv,modv,modm1,shuf));}}
    dit4(a,size,tw1i,tw2i,tw3i,Iinv);
    u32 ninv_scale=modpow((u32)size,MOD-2);
    for(int i=0;i<size;i++)a[i]=mont_s(a[i],ninv_scale);
    int outn=n+m+1;
    for(int i=0;i<outn;i++){if(i)putc(' ');putint((int)a[i]);}
    putc('\n'); fwrite(outbuf,1,outpos,stdout);
    return 0;
}

评测结果