「学习笔记」多项式全家桶

FFT

struct node{
    long double x,y;
    node(){}
    node(long double xx,long double yy){x=xx; y=yy; return ;}
    node operator +(const node &a)const{return node(x+a.x,y+a.y);}
    node operator -(const node &a)const{return node(x-a.x,y-a.y);}
    node operator *(const node &a)const{
        return node(x*a.x-y*a.y,x*a.y+y*a.x);}
    }
    inline node conj(){return node(x,-y);}
};
inline void FFT(node *f,int lim,int opt){
	rep(i,0,lim-1){
        r[i]=r[i>>1]>>1|((i&1)?(lim>>1):0)；
        if(i<r[i]) swap(f[i],f[r[i]]);
    }
    for(int p=2;p<=lim;p<<=1){
	    int len=p>>1; node tt=W[0]=node(1,0); 
        rep(i,1,len-1) W[i]=node(cos(pi/len*i),opt*sin(pi/len*i));
	    for(int k=0;k<lim;k+=p){
            for(int l=k;l<k+len;++l){
                tt=f[l+len]*W[l-k];
                f[len+l]=f[l]-tt,f[l]=f[l]+tt;
            }
        }
    } if(opt==-1) rep(i,0,lim-1) f[i].x/=lim; return ;
}

提升精度的方式是按照上面给出的代码写，而不是整个 buf 每次乘

万径人踪灭

答案应该可以写成：回文子序列数 \(-\) 回文子串个数

后面的随便做一下就行了，比如 \(hash+\) 二分或 \(Manacher\)

然后考虑前面的部分：

令 \(f_{i}=\sum\limits_{j<i} [s_j=s_{i\times 2-j}]\)

那么以 \(i\) 为对称半径的答案就是 \(2^{f_i}-1\)

设 \(a=1,b=0\) 那么 \(1\times 1=1,0\times 0=0\times 1=1 \times 0=0\) 卷积上去即可

再令 \(a=0,b=1\) 做一遍一样的

NTT

inline void NTT(vector<int> &f,int lim,int opt){ f.resize(lim);
    rep(i,0,lim-1){
        r[i]=r[i>>1]>>1|((i&1)?(lim>>1):0)；
        if(i<r[i]) swap(f[i],f[r[i]]);
    }
    for(int p=2;p<=lim;p<<=1){            
        int len=p>>1; W[0]=1; W[1]=ksm(3,(mod-1)/p); 
        if(opt==-1) W[1]=ksm(W[1],mod-2); 
        rep(j,2,len-1) W[j]=mul(W[j-1],W[1]);
        for(int k=0;k<lim;k+=p){
            for(int l=k,tt;l<k+len;++l){
                tt=mul(W[l-k],f[l+len])；
                f[l+len]=del(f[l],tt),ckadd(f[l],tt);
            }
        }
    } 
    if(opt==-1){
        for(int i=0,tmp=ksm(lim,mod-2);i<lim;++i) ckmul(f[i],tmp); 
    }
    return ;
}

写代码的习惯是很重要的，预处理单位根和最后只做一次快速幂还是挺重要的

一个人的高三楼

其实也是 \(Luogu5488\) 差分和前缀和

前缀和就是 $\Delta A=A * 1 $

那么对 \(1\) 多项式快速幂可以得到一个比较慢的做法

那么考虑实际意义：\(\Delta^k A_i\) 就是所有的加和为 \(i\) 的位置的求和

那么系数是 \(\binom {i+k-1}{k-1},NTT\) 即可

差分的初始卷的应该是 \(1-x\)（简单构造可以得到），二项式定理就可以得到系数

（理解起来貌似是很简单的）

或者也可以用生成函数的封闭形式来得到这题……

HAOI2018 染色

\[ans=\sum\limits_{i=0} ^{min(m,\lfloor \frac{n}S\rfloor)} w_i\times method(i) \]

做一些卷积题目就会猜一下这个 \(method(i)\) 可能是要卷积出来的东西

那么推一下这个 \(method(i)\)

