代码

// DUCK_PUBLIC_SURROGATE_BENCHMARK
// Independent public uint32 middle products; complete uint64-Nuss/Kara equality.
#pragma once
#include <bits/stdc++.h>
#include <immintrin.h>
#pragma GCC push_options
#pragma GCC target("avx2")
#pragma GCC optimize("O3")
namespace nuss {
using u32=uint32_t;using u64=uint64_t;
template<class T>struct Storage{T*p;unsigned n;explicit Storage(unsigned count):p(nullptr),n(count){void*q=nullptr;if(posix_memalign(&q,64,size_t(n)*sizeof(T)))std::abort();p=(T*)q;for(unsigned i=0;i<n;++i)new(p+i)T{};}~Storage(){free(p);}Storage(const Storage&)=delete;T&operator[](unsigned i){return p[i];}const T&operator[](unsigned i)const{return p[i];}};
struct V4{__m256i v;V4():v(_mm256_setzero_si256()){}explicit V4(__m256i z):v(z){}friend V4 operator+(V4 a,V4 b){return V4(_mm256_add_epi64(a.v,b.v));}friend V4 operator-(V4 a,V4 b){return V4(_mm256_sub_epi64(a.v,b.v));}};
inline u64 mul(u64 a,u64 b){return a*b;}
inline V4 mul(V4 a,u64 b){const __m256i bv=_mm256_set1_epi64x(b),sw=_mm256_set1_epi64x((b>>32)|(b<<32));__m256i low=_mm256_mul_epu32(a.v,bv),cross=_mm256_mullo_epi32(a.v,sw);cross=_mm256_add_epi64(cross,_mm256_srli_epi64(cross,32));return V4(_mm256_add_epi64(low,_mm256_slli_epi64(cross,32)));}
template<unsigned N,class V>__attribute__((always_inline))inline void small_product(const V*a,const u64*b,V*out){if constexpr(N==1){out[0]=mul(a[0],b[0]);out[1]=V{};}else{constexpr unsigned S=N/2;V sa[S],cross[2*S];u64 sb[S];for(unsigned i=0;i<S;++i){sa[i]=a[i]+a[S+i];sb[i]=b[i]+b[S+i];}small_product<S>(a,b,out);small_product<S>(a+S,b+S,out+2*S);small_product<S>(sa,sb,cross);for(unsigned i=0;i<2*S;++i)cross[i]=cross[i]-out[i]-out[2*S+i];for(unsigned i=0;i<2*S;++i)out[i+S]=out[i+S]+cross[i];}}
template<class V>void full_product(const V*a,const u64*b,V*out,unsigned n,V*scratch,u64*bs){
 if(n<=8){if(n==1)small_product<1>(a,b,out);else if(n==2)small_product<2>(a,b,out);else if(n==4)small_product<4>(a,b,out);else small_product<8>(a,b,out);return;}
 unsigned s=n/2;V*sumA=scratch,*cross=sumA+s,*next=cross+2*s;u64*sumB=bs,*nextB=bs+s;
 for(unsigned i=0;i<s;++i){sumA[i]=a[i]+a[s+i];sumB[i]=b[i]+b[s+i];}
 full_product(a,b,out,s,next,nextB);full_product(a+s,b+s,out+2*s,s,next,nextB);full_product(sumA,sumB,cross,s,next,nextB);
 for(unsigned i=0;i<2*s;++i)cross[i]=cross[i]-out[i]-out[i+2*s];for(unsigned i=0;i<2*s;++i)out[i+s]=out[i+s]+cross[i];
}
template<class V>void rotate(const V*src,V*dst,unsigned len,unsigned shift){shift%=2*len;bool neg=shift>=len;if(neg)shift-=len;for(unsigned i=0;i<len-shift;++i)dst[i+shift]=neg?V{}-src[i]:src[i];for(unsigned i=len-shift;i<len;++i)dst[i+shift-len]=neg?src[i]:V{}-src[i];}
template<class V>void fft(V*a,unsigned count,unsigned len,bool inverse){Storage<V>tmp(len);if(!inverse){for(unsigned rank=count;rank>=2;rank/=2){unsigned half=rank/2,unit=2*len/rank;for(unsigned block=0;block<count;block+=rank)for(unsigned j=0;j<half;++j){V*x=a+(block+j)*len,*y=a+(block+j+half)*len;for(unsigned t=0;t<len;++t){tmp[t]=x[t]-y[t];x[t]=x[t]+y[t];}rotate(tmp.p,y,len,unit*j);}}}else for(unsigned rank=2;rank<=count;rank*=2){unsigned half=rank/2,unit=2*len/rank;for(unsigned block=0;block<count;block+=rank)for(unsigned j=0;j<half;++j){V*x=a+(block+j)*len,*y=a+(block+j+half)*len;rotate(y,tmp.p,len,(2*len-unit*j)%(2*len));for(unsigned t=0;t<len;++t){V old=x[t];x[t]=old+tmp[t];y[t]=old-tmp[t];}}}}
struct Shape{unsigned h,K,M,N,L,logN,logM;explicit Shape(unsigned size):h(size){assert(h>=4&&(h&(h-1))==0);unsigned lg=__builtin_ctz(h);K=1u<<((lg+1)/2);M=h/K;N=2*K;L=2*M;logN=__builtin_ctz(N);logM=__builtin_ctz(M);assert(2*L%N==0);}};
class Plan64{
public:
 Shape s;Storage<u64>bhat;explicit Plan64(const u32*b,unsigned h):s(h),bhat(s.N*s.L){for(unsigned j=0;j<s.K;++j)for(unsigned t=0;t<s.M;++t)bhat[j*s.L+t]=b[s.h-1-j*s.M-t];fft(bhat.p,s.N,s.L,false);}
 bool apply(const u32*a,u32*out,bool check_divisibility=false)const{Storage<u64>x(s.N*s.L),prod(2*s.L),scratch(4*s.L),bs(2*s.L);for(unsigned j=0;j<s.N;++j)for(unsigned t=0;t<s.M;++t)x[j*s.L+t]=a[j*s.M+t];fft(x.p,s.N,s.L,false);for(unsigned j=0;j<s.N;++j){full_product(x.p+j*s.L,bhat.p+j*s.L,prod.p,s.L,scratch.p,bs.p);for(unsigned t=0;t<s.L;++t)x[j*s.L+t]=prod[t]-prod[t+s.L];}fft(x.p,s.N,s.L,true);if(check_divisibility)for(unsigned i=0;i<s.N*s.L;++i)if(x[i]&(s.N-1))return false;for(unsigned t=0;t<s.h;++t){unsigned k=s.h-1+t,j=k>>s.logM,r=k&(s.M-1),prev=(j+s.N-1)&(s.N-1);out[t]=u32(x[j*s.L+r]>>s.logN)+u32(x[prev*s.L+r+s.M]>>s.logN);}return true;}
};
class Plan4{
 const Plan64&plan;Storage<V4>x,prod,scratch;Storage<u64>bs;
public:
 explicit Plan4(const Plan64&p):plan(p),x(p.s.N*p.s.L),prod(2*p.s.L),scratch(4*p.s.L),bs(2*p.s.L){}
 bool apply(const u32*a,u32*out,unsigned stride=4,bool check_divisibility=false){const auto&s=plan.s;for(unsigned i=0;i<s.N*s.L;++i)x[i]=V4{};for(unsigned j=0;j<s.N;++j)for(unsigned t=0;t<s.M;++t)x[j*s.L+t]=V4(_mm256_cvtepu32_epi64(_mm_loadu_si128((const __m128i*)(a+stride*(j*s.M+t)))));fft(x.p,s.N,s.L,false);for(unsigned j=0;j<s.N;++j){full_product(x.p+j*s.L,plan.bhat.p+j*s.L,prod.p,s.L,scratch.p,bs.p);for(unsigned t=0;t<s.L;++t)x[j*s.L+t]=prod[t]-prod[t+s.L];}fft(x.p,s.N,s.L,true);if(check_divisibility){__m256i mask=_mm256_set1_epi64x(s.N-1);for(unsigned i=0;i<s.N*s.L;++i)if(!_mm256_testz_si256(x[i].v,mask))return false;}const __m128i sh=_mm_cvtsi32_si128(s.logN);const __m256i perm=_mm256_setr_epi32(0,2,4,6,1,3,5,7);for(unsigned t=0;t<s.h;++t){unsigned k=s.h-1+t,j=k>>s.logM,r=k&(s.M-1),prev=(j+s.N-1)&(s.N-1);__m256i v=_mm256_srl_epi64(_mm256_add_epi64(x[j*s.L+r].v,x[prev*s.L+r+s.M].v),sh);_mm_storeu_si128((__m128i*)(out+stride*t),_mm256_castsi256_si128(_mm256_permutevar8x32_epi32(v,perm)));}return true;}
};
}
#pragma GCC pop_options

