代码

// 1002 in-place single-table radix-8 NTT (mont8 twiddles, shoup internal).
// Layout: A in-place in tasklib a[0..2^21) (extends into b's old region + part of c),
// B in static buffer[2^21], ONE twiddle table (forward used for 2 fwd NTTs, then
// overwritten with inverse twiddles for the inverse NTT). Output copy buffer->c.
#pragma GCC target("avx2")
#pragma GCC optimize("O3,unroll-loops")
#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<<21;
alignas(32) static u32 buffer[MAXN];
alignas(32) static u32 tw[MAXN+64];
alignas(32) static u32 tmp[MAXN/8+8];

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 u32 add_s(u32 a,u32 b){u32 s=a+b; return s>=MOD? s-MOD:s;}
static inline u32 sub_s(u32 a,u32 b){u32 d=a-b; return (int)d<0? d+MOD:d;}

static inline __m256i mont8(__m256i x,__m256i y,const __m256i&ninv,const __m256i&modv,const __m256i&shuf){
    __m256i ao=_mm256_srli_epi64(x,32);
    __m256i bo=_mm256_srli_epi64(y,32);
    __m256i te=_mm256_mul_epu32(x,y);
    __m256i to=_mm256_mul_epu32(ao,bo);
    __m256i me=_mm256_mul_epu32(te,ninv);
    __m256i mo=_mm256_mul_epu32(to,ninv);
    __m256i ue=_mm256_srli_epi64(_mm256_add_epi64(te,_mm256_mul_epu32(me,modv)),32);
    __m256i uo=_mm256_srli_epi64(_mm256_add_epi64(to,_mm256_mul_epu32(mo,modv)),32);
    __m256i u=_mm256_or_si256(ue,_mm256_slli_epi64(uo,32));
    u=_mm256_min_epu32(u,_mm256_sub_epi32(u,modv));
    return u;
}
const u64 MAGIC=18479187002ULL;
static inline u32 shp_quot(u32 w){ return (u32)(((u64)w*MAGIC)>>32); }
static inline u32 shoup_mul(u32 a,u32 w,u32 wp){u32 q=(u32)(((u64)a*wp)>>32);u32 r=(u32)((u64)a*w-(u64)q*MOD);if(r>=MOD)r-=MOD;return r;}
static inline __m256i shoup_mul8(__m256i a,__m256i w,__m256i wp,const __m256i&modv,const __m256i&shuf){
    __m256i lo=_mm256_mullo_epi32(a,w);
    __m256i qe=_mm256_mul_epu32(a,wp);
    __m256i qo=_mm256_mul_epu32(_mm256_srli_epi64(a,32),_mm256_srli_epi64(wp,32));
    __m256i qe2=_mm256_srli_epi64(qe,32);
    __m256i qo2=_mm256_srli_epi64(qo,32);
    __m256i q=_mm256_or_si256(qe2,_mm256_slli_epi64(qo2,32));
    __m256i r=_mm256_sub_epi32(lo,_mm256_mullo_epi32(q,modv));
    return _mm256_min_epu32(r,_mm256_sub_epi32(r,modv));
}
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 ntt8_vec(u32*x,int n,u32 Z1,u32 Z3,u32 Iv,u32 Z1p,u32 Z3p,u32 Ivp){
    const __m256i ninv=_mm256_set1_epi32(NINV),modv=_mm256_set1_epi32(MOD);
    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 z1=_mm256_set1_epi32(Z1),z3=_mm256_set1_epi32(Z3),iv=_mm256_set1_epi32(Iv);
    const __m256i z1m=_mm256_set1_epi32(mont_s(Z1,R2)),z3m=_mm256_set1_epi32(mont_s(Z3,R2)),ivm=_mm256_set1_epi32(mont_s(Iv,R2));
    for(int i=0;i<n;i+=64){
        __m256i a0=_mm256_loadu_si256((__m256i*)(x+i+0)), a1=_mm256_loadu_si256((__m256i*)(x+i+8)), a2=_mm256_loadu_si256((__m256i*)(x+i+16)), a3=_mm256_loadu_si256((__m256i*)(x+i+24));
        __m256i a4=_mm256_loadu_si256((__m256i*)(x+i+32)), a5=_mm256_loadu_si256((__m256i*)(x+i+40)), a6=_mm256_loadu_si256((__m256i*)(x+i+48)), a7=_mm256_loadu_si256((__m256i*)(x+i+56));
        __m256i b0=_mm256_unpacklo_epi32(a0,a1), b1=_mm256_unpackhi_epi32(a0,a1);
        __m256i b2=_mm256_unpacklo_epi32(a2,a3), b3=_mm256_unpackhi_epi32(a2,a3);
        __m256i b4=_mm256_unpacklo_epi32(a4,a5), b5=_mm256_unpackhi_epi32(a4,a5);
        __m256i b6=_mm256_unpacklo_epi32(a6,a7), b7=_mm256_unpackhi_epi32(a6,a7);
        __m256i c0=_mm256_unpacklo_epi64(b0,b2), c1=_mm256_unpackhi_epi64(b0,b2);
        __m256i c2=_mm256_unpacklo_epi64(b1,b3), c3=_mm256_unpackhi_epi64(b1,b3);
        __m256i c4=_mm256_unpacklo_epi64(b4,b6), c5=_mm256_unpackhi_epi64(b4,b6);
        __m256i c6=_mm256_unpacklo_epi64(b5,b7), c7=_mm256_unpackhi_epi64(b5,b7);
        __m256i t0=_mm256_permute2f128_si256(c0,c4,0x20), t1=_mm256_permute2f128_si256(c1,c5,0x20);
        __m256i t2=_mm256_permute2f128_si256(c2,c6,0x20), t3=_mm256_permute2f128_si256(c3,c7,0x20);
        __m256i t4=_mm256_permute2f128_si256(c0,c4,0x31), t5=_mm256_permute2f128_si256(c1,c5,0x31);
        __m256i t6=_mm256_permute2f128_si256(c2,c6,0x31), t7=_mm256_permute2f128_si256(c3,c7,0x31);
        __m256i d0=addm(t0,t4,modv), d1=addm(t1,t5,modv), d2=addm(t2,t6,modv), d3=addm(t3,t7,modv);
        __m256i e0=subm(t0,t4,modv);
        __m256i e1=mont8(subm(t1,t5,modv),z1m,ninv,modv,shuf);
        __m256i e2=mont8(subm(t2,t6,modv),ivm,ninv,modv,shuf);
        __m256i e3=mont8(subm(t3,t7,modv),z3m,ninv,modv,shuf);
        __m256i s02=addm(d0,d2,modv), t02=subm(d0,d2,modv);
        __m256i s13=addm(d1,d3,modv), t13=mont8(subm(d1,d3,modv),ivm,ninv,modv,shuf);
        __m256i D0=addm(s02,s13,modv), D2=subm(s02,s13,modv);
        __m256i D1=addm(t02,t13,modv), D3=subm(t02,t13,modv);
        __m256i u02=addm(e0,e2,modv), v02=subm(e0,e2,modv);
        __m256i u13=addm(e1,e3,modv), v13=mont8(subm(e1,e3,modv),ivm,ninv,modv,shuf);
        __m256i E0=addm(u02,u13,modv), E2=subm(u02,u13,modv);
        __m256i E1=addm(v02,v13,modv), E3=subm(v02,v13,modv);
        __m256i rb0=_mm256_unpacklo_epi32(D0,E0), rb1=_mm256_unpackhi_epi32(D0,E0);
        __m256i rb2=_mm256_unpacklo_epi32(D1,E1), rb3=_mm256_unpackhi_epi32(D1,E1);
        __m256i rb4=_mm256_unpacklo_epi32(D2,E2), rb5=_mm256_unpackhi_epi32(D2,E2);
        __m256i rb6=_mm256_unpacklo_epi32(D3,E3), rb7=_mm256_unpackhi_epi32(D3,E3);
        __m256i rc0=_mm256_unpacklo_epi64(rb0,rb2), rc1=_mm256_unpackhi_epi64(rb0,rb2);
        __m256i rc2=_mm256_unpacklo_epi64(rb1,rb3), rc3=_mm256_unpackhi_epi64(rb1,rb3);
        __m256i rc4=_mm256_unpacklo_epi64(rb4,rb6), rc5=_mm256_unpackhi_epi64(rb4,rb6);
        __m256i rc6=_mm256_unpacklo_epi64(rb5,rb7), rc7=_mm256_unpackhi_epi64(rb5,rb7);
        __m256i r0=_mm256_permute2f128_si256(rc0,rc4,0x20), r1=_mm256_permute2f128_si256(rc1,rc5,0x20);
        __m256i r2=_mm256_permute2f128_si256(rc2,rc6,0x20), r3=_mm256_permute2f128_si256(rc3,rc7,0x20);
        __m256i r4=_mm256_permute2f128_si256(rc0,rc4,0x31), r5=_mm256_permute2f128_si256(rc1,rc5,0x31);
        __m256i r6=_mm256_permute2f128_si256(rc2,rc6,0x31), r7=_mm256_permute2f128_si256(rc3,rc7,0x31);
        _mm256_storeu_si256((__m256i*)(x+i+0),r0); _mm256_storeu_si256((__m256i*)(x+i+8),r1);
        _mm256_storeu_si256((__m256i*)(x+i+16),r2); _mm256_storeu_si256((__m256i*)(x+i+24),r3);
        _mm256_storeu_si256((__m256i*)(x+i+32),r4); _mm256_storeu_si256((__m256i*)(x+i+40),r5);
        _mm256_storeu_si256((__m256i*)(x+i+48),r6); _mm256_storeu_si256((__m256i*)(x+i+56),r7);
    }
}

