『组合数学总结3：多项式算法』
组合数学知识点的总结第三篇
Preface
前排提示：本文数学公式较多，加载\(\LaTeX\)需要一定时间，可能会导致浏览器暂时卡顿，请耐心等待数学公式正常显示.
前两篇：『组合数学总结1：基础组合数学和组合原理』，『组合数学总结2：生成函数和特殊计数数列』
第二篇特殊计数数列部分暂时停更，本文完结后会继续更新.
应教练安排，本文短期内不会再更新了.
\(\mathrm{Update}\)：\(2020.5.15\) 远古博客重新开更了.
\(\mathrm{Update}\)：\(2020.6.8\) 本博客完结.
多项式基础操作
Fast Fourier Transform
用来计算两个多项式的卷积，即给定多项式\(A(x)=\sum_{i=0}^{n-1}a_ix^i,B(x)=\sum_{i=0}^{m-1}b_ix_i\)，求多项式\(C(x)=\sum_{i=0}^{n+m-1}(\sum_{j=0}^ia_jb_{i-j})x^i\)的各项系数.
思路：系数表示直接求卷积是\(\mathcal{O}(n^2)\)的，但是点值表示求卷积是\(\mathcal{O}(n)\)的，考虑选一些合适的点，然后快速插值.
对于多项式\(F(x)\)，我们可以定义序列算子\(\mathrm{DFT}\)为其在\(n\)次单位根下各个\(n\)点的点值，即：
注意到单位根的定义为\(\omega_n^k=\cos\left(\frac{2\pi k}{n}\right)+i\sin\left(\frac{2\pi k}{n}\right)=e^{\frac{2\pi k}{n}i}\)，具有如下的几个性质：
-
\(\omega_n^n=1\)，这是单位根的原始定义，也可以根据复数乘法的几何意义（辐角相加，模相乘）简单地推出.
-
\(\omega_{dn}^{kn}=\omega_n^k\)，根据欧拉公式的展开形式可以推出.
-
\(\omega_n^{k+\frac{n}{2}}=-\omega_n^k\)，根据几何意义可以推出，也可以根据三角函数公式推出.
-
\(\omega_n^k=\omega_n^{k\bmod n}\)，根据单位根的原始定义可以推出.
\(\mathrm{FFT}\)算法就是在\(O(n\log_2 n)\)的时间内对一个多项式进行\(\mathrm{DFT}\)运算，这样就实现了插值的效果.
以下推导要求\(n\)为\(2\)的整次幂，假设我们要操作的多项式为\(F(x)=\sum_{i=0}^{n-1}f_ix^i\)，令：
那么有\(F(x)=F_0(x^2)+xF_1(x^2)\)，现在根据我们的目标，带入\(x=\omega_n^{k}\)，得到\(F(\omega_n^k)=F_0(\omega_\frac{n}{2}^k)+\omega_n^kF_1(\omega_\frac{n}{2}^k)\)，带入\(x=\omega_n^{k+\frac{n}{2}}\)，得到：\(F(\omega_n^k)=F_0(\omega_\frac{n}{2}^k)-\omega_n^kF_1(\omega_\frac{n}{2}^k)\)，那么我们只需求\(F_0(x),F_1(x)\)的\(\mathrm{DFT}\)，于是可以分治求解，时间复杂度\(\mathcal{O}(n\log_2 n)\).
到这里还没有结束，我们知道一个多项式的\(\mathrm{DFT}\)后还要有对应的插值算法求出它的系数表示，我们设\(G(x)=\mathrm{DFT}\ F(x)=\sum_{i=0}^{n-1}F(\omega_n^{i})x^i\)，不妨求\(\mathrm{DFT}\ G(x)\)，得到：
当且仅当\(j+k\equiv 0\pmod n\)的时候最后一步等比数列求和不成立，有\(\sum_{i=0}^{n-1}\left(\omega_n^{j+k}\right)^i=n\)，综上有：
那么我们只要对\(\mathrm{DFT}\ F(x)\)再做一遍\(\mathrm{DFT}\)，再稍加处理就可以得到系数表示了，时间复杂度仍然是\(\mathcal{O}(n\log_2n)\).
观察分治算法的流程，我们发现分治树底层的标号是原序列的二进制翻转，所以我们可以事先进行二进制翻转，然后两两合并就可以实现非递归版的算法，代码如下：
#include <bits/stdc++.h>
using namespace std;
struct Complex
{
double Re,Im; Complex (double R = 0,double I = 0) { Re = R , Im = I; }
friend Complex operator + (Complex a,Complex b) { return Complex( a.Re + b.Re , a.Im + b.Im ); }
friend Complex operator - (Complex a,Complex b) { return Complex( a.Re - b.Re , a.Im - b.Im ); }
friend Complex operator * (Complex a,Complex b) { return Complex( a.Re*b.Re - a.Im*b.Im , a.Re*b.Im + a.Im*b.Re ); }
};
typedef vector<Complex> Poly; vector<int> rev;
const int N = 1<<21;
const long double Pi = acos(-1.0);
inline int Read(void)
{
int x = 0 , w = 0; char ch = ' ';
for (; !isdigit(ch); ch = getchar()) w |= ch == '-';
for (; isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
return w ? -x : x;
}
inline void FFT(Poly &p,int n,int op)
{
p.resize(n); static Complex w,omega,x,y;
for (int i = 0; i < n; i++) if ( rev[i] > i ) swap( p[rev[i]] , p[i] );
for (int i = 1; w = Complex( cos(Pi/i) , sin(Pi/i) ) , i < n; i <<= 1)
for (int j = 0; omega = Complex(1,0) , j < n; j += i<<1)
for (int k = 0; k < i; k++ , omega = omega * w)
x = p[j+k] , y = omega * p[i+j+k] , p[j+k] = x+y , p[i+j+k] = x-y;
if ( op == 1 ) return void(); reverse( p.begin() + 1 , p.end() );
for (int i = 0; i < n; i++) p[i].Re /= n , p[i].Im /= n;
}
inline int Calc(int n) { int res = 1; while ( res < n ) res <<= 1; return res; }
inline void Flip(int n)
{
int k = 0; while ( 1<<k < n ) ++k;
rev.resize(n) , rev[0] = 0 , --k;
for (int i = 1; i < n; i++)
rev[i] = ( rev[i>>1] >> 1 ) | ( (i&1) << k );
}
inline Poly Multiply(Poly a,Poly b)
{
if ( !a.size() || !b.size() ) return {};
int _n = a.size() + b.size() - 1 , n = Calc(_n);
Flip(n) , FFT(a,n,1) , FFT(b,n,1);
for (int i = 0; i < n; i++) a[i] = a[i] * b[i];
return FFT(a,n,-1) , a.resize(_n) , a;
}
int main(void)
{
int n = Read() , m = Read(); Poly a,b;
a.resize(n+1) , b.resize(m+1);
for (int i = 0; i <= n; i++) a[i].Re = Read();
for (int i = 0; i <= m; i++) b[i].Re = Read();
Poly c = Multiply( a , b );
for (int i = 0; i <= n + m; i++)
printf( "%d%c" , int(c[i].Re+0.5) , " \n"[i==n+m] );
return 0;
}
Fast Number-Theoretic Transform
当\(n=10^5\)，模数达到\(10^9\)级别的时，\(\mathrm{double}\)不能胜任\(10^{23}\)级别的精度，因此需要一种可以处理模意义下多项式卷积的算法，即\(\mathrm{NTT}\).
定义\(a\)模\(p\)的阶\(\delta_p(a)\)为最小的正整数\(x\)，满足\(a^x\equiv 1\pmod p\)，定义\(a\)为模\(p\)意义下的一个原根，当且仅当\(\delta_p(a)=\varphi(p)\).
原根具有如下三个性质：
- 若模\(p\)意义下存在原根，则一定存在\(\varphi(\varphi(p))\)个原根.
- 模\(p\)意义下原根存在，当且仅当\(p=2,4,k^a,2k^a\)，其中\(k\)为奇素数且\(a\geq 1\).
- 若\(p\)为素数，\(g\)为模\(p\)意义下的一个原根，则\(g^i\bmod p\ (i\in[0,p)\cap\Z)\)各不相同.
我们可以令\(\omega_n=g^{\frac{p-1}{n}}\bmod p\)，考察原本单位根具有的性质都还存在的条件：\(1.\) \(p\)是质数 \(2.\) \(p\)有原根 \(3.\) \(p-1\)是\(2^m\)的倍数，\(2^m\geq n\).
当模数\(p\)满足这些条件的时候我们就可以执行这样的\(\mathrm{DFT}\)算法，进行模意义下的运算，这样的模数又称为\(\mathrm{NTT}\)模数，常见的有：\(998244353=2^{23}\times7\times 17+1,1004535809=479\times 2^{21}+1\)，它们的原根都是\(3\).
#include <bits/stdc++.h>
using namespace std;
const int N = 1<<21 , Mod = 998244353 , G = 3;
inline int inc(int a,int b) { return a + b >= Mod ? a + b - Mod : a + b; }
inline int dec(int a,int b) { return a - b < 0 ? a - b + Mod : a - b; }
inline int mul(int a,int b) { return 1LL * a * b % Mod; }
inline void Inc(int &a,int b) { return void( a = inc( a , b ) ); }
inline void Dec(int &a,int b) { return void( a = dec( a , b ) ); }
inline void Mul(int &a,int b) { return void( a = mul( a , b ) ); }
inline int fastpow(int a,int b) { int c = 1; for ( ; b; b >>= 1 , Mul(a,a)) if (1&b) Mul(c,a); return c; }
inline int Inv(int x) { return fastpow( x , Mod - 2 ); }
typedef vector<int> Poly; Poly rev;
inline int Read(void)
{
int x = 0 , w = 0; char ch = ' ';
while ( !isdigit(ch) ) w |= ch == '-' , ch = getchar();
while ( isdigit(ch) ) x = x * 10 + ch - 48 , ch = getchar();
return w ? -x : x;
}
inline int Calc(int n) { int p = 1; while ( p < n ) p <<= 1; return p; }
inline void Flip(int n)
{
int k = 1; while ( 1<<k < n ) k++;
rev.resize(n) , rev[0] = 0 , k--;
for (int i = 1; i < n; i++)
rev[i] = ( rev[i>>1] >> 1 ) | ( (i&1) << k );
}
inline void NTT(Poly &p,int n,int op)
{
p.resize(n); static int w,omega,x,y;
for (int i = 0; i < n; i++) if ( rev[i] > i ) swap( p[rev[i]] , p[i] );
for (int i = 1; w = fastpow( G , (Mod-1) / (i<<1) ) , i < n; i <<= 1)
for (int j = 0; omega = 1 , j < n; j += i<<1)
for (int k = 0; k < i; k++ , Mul( omega , w ))
x = p[j+k] , y = mul( omega , p[i+j+k] ) , p[j+k] = inc(x,y) , p[i+j+k] = dec(x,y);
if ( op == 1 ) return void(); reverse( p.begin() + 1 , p.end() );
for (int i = 0 , Invn = Inv(n); i < n; i++) Mul( p[i] , Invn );
}
inline Poly Multiply(Poly a,Poly b)
{
if ( !a.size() || !b.size() ) return {};
int _n = a.size() + b.size() - 1 , n = Calc(_n);
Flip(n) , NTT(a,n,1) , NTT(b,n,1);
for (int i = 0; i < n; i++) Mul( a[i] , b[i] );
return NTT(a,n,-1) , a.resize(_n) , a;
}
int main(void)
{
int n = Read() , m = Read(); Poly a,b;
a.resize(n+1) , b.resize(m+1);
for (int i = 0; i <= n; i++) a[i] = Read();
for (int i = 0; i <= m; i++) b[i] = Read();
a = Multiply( a , b );
for (int i = 0; i < a.size(); i++)
printf( "%d%c" , a[i] , " \n"[i==a.size()] );
return 0;
}
Split Fast Fourier Transform
拆系数\(\mathrm{FFT}\)，俗称\(\mathrm{MTT}\)，用于计算任意模数但是值域又比较大的多项式卷积。
对于\(A(x)\)和\(B(x)\)的卷积，令\(A(x)=2^{15}A_0(x)+A_1(x),B(x)=2^{15}B_0(x)+B_1(x)\)，那么：
直接算的话要序列\(\mathrm{DFT}\)和\(\mathrm{IDFT}\)总共\(7\)次，常数太大，考虑优化。
结论：对于对应系数互为共轭复数的多项式\(P(x),Q(x)\)，有\((\mathrm{DFT}\ P(x))_i=\overline{(\mathrm{DFT}\ Q(x))_{(n-i)\bmod n}}\).
推导过程略繁琐，利用上述结论，直接将\(\mathrm{DFT}\)的过程两两合并，就可以优化到\(2\)次\(\mathrm{DFT}\).
\(\mathrm{DFT}\)后的序列数字都是复数，实部虚部都有，所以不能用上述优化合并\(\mathrm{IDFT}\)，但是鉴于\(\mathrm{IDFT}\)之后只剩下实部，所以可以直接一个放在虚部一个放在实部进行合并，同样只需\(2\)次\(\mathrm{IDFT}\)。
#include <bits/stdc++.h>
using namespace std;
struct Complex
{
double Re,Im; Complex (double R = 0,double I = 0) { Re = R , Im = I; }
friend Complex operator + (Complex a,Complex b) { return Complex( a.Re + b.Re , a.Im + b.Im ); }
friend Complex operator - (Complex a,Complex b) { return Complex( a.Re - b.Re , a.Im - b.Im ); }
friend Complex operator * (Complex a,Complex b) { return Complex( a.Re*b.Re - a.Im*b.Im , a.Re*b.Im + a.Im*b.Re ); }
};
const long double Pi = acos(-1);
const int N = 1<<18;
typedef vector<int> Poly;
typedef long long ll;
Poly rev; Complex omega[N<<1]; int Mod;
inline int Read(void)
{
int x = 0 , w = 0; char ch = ' ';
while ( !isdigit(ch) ) w |= ch == '-' , ch = getchar();
while ( isdigit(ch) ) x = x * 10 + ch - 48 , ch = getchar();
return w ? -x : x;
}
inline Complex Conj(Complex z) { return Complex( z.Re , -z.Im ); }
inline int Calc(int n) { int k = 1; while ( k < n ) k <<= 1; return k; }
inline void Flip(int n)
{
int k = 0; while ( 1<<k < n ) ++k;
rev.resize(n) , rev[0] = 0 , --k;
for (int i = 1; i < n; i++)
rev[i] = ( rev[i>>1] >> 1 ) | ( (i&1) << k );
}
inline void Init(void)
{
for (int i = 1 , j; i < N; i <<= 1)
for (j = 1 , omega[i] = Complex(1,0); j < i; j++)
if ( (j&31) == 1 ) omega[i+j] = Complex( cos(Pi*j/i) , sin(Pi*j/i) );
else omega[i+j] = omega[i+j-1] * omega[i+1];
}
inline void FFT(vector<Complex> &p,int n,int op)
{
p.resize(n); static Complex x,y;
for (int i = 0; i < n; i++) if ( rev[i] > i ) swap( p[rev[i]] , p[i] );
for (int i = 1; i < n; i <<= 1) for (int j = 0; j < n; j += i<<1)
for (int k = 0; x = p[j+k] , y = omega[i+k] * p[i+j+k] , k < i; k++)
p[j+k] = x + y , p[i+j+k] = x - y;
if ( op == 1 ) return void(); reverse( p.begin() + 1 , p.end() );
for (int i = 0; i < n; i++) p[i].Re /= n , p[i].Im /= n;
}
inline Poly MTT(Poly a,Poly b)
{
if ( !a.size() || !b.size() ) return {};
int _n = a.size() + b.size() - 1 , n = Calc(_n); Flip(n);
static vector<Complex> A,B,P,Q; static Complex a0,a1,b0,b1;
A.resize(n) , B.resize(n) , P.resize(n) , Q.resize(n) , a.resize(n) , b.resize(n);
for (int i = 0; i < n; i++)
A[i] = Complex( a[i]>>15 , a[i] & ( (1<<15) - 1 ) ),
B[i] = Complex( b[i]>>15 , b[i] & ( (1<<15) - 1 ) );
FFT(A,n,1) , FFT(B,n,1) , P[0] = Conj(A[0]) , Q[0] = Conj(B[0]);
for (int i = 1; i < n; i++) P[i] = Conj(A[n-i]) , Q[i] = Conj(B[n-i]);
for (int i = 0; i < n; i++)
a0 = ( A[i] + P[i] ) * Complex(0.5,0) , a1 = ( P[i] - A[i] ) * Complex(0,0.5),
b0 = ( B[i] + Q[i] ) * Complex(0.5,0) , b1 = ( Q[i] - B[i] ) * Complex(0,0.5),
A[i] = a0 * b0 + Complex(0,1) * a0 * b1 , B[i] = b0 * a1 + Complex(0,1) * a1 * b1;
FFT(A,n,-1) , FFT(B,n,-1);
static Poly c; static long long x,y,z; c.resize(_n);
for (int i = 0; i < _n; i++)
x = (ll)( A[i].Re + 0.5 ) % Mod , z = (ll)( B[i].Im + 0.5 ) % Mod ,
y = (ll)( A[i].Im + B[i].Re + 0.5 ) % Mod , c[i] = ( (x<<30) + (y<<15) + z ) % Mod;
return c;
}
int main(void)
{
int n,m; Poly a,b,c;
n = Read() , m = Read() , Mod = Read() , Init();
a.resize(n+1) , b.resize(m+1);
for (int i = 0; i <= n; i++) a[i] = Read();
for (int i = 0; i <= m; i++) b[i] = Read();
c = MTT( a , b );
for (int i = 0; i <= n+m; i++)
printf( "%d%c" , c[i] , " \n"[i==n+m] );
return 0;
}
Poly Inverse
已知多项式\(F(x)\)，求一个多项式\(G(x)\)满足\(F(x)G(x)\equiv 1\pmod {x^n}\).
假设我们求出了\(H(x)\)满足\(F(x)H(x)\equiv 1\pmod{x^{\frac{n}{2}}}\)，那么有：
可以倍增求解，时间复杂度\(T(n)=T(\frac{n}{2})+\mathcal{O}(n\log_2 n)=\mathcal{O}(n\log_2 n)\).
inline void Init(void)
{
for (int i = 1 , j; i < N; i <<= 1)
for (j = 1 , omega[i] = 1; j < i; j++)
if ( j == 1 ) omega[i+j] = fastpow( G , (Mod-1) / (i<<1) );
else omega[i+j] = mul( omega[i+j-1] , omega[i+1] );
}
inline void NTT(Poly &p,int n,int op)
{
p.resize(n); static int x,y;
for (int i = 0; i < n; i++) if ( rev[i] > i ) swap( p[rev[i]] , p[i] );
for (int i = 1; i < n; i <<= 1) for (int j = 0; j < n; j += i<<1)
for (int k = 0; x = p[j+k] , y = mul(p[i+j+k],omega[i+k]) , k < i; ++k)
p[j+k] = inc( x , y ) , p[i+j+k] = dec( x , y );
if ( op == 1 ) return void(); reverse( p.begin() + 1 , p.end() );
for (int i = 0 , Invn = Inv(n); i < n; i++) Mul( p[i] , Invn );
}
inline Poly Inverse(Poly f,int n)
{
int _n = Calc(n) , Len; static Poly g,h;
g.resize(1) , f.resize(_n) , g[0] = Inv(f[0]);
for (int i = 2; i <= _n; i <<= 1)
{
h.resize(i) , Len = i << 1 , Flip(Len);
for (int j = 0; j < i; j++) h[j] = f[j];
NTT( h , Len , 1 ) , NTT( g , Len , 1 );
for (int j = 0; j < Len; j++)
g[j] = mul( g[j] , dec( 2 , mul( h[j] , g[j] ) ) );
NTT( g , Len , -1 ) , g.resize(i);
}
return g.resize(n) , g;
}
Divide FFT in Poly Inverse
给定序列 \(g_{1\cdots n - 1}\)，求序列 \(f_{0\cdots n - 1}\)，其中 \(f_i=\sum_{j=1}^if_{i-j}g_j\)，边界为 \(f_0=1\).
此类问题也称为半在线卷积，可以使用多种方法解决.
令\(g_0=0\)，设序列\(f,g\)的生成函数分别为\(F(x),G(x)\)，则：
在\(\pmod{x^n}\)的意义下多项式求逆即可，时间复杂度\(\mathcal{O}(n \log_2n)\).
Poly Division
已知\(n\)次多项式\(F(x)\)和\(m\)次多项式\(G(x)\)，求多项式\(Q(x),R(x)\)，满足：
\(1.\) \(Q(x)\)为\(n-m\)次多项式，\(R(x)\)的次数小于\(m\)
\(2.\) \(F(x)=G(x)\times Q(x)+R(x)\)
定义\(F^{rev}(x)\)为多项式\(F(x)\)系数翻转后的表达式，即：\(F^{rev}(x)=x^nF\left(\frac{1}{x}\right)\).
考虑对\(G^{rev}\)做一个多项式求逆，然后再把它和\(F^{rev}\)乘起来就是\(Q^{rev}\)了.
然后\(R\)直接算就可以了：\(R(x)=F(x)-Q(x)\times G(x)\)，时间复杂度\(\mathcal{O}(n\log_2n)\).
inline Poly Minus(Poly a,Poly b)
{
int n = max( a.size() , b.size() );
a.resize(n) , b.resize(n); Poly c = a;
for (int i = 0; i < c.size(); i++) Dec( c[i] , b[i] );
return c;
}
inline pair<Poly,Poly> Division(Poly a,Poly b)
{
int n = a.size() , m = b.size(); static Poly Q,R,A,B;
if ( n < m ) return R = a , R.resize(m-1) , make_pair(Poly{},R);
A = a , reverse( A.begin() , A.end() ) , A.resize( n - m + 1 );
B = b , reverse( B.begin() , B.end() ) , B.resize( n - m + 1 );
Q = Multiply( A , Inverse( B , n - m + 1 ) );
Q.resize( n - m + 1 ) , reverse( Q.begin() , Q.end() );
R = Minus( a , Multiply( b , Q ) ) , R.resize( m - 1 );
return make_pair( Q , R );
}
int main(void)
{
int n = Read() , m = Read(); Poly A(n+1),B(m+1); Init();
for (int i = 0; i <= n; i++) A[i] = Read();
for (int i = 0; i <= m; i++) B[i] = Read();
Poly Q,R; tie( Q , R ) = Division( A , B );
for (int i = 0; i < Q.size(); i++) printf( "%d%c" , Q[i] , " \n"[i==Q.size()-1] );
for (int i = 0; i < R.size(); i++) printf( "%d%c" , R[i] , " \n"[i==R.size()-1] );
return 0;
}
Poly Sqrt
已知多项式\(A(x)\)，求多项式\(B(x)\)使得\(B^2(x) \equiv A(x) \pmod{x^n}\)
假设我们求出了\(B'(x)\)满足\(B'^2(x)\equiv A(x)\pmod{x^{\frac{n}{2}}}\)，那么有：
可以倍增求解，时间复杂度\(T(n)=T(\frac{n}{2})+\mathcal{O}(n\log_2 n)=\mathcal{O}(n\log_2 n)\).
inline Poly Pow2(Poly p)
{
int n = p.size() * 2 - 1 , _n = Calc(n);
Flip(_n) , NTT( p , _n , 1 );
for (int i = 0; i < _n; i++) p[i] = mul( p[i] , p[i] );
return NTT( p , _n , -1 ) , p.resize(n) , p;
}
inline Poly Sqrt(Poly f,int n)
{
int _n = Calc(n); static Poly g,h;
g.resize(1) , f.resize(_n) , g[0] = 1;
for (int i = 2; i <= _n; i <<= 1)
{
h = Pow2(g) , h.resize(i);
for (int j = 0; j < i; j++) Inc( h[j] , f[j] );
g = Multiply( h , Inverse( g , i ) ) , g.resize(i);
for (int j = 0; j < i; j++) Mul( g[j] , inv[2] );
}
return g.resize(n) , g;
}
Poly Ln
已知多项式\(A(x)\)，求多项式\(B(x)\)使得\(B(x) \equiv \ln A(x) \pmod{x^n}\)
对等式两边同时求导，得到：
多项式的积分和微分不难\(\mathcal{O}(n)\)实现，所以时间复杂度为\(\mathcal{O}(n\log_2n)\).
inline Poly Derivative(Poly p)
{
int n = p.size();
for (int i = 1; i < n; i++) p[i-1] = mul( p[i] , i );
return p.resize(n-1) , p;
}
inline Poly Integral(Poly p)
{
p.push_back(0); int n = p.size();
for (int i = n - 1; i >= 1; i--) p[i] = mul( p[i-1] , inv[i] );
return p[0] = 0 , p;
}
inline Poly Ln(Poly p,int n)
{
Poly f = Integral( Multiply( Derivative(p) , Inverse(p,n) ) );
return f.resize(n) , f;
}
Newton-Raphson Method
首先，一般的牛顿迭代是一种用于求解非线性方程\(F(x)=0\)近似解的方法，现在我们用它来求解多项式方程在模意义下的解，所以我们可以先定义\(F:P\rightarrow P\)表示一个多项式到多项式的函数。
考虑上述很多个算法都用到了的倍增思想，假设我们已经求出了\(\pmod {x^{\frac{n}{2}}}\)意义下的答案\(G_0(x)\)，满足\(F(G_0(x))\equiv0\pmod{x^{\frac{n}{2}}}\)，我们假设现在要的答案是\(G(x)\)，满足\(F(G(x))\equiv 0\pmod{x^n}\). 不妨将上式在\(G_0(x)\)处泰勒展开：
因为\(G(x)-G_0(x)\equiv 0\pmod{x^{\frac{n}{2}}}\)，所以有\(\left(G(x)-G_0(x)\right)^2\equiv 0\pmod{x^n}\)，那么上式就可以直接化简成：
整理得到如下递推公式：
值得注意的是上式中\(F'(x)=\frac{\mathrm{d}F(x)}{\mathrm{d}x}\)，而非\(\frac{\mathrm{d}F(G(x))}{\mathrm{d}x}\)，所以牛顿迭代法是我们针对某个具体的函数\(F\)来找递推公式的手段，\(F'(x)\)也要我们提前手算好.
不妨使用\(\mathrm{Poly\ Sqrt}\)的例子，令多项式函数\(F(t)=t^2-A(x)\)，则求解\(B^2(x)\equiv A(x)\pmod{x^n}\)等价于求解方程\(F(t)\equiv 0\pmod{x^n}\)的一个零点\(B(x)\)，带入牛顿迭代递推公式，得到：
与之前我们直接求出的递推方程一致.
Poly Exp
已知多项式\(A(x)\)，求多项式\(B(x)\)使得\(B(x) \equiv \exp A(x) \pmod{x^n}\)
对等式两边同时取对数，得到\(\ln B(x)\equiv A(x)\pmod{x^n}\)，令\(F(t)=\ln t-A(x)\)，则要求多项式函数\(F(t)\)的零点\(B(x)\)，带入牛顿迭代递推公式：
\(\mathcal{O}(n\log_2n)\) 递推求解即可.
inline Poly Exp(Poly f,int n)
{
int _n = Calc(n) , Len; static Poly g,h;
f.resize(_n) , g.resize(1) , g[0] = 1;
for (int i = 2; i <= _n; i <<= 1)
{
h = Ln( g , i ) , Len = i << 1 , Flip(Len);
for (int j = 0; j < i; j++) h[j] = dec( f[j] , h[j] );
Inc( h[0] , 1 ) , NTT( h , Len , 1 ) , NTT( g , Len , 1 );
for (int j = 0; j < Len; j++) g[j] = mul( g[j] , h[j] );
NTT( g , Len , -1 ) , g.resize(i);
}
return g.resize(n) , g;
}
Poly Pow
已知多项式\(A(x)\)，且\([x^0]A(x)=1\)，求多项式\(B(x)\)使得\(B(x) \equiv A^k(x) \pmod{x^n}\)
这个问题很简单，只要利用幂指恒等式转化一下就可以求解：
当然前提条件是\([x^0]A(x)=1\)，否则\(\ln A(x)\)无意义，需要用一些特殊的处理技巧，据代码实践，不是很\(\mathrm{practical}\)，意义也不大.
inline Poly Fastpow(Poly f,int n,int k)
{
f = Ln( f , n );
for (int i = 0; i < n; i++) Mul( f[i] , k );
return f = Exp( f , n ) , f;
}
Lagrange Inversion
对于两个多项式\(F(x)\)和\(G(x)\)常数项为\(0\)且一次项不为\(0\)，若有\(F(G(x))=x\)或\(G(F(x))=x\)，则称\(F(x)\)和\(G(x)\)互为复合逆.
对于一个给定的多项式\(F(x)\)，求其复合逆还没有很好的办法，据说最快的时间复杂度也只做到了\(\mathcal{O}\left((n\log_2n)^{1.5}\right)\)，而且代码复杂，实际意义不大，但是对于求每一个多项式复合逆一项的系数，根据拉格朗日反演可以很快地求解：
这其实是一个抽象代数里面的结论，所以会有\([x^{-1}]\)这种奇怪的东西，至于证明，可以参考这里.
在多项式运算里面，如果\(F(x)\)常数项为\(0\)且一次项不为\(0\)，可以令\(F'(x)=F(x)/x\)，于是就有网上很多博客所写的结论了：
所以写代码时，对于\(F'(x)=F(x)/x\)这种多项式，要额外注意细节.
inline int Lagrange(Poly f,int n)
{
f = Inverse( f , n ) , f = Fastpow( f , n , n );
return mul( f[n-1] , inv[n] );
}
Divide FFT in cdq's Divide Algorithm
给定序列 \(g_{1\cdots n - 1}\)，求序列 \(f_{0\cdots n - 1}\)，其中 \(f_i=\sum_{j=1}^if_{i-j}g_j\)，边界为 \(f_0=1\).
这个问题我们在上文中已经使用生成函数方法解决过了，但是更通用的方法实际上是\(\mathrm{cdq}\)分治，对于一般形式的半在线卷积，生成函数方法可能无法奏效.
首先我们\(\mathrm{cdq}\)分治的思想来设计分治函数，由于贡献可以拆分，所以我们定义\(\mathrm{d}(l,r)\)为计算区间\([l,r)\)的\(f\)值（左闭右开区间），先计算\(\mathrm{d}(l,\lfloor\frac{r+l+1}{2}\rfloor)\)，求出左边区间对右边区间的贡献，再计算\(\mathrm{d}(\lfloor\frac{r+l+1}{2}\rfloor,r)\)，计算贡献的部分就可以用\(\mathrm{FFT}\)优化，于是达到求解原问题的目的.
我们详细来看看，首先令\(\mathrm{mid}=\lfloor\frac{l+r+1}{2}\rfloor\)，则对于右半区间的某个点\(x\)，左半区间对它的贡献为：
这时候就有一个问题，我们明明只知道\([l,\mathrm{mid})\)内的\(f\)值，可是求和式里却出现了\(i\geq \mathrm{mid}\)的情况，事实上，这个求和上下界是为了方便后面卷积的推导，尽管这个时候\(f_i(i\geq \mathrm{mid})\)的值还不知道，当做\(0\)处理也不影响答案.
上面的上限界转化很显然，但是我们要的卷积式还不是这样的，我们更希望有这样的卷积式：
这时候我们就要找\(f\)和\(f'\)，\(g\)和\(g'\)的关系，以解决卷积：显然\(f'_i=f_{i+l}\)，\(g'_i=g_{i+1}\). 根据\(x\in [\mathrm{mid},r)\)以及求和号上下界\(i\in[0,x-l-1]\subseteq[0,r-l-1)\)，只需处理出\(f'_{0\sim \mathrm{mid}-l-1},g'_{0\sim r-l-2}\)，求其卷积，把第\(x-l-1\)项系数贡献到\(f_x\)就好了.
鉴于\(\mathcal{O}(n\log^2n)\)的时间复杂度过大，我们通常还会加上这样一个优化：
注意到\(x\in[\mathrm{mid},r)\)，所以\(x-l-1\in[\mathrm{mid}-l-1,r-l-1)\)，那么我们不在乎序列卷积第\(r-l-1\)项及以后那些项的系数是多少，可以不做计算.
这样对常数优化很有帮助，可以减少近乎一半的运行时间.
inline Poly Multiply(Poly a,Poly b,int m = 0)
{
if ( !a.size() || !b.size() ) return {};
int n = !m ? a.size() + b.size() - 1 : m , _n = Calc(n);
Flip(_n) , NTT( a , _n , 1 ) , NTT( b , _n , 1 );
for (int i = 0; i < _n; i++) Mul( a[i] , b[i] );
return NTT( a , _n , -1 ) , a.resize(n) , a;
}
inline void cdqDivide(int l,int r,Poly &f,const Poly &g)
{
if ( l + 1 == r ) return void( l?:f[0] = 1 );
int mid = ( l + r + 1 ) >> 1; Poly h;
cdqDivide( l , mid , f , g );
// According to relationships between g and g' , f and f'
// The sections are f[l,mid-1] and g[1,r-l-1]
// Remember that the section in C++ don't include the right endpoint
// So the sections are f[l,mid) and g[1,r-l)
h = Multiply( Poly( f.begin() + l , f.begin() + mid ) ,
Poly( g.begin() + 1 , g.begin() + r - l ) , r - l - 1 );
for (int i = mid; i < r; i++) Inc( f[i] , h[i-l-1] );
return cdqDivide( mid , r , f , g );
}
Poly Exp in cdq's Divide Algorithm
上述分治乘法在求解多项式指数函数上也有很好的成效，我们可以进行如下的推导：
同时求左右式\(x^n\)项系数：
其中\(a'_i=[x^i]A'(x)\)，那么问题又转化到了半在线卷积问题，使用类似的方法推导也可以在\(\mathcal{O}(n\log^2 n)\)内求解，这时候代码就比牛顿迭代法求\(\exp\)要好写很多，运行时间也不劣于牛顿迭代法.
inline void cdqDivide(int l,int r,Poly &f,const Poly &g)
{
if ( l + 1 == r ) return void( f[l] = !l ? 1 : mul( f[l] , inv[l] ) );
int mid = ( l + r + 1 ) >> 1; Poly h;
cdqDivide( l , mid , f , g );
h = Multiply( Poly( f.begin() + l , f.begin() + mid ) ,
Poly( g.begin() , g.begin() + r - l - 1 ) , r - l - 1 );
for (int i = mid; i < r; i++) Inc( f[i] , h[i-l-1] );
return cdqDivide( mid , r , f , g );
}
inline Poly Exp(Poly f,int n)
{
f.resize(n) , f = Derivative(f); Poly g(n);
return cdqDivide( 0 , n , g , f ) , g;
}
例题
[Luogu 5748] 集合划分计数
题意
求\(\mathrm{Bell}_n\)，\(n\leq 10^5\)，多组询问，\(T\leq 10^5\).
题解
在『组合数学总结2：生成函数和特殊计数数列』一文中推导过\(\mathrm{Bell}_n\)的指数型生成函数\(\hat W(x)=e^{e^x-1}\)，预处理\([x^n](e^x-1)=\frac{[n\not =0]}{n!}\)，直接用多项式\(\exp\)算出来就好了，最后乘上\(n!\)即可.
这道题也帮我们复习了指数型生成函数的意义：对于拼接是二项卷积的组合对象，我们先规定其生成函数格式为\(F(x)=\sum_{i=0}^\infty f_i\frac{x^i}{i!}\)，然后用普通卷积实现二项卷积，最后再乘上\(n!\)还原成原来的答案.
[CF623E] Transforming Sequence
题意
对于一个整数序列 \(\{a_1, a_2, \ldots, a_n\}\)，我们定义它的变换为 \(\{b_1, b_2, \ldots, b_n\}\)，其中 \(b_i = a_1 | a_2 | \ldots | a_i\)，其中 \(|\) 表示二进制按位或运算.
现在求有多少个长为 \(n\) 的序列 \(a\)，满足 \(1\leq a_i \leq 2^k - 1\)，使得它的变换 \(b\) 是严格单调递增的，对 \(10^9+7\) 取模.
\(1\leq n \leq 10^{18}\)，\(1\leq k \leq 3 \times 10^4\).
题解
首先应该想到每次在原数列后面加一个数字，应该至少出现一个数位上的\(1\)使得前面所有的数字在这一位上没有\(1\)，这样才能保证\(b\)序列递增. 于是可以想到这样一个朴素的\(\mathrm{dp}\)：设\(f_{i,j}\)代表序列长度为\(i\)，有\(j\)个位置占用了\(1\)（不考虑是哪\(j\)个，只在乎有\(j\)个\(1\)），状态转移方程如下：
组合意义就是考虑新增\(j-p\)个不同位的\(1\)，那么之前的\(p\)个位置就可以随便选或不选.
组合数可以拆掉用\(\mathrm{FFT}\)加速卷积，时间复杂度\(\mathcal{O}(nk\log_2 k)\)，当然这样还不够快.
考虑到有一个\(2^p\)的系数，不便于转成生成函数求解，但是我们仍然能够发现每一个位的转移的平等的，可以考虑推式子，能否一次转移多个位：
组合意义是一样的，同样可以用卷积加速，那么对于\(f_n\)的求解，类似于快速幂的方法二进制分解一下\(n\)，然后求卷积递推即可，这里还需要用到\(\mathrm{MTT}\)，最后的答案是\(\sum_{i=1}^k \binom{k}{i}f_{n,i}\).
[ZJOI2014] 力
题意
给定序列\(q\)，对于\(i=1\sim n\)，求\(E_i=\sum_{j=1}^{i-1}\frac{q_j}{(i-j)^2}-\sum_{j=i+1}^n\frac{q_j}{(i-j)^2}\)，\(n\leq 10^5\).
题解
首先应该想到化简式子，令\(g_i=i^{-2}\)，\(q_0=0\)，则
式子前半部分已经是卷积的形式了，我们可以用\(\mathrm{FFT}\)算，后半部分还不是，可以把它像卷积的形式靠拢：
不妨令\(n-i=t\)，则\(q_{i+j}g_j=q_{j+n-t}g_j\)，考虑到两个序列的卷积要使其下标和为\(t\)，那么仅需让\(q_{x}=q'_{n-x}\)，就有\(q_{j+n-t}g_j=q'_{t-j}g_j\)，那么预处理\(q'\)和\(g\)两个数组，然后就可以用卷积算了。
[Luogu 4233] 射命丸文的笔记
题意
求所有\(n\)个点带标号有哈密顿回路的竞赛图中哈密顿回路数量的平均值，\(n\leq 10^5\).
题解
首先所有竞赛图里哈密顿回路的总数量很好求，\(s=\frac{n!}{n}\times 2^{\binom{n}{2}-n}\)，组合意义就是枚举\(n\)个点在哈密顿回路上环排列的方案数，再将剩下\(\binom{n}{2}-n\)条边随意定向即可.
那么题意转化为求有哈密顿回路的竞赛图总数，显然，有哈密顿回路的竞赛图是强连通图，容易证明，一张强联通竞赛图必有哈密顿回路，那么问题转化为强联通竞赛图计数.
考虑到竞赛图里每两个点的相对关系都是确定的，所以其\(\mathrm{SCC}\)缩点后形成的\(\mathrm{DAG}\)是一条链，那么我们很容易用\(n\)个点强联通竞赛图的数量（设为\(f_n\)）来表示\(n\)个点竞赛图总数（设为\(g_n\)），可以枚举拓扑序最小的\(\mathrm{SCC}\)大小，得到：
不妨设\(\{f_n\}\)的指数型生成函数为\(\hat F(x)\)，\(\{g_n\}\)的指数型生成函数为\(\hat G(x)\)，那么我们很容易想到要写出\(\hat F(x)\times \hat G(x)\)和递推式\((1)\)取比对，即在\(n\geq 1\)时，有：
所以我们有：
至于为什么左边要减掉\(g_0\)，右边要减掉\(f_0g_0\)，那是因为递推式\((1)\)仅在\(n\geq 1\)时成立，所以把等式两边\(x^0\)项系数先扣掉. 根据递推式\((1)\)，\(g_0\)有意义，要求\(g_0=1\)，所以可以化简得到：
回顾可知，我们推出式\((2)\)使用的唯一条件就是递推式\((1)\)成立，而在原题中\(f_0\)没有意义，且\(f_0\)的意义在递推式中也得不到体现（即\(f_0\)取什么值都不会影响到递推式正确性），为了计算方便，我们完全可以令\(f_0=0\)，所以有：
预处理出\(g_n=\frac{2^{n\times (n-1)/2}}{n!},g_0=g_1=1\)，在\(\pmod {x^n}\)意义下多项式求逆即可.
[BZOJ 3625] 小朋友与二叉树
题意
定义一棵二叉树是合法的，当且仅当每个节点的权值都在给定正整数集合\(\mathrm{C}\)中. 对于每个\(S ∈ [1,m]\)，求有多少棵不同的，权值和为\(S\)的合法二叉树. \(|\mathrm{C}|,m ≤ 10^5\)，\(\mathrm{C}\)中的元素不超过\(10^5\).
题解
设权值和为\(i\)的二叉树有\(f_i\)种，则可以写出如下状态转移方程：
右边的求和显然是一个卷积，设\(\{f_n\}\)的生成函数为\(F(x)\)，集合\(\mathrm{C}\)的\(01\)生成函数为\(C(x)\)，那么根据状态转移方程，有\(F(x)=F^2(x)C(x)+1\)，解得\(F(x)=\frac{1\pm \sqrt{1-4C(x)}}{2C(x)}\)，这个生成函数和卡特兰数的生成函数是基本一样的，求一下\(x\rightarrow0\)的极限发现应该取\(+\)号.
那么处理出\(C(x)\)，在\(\pmod {x^{m+1}}\)意义下多项式开根，多项式求逆即可.
[BZOJ 3684] 大朋友和多叉树
题意
定义一棵树是合法的，当且仅当每个节点的孩子个数都在给定集合\(\mathrm{C}\)中. 求有多少棵不同的，非叶节点个数为\(s\)的合法树. 定义只有一个点的树一定合法，非叶节点个数视为为\(1\). \(|\mathrm{C}|, s ≤ 10^5\)，\(\mathrm{C}\)中的元素不超过\(10^5\).
题解
设非叶节点个数为\(i\)的树数量为\(f_i\)，\(f_0=0\)，则可以写出如下状态转移方程：
设\(\{f_n\}\)的生成函数为\(F(x)\)，则\(F(x)=x+\sum_{k\in \mathrm{C}}F^k(n)\)，移项可得\(F(x)-\sum_{k\in \mathrm{C}}F^k(n)=x\)，不妨设\(G(t)=t-\sum_{k\in\mathrm{C}}t^k\)，则要求多项式\(G(x)\)的复合逆\(F(x)\)，使其满足\(G(F(x))=x\)，根据拉格朗日反演公式：
即可求解\([x^s]F(x)\)，在\(\pmod{x^s}\)意义下多项式求逆，多项式快速幂即可.
[Lougu 4389] 付公主的背包
题意
有\(n\)种商品，每种商品体积为\(v_i\)，给定\(m\)，对于\(s∈[1,m]\)，请你回答用这些商品恰好装\(s\)体积的完全背包的方案数，\(n,m\leq 10^5\).
题解
设\(f_{i,j}\)表示考虑了前\(i\)个物品，装满体积为\(j\)的背包的方案数，则状态转移方程如下：
设\(F_i(x)\)为\(\{f_{i,n}\}\)的生成函数，则状态转移方程在生成函数上表达如下：
所以最终答案的生成函数\(F(x)=F_n(x)=\prod_{i=1}^n\frac{1}{1-x^{v_i}}\)，使用幂指恒等式转化：
最后一步用到了\(\ln \frac{1}{1-x}\)的展开式，这个结论在生成函数那一篇博客里推导过. 那么转化到这个形式就可以预处理\(\pmod{x^m}\)意义下括号内生成函数的值，然后多项式\(\exp\)求一下就好了.
[Luogu 4841] 城市规划
题意
求出\(n\)个点的简单有标号无向连通图数目，\(n\leq 10^5\).
题解
根据生成函数篇的讨论，设\(\hat G(x)\)为简单无向图数量的生成函数，则无向联通图的数量生成函数\(\hat F(x)=\ln \hat G(x)\)，那么预处理\(\hat G(x)\)，在\(\pmod{x^{n}}\)意义下求\(\ln\)即可.
[HAOI 2018] 染色
题意
给定\(n\)个格子，每个格子可以染\(m\)种不同的颜色，令\(w_k\)表示恰好出现\(s\)次的颜色有\(k\)种时对应方案的权值，求所有方案权值和，答案对\(98244353\)取模.
\(n\leq 10^7,m\leq 10^5,s\leq 150\).
题解
令\(f_k\)代表恰好出现\(s\)次的颜色有\(k\)种的方案数，显然答案为\(\sum_{i=1}^m f_iw_i\).
考虑对\(k\)进行容斥，设\(g_k\)代表恰好出现\(s\)次的颜色至少有\(k\)种的方案数，显然有：
组合意义就是强选\(k\)中颜色恰好出现了\(s\)次，然后剩下随便填，并乘上多重集排列数.
根据二项式反演公式（组合数学总结第一篇当中有介绍），有：
对于这种卷积，可以用类似于 \([\mathrm{ZJOI}2014]\)力 那题的方法求解，那么原问题解决.
这道题启示我们熟悉的二项式反演的结论还可以用\(\mathrm{NTT}\)优化.
[Luogu 4705] 玩游戏
题意
给定\(\{a_i\}\)，\(\{b_i\}\)，求\(\frac{1}{nm}\sum_{i=1}^n\sum_{j=1}^m(a_i+b_j)^k\).
题解
\(\frac{1}{nm}\)先不考虑，将后式用二项式定理展开：
不妨设
显然所求式的生成函数\(\hat F(x)=\hat A(x)\times \hat B(x)\)，那么问题转化为求解\(A(x), B(x)\).
根据生成函数公式，显然有：
这样直接算肯定是不行的，那么我们的思路就是想方设法把分式给消掉，便于处理出\(A(x)\). 于是可以想到\(\ln'(x)=\frac{1}{x}\)，事实上，用到了这样一个公式来转化：
也就是说
将之前的式子向上式转化：
那么显然有：
那么一切都好办了起来，对于\(\prod (1-a_ix)\)，我们可以使用分治乘法：
设多项式\(A_i=1-a_ix\)，令\(T_{l,r}=\prod_{i=l}^rA_i\)，那么\(T_{l,r}=T_{l,\lfloor\frac{l+r}{2}\rfloor}\times T_{\lfloor\frac{l+r}{2}\rfloor,r}\).
时间复杂度为\(T(n)=2T(\frac{n}{2})\times \mathcal{O}(n\log_2 n)=\mathcal{O}(n\log^2 n)\).
剩下的只要写写多项式求\(\ln\)和求导即可.
[LOJ 6183] 看无可看
题意
已知\(f(0),f(1)\)为常数，令\(f(i)=3f(i-1)+2f(i-2)(i\geq 2)\)，给定序列\(a_{1\sim n}\)，求：
\(n\leq 10^5,\forall i\in[1,n]\ a_i\leq 10^9\).
题解
令\(f_n=f(n)\)，根据常系数线性递推数列解的规律，可以列出特征方程\(x^2=3x+2\)，解得\(x=3/-1\).
那么可以设\(f_n=\alpha\times 3^n+\beta\times (-1)^n\)，带入\(f_0=f(0),f_1=f(1)\)可得\(\alpha=\frac{1}{4}(f_0+f_1),\beta=\frac{1}{4}(3f_0-f_1)\).
于是原题就转化成了求这个式子：
那么我们要求形如这样的式子：
解决方法是构造多项式如下
显然有\(\alpha\sum_{S\in [n],|S|=k}\left( \prod_{x\in S}3^{a_x} \right)=\alpha[x^k]F(x)\)，鉴于模数较小，直接用\(\mathrm{FFT}\)实现分治乘法即可.
[CF553E] Kyoya and Train
题意
给定一张 \(n\) 个点 \(m\) 条边的无重边无自环的有向图，你要从 \(1\) 号点到 \(n\) 号点去.
- 如果你在 \(t\) 时刻之后到达 \(n\) 号点，你要交 \(x\) 元的罚款.
- 每条边从 \(a_i\) 到 \(b_i\)，走过它需要花费 \(c_i\) 元，多次走过同一条边需要多次花费.
- 走过每条边所需的时间是随机的，对于 \(k \in [1,t]\)，\(\frac{p_{i,k}}{10^5}\) 表示走过第 \(i\) 条边需要时间 \(k\) 的概率. 因此如果多次走过同一条边，所需的时间也可能不同.
- 你希望花费尽可能少的钱（花费与罚款之和）到达 \(n\) 号点，因此每到达一个点，你可能会更改原有的计划.
- 求在最优决策下，你期望花费的钱数.
\(n \le 50\)，\(m \le 100\)，\(t \le 2 \times 10^4\)，\(x,c_i \le 10^6\)，\(\sum_{k=1}^t p_{i,k} = 10^5\)，答案精度误差 \(\le 10^{-6}\).
题解
首先应该要想到一个朴素\(\mathrm{dp}\)：\(f_{i,j}\)表示在\(n\)时刻\(j\)到达\(i\)号点前提下，前往终点的最小期望花费，则\(f_{n,j}=\begin{cases}0,j\leq t\\ x,j> t\end{cases}\)，剩下的点可以列出状态转移方程：
用于每一次转移\(j\)单增，所以这个\(\mathrm{dp}\)可以看做在以时间为层数的分层图上\(\mathrm{dp}\)，这个分层图显然是\(\mathrm{DAG}\)，所以这样做是可以的，时间复杂度\(\mathcal{O}((n+m)\times t^2)\).
这样做显然还会超时，考虑到\(f\)以时间为阶段，对第二维有依赖，可以尝试用时间分治来优化. 定义\(\mathrm{d}(l,r)\)为计算时间在\([l,r]\)范围内所有\(f\)值的函数，则我们可以先计算\([\mathrm{mid}+1,r]\)的\(f\)值，然后执行转移，再计算\([l,\mathrm{mid}]\)的\(f\)值，其中\(\mathrm{mid}=\left\lfloor\frac{l+r}{2}\right\rfloor\).
由于状态转移方程类似于卷积，所以我们不妨考虑每一条边的转移，然后用\(\mathrm{FFT}\)计算期望. 对于时间点\(x\in [l,\mathrm{mid}]\)，图上的点\(k\)，分治区间内其他点对其贡献为：
考虑到为了方便推卷积，求和式的下界从\(x+1\)开始，但是一部分\(f\)仍然未被计算，所以当\(0\)处理不影响答案，这点与分治乘法中的思想相同. 但是不同点在于这题里我们要把卷积求出的\(\mathrm{dp}\)值存在一个数组\(g\)里面，不能直接更新到\(f\)上，否则有一部分\(f\)的值就不为零了.
那么不妨设\(g_{i,j}\)表示在时间点\(j\)，走第\(i\)条边到终点的最小期望花费（当然设点为状态也可以），并在分治边界用\(g\)更新\(f\)，那么转移就可以改写成：
\(\mathrm{end}_i\)表示边\(i\)的终点，相当于原式枚举边的终点\(y\). 那么如果可以在\(\mathcal{O}(t\log_2t)\)的时间内计算转移，总的时间复杂度就优化到了\(\mathcal{O}((n+m)t\log^2 t)\)，接下来只要想办法把这个转移式推成卷积即可，这样就是常规套路了.
我们枚举每一条边计算，所以第一维暂时省去，可以推导如下：
令\(f'_{i}=f_{i+x+1}\)，\(p'_{i}=p'_{i+1}\)，其中\(i\in[0,r-x-1]\)，那么：
这样就可以直接卷积计算了，虽然这题本质上是结合了时间分治算法优化\(\mathrm{dp}\)，但是推导卷积的思路和分治\(\mathrm{FFT}\)是一样的，所以也可以说是一道分治\(\mathrm{FFT}\)的题.
[LOJ 575] 不等关系
题意
给定一个符号序列\(\mathrm{op}_{1\sim n}\)，求有多少个\(1\sim n+1\)的排列\(\{a\}\)满足：
其中\(\forall i\in[1,n],\mathrm{op}_i\in \{<,>\}\)，\(n\leq 10^5\).
题解
考虑到一串只有\(<\)的符号序列的方案数是唯一的，也就是递增序列，那么我们可以针对\(>\)计数，使用容斥原理.
设\(f_i\)表示前缀\(i\)合法的方案数，那么我们不妨考虑离\(i\)最近的一个大于号，设为\(\mathrm{op}_j\)，那么理应有\(f_i=f_j\binom{i}{j}\)，也就是\(i\)个数里选\(j\)个放在前\(j\)个位置，方案数为\(f_j\)，区间\([j+1,i]\)都是小于号，也就是单增序列，所以方案数唯一. 但是这样可能会算重，\(\mathrm{op}_j\)表示\(a_j>a_{j+1}\)，若随意放置的时候\(a_j<a_{j+1}\)，就不合法了. 一旦\(a_j<a_{j+1}\)，那么也就是说最后那一段单增区间和倒数第二段单增区间连在了一起，再考虑倒数第二个大于号的位置，减掉连在一起的贡献，发现会减多，和倒数第三段又连在了一起，以此类推，可以容斥：
根据常规套路，我们把组合数拆掉：
令\(p_j=[(\mathrm{op}_j\ is\ >)\ \mathrm{or}\ (j=0)]\)，整理：
令\(F_{i}=\frac{p_jf_j(-1)^{-c_j}}{j!}\)，整理：
这样就可以分治\(\mathrm{FFT}\)求解了.
[HDU 5279] YJC plays Minecraft
题意
给定\(n\)个完全图，第\(i\)个完全图的大小为\(a_i\)，第\(i\)个完全图的\(a_i\)号点与第\(i+1\)个完全图的\(1\)号点有一条边相连，第\(n\)个完全图的\(a_n\)号点也与第\(1\)个完全图的\(1\)号点相连. 求有多少种删边的方案，使得该图删去若干条边以后是一个没有环的无向图. 方案数对\(998244353\)取模.
\(n\leq 10^5,\forall i\in[1,n]\ a_i\leq 10^5\).
题解
首先对每一个子图考虑，由于是完全图，所以删边的方案也可以看做是选择一些边的方案. 将完全图保留一些边形成森林，那么每一个子图内也就没有环了. 假设\(n\)个点带标号构成森林的方案数为\(f_n\)，那么总的方案数应该可以想到是\(2^n\prod_{i=1}^nf_{a_i}\)，\(2^n\)表示每两个完全图之间的边都可以选或不选.
但是考虑到原图中没两个子图之间的边也可能最终成环，所以准确的答案是\(2^n\prod_{i=1}^nf_{a_i}-\prod_{i=1}^ng_{a_i}\)，\(g_{n}\)表示\(n\)个点带标号森林且\(1\)号点和\(n\)号点联通的方案数.
根据生成函数的知识我们首先应该想到求\(n\)个点带标号树的方案数，根据\(\mathrm{Prufer}\)序列和无根树一一对应的关系，有\(h_n=n^{n-2}\)，\(h_n\)表示\(n\)个点带标号无根树的数量. 那么可以设\(\{h_n\}\)的指数型生成函数为\(\hat H(x)\)，\(\{f_n\}\)的指数型生成函数为\(\hat F(x)\)，显然有\(\hat F(x)=e^{\hat H(x)}\)，那么这两个都是可以求解的.
现在的问题就转化为求解\(g_n\)，我们考虑补集转化，令\(g_n=f_n-p_n\)，则\(p_n\)表示表示\(n\)个点带标号森林且\(1\)号点和\(n\)号点不联通的方案数，然后列出如下状态转移方程：
其中\(i\)枚举\(1\)号点所在联通块的大小，由于不与\(n\)号点联通，所以大小至多为\(n-1\)，选择剩下\(i-1\)个点的方案数为\(\binom{n-2}{i-1}\)，他们共同构成一棵树，方案数为\(h_i\)，剩下的点形成森林，方案数为\(f_{n-i}\).
整理状态转移方程为\(\mathrm{EGF}\)的形式：
应该不难发现\(\frac{p_n}{(n-2)!}=[x^{n-2}]\hat P''(x)\)，\(\frac{h_i}{(i-1)!}=[x^{i-1}]\hat H'(x)\)，\(\frac{f_{n-i}}{(n-i-1)!}=[x^{n-i}]\hat F'(x)\)，所以调整一下求和号上下界可以得到：
于是可以得出
那么一切都容易解决了.
[SDOI 2017] 遗忘的集合
题意
给定一个集合\(\mathrm{S}\)，其中每个元素大小都不超过\(n\)，令\(f_i\)表示取集合\(\mathrm{S}\)中若干元素（同一元素可以取多次，也可以一次都不取），使其和为\(i\)的方案数对\(p\)取模的值.
现在给定序列\(f_{1\sim n}\)，求一个字典序最小的集合\(\mathrm{S}\)符合题意，保证有解.
\(n< 2^{18},p<2^{30},\forall i\in[1,n]\ f_i<p\).
题解
很有意思的一道题. 首先应该看出来序列\(f\)的构造过程类似于完全背包方案数，集合\(\mathrm{S}\)有元素\(i\)则说明有一个体积为\(i\)的物品，那么我们可以令\(s_i=0/1\)表示集合\(\mathrm{S}\)是否有大小为\(i\)的元素. 设\(F(x)\)表示序列\(\{f_n\}\)的生成函数，根据付公主的背包一题中推出来的完全背包生成函数解的表示，有：
同样的套路，我们对两边同时取\(\ln\)，则有：
我们早就推过序列\(\{0,1,\frac{1}{2},\frac{1}{3},\cdots\}\)的生成函数是\(\ln \frac{1}{1-x}\)，所以可以直接带进去：
上面的式子只要交换求和符号就可以得到，那么接下来的思路就应该很显然了. 令\(g_n=n[x^n]\ln F(x)\)，\(a_n=n\times s_n\)，则有：
调和级数地枚举算贡献即可求出每一个\(a_i\)，进而求出\(s_i\)，由构造过程可知解是唯一的.
求\(g\)需要用到多项式\(\ln\)，鉴于模数不定，还需要\(\mathrm{MTT}\).
[UOJ 50] 链式反应
题意
有一棵\(n\)个点的带标号的树，其任意一个节点都满足一下两个条件之一：
- 该节点是叶节点，没有儿子
- 该节点有且仅有\(c+2\)个儿子，其中\(c\)个儿子必须是叶节点，以黑色边相连，且满足\(c\in A\)，剩下两个儿子以白色边相连，不作限制
求不同二叉树的方案数，其中\(A\)是给定数集，\(|A|\leq n\leq 2\times 10^5\).
题解
真是一道不可多得的分治\(\mathrm{FFT}\)好题.
首先设答案的指数型生成函数是\(\hat F(x)\)，集合\(A\)的指数型生成函数是\(\hat A(x)\)，那么应该容易想到：
首先子树大小的分配类似于背包，所以\(\hat F^2(x)\hat A(x)\)不难解释，因为根节点也占总大小中的一个，所以要\(\int\)向右平移一下，又因为链接两个不限制儿子的边不作区分，所以方案数乘\(\frac{1}{2}\). 至于加了一个\(x\)则是因为只有一个点的树（根节点就是叶子节点）也是可行的方案.
设\(f_n\)表示\(n\)的节点的树的方案数，写出\(\mathrm{dp}\)方程也许更好理解上面的生成函数表达式：
那么接下来就有两种思路了，其一是使用生成函数方法，根据牛顿迭代的思想解常微分方程，然后递推，据说时间复杂度是\(\mathcal{O}(n\log_2 n)\)，难度比较大还要写多项式全家桶.
还有一种思路是卷积优化\(\mathrm{dp}\)，接下来重点讲.
首先根据\(\mathrm{EGF}\)的套路整理式子，令\(F_i=\frac{f_i}{i!}\)，\(A_{i}=\frac{[i\in A]}{i!}\)，那么有：
其中边界是\(F_1=f_1=1\). 考虑化简：
令\(F_0=0\)，\(G_p=\sum_{i=0}^pF_iF_{p-i}\)，则有：
我们可以对\(G\)做分治\(\mathrm{FFT}\)，具体如下：
假设分治区间为\([l,r)\)，右部点\(x\in[\mathrm{mid},r)\)的贡献为：
其中\(G'_i=G_{i+l}\)，只要对\(A[0,r-l-1)\)和\(G'[0,\mathrm{mid}-l)=G[l,\mathrm{mid})\)做卷积，取第\(x-l-1\)项系数贡献到\(F_x\)上就可以了.
但是现在最大的问题就是我们如何在递归完分治区间\([l,\mathrm{mid})\)之后计算出\(G[l,\mathrm{mid})\)的值. 仔细一想的话并没有那么容易，因为\(G\)是受到\(F\)约束的.
这里我们不妨宕开一笔讲一下另一个\(\mathrm{cdq}\)分治优化卷积的问题：
\(\mathrm{Preblem\ ex.}\) 已知\(f_1\)，求\(f_{2\sim n}\)，其中\(f_i=\sum_{j=1}^{i-1}f_jf_{i-j}(i\geq 2)\)
我们不妨先按照之前一般分治\(\mathrm{FFT}\)的思路考虑，把求和式里的第一个\(f\)当成\(f\)，把第二个\(f\)当成常数序列\(g\)，然后考虑贡献：
对于右部点\(x\in[\mathrm{mid},r)\)，其贡献为：
其中\(f'_i=f_{i+l},f''_{i}=f_{i+1}\)，那么我们就要对\(f'[0,\mathrm{mid}-l)=f[l,\mathrm{mid})\)和\(f''[0,r-l-1)=f[1,r-l)\)这两段做卷积.
显然\(f[l,\mathrm{mid})\)我们已经递归求出来了，但是很可惜我们无法保证\(r-l\leq \mathrm{mid}\)，也就是说\(f[1,r-l)\)可能还有一部分值不知道.
此时就没有好的直接处理方案，我们只好考察什么时候会出现\(r-l>\mathrm{mid}\).
\(\mathrm{Cases1:}\)
注意到第一个左端点非\(1\)的分治区间\([\mathrm{mid},r)\)满足\(r-l\leq l\leq \mathrm{mid}\)，并且继续往下递归时要么\(r\)减小，要么\(l\)增大，都有\(r-l\leq l\leq \mathrm{mid}\)成立，即满足\(r-l\leq \mathrm{mid}\)，直接正常地分治\(\mathrm{FFT}\)即可.
\(\mathrm{Case2:}\)
剩下的也就是形容\([1,r)\)，此时不一定满足\(r-1\leq \mathrm{mid}\).
回到上面的推导，我们此时要求\(f[1,\mathrm{mid})\)和\(f[1,r-1)\)的卷积，然后贡献到\(f[\mathrm{mid},r)\)上，我们可以把它看做两次贡献：把\(f[1,\mathrm{mid})\)和\(f[1,\mathrm{mid})\)做卷积的贡献，把\(f[1,\mathrm{mid})\)和\(f[\mathrm{mid},r-1)\)做卷积的贡献.
显然第一部分贡献我们可以直接卷积求出，第二部分贡献还不知道. 但是，区间\([\mathrm{mid},r-1)\)应该归类在\(\mathrm{Case1}\)里面，可以求解，那么我们不妨把这次贡献延迟到递归完右半区间在计算.
具体来说，对于处理形如\(\mathrm{Case1}\)的区间时，分治函数定义的任务就多了一个，弥补前缀区间漏算的贡献. 而前缀区间对于某个右半区间\([\mathrm{mid},r)\)的贡献，恰好可以用\(f[1,l)\times f[l,\mathrm{mid})\)的某一项表示，与\(\mathrm{Case1}\)本来就要算的贡献有包含关系，所以在处理\(\mathrm{Case1}\)时，只需把\(f[1,r-l)\)前\(l-1\)项系数乘\(2\)，再直接计算\(f[1,r-l)\times f[l,\mathrm{mid})\)，贡献到右半区间即可.
总的来说，我们通过讨论完善了分治函数：对于\(\mathrm{Case2}\)的区间，我们先计算可以计算的部分，对于\(\mathrm{Case1}\)的区间，我们计算了左半部分对右半部分的贡献，并且顺带弥补了之前前缀区间对其漏算的贡献，那么这样计算出的前缀贡献就是真实有效的了，根据\(\mathrm{cdq}\)分治算法的性质，最终求出的答案也包含了所有应该计算的贡献.
对于代码实现，只要卷积时对\(f[1,r-l)\)前\(l-1\)项系数乘\(2\)再卷积就好.
回到原来的问题，我们观察\(G\)数组的定义：
既然通过\(G\)可以轻松计算\(F\)左半区间对右半区间的贡献，那么我们在分治函数里，利用刚刚讨论的方法，顺带地用\(F\)求解右半区间的\(G\)即可.
此处必须有参考代码：
inline void cdqDivide(int l,int r,Poly &f,Poly &g,const Poly &A)
{
if ( l + 1 == r ) return void( l == 1 ? f[l] = 1 : f[l] = mul( f[l] , mul( inv[2] , inv[l] ) ) );
int mid = ( l + r + 1 ) >> 1; Poly h;
cdqDivide( l , mid , f , g , A );
h = Multiply( Poly( A.begin() , A.begin() + r - l - 1 ) ,
Poly( g.begin() + l , g.begin() + mid ) , r - l - 1 );
for (int i = mid; i < r; i++) Inc( f[i] , h[i-l-1] ); h.resize(r-l);
for (int i = 0; i < r - l - 1; i++)
if ( i + 1 < mid ) h[i] = mul( f[i+1] , i + 1 < l ? 2 : 1 );
else { h.resize(i); break; }
h = Multiply( h , Poly( f.begin() + l , f.begin() + mid ) , r - l - 1 );
for (int i = mid; i < r; i++) Inc( g[i] , h[i-l-1] );
return cdqDivide( mid , r , f , g , A );
}
其他习题
- [Luogu 5900] 无标号无根树计数
- [HNOI2019] 白兔之舞
- [LOJ 6569] 仙人掌计数
- [Lougu 5434] 有标号荒漠计数
- [LOJ 6363] 地底蔷薇
- [CTS 2018] 青蕈领主
- [Luogu 5644] 猎人杀
- [UOJ 269] 如何优雅的求和
- [CTS 2019] 珍珠
以上习题不按照难度排序.
Epilogue
生成函数的题还有很多，但是入门的知识到这里就差不多结束了. 本博客到此结束，以后可能会补充习题的题解.