#pragma once
#pragma once
#pragma once
#pragma once
#pragma GCC target("avx2")
#pragma GCC optimize("O3")
#include <bits/stdc++.h>
#include <immintrin.h>
namespace adaptive_mod {using u32=uint32_t;
inline void direct_range(unsigned n,const unsigned* __restrict a,const unsigned* __restrict b,unsigned* __restrict c,unsigned l,unsigned r){
 if(n<64){for(unsigned j=l+1;j<std::min(r,n-l);++j)for(unsigned start=j;start+l<n;start+=j)for(unsigned k=l;k<std::min(r,std::min(j,n-start));++k)c[k]+=a[start+k]*b[j];return;}
 for(unsigned j=l+1;j<std::min(r,65u);++j){u32 sum[64]={};for(unsigned start=j;start+l<n;start+=j){unsigned end=std::min(j,n-start);for(unsigned k=l;k<end;++k)sum[k]+=a[start+k];}for(unsigned k=l;k<j;++k)c[k]+=sum[k]*b[j];}
 const unsigned limit=std::min(r,n/2+1),jlo=std::max(l+1,65u);if(jlo>=r)return;
 for(unsigned k0=l;k0<limit;k0+=32){__m256i z[4],idx[4],endmask[4];for(unsigned v=0;v<4;++v){z[v]=_mm256_loadu_si256((const __m256i*)(c+k0+8*v));idx[v]=_mm256_setr_epi32(k0+8*v,k0+8*v+1,k0+8*v+2,k0+8*v+3,k0+8*v+4,k0+8*v+5,k0+8*v+6,k0+8*v+7);endmask[v]=_mm256_cmpgt_epi32(_mm256_set1_epi32(r),idx[v]);}
  unsigned first=std::max(jlo,k0+1),first_full=std::max(jlo,k0+32),qmax=(n-1-k0)/first;
  for(unsigned q=1;q<=qmax;++q){unsigned last=std::min(r-1,(n-1-k0)/q),full_last=std::min(r-1,(n-k0-32)/q),j=first;
   for(;j<=last&&j<first_full;++j){__m256i bj=_mm256_set1_epi32(b[j]);unsigned start=q*j+k0;for(unsigned v=0;v<4;++v){__m256i mask=_mm256_and_si256(endmask[v],_mm256_and_si256(_mm256_cmpgt_epi32(_mm256_set1_epi32(j),idx[v]),_mm256_cmpgt_epi32(_mm256_set1_epi32(n-q*j),idx[v])));__m256i x=start+8*v<n?_mm256_maskload_epi32((const int*)(a+start+8*v),mask):_mm256_setzero_si256();z[v]=_mm256_add_epi32(z[v],_mm256_mullo_epi32(x,bj));}}
   for(;j<=full_last;++j){__m256i bj=_mm256_set1_epi32(b[j]);const u32*x=a+q*j+k0;for(unsigned v=0;v<4;++v){__m256i value=_mm256_loadu_si256((const __m256i*)(x+8*v));if(k0+32>r)value=_mm256_and_si256(value,endmask[v]);z[v]=_mm256_add_epi32(z[v],_mm256_mullo_epi32(value,bj));}}
   for(;j<=last;++j){__m256i bj=_mm256_set1_epi32(b[j]);unsigned start=q*j+k0;for(unsigned v=0;v<4;++v){__m256i mask=_mm256_and_si256(endmask[v],_mm256_and_si256(_mm256_cmpgt_epi32(_mm256_set1_epi32(j),idx[v]),_mm256_cmpgt_epi32(_mm256_set1_epi32(n-q*j),idx[v])));__m256i x=start+8*v<n?_mm256_maskload_epi32((const int*)(a+start+8*v),mask):_mm256_setzero_si256();z[v]=_mm256_add_epi32(z[v],_mm256_mullo_epi32(x,bj));}}
  }
  for(unsigned v=0;v<4;++v)_mm256_storeu_si256((__m256i*)(c+k0+8*v),z[v]);
 }
}
}

