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

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;
}