\[method(i)=\binom M i\binom N {iS} (m-i)^{n-iS}\times \frac{(iS)!}{(S!)^i} \]

然后就错了，因为最后的乱放可能再次出现有 \(S\) 个颜色

那么容斥这个过程，设最终出现了 \(S\) 次的颜色种类数是 \(j\) 那么得到如下：

\[method(i)=\binom M i\binom N{iS}\frac{(iS)!}{(S!)^i} \sum\limits _ {j=i} ^M (-1)^{j-i} \binom {M-i}{j-i} \binom {n-iS}{(j-i)S} \frac{((j-i)S)!}{(S!)^{j-i}}(m-j)^{N-jS} \]

这式子貌似不短

那么化简完了 \(NTT\) 即可得到答案

把所有的组合数打开得到

\[method(i)=\frac{N!M!}{i!} \sum\limits_{j=i}^M \frac{(-1)^{j-i}}{(j-i)!}\frac{(M-j)^{n-jS}}{(M-j)!(N-jS)(S!)^j} \]

减法卷积可以翻转完了直接做，但是如果有梦想也可以推成加法卷积

\[\sum_{j=0}^N\frac{N!M!(m-j)^{N-jS}}{(S!)^j(N-jS)!(M-j)!} \sum_{i=0}^j \frac{w_i}{i!} \times\frac{(-1)^{j-i}}{(j-i)!} \]

\(\Theta(n\log n+n)\) 即可

分治FFT

还是 \(cdq\) 的那一套

每次对于当前区间处理当前左边和 \(g\) 的卷积，把答案累加到右侧

按照左中右的顺序进行处理

有些细节：开始的时候把序列长度补成 \(2^x\) 的形式，给右边添加的时候的起点需要思考

这样就会做快速求第二类斯特林数的一列的形式

inline void NTT(int *f,int lim,int fl){
    for(int i=0;i<lim;++i) if(i<R[i]) swap(f[i],f[R[i]]);
    for(int p=2;p<=lim;p<<=1){
        int len=p>>1,tmp=ksm(3,(mod-1)/p); if(fl==-1) tmp=ksm(tmp,mod-2);
        for(int k=0;k<lim;k+=p){
            int buf=1; for(int l=k;l<k+len;++l){
                int tt=mul(buf,f[l+len]); 
                f[l+len]=del(f[l],tt); f[l]=add(f[l],tt);
                buf=mul(buf,tmp);
            }
        }
    } 
    if(fl==-1) for(int i=0,tp=ksm(lim,mod-2);i<lim;++i) f[i]=mul(f[i],tp); 
    return ;
}
inline void solve(int l,int r){
    if(l+1==r) return ; int mid=l+((r-l)>>1); solve(l,mid); 
    int len=r-l; 
    for(int i=0;i<len;++i) R[i]=R[i>>1]>>1|((i&1)?(len>>1):0);
    for(int i=0;i<len;++i) G[i]=g[i]; 
    for(int i=l;i<mid;++i) F[i-l]=f[i]; 
    for(int i=mid;i<r;++i) F[i-l]=0;
    NTT(F,len,1); NTT(G,len,1); 
    for(int i=0;i<len;++i) F[i]=mul(F[i],G[i]); 
    NTT(F,len,-1);
    for(int i=mid;i<r;++i) f[i]=add(f[i],F[i-l]);
    return solve(mid,r);
}

任意模数NTT

考虑 \(FFT\) 的精度如果在系数是 \(10^9\) 的情况会爆炸，那么考虑优化这个直接 \(FFT\) 再取模的做法

把当前多项式的每个 \(A_i\) 拆成 \(A_0[i]\times k+A_1[i]\) 这里 \(A_0[i]=A[i]/k,A_1[i]=A[i]\%k\)，\(B\) 也拆一下系数

\(DFT\) 求点值的时候 考虑把两个并成一个，因为每个数一开始的虚部均为 \(0\)

构造 \(P_i=A_i+iB_i,Q_i=A_i-iB_i\)，并观察 \(P,Q\) 的 \(DFT\) 定义式子 \(P[k],Q[n-k]\) 发现是