#pragma once
#include <bits/stdc++.h>
#include <immintrin.h>
#pragma GCC optimize("O3")
#pragma GCC target("avx2")
namespace mod_karatsuba {
using u32=uint32_t;using u64=uint64_t;
// out[t] = sum_i a[i+t]*b[i], 0<=t<h. a has at least 2h words.
__attribute__((noinline)) inline void middle(const u32* __restrict a,const u32* __restrict b,u32* __restrict out,unsigned h,u32*scratch){
 if(h<=32){
  if(h<8){for(unsigned t=0;t<h;++t){u32 z=0;for(unsigned i=0;i<h;++i)z+=a[i+t]*b[i];out[t]=z;}return;}
  for(unsigned t=0;t<h;t+=8){__m256i z=_mm256_setzero_si256();for(unsigned i=0;i<h;++i)z=_mm256_add_epi32(z,_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(a+i+t)),_mm256_set1_epi32(b[i])));_mm256_storeu_si256((__m256i*)(out+t),z);}return;
 }
 const unsigned s=h/2;u32*x=scratch,*y=x+h,*bs=y+h,*tmp=bs+s,*next=tmp+s;
 for(unsigned i=0;i<h;++i){x[i]=a[i]-a[s+i];y[i]=a[h+i]-a[s+i];}
 for(unsigned i=0;i<s;++i)bs[i]=b[i]+b[s+i];
 middle(x,b,out,s,next);middle(y,b+s,out+s,s,next);middle(a+s,bs,tmp,s,next);
 for(unsigned i=0;i<s;++i){out[i]+=tmp[i];out[s+i]+=tmp[i];}
}
__attribute__((noinline)) inline void cross(unsigned n,const u32*a,const u32*b,u32*c,unsigned l,unsigned r,unsigned qlo=1,unsigned qhi=~0u){
 const unsigned w=(r-l)/2,m=l+w,jend=std::min(r,n);if(!w||m>=jend||l+m>=n)return;
 qhi=std::min(qhi,(n-1-l)/m);
 for(unsigned q=qlo;q<=qhi;++q){unsigned wanted=std::min(w,n-q*m-l),h=1;while(h<(wanted+q-1)/q)h*=2;
  std::vector<u32>bb(h),aa(2*h+8),z(h),scratch(6*h+32);
  for(unsigned jb=m;jb<jend;jb+=h){unsigned count=std::min(h,jend-jb);std::fill(bb.begin(),bb.end(),0);std::copy(b+jb,b+jb+count,bb.data());
   for(unsigned rho=0;rho<std::min(q,wanted);++rho){unsigned len=(wanted-1-rho)/q+1;u64 start=u64(q)*jb+l+rho;
    for(unsigned i=0;i<2*h;++i){u64 at=start+u64(q)*i;aa[i]=at<n?a[at]:0;}
    middle(aa.data(),bb.data(),z.data(),h,scratch.data());
    for(unsigned t=0;t<len;++t)c[l+rho+q*t]+=z[t];
   }
  }
 }
}
__attribute__((noinline)) inline void triangle(unsigned n,const u32*a,const u32*b,u32*c,unsigned l,unsigned r){
 if(r-l<=32){for(unsigned j=l+1;j<std::min(r,n-l);++j)for(unsigned st=j;st+l<n;st+=j)for(unsigned k=l;k<std::min(r,std::min(j,n-st));++k)c[k]+=a[st+k]*b[j];return;}
 if(2ull*l+1>=n)return;unsigned m=l+(r-l)/2;triangle(n,a,b,c,l,m);triangle(n,a,b,c,m,r);cross(n,a,b,c,l,r);
}
}