static void dif8(u32*x,int n,const u32*tw,u32 Z1,u32 Z3,u32 Iv,u32 Z1p,u32 Z3p,u32 Ivp){
    const __m256i ninv=_mm256_set1_epi32(NINV),modv=_mm256_set1_epi32(MOD);
    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 z1=_mm256_set1_epi32(Z1),z3=_mm256_set1_epi32(Z3),iv=_mm256_set1_epi32(Iv);
    const __m256i z1m=_mm256_set1_epi32(mont_s(Z1,R2)),z3m=_mm256_set1_epi32(mont_s(Z3,R2)),ivm=_mm256_set1_epi32(mont_s(Iv,R2));
    int off=0;
    for(int len=n;len>8;len>>=3){int s=len>>3;
        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 x4=_mm256_loadu_si256((__m256i*)(y+p+4*s));
                __m256i x5=_mm256_loadu_si256((__m256i*)(y+p+5*s));
                __m256i x6=_mm256_loadu_si256((__m256i*)(y+p+6*s));
                __m256i x7=_mm256_loadu_si256((__m256i*)(y+p+7*s));
                const u32* wb=tw+off+p*7;
                __m256i w1=_mm256_loadu_si256((__m256i*)(wb+0));
                __m256i w2=_mm256_loadu_si256((__m256i*)(wb+8));
                __m256i w3=_mm256_loadu_si256((__m256i*)(wb+16));
                __m256i w4=_mm256_loadu_si256((__m256i*)(wb+24));
                __m256i w5=_mm256_loadu_si256((__m256i*)(wb+32));
                __m256i w6=_mm256_loadu_si256((__m256i*)(wb+40));
                __m256i w7=_mm256_loadu_si256((__m256i*)(wb+48));
                __m256i d0=addm(x0,x4,modv),d1=addm(x1,x5,modv),d2=addm(x2,x6,modv),d3=addm(x3,x7,modv);
                __m256i e0=subm(x0,x4,modv);
                __m256i e1=mont8(subm(x1,x5,modv),z1m,ninv,modv,shuf);
                __m256i e2=mont8(subm(x2,x6,modv),ivm,ninv,modv,shuf);
                __m256i e3=mont8(subm(x3,x7,modv),z3m,ninv,modv,shuf);
                __m256i s02=addm(d0,d2,modv),t02=subm(d0,d2,modv);
                __m256i s13=addm(d1,d3,modv),t13=mont8(subm(d1,d3,modv),ivm,ninv,modv,shuf);
                __m256i D0=addm(s02,s13,modv),D2=subm(s02,s13,modv);
                __m256i D1=addm(t02,t13,modv),D3=subm(t02,t13,modv);
                __m256i u02=addm(e0,e2,modv),v02=subm(e0,e2,modv);
                __m256i u13=addm(e1,e3,modv),v13=mont8(subm(e1,e3,modv),ivm,ninv,modv,shuf);
                __m256i E0=addm(u02,u13,modv),E2=subm(u02,u13,modv);
                __m256i E1=addm(v02,v13,modv),E3=subm(v02,v13,modv);
                _mm256_storeu_si256((__m256i*)(y+p),D0);
                _mm256_storeu_si256((__m256i*)(y+p+s),mont8(E0,w1,ninv,modv,shuf));
                _mm256_storeu_si256((__m256i*)(y+p+2*s),mont8(D1,w2,ninv,modv,shuf));
                _mm256_storeu_si256((__m256i*)(y+p+3*s),mont8(E1,w3,ninv,modv,shuf));
                _mm256_storeu_si256((__m256i*)(y+p+4*s),mont8(D2,w4,ninv,modv,shuf));
                _mm256_storeu_si256((__m256i*)(y+p+5*s),mont8(E2,w5,ninv,modv,shuf));
                _mm256_storeu_si256((__m256i*)(y+p+6*s),mont8(D3,w6,ninv,modv,shuf));
                _mm256_storeu_si256((__m256i*)(y+p+7*s),mont8(E3,w7,ninv,modv,shuf));
            }
            for(;p<s;p++){
                u32 x0=y[p],x1=y[p+s],x2=y[p+2*s],x3=y[p+3*s],x4=y[p+4*s],x5=y[p+5*s],x6=y[p+6*s],x7=y[p+7*s];
                const u32* wb=tw+off+p*7;
                u32 d0=add_s(x0,x4),d1=add_s(x1,x5),d2=add_s(x2,x6),d3=add_s(x3,x7);
                u32 e0=sub_s(x0,x4),e1=shoup_mul(sub_s(x1,x5),Z1,Z1p),e2=shoup_mul(sub_s(x2,x6),Iv,Ivp),e3=shoup_mul(sub_s(x3,x7),Z3,Z3p);
                u32 s02=add_s(d0,d2),t02=sub_s(d0,d2);
                u32 s13=add_s(d1,d3),t13=shoup_mul(sub_s(d1,d3),Iv,Ivp);
                u32 D0=add_s(s02,s13),D2=sub_s(s02,s13);
                u32 D1=add_s(t02,t13),D3=sub_s(t02,t13);
                u32 u02=add_s(e0,e2),v02=sub_s(e0,e2);
                u32 u13=add_s(e1,e3),v13=shoup_mul(sub_s(e1,e3),Iv,Ivp);
                u32 E0=add_s(u02,u13),E2=sub_s(u02,u13);
                u32 E1=add_s(v02,v13),E3=sub_s(v02,v13);
                y[p]=D0;
                y[p+s]=mont_s(E0,wb[0]);
                y[p+2*s]=mont_s(D1,wb[8]);
                y[p+3*s]=mont_s(E1,wb[16]);
                y[p+4*s]=mont_s(D2,wb[24]);
                y[p+5*s]=mont_s(E2,wb[32]);
                y[p+6*s]=mont_s(D3,wb[40]);
                y[p+7*s]=mont_s(E3,wb[48]);
            }
        }
        off+=7*s;
    }
    ntt8_vec(x,n,Z1,Z3,Iv,Z1p,Z3p,Ivp);
}
static void dit8(u32*x,int n,const u32*tw,u32 Z1,u32 Z3,u32 Iv,u32 Z1p,u32 Z3p,u32 Ivp){
    const __m256i ninv=_mm256_set1_epi32(NINV),modv=_mm256_set1_epi32(MOD);
    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 z1=_mm256_set1_epi32(Z1),z3=_mm256_set1_epi32(Z3),iv=_mm256_set1_epi32(Iv);
    const __m256i z1m=_mm256_set1_epi32(mont_s(Z1,R2)),z3m=_mm256_set1_epi32(mont_s(Z3,R2)),ivm=_mm256_set1_epi32(mont_s(Iv,R2));
    ntt8_vec(x,n,Z1,Z3,Iv,Z1p,Z3p,Ivp);
    int off=7;
    for(int len=64;len<=n;len<<=3){int s=len>>3;
        for(int i=0;i<n;i+=len){u32*y=x+i;int p=0;
            for(;p+8<=s;p+=8){
                const u32* wb=tw+off+p*7;
                __m256i x0=_mm256_loadu_si256((__m256i*)(y+p));
                __m256i x1=mont8(_mm256_loadu_si256((__m256i*)(y+p+s)),_mm256_loadu_si256((__m256i*)(wb+0)),ninv,modv,shuf);
                __m256i x2=mont8(_mm256_loadu_si256((__m256i*)(y+p+2*s)),_mm256_loadu_si256((__m256i*)(wb+8)),ninv,modv,shuf);
                __m256i x3=mont8(_mm256_loadu_si256((__m256i*)(y+p+3*s)),_mm256_loadu_si256((__m256i*)(wb+16)),ninv,modv,shuf);
                __m256i x4=mont8(_mm256_loadu_si256((__m256i*)(y+p+4*s)),_mm256_loadu_si256((__m256i*)(wb+24)),ninv,modv,shuf);
                __m256i x5=mont8(_mm256_loadu_si256((__m256i*)(y+p+5*s)),_mm256_loadu_si256((__m256i*)(wb+32)),ninv,modv,shuf);
                __m256i x6=mont8(_mm256_loadu_si256((__m256i*)(y+p+6*s)),_mm256_loadu_si256((__m256i*)(wb+40)),ninv,modv,shuf);
                __m256i x7=mont8(_mm256_loadu_si256((__m256i*)(y+p+7*s)),_mm256_loadu_si256((__m256i*)(wb+48)),ninv,modv,shuf);
                __m256i d0=addm(x0,x4,modv),d1=addm(x1,x5,modv),d2=addm(x2,x6,modv),d3=addm(x3,x7,modv);
                __m256i e0=subm(x0,x4,modv);
                __m256i e1=mont8(subm(x1,x5,modv),z1m,ninv,modv,shuf);
                __m256i e2=mont8(subm(x2,x6,modv),ivm,ninv,modv,shuf);
                __m256i e3=mont8(subm(x3,x7,modv),z3m,ninv,modv,shuf);
                __m256i s02=addm(d0,d2,modv),t02=subm(d0,d2,modv);
                __m256i s13=addm(d1,d3,modv),t13=mont8(subm(d1,d3,modv),ivm,ninv,modv,shuf);
                __m256i D0=addm(s02,s13,modv),D2=subm(s02,s13,modv);
                __m256i D1=addm(t02,t13,modv),D3=subm(t02,t13,modv);
                __m256i u02=addm(e0,e2,modv),v02=subm(e0,e2,modv);
                __m256i u13=addm(e1,e3,modv),v13=mont8(subm(e1,e3,modv),ivm,ninv,modv,shuf);
                __m256i E0=addm(u02,u13,modv),E2=subm(u02,u13,modv);
                __m256i E1=addm(v02,v13,modv),E3=subm(v02,v13,modv);
                _mm256_storeu_si256((__m256i*)(y+p),D0);
                _mm256_storeu_si256((__m256i*)(y+p+s),E0);
                _mm256_storeu_si256((__m256i*)(y+p+2*s),D1);
                _mm256_storeu_si256((__m256i*)(y+p+3*s),E1);
                _mm256_storeu_si256((__m256i*)(y+p+4*s),D2);
                _mm256_storeu_si256((__m256i*)(y+p+5*s),E2);
                _mm256_storeu_si256((__m256i*)(y+p+6*s),D3);
                _mm256_storeu_si256((__m256i*)(y+p+7*s),E3);
            }
            for(;p<s;p++){
                u32 x0=y[p],x1=mont_s(y[p+s],tw[off+p*7+0]),x2=mont_s(y[p+2*s],tw[off+p*7+8]),x3=mont_s(y[p+3*s],tw[off+p*7+16]);
                u32 x4=mont_s(y[p+4*s],tw[off+p*7+24]),x5=mont_s(y[p+5*s],tw[off+p*7+32]),x6=mont_s(y[p+6*s],tw[off+p*7+40]),x7=mont_s(y[p+7*s],tw[off+p*7+48]);
                u32 d0=add_s(x0,x4),d1=add_s(x1,x5),d2=add_s(x2,x6),d3=add_s(x3,x7);
                u32 e0=sub_s(x0,x4),e1=shoup_mul(sub_s(x1,x5),Z1,Z1p),e2=shoup_mul(sub_s(x2,x6),Iv,Ivp),e3=shoup_mul(sub_s(x3,x7),Z3,Z3p);
                u32 s02=add_s(d0,d2),t02=sub_s(d0,d2);
                u32 s13=add_s(d1,d3),t13=shoup_mul(sub_s(d1,d3),Iv,Ivp);
                u32 D0=add_s(s02,s13),D2=sub_s(s02,s13);
                u32 D1=add_s(t02,t13),D3=sub_s(t02,t13);
                u32 u02=add_s(e0,e2),v02=sub_s(e0,e2);
                u32 u13=add_s(e1,e3),v13=shoup_mul(sub_s(e1,e3),Iv,Ivp);
                u32 E0=add_s(u02,u13),E2=sub_s(u02,u13);
                u32 E1=add_s(v02,v13),E3=sub_s(v02,v13);
                y[p]=D0;y[p+s]=E0;y[p+2*s]=D1;y[p+3*s]=E1;y[p+4*s]=D2;y[p+5*s]=E2;y[p+6*s]=D3;y[p+7*s]=E3;
            }
        }
        off+=7*s;
    }
}
static void dit8_outer(u32*x,int n,const u32*tw,u32 Z1,u32 Z3,u32 Iv,u32 Z1p,u32 Z3p,u32 Ivp){
    const __m256i ninv=_mm256_set1_epi32(NINV),modv=_mm256_set1_epi32(MOD);
    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 z1=_mm256_set1_epi32(Z1),z3=_mm256_set1_epi32(Z3),iv=_mm256_set1_epi32(Iv);
    const __m256i z1m=_mm256_set1_epi32(mont_s(Z1,R2)),z3m=_mm256_set1_epi32(mont_s(Z3,R2)),ivm=_mm256_set1_epi32(mont_s(Iv,R2));
    int off=7;
    for(int len=64;len<=n;len<<=3){int s=len>>3;
        for(int i=0;i<n;i+=len){u32*y=x+i;int p=0;
            for(;p+8<=s;p+=8){
                const u32* wb=tw+off+p*7;
                __m256i x0=_mm256_loadu_si256((__m256i*)(y+p));
                __m256i x1=mont8(_mm256_loadu_si256((__m256i*)(y+p+s)),_mm256_loadu_si256((__m256i*)(wb+0)),ninv,modv,shuf);
                __m256i x2=mont8(_mm256_loadu_si256((__m256i*)(y+p+2*s)),_mm256_loadu_si256((__m256i*)(wb+8)),ninv,modv,shuf);
                __m256i x3=mont8(_mm256_loadu_si256((__m256i*)(y+p+3*s)),_mm256_loadu_si256((__m256i*)(wb+16)),ninv,modv,shuf);
                __m256i x4=mont8(_mm256_loadu_si256((__m256i*)(y+p+4*s)),_mm256_loadu_si256((__m256i*)(wb+24)),ninv,modv,shuf);
                __m256i x5=mont8(_mm256_loadu_si256((__m256i*)(y+p+5*s)),_mm256_loadu_si256((__m256i*)(wb+32)),ninv,modv,shuf);
                __m256i x6=mont8(_mm256_loadu_si256((__m256i*)(y+p+6*s)),_mm256_loadu_si256((__m256i*)(wb+40)),ninv,modv,shuf);
                __m256i x7=mont8(_mm256_loadu_si256((__m256i*)(y+p+7*s)),_mm256_loadu_si256((__m256i*)(wb+48)),ninv,modv,shuf);
                __m256i d0=addm(x0,x4,modv),d1=addm(x1,x5,modv),d2=addm(x2,x6,modv),d3=addm(x3,x7,modv);
                __m256i e0=subm(x0,x4,modv);
                __m256i e1=mont8(subm(x1,x5,modv),z1m,ninv,modv,shuf);
                __m256i e2=mont8(subm(x2,x6,modv),ivm,ninv,modv,shuf);
                __m256i e3=mont8(subm(x3,x7,modv),z3m,ninv,modv,shuf);
                __m256i s02=addm(d0,d2,modv),t02=subm(d0,d2,modv);
                __m256i s13=addm(d1,d3,modv),t13=mont8(subm(d1,d3,modv),ivm,ninv,modv,shuf);
                __m256i D0=addm(s02,s13,modv),D2=subm(s02,s13,modv);
                __m256i D1=addm(t02,t13,modv),D3=subm(t02,t13,modv);
                __m256i u02=addm(e0,e2,modv),v02=subm(e0,e2,modv);
                __m256i u13=addm(e1,e3,modv),v13=mont8(subm(e1,e3,modv),ivm,ninv,modv,shuf);
                __m256i E0=addm(u02,u13,modv),E2=subm(u02,u13,modv);
                __m256i E1=addm(v02,v13,modv),E3=subm(v02,v13,modv);
                _mm256_storeu_si256((__m256i*)(y+p),D0);
                _mm256_storeu_si256((__m256i*)(y+p+s),E0);
                _mm256_storeu_si256((__m256i*)(y+p+2*s),D1);
                _mm256_storeu_si256((__m256i*)(y+p+3*s),E1);
                _mm256_storeu_si256((__m256i*)(y+p+4*s),D2);
                _mm256_storeu_si256((__m256i*)(y+p+5*s),E2);
                _mm256_storeu_si256((__m256i*)(y+p+6*s),D3);
                _mm256_storeu_si256((__m256i*)(y+p+7*s),E3);
            }
            for(;p<s;p++){
                u32 x0=y[p],x1=mont_s(y[p+s],tw[off+p*7+0]),x2=mont_s(y[p+2*s],tw[off+p*7+8]),x3=mont_s(y[p+3*s],tw[off+p*7+16]);
                u32 x4=mont_s(y[p+4*s],tw[off+p*7+24]),x5=mont_s(y[p+5*s],tw[off+p*7+32]),x6=mont_s(y[p+6*s],tw[off+p*7+40]),x7=mont_s(y[p+7*s],tw[off+p*7+48]);
                u32 d0=add_s(x0,x4),d1=add_s(x1,x5),d2=add_s(x2,x6),d3=add_s(x3,x7);
                u32 e0=sub_s(x0,x4),e1=shoup_mul(sub_s(x1,x5),Z1,Z1p),e2=shoup_mul(sub_s(x2,x6),Iv,Ivp),e3=shoup_mul(sub_s(x3,x7),Z3,Z3p);
                u32 s02=add_s(d0,d2),t02=sub_s(d0,d2);
                u32 s13=add_s(d1,d3),t13=shoup_mul(sub_s(d1,d3),Iv,Ivp);
                u32 D0=add_s(s02,s13),D2=sub_s(s02,s13);
                u32 D1=add_s(t02,t13),D3=sub_s(t02,t13);
                u32 u02=add_s(e0,e2),v02=sub_s(e0,e2);
                u32 u13=add_s(e1,e3),v13=shoup_mul(sub_s(e1,e3),Iv,Ivp);
                u32 E0=add_s(u02,u13),E2=sub_s(u02,u13);
                u32 E1=add_s(v02,v13),E3=sub_s(v02,v13);
                y[p]=D0;y[p+s]=E0;y[p+2*s]=D1;y[p+3*s]=E1;y[p+4*s]=D2;y[p+5*s]=E2;y[p+6*s]=D3;y[p+7*s]=E3;
            }
        }
        off+=7*s;
    }
}