\[P[k]=\sum_{i=0}^{n-1} (A[i]+ib[i])(\cos \frac{2\pi ik}n + isin \frac{2\pi ik}n) \]

\[Q[n-k]=\sum_{i=0}^{n-1} (A[i]-ib[i])(\cos \frac{2\pi ik}n - isin \frac{2\pi ik}n) \]

由于 \(\omega^{-jk}_n=\cos \frac{2\pi ik}n - isin \frac{2\pi ik}n\) 即三角函数里面变号，所以就解释得通了

对着展开就能发现是共轭复数，所以得到其一，翻转取共轭就得到了另一个

那么解方程组可以的到 \(DFT\) 的结果（复数科技妙）

得到 \(DFT\) 的结果之后，考虑怎么 \(IDFT\) 四个多项式得到最终答案

接着构造 \(P_i=A_0[i]B_0[i]+iA_1[i]B_0[i],Q_i=iA_1[i]B_1[i]+A_0[i]B_1[i]\)

因为这些个多项式系数都是实数，那么虚不会变成实，实不会变成虚，所以直接 \(IDFT\) 回去得到答案

精度不足就可以做三组 \(\rm FFT\) 来完成

很厉害的东西

inline void MTT(){
    for(int i=0;i<=n;++i) A0[i].x=a[i]>>15,A0[i].y=a[i]&32767;
    for(int i=0;i<=m;++i) B0[i].x=b[i]>>15,B0[i].y=b[i]&32767;
    int lim=1; while(lim<=(n+m)) lim<<=1; 
    for(int i=0;i<lim;++i) r[i]=r[i>>1]>>1|((i&1)?(lim>>1):0);
    FFT(A0,lim,1); FFT(B0,lim,1);
    for(int i=0;i<lim;++i) A1[i]=A0[i?lim-i:i].conj(),B1[i]=B0[i?lim-i:i].conj();
    node p,q;
    for(int i=0;i<lim;++i){
        p=A0[i],q=A1[i];
        A0[i]=(p+q)*node(0.5,0),A1[i]=(q-p)*node(0,0.5);
    }
    for(int i=0;i<lim;++i){
        p=B0[i],q=B1[i];
        B0[i]=(p+q)*node(0.5,0),B1[i]=(q-p)*node(0,0.5);
    }
    /*
    你可能对于是 (q-p)*node(0,0.5) 产生疑惑，复数的三角表示的运算规则，辐角相加
    原来就有一个90度，再转90度就到x轴负半轴了，所以要乘 -1 转回来
    */
    for(int i=0;i<lim;++i){
        P[i]=A0[i]*B0[i]*node(1,0)+A1[i]*B0[i]*node(0,1);
        Q[i]=A1[i]*B1[i]*node(0,1)+A0[i]*B1[i]*node(1,0);
    }
    FFT(P,lim,-1); FFT(Q,lim,-1);
    for(int i=0;i<=n+m;++i){
        P[i].x=P[i].x/lim+0.5;
        P[i].y=P[i].y/lim+0.5;
        Q[i].x=Q[i].x/lim+0.5;
        Q[i].y=Q[i].y/lim+0.5;
        int r1=(int)(P[i].x)%mod*(1ll<<30)%mod;
        int r2=((int)(P[i].y)%mod+(int)(Q[i].x)%mod)%mod*(1ll<<15)%mod;
        int r3=(int)(Q[i].y)%mod;
        print(add(r1,add(r2,r3));
    }
    puts("");
    return ;    
}

多项式求逆

已知 \(H(x)F(x)\equiv 1\mod x^{\frac n 2}\)，求

\[G(x)F(x)=1 \mod x^n \]

必有

\[\begin{aligned} G(x)F(x)\equiv1 \mod x^{\frac n 2}\\ G(x)-H(x))F(x)\equiv 0\mod x^{\frac n2}\\ (G(x)-H(x))^2\equiv 0\mod x^n\\ G^2(x)-2H(x)G(x)+H^2(x)\equiv 0\mod x^n\\ G^2(x)F(x)-2F(x)G(x)H(x)+H^2(x)\equiv 0 \mod x^n\\ \end{aligned}\]