namespace mod_karatsuba {
__attribute__((noinline)) inline void direct_cross(unsigned n,const u32*a,const u32*b,u32*c,unsigned l,unsigned r,unsigned qlo,unsigned qhi){
 unsigned m=l+(r-l)/2;for(unsigned k0=l;k0<m;k0+=32){__m256i z[4],mask[4];for(unsigned v=0;v<4;++v){z[v]=_mm256_loadu_si256((const __m256i*)(c+k0+8*v));mask[v]=_mm256_setr_epi32(k0+v*8<m?~0u:0u,k0+v*8+1<m?~0u:0u,k0+v*8+2<m?~0u:0u,k0+v*8+3<m?~0u:0u,k0+v*8+4<m?~0u:0u,k0+v*8+5<m?~0u:0u,k0+v*8+6<m?~0u:0u,k0+v*8+7<m?~0u:0u);}
  unsigned lastq=std::min(qhi,(n-1-k0)/m);for(unsigned q=qlo;q<=lastq;++q){unsigned last=std::min(r-1,(n-1-k0)/q),full=std::min(r-1,(n-k0-32)/q),j=m;for(;j<=full;++j){__m256i bj=_mm256_set1_epi32(b[j]);for(unsigned v=0;v<4;++v){__m256i x=_mm256_loadu_si256((const __m256i*)(a+q*j+k0+8*v));if(k0+32>m)x=_mm256_and_si256(x,mask[v]);z[v]=_mm256_add_epi32(z[v],_mm256_mullo_epi32(x,bj));}}for(;j<=last;++j){__m256i bj=_mm256_set1_epi32(b[j]);for(unsigned v=0;v<4;++v){unsigned base=q*j+k0+8*v;__m256i valid=_mm256_and_si256(mask[v],_mm256_setr_epi32(base<n?~0u:0u,base+1<n?~0u:0u,base+2<n?~0u:0u,base+3<n?~0u:0u,base+4<n?~0u:0u,base+5<n?~0u:0u,base+6<n?~0u:0u,base+7<n?~0u:0u));__m256i x=base<n?_mm256_maskload_epi32((const int*)(a+base),valid):_mm256_setzero_si256();z[v]=_mm256_add_epi32(z[v],_mm256_mullo_epi32(x,bj));}}}
  for(unsigned v=0;v<4;++v)_mm256_storeu_si256((__m256i*)(c+k0+8*v),z[v]);
 }
}
__attribute__((noinline)) inline void mixed_triangle(unsigned n,const u32*a,const u32*b,u32*c,unsigned l,unsigned r,unsigned threshold=32){
 if(r-l<=256||2ull*l+1>=n){adaptive_mod::direct_range(n,a,b,c,l,r);return;}
 unsigned w=(r-l)/2,m=l+w,qmax=(n-1-l)/m;
 mixed_triangle(n,a,b,c,l,m,threshold);mixed_triangle(n,a,b,c,m,r,threshold);
 if(qmax==0)return;
 unsigned karatsuba_max=std::min(qmax,w/threshold);
 if(karatsuba_max)cross(n,a,b,c,l,r,1,karatsuba_max);
 if(karatsuba_max<qmax)direct_cross(n,a,b,c,l,r,karatsuba_max+1,qmax);
}
}