// Fused: pointwise-multiply B by A (mont8) then the 8-point leaf DFT, in ONE pass.
// Replaces the separate pointwise pass (saves its 16MB write+re-read traffic).
static void pointwise_and_leaf8(u32*B,u32*A,int n,u32 Z1,u32 Z3,u32 Iv,u32 Z1p,u32 Z3p,u32 Ivp){
    const __m256i ninv=_mm256_set1_epi32(NINV),modv=_mm256_set1_epi32(MOD);
    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 z1=_mm256_set1_epi32(Z1),z3=_mm256_set1_epi32(Z3),iv=_mm256_set1_epi32(Iv);
    const __m256i z1m=_mm256_set1_epi32(mont_s(Z1,R2)),z3m=_mm256_set1_epi32(mont_s(Z3,R2)),ivm=_mm256_set1_epi32(mont_s(Iv,R2));
    for(int i=0;i<n;i+=64){
        __m256i y0=_mm256_loadu_si256((__m256i*)(A+i+0)), y1=_mm256_loadu_si256((__m256i*)(A+i+8)), y2=_mm256_loadu_si256((__m256i*)(A+i+16)), y3=_mm256_loadu_si256((__m256i*)(A+i+24));
        __m256i y4=_mm256_loadu_si256((__m256i*)(A+i+32)), y5=_mm256_loadu_si256((__m256i*)(A+i+40)), y6=_mm256_loadu_si256((__m256i*)(A+i+48)), y7=_mm256_loadu_si256((__m256i*)(A+i+56));
        __m256i a0=mont8(_mm256_loadu_si256((__m256i*)(B+i+0)),y0,ninv,modv,shuf);
        __m256i a1=mont8(_mm256_loadu_si256((__m256i*)(B+i+8)),y1,ninv,modv,shuf);
        __m256i a2=mont8(_mm256_loadu_si256((__m256i*)(B+i+16)),y2,ninv,modv,shuf);
        __m256i a3=mont8(_mm256_loadu_si256((__m256i*)(B+i+24)),y3,ninv,modv,shuf);
        __m256i a4=mont8(_mm256_loadu_si256((__m256i*)(B+i+32)),y4,ninv,modv,shuf);
        __m256i a5=mont8(_mm256_loadu_si256((__m256i*)(B+i+40)),y5,ninv,modv,shuf);
        __m256i a6=mont8(_mm256_loadu_si256((__m256i*)(B+i+48)),y6,ninv,modv,shuf);
        __m256i a7=mont8(_mm256_loadu_si256((__m256i*)(B+i+56)),y7,ninv,modv,shuf);
        __m256i b0=_mm256_unpacklo_epi32(a0,a1), b1=_mm256_unpackhi_epi32(a0,a1);
        __m256i b2=_mm256_unpacklo_epi32(a2,a3), b3=_mm256_unpackhi_epi32(a2,a3);
        __m256i b4=_mm256_unpacklo_epi32(a4,a5), b5=_mm256_unpackhi_epi32(a4,a5);
        __m256i b6=_mm256_unpacklo_epi32(a6,a7), b7=_mm256_unpackhi_epi32(a6,a7);
        __m256i c0=_mm256_unpacklo_epi64(b0,b2), c1=_mm256_unpackhi_epi64(b0,b2);
        __m256i c2=_mm256_unpacklo_epi64(b1,b3), c3=_mm256_unpackhi_epi64(b1,b3);
        __m256i c4=_mm256_unpacklo_epi64(b4,b6), c5=_mm256_unpackhi_epi64(b4,b6);
        __m256i c6=_mm256_unpacklo_epi64(b5,b7), c7=_mm256_unpackhi_epi64(b5,b7);
        __m256i t0=_mm256_permute2f128_si256(c0,c4,0x20), t1=_mm256_permute2f128_si256(c1,c5,0x20);
        __m256i t2=_mm256_permute2f128_si256(c2,c6,0x20), t3=_mm256_permute2f128_si256(c3,c7,0x20);
        __m256i t4=_mm256_permute2f128_si256(c0,c4,0x31), t5=_mm256_permute2f128_si256(c1,c5,0x31);
        __m256i t6=_mm256_permute2f128_si256(c2,c6,0x31), t7=_mm256_permute2f128_si256(c3,c7,0x31);
        __m256i d0=addm(t0,t4,modv), d1=addm(t1,t5,modv), d2=addm(t2,t6,modv), d3=addm(t3,t7,modv);
        __m256i e0=subm(t0,t4,modv);
        __m256i e1=mont8(subm(t1,t5,modv),z1m,ninv,modv,shuf);
        __m256i e2=mont8(subm(t2,t6,modv),ivm,ninv,modv,shuf);
        __m256i e3=mont8(subm(t3,t7,modv),z3m,ninv,modv,shuf);
        __m256i s02=addm(d0,d2,modv), t02=subm(d0,d2,modv);
        __m256i s13=addm(d1,d3,modv), t13=mont8(subm(d1,d3,modv),ivm,ninv,modv,shuf);
        __m256i D0=addm(s02,s13,modv), D2=subm(s02,s13,modv);
        __m256i D1=addm(t02,t13,modv), D3=subm(t02,t13,modv);
        __m256i u02=addm(e0,e2,modv), v02=subm(e0,e2,modv);
        __m256i u13=addm(e1,e3,modv), v13=mont8(subm(e1,e3,modv),ivm,ninv,modv,shuf);
        __m256i E0=addm(u02,u13,modv), E2=subm(u02,u13,modv);
        __m256i E1=addm(v02,v13,modv), E3=subm(v02,v13,modv);
        __m256i rb0=_mm256_unpacklo_epi32(D0,E0), rb1=_mm256_unpackhi_epi32(D0,E0);
        __m256i rb2=_mm256_unpacklo_epi32(D1,E1), rb3=_mm256_unpackhi_epi32(D1,E1);
        __m256i rb4=_mm256_unpacklo_epi32(D2,E2), rb5=_mm256_unpackhi_epi32(D2,E2);
        __m256i rb6=_mm256_unpacklo_epi32(D3,E3), rb7=_mm256_unpackhi_epi32(D3,E3);
        __m256i rc0=_mm256_unpacklo_epi64(rb0,rb2), rc1=_mm256_unpackhi_epi64(rb0,rb2);
        __m256i rc2=_mm256_unpacklo_epi64(rb1,rb3), rc3=_mm256_unpackhi_epi64(rb1,rb3);
        __m256i rc4=_mm256_unpacklo_epi64(rb4,rb6), rc5=_mm256_unpackhi_epi64(rb4,rb6);
        __m256i rc6=_mm256_unpacklo_epi64(rb5,rb7), rc7=_mm256_unpackhi_epi64(rb5,rb7);
        __m256i r0=_mm256_permute2f128_si256(rc0,rc4,0x20), r1=_mm256_permute2f128_si256(rc1,rc5,0x20);
        __m256i r2=_mm256_permute2f128_si256(rc2,rc6,0x20), r3=_mm256_permute2f128_si256(rc3,rc7,0x20);
        __m256i r4=_mm256_permute2f128_si256(rc0,rc4,0x31), r5=_mm256_permute2f128_si256(rc1,rc5,0x31);
        __m256i r6=_mm256_permute2f128_si256(rc2,rc6,0x31), r7=_mm256_permute2f128_si256(rc3,rc7,0x31);
        _mm256_storeu_si256((__m256i*)(B+i+0),r0); _mm256_storeu_si256((__m256i*)(B+i+8),r1);
        _mm256_storeu_si256((__m256i*)(B+i+16),r2); _mm256_storeu_si256((__m256i*)(B+i+24),r3);
        _mm256_storeu_si256((__m256i*)(B+i+32),r4); _mm256_storeu_si256((__m256i*)(B+i+40),r5);
        _mm256_storeu_si256((__m256i*)(B+i+48),r6); _mm256_storeu_si256((__m256i*)(B+i+56),r7);
    }
}