这时候因为 \(G(x)F(x)\equiv 0\mod x^n\) 所以就得到了答案式子：

\[G(x)\equiv 2H(x)-F(x)H(x)^2\mod x^n \]

递归即可，边界为 \(G(0)\equiv \frac 1{F(0)} \mod 998244353\)，复杂度 \(\Theta(n\log n)\)

注意递归时是 \(\lceil\frac n2\rceil\) 再合并能得到 \(\bmod x^n\) 意义下的结果，下取整就错了

inline void solve(int Len){
    if(Len==1) return G[0]=ksm(F[0],mod-2),void(); solve((Len+1)>>1);
    int lim=1; while(lim<=(Len<<1)) lim<<=1;  
    for(int i=0;i<Len;++i) H[i]=F[i]; 
    for(int i=Len;i<lim;++i) H[i]=0;
    P.MTT(H,G,Len,lim); P.MTT(H,G,Len,lim);
    for(int i=0;i<Len;++i) G[i]=del(mul(G[i],2),H[i]);
    for(int i=Len;i<lim;++i) G[i]=0; 
    return ;
}

城市规划

设 \(dp[i]\) 表示当前 \(i\) 个点构成的无项联通图的数目，\(dp[1]=1\)

转移考虑 \(i\) 个点每个点都可以向下一个点连边，所以是 \(dp[i+1]=dp[i] \times (2^i-1)\)

但是 \((5,7),(4,7)\) 的连边方式就可以保证 \((4,5)\)联通，并不需要 \(1\to 4\) 的哪个点和 \(5\) 连

设 \(n\) 个点的不要求联通的无向联通图数目为 \(g(n)=2^{\frac{n(n-1)}2-1}\)

设 \(f(n)\) 为联通的无向图，那么就有如下的式子：

\[g(n)=\sum_{i=1}^n \binom {n-1}{i-1} g(n-i) f(i) \]

具体就是枚举第一个点所在的联通块里面有多少个点

设 \(F(x)=\frac{f(x)}{(x-1)!},G(x)=\frac{g(x)}{(x-1)!},H(x)=\frac{g(x)}{x!}\)

写成生成函数的式子发现就变成了卷积，多项式求逆即可

• 吃一堑涨一智：多项式求逆中的数组别用混

多项式 ln

\[G(x)=\ln(F(x)) \mod x^n \]

两边同时对复合函数求导得到

\[G'(x)=\ln'(F(x))F'(x) \mod x^n \]

\[G'(x)=\frac{F'(x)}{F(x)} \mod x^n \]

多项式求逆完了积分回去即可

注意在进行多项式 \(\ln\) 的过程中需要保证 \(F[0]=1\) 否则求导积分的过程会出问题，具体出什么问题我也不知道

inline void Dao(int *G,int *F,int len){
    for(int i=0;i<len-1;++i) G[i]=mul(F[i+1],i+1);
    return ;  
}
inline void Ji(int *G,int *F,int len){
    for(int i=0;i<len;++i) G[i+1]=mul(F[i],ksm(i+1,mod-2));
    return ;
}
inline void Ln(){
    P.Dao(G,F,n); P.Inv(F,H,n);
    int lim=1; while(lim<(n<<1)) lim<<=1; 
    P.NTT(G,lim,1); P.NTT(H,lim,1); 
    for(int i=0;i<lim;++i) G[i]=mul(G[i],H[i]); P.NTT(G,lim,-1);
    P.Ji(ans,G,n);
    return ;
}

多项式 exp

先推式子：\(F(x)=e^{A(x)}\mod x^n\)

\[\ln F(x)-A(x)=0\mod x^n \]

那么我们需要求出来的就是 \(G(F(x))=\ln F(x) -A(x)=0\)

这时候有个科技叫牛顿迭代，就是快速求函数零点的一个方法，本质上利用了导函数是对应值点处函数的斜率这一性质