namespace mod_karatsuba {
template<unsigned H>__attribute__((always_inline))inline void tiny_vec(const u32*a,const u32*b,u32*out){
 __m256i z[H];for(unsigned t=0;t<H;++t)z[t]=_mm256_setzero_si256();for(unsigned i=0;i<H;++i){__m256i bb=_mm256_set1_epi32(b[i]);for(unsigned t=0;t<H;++t)z[t]=_mm256_add_epi32(z[t],_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(a+8*(i+t))),bb));}for(unsigned t=0;t<H;++t)_mm256_storeu_si256((__m256i*)(out+8*t),z[t]);
}
__attribute__((noinline)) inline void middle_vec(const u32*a,const u32*b,u32*out,unsigned h,u32*scratch){
 if(h<=4){if(h==1)tiny_vec<1>(a,b,out);else if(h==2)tiny_vec<2>(a,b,out);else tiny_vec<4>(a,b,out);return;}
 unsigned s=h/2;u32*x=scratch,*y=x+8*h,*tmp=y+8*h,*bs=tmp+8*s,*next=bs+s;
 for(unsigned i=0;i<h;++i){__m256i mid=_mm256_loadu_si256((const __m256i*)(a+8*(s+i)));_mm256_storeu_si256((__m256i*)(x+8*i),_mm256_sub_epi32(_mm256_loadu_si256((const __m256i*)(a+8*i)),mid));_mm256_storeu_si256((__m256i*)(y+8*i),_mm256_sub_epi32(_mm256_loadu_si256((const __m256i*)(a+8*(h+i))),mid));}
 for(unsigned i=0;i<s;++i)bs[i]=b[i]+b[s+i];
 middle_vec(x,b,out,s,next);middle_vec(y,b+s,out+8*s,s,next);middle_vec(a+8*s,bs,tmp,s,next);
 for(unsigned i=0;i<s;++i){__m256i t=_mm256_loadu_si256((const __m256i*)(tmp+8*i));_mm256_storeu_si256((__m256i*)(out+8*i),_mm256_add_epi32(t,_mm256_loadu_si256((const __m256i*)(out+8*i))));_mm256_storeu_si256((__m256i*)(out+8*(s+i)),_mm256_add_epi32(t,_mm256_loadu_si256((const __m256i*)(out+8*(s+i)))));}
}
__attribute__((noinline)) inline void cross_vec(unsigned n,const u32*a,const u32*b,u32*c,unsigned l,unsigned r,unsigned qlo=1,unsigned qhi=~0u){
 const unsigned w=(r-l)/2,m=l+w,jend=std::min(r,n);if(!w||m>=jend||l+m>=n)return;
 qhi=std::min(qhi,(n-1-l)/m);if(qlo<=7)cross(n,a,b,c,l,r,qlo,std::min(qhi,7u));qlo=std::max(8u,qlo);
 for(unsigned q=qlo;q<=qhi;++q){unsigned wanted=std::min(w,n-q*m-l),h=1;while(h<(wanted+q-1)/q)h*=2;
  std::vector<u32>bb(h),aa(16*h+8),z(8*h),scratch(48*h+32);unsigned jhi=std::min(jend,(n-1-l)/q+1);
  for(unsigned jb=m;jb<jhi;jb+=h){unsigned count=std::min(h,jhi-jb);std::fill(bb.begin(),bb.end(),0);std::copy(b+jb,b+jb+count,bb.data());
   for(unsigned rho=0;rho<std::min(q,wanted);rho+=8){u64 start=u64(q)*jb+l+rho;
    for(unsigned i=0;i<2*h;++i){u64 at=start+u64(q)*i;__m256i v;if(at+8<=n)v=_mm256_loadu_si256((const __m256i*)(a+at));else if(at<n){unsigned rem=n-at;__m256i mask=_mm256_cmpgt_epi32(_mm256_set1_epi32(rem),_mm256_setr_epi32(0,1,2,3,4,5,6,7));v=_mm256_maskload_epi32((const int*)(a+at),mask);}else v=_mm256_setzero_si256();_mm256_storeu_si256((__m256i*)(aa.data()+8*i),v);}
    middle_vec(aa.data(),bb.data(),z.data(),h,scratch.data());
    for(unsigned t=0;t<h;++t){unsigned k=rho+q*t;if(k>=wanted)break;unsigned rem=std::min(8u,std::min(q-rho,wanted-k));__m256i value=_mm256_add_epi32(_mm256_loadu_si256((const __m256i*)(c+l+k)),_mm256_loadu_si256((const __m256i*)(z.data()+8*t)));if(rem==8)_mm256_storeu_si256((__m256i*)(c+l+k),value);else{__m256i mask=_mm256_cmpgt_epi32(_mm256_set1_epi32(rem),_mm256_setr_epi32(0,1,2,3,4,5,6,7));_mm256_maskstore_epi32((int*)(c+l+k),mask,value);}}
   }
  }
 }
}
__attribute__((noinline)) inline void vec_triangle(unsigned n,const u32*a,const u32*b,u32*c,unsigned l,unsigned r,unsigned cutoff=256){
 if(r-l<=cutoff||2ull*l+1>=n){adaptive_mod::direct_range(n,a,b,c,l,r);return;}
 unsigned w=(r-l)/2,m=l+w;vec_triangle(n,a,b,c,l,m,cutoff);vec_triangle(n,a,b,c,m,r,cutoff);cross_vec(n,a,b,c,l,r);
}
}