void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c) {
    int size=1; while(size<n+m+1)size<<=1;
    {int s=size,lg=0;while(s>1){s>>=1;lg++;}while(lg%3)size<<=1,lg++;} // pad to power of 8
    u32 ninv_scale=modpow((u32)size,MOD-2);
    u32 md[10], mda[10];
    for(int i=0;i<10;i++){md[i]=mont_s((u32)i,R2); mda[i]=(u32)(((u64)i*ninv_scale)%MOD);}

    // 1. copy b -> buffer (mont form). buffer is static BSS (zero), single call, so
    //    no explicit pad needed. Must happen BEFORE padding a (overwrites b's region).
    // vectorized: md[b[i]] = b[i]*R mod p, b[i]<10 so b[i]*R < 10*R < 3p (fits u32).
    { u32 R=md[1];
      __m256i Rv=_mm256_set1_epi32(R), Pv=_mm256_set1_epi32(MOD);
      int i=0;
      for(;i+8<=m+1;i+=8){
        __m256i b8=_mm256_loadu_si256((__m256i*)(b+i));
        __m256i r=_mm256_mullo_epi32(b8,Rv);
        r=_mm256_min_epu32(r,_mm256_sub_epi32(r,Pv));
        r=_mm256_min_epu32(r,_mm256_sub_epi32(r,Pv));
        _mm256_storeu_si256((__m256i*)(buffer+i),r);
      }
      for(;i<=m;i++)buffer[i]=md[b[i]]; }
    // 2. scale a in place + pad. Only b's old region (a[n+1..n+m+1]) needs zeroing;
    //    c's region (a[n+m+2..size-1]) is already 0 (c initialized to 0 by tasklib).
    { u32 wsc=ninv_scale, wpsc=shp_quot(wsc);
      __m256i wv=_mm256_set1_epi32(wsc), wpv=_mm256_set1_epi32(wpsc), Pv=_mm256_set1_epi32(MOD);
      int i=0;
      for(;i+8<=n+1;i+=8){
        __m256i a8=_mm256_loadu_si256((__m256i*)(a+i));
        __m256i lo=_mm256_mullo_epi32(a8,wv);
        __m256i qe=_mm256_mul_epu32(a8,wpv);
        __m256i qo=_mm256_mul_epu32(_mm256_srli_epi64(a8,32),_mm256_srli_epi64(wpv,32));
        __m256i qe2=_mm256_srli_epi64(qe,32);
        __m256i qo2=_mm256_srli_epi64(qo,32);
        __m256i q=_mm256_or_si256(qe2,_mm256_slli_epi64(qo2,32));
        __m256i r=_mm256_sub_epi32(lo,_mm256_mullo_epi32(q,Pv));
        r=_mm256_min_epu32(r,_mm256_sub_epi32(r,Pv));
        _mm256_storeu_si256((__m256i*)(a+i),r);
      }
      for(;i<=n;i++)a[i]=mda[a[i]]; }
    { __m256i z=_mm256_setzero_si256(); int i=n+1;
      for(;i+8<=n+m+2;i+=8)_mm256_storeu_si256((__m256i*)(a+i),z);
      for(;i<=n+m+1;i++)a[i]=0; }

    u32 w=modpow(G,(MOD-1)/(u32)size), iw=modpow(w,MOD-2);
    u32 Z1=modpow(w,(u32)(size/8)), Z3=modpow(w,3u*(u32)(size/8)), Iv=IM;
    u32 Z1i=modpow(iw,(u32)(size/8)), Z3i=modpow(iw,3u*(u32)(size/8)), Ivi=IINV;
    u32 Z1p=shp_quot(Z1), Z3p=shp_quot(Z3), Ivp=shp_quot(Iv);
    u32 Z1ip=shp_quot(Z1i), Z3ip=shp_quot(Z3i), Ivip=shp_quot(Ivi);

    // forward twiddles -> tw
    {const __m256i ninv=_mm256_set1_epi32(NINV),modv=_mm256_set1_epi32(MOD);
     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);
     int off=0;
     for(int len=size;len>=8;len>>=3){int s=len>>3,step=size/len;
       u32 ws=modpow(w,(u64)step); u32 wm=mont_s(ws,R2);
       tmp[0]=mont_s(1,R2);
       for(int p=1;p<8 && p<s;p++) tmp[p]=mont_s(tmp[p-1],wm);
       if(s>8){
         __m256i w8=_mm256_set1_epi32(mont_s(tmp[7],wm));
         for(int p=8;p<s;p+=8) _mm256_storeu_si256((__m256i*)(tmp+p),mont8(_mm256_loadu_si256((__m256i*)(tmp+p-8)),w8,ninv,modv,shuf));
       }
       for(int p=0;p<s;p+=8){
         __m256i t1=_mm256_loadu_si256((__m256i*)(tmp+p));
         __m256i t2=mont8(t1,t1,ninv,modv,shuf);
         __m256i t4=mont8(t2,t2,ninv,modv,shuf);
         __m256i t3=mont8(t2,t1,ninv,modv,shuf);
         __m256i t6=mont8(t4,t2,ninv,modv,shuf);
         __m256i t5=mont8(t4,t1,ninv,modv,shuf);
         __m256i t7=mont8(t6,t1,ninv,modv,shuf);
         u32* wb=tw+off+p*7;
         _mm256_storeu_si256((__m256i*)(wb+0),t1);
         _mm256_storeu_si256((__m256i*)(wb+8),t2);
         _mm256_storeu_si256((__m256i*)(wb+16),t3);
         _mm256_storeu_si256((__m256i*)(wb+24),t4);
         _mm256_storeu_si256((__m256i*)(wb+32),t5);
         _mm256_storeu_si256((__m256i*)(wb+40),t6);
         _mm256_storeu_si256((__m256i*)(wb+48),t7);}
       off+=7*s;}}

    dif8(a,size,tw,Z1,Z3,Iv,Z1p,Z3p,Ivp);
    dif8(buffer,size,tw,Z1,Z3,Iv,Z1p,Z3p,Ivp);

    // inverse twiddles -> tw (overwrite forward table)
    {const __m256i ninv=_mm256_set1_epi32(NINV),modv=_mm256_set1_epi32(MOD);
     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);
     int off=0;
     for(int len=8;len<=size;len<<=3){int s=len>>3,step=size/len;
       u32 ws=modpow(iw,(u64)step); u32 wm=mont_s(ws,R2);
       tmp[0]=mont_s(1,R2);
       for(int p=1;p<8 && p<s;p++) tmp[p]=mont_s(tmp[p-1],wm);
       if(s>8){
         __m256i w8=_mm256_set1_epi32(mont_s(tmp[7],wm));
         for(int p=8;p<s;p+=8) _mm256_storeu_si256((__m256i*)(tmp+p),mont8(_mm256_loadu_si256((__m256i*)(tmp+p-8)),w8,ninv,modv,shuf));
       }
       for(int p=0;p<s;p+=8){
         __m256i t1=_mm256_loadu_si256((__m256i*)(tmp+p));
         __m256i t2=mont8(t1,t1,ninv,modv,shuf);
         __m256i t4=mont8(t2,t2,ninv,modv,shuf);
         __m256i t3=mont8(t2,t1,ninv,modv,shuf);
         __m256i t6=mont8(t4,t2,ninv,modv,shuf);
         __m256i t5=mont8(t4,t1,ninv,modv,shuf);
         __m256i t7=mont8(t6,t1,ninv,modv,shuf);
         u32* wb=tw+off+p*7;
         _mm256_storeu_si256((__m256i*)(wb+0),t1);
         _mm256_storeu_si256((__m256i*)(wb+8),t2);
         _mm256_storeu_si256((__m256i*)(wb+16),t3);
         _mm256_storeu_si256((__m256i*)(wb+24),t4);
         _mm256_storeu_si256((__m256i*)(wb+32),t5);
         _mm256_storeu_si256((__m256i*)(wb+40),t6);
         _mm256_storeu_si256((__m256i*)(wb+48),t7);}
       off+=7*s;}}

    pointwise_and_leaf8(buffer,a,size,Z1i,Z3i,Ivi,Z1ip,Z3ip,Ivip);
    dit8_outer(buffer,size,tw,Z1i,Z3i,Ivi,Z1ip,Z3ip,Ivip);
    int outn=n+m+1;
    { int i=0;
      while(((unsigned long long)(c+i)&31ULL) && i<outn) c[i]=buffer[i], i++;
      for(;i+8<=outn;i+=8){__m256i v=_mm256_loadu_si256((__m256i*)(buffer+i));_mm256_stream_si256((__m256i*)(c+i),v);}
      for(;i<outn;i++) c[i]=buffer[i]; }
}

评测结果