问题是 \(F(G(x))=0\) 求 \(G(x)\)，牛顿迭代的方法是利用上次求出来的精度较小的 \(G_0(x)\) 得到 \(G(x)=G_0(x)-\frac{F(G_0(x))}{F'(G_0(x))}\)

注意！这里的分母是 \(F'(G_0(x))\) 而不是复合函数求导！！

这个式子比较好的地方是精度翻倍，也就是说对 \(\mod x^{\frac n2}\) 意义下做一次，现在精度就是 \(\mod x^n\) 意义下的了，证明可以用逆多项式同余意义下也是 \(\frac n2\) 次乘起来就翻倍了理解

正经得到证明需要泰勒展开一下

套回原来的式子，大概推导过程如下：

由于牛顿迭代的的精度翻倍，那么可以类似多项式求逆来进行递归，设已经得到了 \(F_0(x)\equiv e^{A(x)}\mod x^{\frac n2}\)，现尝试推 \(F(x)\equiv e^{A(x)}\mod x^n\)

带入上文中提及的 \(G(F(x))=\ln F(x)-A(x)=0\)

带入牛顿迭代表达式：

\[\begin{aligned} F(x)&=F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\\ &=F_0(x)-\frac{\ln F_0(x)-A(x)}{\frac{1}{F_0(x)}} \end{aligned}\]

注意这里只对 \(G(x)\) 求导，所以并不需要复合函数求导法则

inline void exp(int *a,int *b,int n){
    if(n==1) return b[0]=1,void(); exp(a,b,(n+1)>>1); 
    int lim=1; while(lim<=(n<<1)) lim<<=1; Ln(b,Tln,n); 
    for(int i=0;i<n;++i) Tln[i]=del(a[i],Tln[i]); 
    Tln[0]=add(Tln[0],1); 
    for(int i=n;i<lim;++i) Tln[i]=0;
    NTT(Tln,lim,1); NTT(b,lim,1); 
    for(int i=0;i<lim;++i) b[i]=mul(b[i],Tln[i]); 
    NTT(b,lim,-1); 
    for(int i=n;i<lim;++i) b[i]=0; return ;
}

---

\(\ln,\exp\) 在 \(\rm EGF\) 的表示下的一些模型中是为互逆运算：

• 在错排问题中，没有置换环的 \(\rm EGF\) 是 \(\sum_{i\ge 2}(i-1)! \frac{x^i}{i!}=-\ln(1-x)-x\)

  那么错排问题的 \(\rm EGF\) 就是 \(\exp(-\ln(1-x)-x)\)

• 无向连通图的 \(\rm EGF\) 为 \(F(x)\)，无向图（不保证连通） 的 \(\rm EGF\) 是 \(G(x)\) 那么有 \(G(x)=\exp(F(x)),F(x)=\ln G(x)\)

  这个也可以同理到森林 \(G(x)\) 和树 \(F(x)\) 上，也是 \(F(x)=\ln G(x)\) 的关系

  如果限制深度的话可以尝试递推：\(F_k(x)=xe^{F_{k-1}(x)},G_k(x)=\exp [xG_{k-1}(x)]\)

付公主的背包

关键题：乘法取对数变成加法

考虑构造出来生成函数封闭形式就是这个式子了

\[\frac{1}{e^{-\sum_{i=1}\ln (1-x^v)}} \]

设 \(F(x)=1-x^{V},G(x)=\ln F(x)\) 那么：

\[\begin{aligned} &G'(x)=\frac{F'(x)}{F(x)}\\ &G'(x)=-V\sum_{i\ge 0}x^{V-1-Vi}\\ &G(x)=-\sum_{i\ge 0}\frac{x^{V+Vi}}{1+i}\\ &G(x)=-\sum_{i\ge 1}\frac{x^{Vi}}i \end{aligned}\]

调和级数把多项式系数暴力得到之后写一份 \(exp\) 和 多项式求逆即可

多项式开根

设 \(H(x)^2\equiv F(x)\mod x^{\lceil \frac n2\rceil}\)，\(G(x)\equiv H(x) \mod x^{\lceil \frac n2\rceil}\)

剩下的推导和多项式求逆的过程一致，最后得到

\[G(x)=\frac{H^2(x)+F(x)}{2H(x)} \]

实现类似多项式求逆的递归即可

inline void sqr(int *A,int *B,int n){
    if(n==1) return B[0]=1,void(); sqr(A,B,(n+1)>>1); 
    int lim=1; while(lim<(n<<1)) lim<<=1; 
    for(int i=0;i<lim;++i) Tp[i]=0; 
    for(int i=0;i<n;++i) T2[i]=A[i]; 
    for(int i=n;i<lim;++i) T2[i]=0; 
    Inv(B,Tp,n); NTT(Tp,lim,1); NTT(B,lim,1); NTT(T2,lim,1); 
    for(int i=0;i<lim;++i) B[i]=mul(i2,add(B[i],mul(Tp[i],T2[i]))); 
    NTT(B,lim,-1); 
    for(int i=n;i<lim;++i) B[i]=0; 
    return ;
}

多项式快速幂

求 \(F(x)=G^k(x) \mod x^n\)，两边同时 \(\ln\) 则

\[k \ln G(x)=\ln F(x) \mod x^n \]

再来一把 \(exp\) 则答案就有了：

\[e^{k \ln G(x)}=F(x)\mod x^n \]

这个做法依赖于多项式 \(\ln\)，所以需要保证 \(F[0]=0\)

对于 \(F[0]\neq 0\) 的情况，平移即可，注意特判快速幂指数 \(k\) 过大导致的全是 \(0\) 的情况

多项式除法

设 \(F_R(x)=x^nF(\frac{1}x)\)，其实也就是系数翻转（没错，\(R\) 就是 \(Reverse\)）

\[F(x)=G(x)Q(x)+R(x) \]

\[F(\frac 1x)=G(\frac{1}x)Q(\frac 1 x )+R(\frac 1 x ) \]

配上 \(x^n\) 可以得到：

\[x^nF(\frac{1}x)=x^mG(\frac1 x)x^{n-m} Q(\frac{1}x)+x^{n-(m-1)} x^{m-1} R(\frac 1 x) \]

换式子：

\[F_R(x)=G_R(x)Q_R(x)+x^{n-(m-1)} R_R(x) \]

\[F_R(x)\equiv G_R(x)Q_R(x)\mod x^{n-(m-1)} \]

那么求逆得到 \(Q_R(x)\) 之后减法做出来 \(R(x)\) 就行了

写得时候记得清空一些数组

取模和除法是并行的

常系数齐次线性递推

建议阅读 https://www.luogu.com.cn/blog/BJpers2/solution-p4723

做法就是将用 \(x^n\) 模特征多项式（定义为递推数列取负翻转后再添 \(1\) 作为 \([x^k]P\)）

使用倍增求 \(x^n\) 来加速这个过程

Code Display

inline poly Mod(vector<int> F,vector<int> G){
    int n=F.size(),m=G.size(); if(n<m) return F;
    poly Fr=F,Gr=G,R;
    reverse(Fr.begin(),Fr.end());
    reverse(Gr.begin(),Gr.end()); Gr.resize(n-m+1);
    Gr=Inv(Gr,n-m+1);
    poly Q=Mul(Gr,Fr); Q.resize(n-m+1);
    reverse(Q.begin(),Q.end());
    Q=Mul(Q,G);
    rep(i,0,m-2) R.pb(del(F[i],Q[i]));
    return R;
}
int f[N],a[N],n,k;
signed main(){
    n=read(); k=read();
    rep(i,1,k) f[i]=(read()%mod+mod)%mod; 
    rep(i,0,k-1) a[i]=(read()%mod+mod)%mod;
    vector<int> p={1},x={0,1},pmod;
    rep(i,0,k-1) pmod.pb(del(0,f[k-i])); pmod.pb(1);
    for(;n;n>>=1,x=Mod(Mul(x,x),pmod)) if(n&1){
        p=Mod(Mul(p,x),pmod);
    }
    p.resize(k);
    int ans=0;
    rep(i,0,k-1) ckadd(ans,mul(a[i],p[i]));
    print(ans);
    return 0;
}