namespace mod_fused {
using mod_karatsuba::u32;using mod_karatsuba::u64;using mod_karatsuba::cross;
template<unsigned H>__attribute__((always_inline))inline void tiny_vec(const u32*a,const u32*b,u32*out){
 __m256i z[H];for(unsigned t=0;t<H;++t)z[t]=_mm256_setzero_si256();for(unsigned i=0;i<H;++i){__m256i bb=_mm256_set1_epi32(b[i]);for(unsigned t=0;t<H;++t)z[t]=_mm256_add_epi32(z[t],_mm256_mullo_epi32(_mm256_loadu_si256((const __m256i*)(a+8*(i+t))),bb));}for(unsigned t=0;t<H;++t)_mm256_storeu_si256((__m256i*)(out+8*t),z[t]);
}
__attribute__((always_inline))inline void fused8(const u32*a,const u32*b,u32*out){
 __m256i z0=_mm256_setzero_si256(), z1=_mm256_setzero_si256(), z2=_mm256_setzero_si256(), z3=_mm256_setzero_si256(), z4=_mm256_setzero_si256(), z5=_mm256_setzero_si256(), z6=_mm256_setzero_si256(), z7=_mm256_setzero_si256();
 for(unsigned i=0;i<4;++i){__m256i b0=_mm256_set1_epi32(b[i]),b1=_mm256_set1_epi32(b[4+i]),bs=_mm256_add_epi32(b0,b1);
 {__m256i am=_mm256_loadu_si256((const __m256i*)(a+8*(4+i+0))),prod=_mm256_mullo_epi32(am,bs);
 z0=_mm256_add_epi32(z0,prod);z4=_mm256_add_epi32(z4,prod);
 z0=_mm256_add_epi32(z0,_mm256_mullo_epi32(_mm256_sub_epi32(_mm256_loadu_si256((const __m256i*)(a+8*(i+0))),am),b0));
 z4=_mm256_add_epi32(z4,_mm256_mullo_epi32(_mm256_sub_epi32(_mm256_loadu_si256((const __m256i*)(a+8*(8+i+0))),am),b1));}
 {__m256i am=_mm256_loadu_si256((const __m256i*)(a+8*(4+i+1))),prod=_mm256_mullo_epi32(am,bs);
 z1=_mm256_add_epi32(z1,prod);z5=_mm256_add_epi32(z5,prod);
 z1=_mm256_add_epi32(z1,_mm256_mullo_epi32(_mm256_sub_epi32(_mm256_loadu_si256((const __m256i*)(a+8*(i+1))),am),b0));
 z5=_mm256_add_epi32(z5,_mm256_mullo_epi32(_mm256_sub_epi32(_mm256_loadu_si256((const __m256i*)(a+8*(8+i+1))),am),b1));}
 {__m256i am=_mm256_loadu_si256((const __m256i*)(a+8*(4+i+2))),prod=_mm256_mullo_epi32(am,bs);
 z2=_mm256_add_epi32(z2,prod);z6=_mm256_add_epi32(z6,prod);
 z2=_mm256_add_epi32(z2,_mm256_mullo_epi32(_mm256_sub_epi32(_mm256_loadu_si256((const __m256i*)(a+8*(i+2))),am),b0));
 z6=_mm256_add_epi32(z6,_mm256_mullo_epi32(_mm256_sub_epi32(_mm256_loadu_si256((const __m256i*)(a+8*(8+i+2))),am),b1));}
 {__m256i am=_mm256_loadu_si256((const __m256i*)(a+8*(4+i+3))),prod=_mm256_mullo_epi32(am,bs);
 z3=_mm256_add_epi32(z3,prod);z7=_mm256_add_epi32(z7,prod);
 z3=_mm256_add_epi32(z3,_mm256_mullo_epi32(_mm256_sub_epi32(_mm256_loadu_si256((const __m256i*)(a+8*(i+3))),am),b0));
 z7=_mm256_add_epi32(z7,_mm256_mullo_epi32(_mm256_sub_epi32(_mm256_loadu_si256((const __m256i*)(a+8*(8+i+3))),am),b1));}
 }
 _mm256_storeu_si256((__m256i*)(out+0),z0);
 _mm256_storeu_si256((__m256i*)(out+8),z1);
 _mm256_storeu_si256((__m256i*)(out+16),z2);
 _mm256_storeu_si256((__m256i*)(out+24),z3);
 _mm256_storeu_si256((__m256i*)(out+32),z4);
 _mm256_storeu_si256((__m256i*)(out+40),z5);
 _mm256_storeu_si256((__m256i*)(out+48),z6);
 _mm256_storeu_si256((__m256i*)(out+56),z7);
}
__attribute__((noinline)) inline void middle_vec(const u32*a,const u32*b,u32*out,unsigned h,u32*scratch){
 if(h==8){fused8(a,b,out);return;}
 if(h<=4){if(h==1)tiny_vec<1>(a,b,out);else if(h==2)tiny_vec<2>(a,b,out);else tiny_vec<4>(a,b,out);return;}
 unsigned s=h/2;u32*x=scratch,*y=x+8*h,*tmp=y+8*h,*bs=tmp+8*s,*next=bs+s;
 for(unsigned i=0;i<h;++i){__m256i mid=_mm256_loadu_si256((const __m256i*)(a+8*(s+i)));_mm256_storeu_si256((__m256i*)(x+8*i),_mm256_sub_epi32(_mm256_loadu_si256((const __m256i*)(a+8*i)),mid));_mm256_storeu_si256((__m256i*)(y+8*i),_mm256_sub_epi32(_mm256_loadu_si256((const __m256i*)(a+8*(h+i))),mid));}
 for(unsigned i=0;i<s;++i)bs[i]=b[i]+b[s+i];
 middle_vec(x,b,out,s,next);middle_vec(y,b+s,out+8*s,s,next);middle_vec(a+8*s,bs,tmp,s,next);
 for(unsigned i=0;i<s;++i){__m256i t=_mm256_loadu_si256((const __m256i*)(tmp+8*i));_mm256_storeu_si256((__m256i*)(out+8*i),_mm256_add_epi32(t,_mm256_loadu_si256((const __m256i*)(out+8*i))));_mm256_storeu_si256((__m256i*)(out+8*(s+i)),_mm256_add_epi32(t,_mm256_loadu_si256((const __m256i*)(out+8*(s+i)))));}
}
__attribute__((noinline)) inline void cross_vec(unsigned n,const u32*a,const u32*b,u32*c,unsigned l,unsigned r,unsigned qlo=1,unsigned qhi=~0u){
 const unsigned w=(r-l)/2,m=l+w,jend=std::min(r,n);if(!w||m>=jend||l+m>=n)return;
 qhi=std::min(qhi,(n-1-l)/m);if(qlo<=7)cross(n,a,b,c,l,r,qlo,std::min(qhi,7u));qlo=std::max(8u,qlo);
 for(unsigned q=qlo;q<=qhi;++q){unsigned wanted=std::min(w,n-q*m-l),h=1;while(h<(wanted+q-1)/q)h*=2;
  std::vector<u32>bb(h),aa(16*h+8),z(8*h),scratch(48*h+32);unsigned jhi=std::min(jend,(n-1-l)/q+1);
  for(unsigned jb=m;jb<jhi;jb+=h){unsigned count=std::min(h,jhi-jb);std::fill(bb.begin(),bb.end(),0);std::copy(b+jb,b+jb+count,bb.data());
   for(unsigned rho=0;rho<std::min(q,wanted);rho+=8){u64 start=u64(q)*jb+l+rho;
    for(unsigned i=0;i<2*h;++i){u64 at=start+u64(q)*i;__m256i v;if(at+8<=n)v=_mm256_loadu_si256((const __m256i*)(a+at));else if(at<n){unsigned rem=n-at;__m256i mask=_mm256_cmpgt_epi32(_mm256_set1_epi32(rem),_mm256_setr_epi32(0,1,2,3,4,5,6,7));v=_mm256_maskload_epi32((const int*)(a+at),mask);}else v=_mm256_setzero_si256();_mm256_storeu_si256((__m256i*)(aa.data()+8*i),v);}
    middle_vec(aa.data(),bb.data(),z.data(),h,scratch.data());
    for(unsigned t=0;t<h;++t){unsigned k=rho+q*t;if(k>=wanted)break;unsigned rem=std::min(8u,std::min(q-rho,wanted-k));__m256i value=_mm256_add_epi32(_mm256_loadu_si256((const __m256i*)(c+l+k)),_mm256_loadu_si256((const __m256i*)(z.data()+8*t)));if(rem==8)_mm256_storeu_si256((__m256i*)(c+l+k),value);else{__m256i mask=_mm256_cmpgt_epi32(_mm256_set1_epi32(rem),_mm256_setr_epi32(0,1,2,3,4,5,6,7));_mm256_maskstore_epi32((int*)(c+l+k),mask,value);}}
   }
  }
 }
}
__attribute__((noinline)) inline void vec_triangle(unsigned n,const u32*a,const u32*b,u32*c,unsigned l,unsigned r,unsigned cutoff=256){
 if(r-l<=cutoff||2ull*l+1>=n){adaptive_mod::direct_range(n,a,b,c,l,r);return;}
 unsigned w=(r-l)/2,m=l+w;vec_triangle(n,a,b,c,l,m,cutoff);vec_triangle(n,a,b,c,m,r,cutoff);cross_vec(n,a,b,c,l,r);
}
}

