fwt

多项式模板集

NTT&&FFT

NTT板⼦

FFT板⼦

typedeflonglongll;

constintP=998244353,g=3;

constintmaxn=1111111;

intinc(inta,intb){return(a+=b)>=P?a-P:a;}

intqpow(inta,intb){

intres=1;

for(;b;a=1ll*a*a%P,b>>=1)

if(b&1)res=1ll*res*a%P;

returnres;

}

intW[maxn<<2],inv[maxn<<2];//4倍空间

voidprework(intn){

for(inti=1;i

W[i]=1;

for(intj=1,Wn=qpow(g,(P-1)/i>>1);j

}

inv[1]=1;

for(inti=2;i<=n;i++)inv[i]=1ll*(P-P/i)*inv[P%i]%P;//取等

}

voidntt(int*a,intn,intopt){

staticintrev[maxn<<2]={0};//{0}赋初值

for(inti=1;i

if((rev[i]=rev[i>>1]>>1|(i&1?n>>1:0))>i)std::swap(a[i],a[rev[i]]);

for(intq=1;q

for(intp=0;p

for(inti=0,t;i

t=1ll*a[p+q+i]*W[q+i]%P,a[p+q+i]=inc(a[p+i],P-t),a[p+i]=inc(a[p+i],t);

if(~opt)return;

std::rever(a+1,a+n);

for(inti=0;i

}

intgetsize(intn){intx=1;while(x

Copy

任意模数NTT（MTT）

MTT事实上就是将⼀个⼤整数拆分成两部分a×Ba+b，分开相乘最后相加就能保证精度了。直接做FFT次数很多常

数很⼤，在%完myy的论⽂和代码后发现有⼀种只做4次FFT的⽅法，需要⼀些trick和引理。

考虑构造A

i

=a+bi和B

i

=a−bi，如果快速求出FFT(A)和FFT(B)，那么就能加减消元

求出FFT(a)和FFT(b)。myy论⽂中讲如果求出了FFT(A)，可以利⽤其直接推出FFT(B)。

原⽂是这样说的：

B(ωk

)=

¯

A(ω−k)

其中上划线表⽰共轭复数

令ω

k

幅⾓为θ，注意到

//⼿动定义comp

structcomp{doublex,y;};

compoperator+(compa,compb){return(comp){a.x+b.x,a.y+b.y};}

compoperator-(compa,compb){return(comp){a.x-b.x,a.y-b.y};}

compoperator*(compa,compb){return(comp){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}

compW[maxn<<2];

voidprework(intn){

for(inti=1;i

for(intj=0;j

W[i+j]=(comp){cos(PI/i*j),sin(PI/i*j)};

//这⾥直接算，防⽌丢精度。这在MTT中很重要

}

voidfft(comp*a,intn,intopt){//与NTT没有区别

staticintrev[maxn<<2]={0};

for(inti=1;i

if((rev[i]=rev[i>>1]>>1|(i&1?n>>1:0))>i)std::swap(a[i],a[rev[i]]);

for(intq=1;q

for(intp=0;p

for(inti=0;i

compt=a[p+q+i]*W[q+i];a[p+q+i]=a[p+i]-t,a[p+i]=a[p+i]+t;

}

if(~opt)return;

std::rever(a+1,a+n);

for(inti=0;i

}

Copy

a×Ba+b

=a+biA

i

=a−biB

i

FFT(A)FFT(B)

FFT(a)FFT(b)FFT(A)FFT(B)

B()=ωkA()ω−k

¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

ωkθ

B()ωk==(−i)(cos(jθ)+isin(jθ))∑

j=0

n

B

j

ωjk∑

j=0

n

a

j

b

j

=((cos(jθ)+sin(jθ))+i(sin(jθ)−cos(jθ)))∑

j=0

n

a

j

b

j

a

j

b

j

=((cos(−jθ)−sin(−jθ))−i(sin(−jθ)+cos(−jθ)))∑

j=0

n

a

j

b

j

a

j

b

j

=((cos(−jθ)−sin(−jθ))+i(sin(−jθ)+cos(−jθ)))∑

j=0

n

a

j

b

j

a

j

b

j

¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

=(+i)(cos(−jθ)+isin(−jθ))∑

i=0

n

a

j

b

j

¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

==∑

i=0

n

A

j

ω−jk

¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

A()ω−k

¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

B(ωk

)=

n

∑

j=0B

j

ω

jk

=

n

∑

j=0(a

j

−b

j

i)(cos(jθ)+isin(jθ))

=

n

∑

j=0(a

j

cos(jθ)+b

j

sin(jθ))+i(a

j

sin(jθ)−b

j

cos(jθ))

=

n

∑

j=0(a

j

cos(−jθ)−b

j

sin(−jθ))−i(a

j

sin(−jθ)+b

j

cos(−jθ))

=

¯

n

∑

j=0(a

j

cos(−jθ)−b

j

sin(−jθ))+i(a

j

sin(−jθ)+b

j

cos(−jθ))

=

¯

n

∑

i=0(a

j

+b

j

i)(cos(−jθ)+isin(−jθ))

=

¯

n

∑

i=0A

j

ω

−jk

=

¯

A(ω−k)

运⽤这个性质优化可以⼤幅减少FFT的次数。

快速沃尔什变换（FWT）

⽤于解决对下标进⾏位运算卷积问题的⽅法。即

c

k

=

∑

i⊕j=ka

i

×b

j

其中⊕分别为|、&、^的情况。

先考虑|。FWT⼲这样的事情：类似于FFT，对于每⼀个i，它要求出fwt(a)

i

=∑

j|i=i

a

j

，在

O(nlogn)的时间复杂度在两者之间快速变换。然后j|i=i,k|i=i⇔(j|k)|i=i。

()

()

()

voidconv(int*x,int*y,int*z,intn){//z=x*y，长度为n

for(inti=0;i

staticcompa[maxn<<2],b[maxn<<2],da[maxn<<2],db[maxn<<2],dc[maxn<<2],dd[maxn<<2];

for(inti=0;i

a[i]=(comp){x[i]>>15,x[i]&32767},b[i]=(comp){y[i]>>15,y[i]&32767};//将整数拆分成两部分

fft(a,n,1),fft(b,n,1);

for(inti=0;i

intj=(n-1)&(n-i);

staticcompa1,a2,b1,b2;//分离出x和y每个部分的插值结果

a1=(a[i]+conj(a[j]))*(comp){0.5,0};

a2=(a[i]-conj(a[j]))*(comp){0,-0.5};

b1=(b[i]+conj(b[j]))*(comp){0.5,0};

b2=(b[i]-conj(b[j]))*(comp){0,-0.5};

da[i]=a1*b1,db[i]=a1*b2,dc[i]=a2*b1,dd[i]=a2*b2;

}

for(inti=0;i

a[i]=da[i]+db[i]*(comp){0,1},b[i]=dc[i]+dd[i]*(comp){0,1};//IDFT(x+yi)=IDFT(x)+iIDFT(y)，这个将⽤于下⽂

fft(a,n,-1),fft(b,n,-1);

for(inti=0;i

intax=(ll)(a[i].x+0.5)%P,ay=(ll)(a[i].y+0.5)%P,bx=(ll)(b[i].x+0.5)%P,by=(ll)(b[i].y+0.5)%P;//⼀定要转化成ll（数值为2^30*n）

z[i]=(((ll)ax<<30)+((ll)(ay+bx)<<15)+by)%P;

}

}

Copy

=×c

k∑

i⊕j=k

a

i

b

j

⊕

ifwt(a=)

i∑

j|i=i

a

j

O(nlogn)

j|i=i,k|i=i⇔(j|k)|i=i

fwt(a)fwt(b)

如果求出了fwt(a)和fwt(b)，发现有

fwt(a)

i

×fwt(b)

i

=

∑

j|i=ia

j

∑

k|i=ib

k

=

∑

j|i=i,k|i=ia

j

b

k

=

∑

(j|k)|i=ia

j

b

k

=

∑

t|i=i

∑

j|k=ia

j

b

k

=fwt(c)

i

很像系数变点值！考虑怎么变换。我们按最⾼位为0或1来分成两组序列a

[0]

和a

[1]

，不难发现

fwt(a)=merge(fwt(a[0]),fwt(a[0])+fwt(a[1]))

最⾼位是0的序列对应最⾼位是1的⼀定是包含关系，所以右边相加。

同理我们不难发现

a=merge(a[0],a[1]−a[0])

对于&，也满⾜fwt(a)

i

×fwt(b)

i

=fwt(c)

i

。同上我们也可以推导出

fwt(a)=merge(fwt(a[0])+fwt(a[1]),fwt(a[1]))

a=merge(a[0]−a[1],a[1])

对于^稍微⿇烦些，定义x⊗y为x&y中⼆进制下1的个数对2取模。有

构造，则

符合要求，所以能推导出

fwt(a)fwt(b)

fwt(a×fwt(b)

i

)

i

=⎛

⎝∑

j|i=i

a

j⎞

⎠

⎛

⎝∑

k|i=i

b

k⎞

⎠

====fwt(c∑

j|i=i,k|i=i

a

j

b

k∑

(j|k)|i=i

a

j

b

k∑

t|i=i

∑

j|k=i

a

j

b

k

)

i()()

a[0]a[1]

fwt(a)=merge(fwt(),fwt()+fwt())a[0]a[0]a[1]

a=merge(,−)a[0]a[1]a[0]

//or

voidfwt_or(int*a,intn,intopt){

for(intq=1;q

for(intp=0;p

for(inti=0;i

a[p+q+i]=inc(a[p+q+i],~opt?a[p+i]:P-a[p+i]);

}

Copy

fwt(a×fwt(b=fwt(c)

i

)

i

)

i

fwt(a)=merge(fwt()+fwt(),fwt())a[0]a[1]a[1]

a=merge(−,)a[0]a[1]a[1]

x⊗yx&y

(i⊗j)xor(i⊗k)=i⊗(jxork)

fwt(a=−)

i∑

i⊗j=0

a

j∑

i⊗j=1

a

j

fwt(a×fwt(b)

i

)

i

=(−)(−)∑

i⊗j=0

a

j∑

i⊗j=1

a

j∑

i⊗j=0

b

j∑

i⊗j=1

b

j

=−−+∑

i⊗j=0,i⊗k=0

a

j

b

k∑

i⊗j=0,i⊗k=1

a

j

b

k∑

i⊗j=1,i⊗k=0

a

j

b

k∑

i⊗j=1,i⊗k=1

a

j

b

k

=−∑

(i⊗j)xor(i⊗k)=0

a

j

b

k∑

(i⊗j)xor(i⊗k)=1

a

j

b

k

=−∑

i⊗(jxork)=0

a

j

b

k∑

i⊗(jxork)=1

a

j

b

k

=fwt(c)

i

fwt(a)=merge(fwt()+fwt(),fwt()−fwt())a[0]a[1]a[0]a[1]

=merge(,)

[0][1][0][1]

多项式

a=merge(,)

+a[0]a[1]

2

−a[0]a[1]

2

//and&&xor

voidfwt_and(int*a,intn,intopt){

for(intq=1;q

for(intp=0;p

for(inti=0;i

a[p+i]=inc(a[p+i],~opt?a[p+q+i]:P-a[p+q+i]);

}

voidfwt_xor(int*a,intn,intopt){

for(intq=1;q

for(intp=0;p

for(inti=0;i

intt=a[p+q+i];a[p+q+i]=inc(a[p+i],P-t);a[p+i]=inc(a[p+i],t);

if(opt==-1)a[p+i]=1ll*a[p+i]*inv2%P,a[p+q+i]=1ll*a[p+q+i]*inv2%P;

}

}

Copy

#include

usingstd::rever;usingstd::vector;usingstd::swap;usingstd::max;

constintN=100005,P=998244353,inv2=P+1>>1;

typedefvectorPoly;

typedeflonglongLL;

intinc(inta,intb){return(a+=b)>=P?a-P:a;}

intpow(inta,intb){

intt=1;

for(;b;b>>=1,a=1LL*a*a%P)

if(b&1)t=1LL*t*a%P;

returnt;

}

intW[N*4],inv[N*4];

voidprework(intn){

for(inti=1;i

for(intj=W[i]=1,Wn=pow(3,(P-1)/i>>1);j

W[i+j]=1LL*W[i+j-1]*Wn%P;

inv[1]=1;

for(inti=2;i<=n;i++)inv[i]=1LL*(P-P/i)*inv[P%i]%P;

}

voidfft(Poly&a,intn,intopt){

(n);

staticintrev[N*4];

for(inti=1;i

if((rev[i]=rev[i>>1]>>1|(i&1?n>>1:0))>i)swap(a[i],a[rev[i]]);

for(intq=1;q

for(intp=0;p

for(inti=0,t;i

t=1LL*W[q+i]*a[p+q+i]%P,a[p+q+i]=inc(a[p+i],P-t),a[p+i]=inc(a[p+i],t);

if(opt)return;

for(inti=0,inv=pow(n,P-2);i

rever(()+1,());

}

Polypoly_inv(PolyA){

PolyB(1,pow(A[0],P-2)),C(2);

for(intL=1;L<();L<<=1){

(C=A).resize(L*2);fft(B,L*4,1),fft(C,L*4,1);

for(inti=0;i

fft(B,L*4,0);(L*2);

}

(()),B;

}

intgetsz(intn){intx=1;while(x

Copy

intgetsz(intn){intx=1;while(x

voidfix(Poly&A){intx=();while(x>1&&!A[x-1])x--;(x);}

Polyoperator+(PolyA,PolyB){

(max((),()));

for(inti=0;i<();i++)A[i]=inc(A[i],B[i]);

returnfix(A),A;

}

Polyoperator-(PolyA,PolyB){

(max((),()));

for(inti=0;i<();i++)A[i]=inc(A[i],P-B[i]);

returnfix(A),A;

}

Polyoperator*(intk,PolyA){

for(inti=0;i<();i++)A[i]=1LL*k*A[i]%P;

returnA;

}

Polyoperator*(PolyA,PolyB){

intL=getsz(()+()-1);

fft(A,L,1),fft(B,L,1);

for(inti=0;i

returnfft(A,L,0),fix(A),A;

}

Polyoperator/(PolyA,PolyB){

intn=()-()+1;

rever((),());(n);

rever((),());(n);

returnA=A*poly_inv(B),(n),rever((),()),fix(A),A;

}

Polyoperator%(PolyA,PolyB){returnA-A/B*B;}

Polypoly_deri(PolyA){

for(inti=0;i<()-1;i++)A[i]=1LL*(i+1)*A[i+1]%P;

(()-1),A;

}

Polypoly_int(PolyA){

for(inti=()-1;i;i--)A[i]=1LL*A[i-1]*inv[i]%P;

returnA[0]=0,A;

}

Polypoly_sqrt(PolyA){

PolyB(1,1),iB,C(2);

for(intL=1;L<();L<<=1){

(C=A).resize(L*2);(L*2);iB=poly_inv(B);

fft(B,L*4,1),fft(iB,L*4,1),fft(C,L*4,1);

for(inti=0;i

fft(B,L*4,0);(L*2);

}

(()),B;

}

Polypoly_ln(PolyA){

PolyB=poly_deri(A)*poly_inv(A);

(()),poly_int(B);

}

Polypoly_exp(PolyA){

PolyB(1,1),C;

for(intL=1;L<();L<<=1)

(L*2),C=A+Poly(1,1)-poly_ln(B),(L*2),B=B*C;

(()),B;

}

Polypoly_pow(PolyA,intk){

returnpoly_exp(k*poly_ln(A));

}

#definelc(o<<1)

#definerc(o<<1|1)

PolyQ[N*4];

voidbuild(into,intl,intr,intx[]){

if(l==r){Q[o].push_back(P-x[l]),Q[o].push_back(1);return;}

intmid=l+r>>1;

build(lc,l,mid,x),build(rc,mid+1,r,x);

Q[o]=Q[lc]*Q[rc];

}

}

voidcalc(PolyA,into,intl,intr,intx[]){

if(l==r){x[l]=A[0];return;}

intmid=l+r>>1;

calc(A%Q[lc],lc,l,mid,x),calc(A%Q[rc],rc,mid+1,r,x);

}

voidpoly_calc(PolyA,intn,intx[]){

build(1,1,n,x);calc(A,1,1,n,x);

}

Polyinter(into,intl,intr,intx[],inty[]){

if(l==r)returnPoly(1,1LL*y[l]*pow(x[l],P-2)%P);

intmid=l+r>>1;

returninter(lc,l,mid,x,y)*Q[rc]+inter(rc,mid+1,r,x,y)*Q[lc];

}

Polypoly_inter(intn,intx[],inty[]){

returnbuild(1,1,n,x),calc(poly_deri(Q[1]),1,1,n,x),inter(1,1,n,x,y);

}

intn;PolyA;

intmain(){

scanf("%d",&n);(n);prework(n*2);

for(inti=0;i

A=poly_exp(A);

for(inti=0;i

return0;

}

Processingmath:86%