#ifndef DUCK_SHA256_HPP
#define DUCK_SHA256_HPP
#include <stddef.h>
#include <stdint.h>
#include <string.h>
// Portable SHA-256; input is the exact byte sequence, with no native-word hashing.
namespace duck_sha256 {
static inline uint32_t rotate(uint32_t x,unsigned n) {return (x>>n)|(x<<(32-n));}
static void compress(uint32_t h[8],const unsigned char block[64]) {
    static const uint32_t k[64]={
      0x428a2f98u,0x71374491u,0xb5c0fbcfu,0xe9b5dba5u,0x3956c25bu,0x59f111f1u,0x923f82a4u,0xab1c5ed5u,
      0xd807aa98u,0x12835b01u,0x243185beu,0x550c7dc3u,0x72be5d74u,0x80deb1feu,0x9bdc06a7u,0xc19bf174u,
      0xe49b69c1u,0xefbe4786u,0x0fc19dc6u,0x240ca1ccu,0x2de92c6fu,0x4a7484aau,0x5cb0a9dcu,0x76f988dau,
      0x983e5152u,0xa831c66du,0xb00327c8u,0xbf597fc7u,0xc6e00bf3u,0xd5a79147u,0x06ca6351u,0x14292967u,
      0x27b70a85u,0x2e1b2138u,0x4d2c6dfcu,0x53380d13u,0x650a7354u,0x766a0abbu,0x81c2c92eu,0x92722c85u,
      0xa2bfe8a1u,0xa81a664bu,0xc24b8b70u,0xc76c51a3u,0xd192e819u,0xd6990624u,0xf40e3585u,0x106aa070u,
      0x19a4c116u,0x1e376c08u,0x2748774cu,0x34b0bcb5u,0x391c0cb3u,0x4ed8aa4au,0x5b9cca4fu,0x682e6ff3u,
      0x748f82eeu,0x78a5636fu,0x84c87814u,0x8cc70208u,0x90befffau,0xa4506cebu,0xbef9a3f7u,0xc67178f2u};
    uint32_t w[64];
    for(unsigned i=0;i<16;++i)w[i]=(uint32_t(block[4*i])<<24)|(uint32_t(block[4*i+1])<<16)|(uint32_t(block[4*i+2])<<8)|block[4*i+3];
    for(unsigned i=16;i<64;++i){uint32_t a=w[i-15],b=w[i-2];w[i]=w[i-16]+(rotate(a,7)^rotate(a,18)^(a>>3))+w[i-7]+(rotate(b,17)^rotate(b,19)^(b>>10));}
    uint32_t a=h[0],b=h[1],c=h[2],d=h[3],e=h[4],f=h[5],g=h[6],v=h[7];
    for(unsigned i=0;i<64;++i){
        uint32_t t1=v+(rotate(e,6)^rotate(e,11)^rotate(e,25))+((e&f)^(~e&g))+k[i]+w[i];
        uint32_t t2=(rotate(a,2)^rotate(a,13)^rotate(a,22))+((a&b)^(a&c)^(b&c));
        v=g;g=f;f=e;e=d+t1;d=c;c=b;b=a;a=t1+t2;
    }
    h[0]+=a;h[1]+=b;h[2]+=c;h[3]+=d;h[4]+=e;h[5]+=f;h[6]+=g;h[7]+=v;
}
inline void digest(const void*data,size_t length,unsigned char out[32]) {
    uint32_t h[8]={0x6a09e667u,0xbb67ae85u,0x3c6ef372u,0xa54ff53au,0x510e527fu,0x9b05688cu,0x1f83d9abu,0x5be0cd19u};
    const unsigned char*p=(const unsigned char*)data;
    size_t whole=length/64;
    for(size_t i=0;i<whole;++i)compress(h,p+64*i);
    unsigned char tail[128]={};size_t remain=length%64;
    if(remain)memcpy(tail,p+64*whole,remain);
    tail[remain]=0x80;size_t padded=remain<56?64:128;
    uint64_t bits=uint64_t(length)*8;
    for(unsigned i=0;i<8;++i)tail[padded-1-i]=(unsigned char)(bits>>(8*i));
    compress(h,tail);if(padded==128)compress(h,tail+64);
    for(unsigned i=0;i<8;++i)for(unsigned j=0;j<4;++j)out[4*i+j]=(unsigned char)(h[i]>>(24-8*j));
}
}
#endif


#include <x86intrin.h>
static uint64_t ticks(){_mm_lfence();return __rdtsc();}
struct PublicRng{uint64_t s;unsigned next(){uint64_t z=(s+=0x9e3779b97f4a7c15ULL);z=(z^(z>>30))*0xbf58476d1ce4e5b9ULL;z=(z^(z>>27))*0x94d049bb133111ebULL;return unsigned(z^(z>>31));}};
int main(){unsigned char hashes[64];unsigned at=0;PublicRng rg{0x4433998866112255ULL};for(unsigned H:{16384u,32768u}){constexpr unsigned R=4;std::vector<unsigned>a(16*H),b(H),k(8*H),n(8*H),scratch(48*H+32);for(auto&v:a)v=rg.next();for(auto&v:b)v=rg.next();nuss::Plan64 bp(b.data(),H);nuss::Plan4 p(bp);uint64_t t0=ticks();for(unsigned r=0;r<R;++r)mod_fused::middle_vec(a.data(),b.data(),k.data(),H,scratch.data());uint64_t t1=ticks();for(unsigned r=0;r<R;++r){if(!p.apply(a.data(),n.data(),8,true)||!p.apply(a.data()+4,n.data()+4,8,true))return 1;}uint64_t t2=ticks();if(k!=n)return 2;duck_sha256::digest(n.data(),n.size()*4,hashes+32*at++);printf("H%u R%u K%llu N%llu\n",H,R,(unsigned long long)(t1-t0),(unsigned long long)(t2-t1));}unsigned char hash[32];duck_sha256::digest(hashes,64,hash);for(auto v:hash)printf("%02x",v);puts(" PASS");}

评测结果