前言
本篇文章为基础的多项式计算的集合。
文章中涉及到的多项式计算，通常会在生成函数中使用到。
另外，模板题可以去洛谷上面找，一搜就有。
代码有一点需要注意：如果你看到上下文中代码的调用方式不一样，不要惊慌。
这是因为，我的写法是，给每个运算开辟一个 namespace 。如果是牛顿迭代法计算的话，我会再重载一个接口函数。
牛顿迭代法的接口长成这个样子：
void 运算名( int *ret, int *A, const int n ) //分别表示“返回值”，“输入多项式”和“多项式长度”
{
清空运算辅助数组;
复制一遍 A ;
运算名( n ); //这个就是我下文将给出的函数，只有一个参数，表示“多项式长度”
将答案复制到 ret 里面;
}
多项式乘法
这是所有多项式计算的基础，由于真的很基础，所以直接贴代码了。
喜闻乐见的FFT：
void FFT( comp *coe, int len, int type )
{
int lg2 = log2( len );
for( int i = 0, rev ; i < len ; i ++ )
{
rev = 0;
for( int j = 0 ; j < lg2 ; j ++ ) rev |= ( ( i >> j ) & 1 ) << ( lg2 - j - 1 );
if( rev < i ) swapp( coe[rev], coe[i] );
}
comp wp, w, wo, we;
for( int s = 2, p ; s <= len ; s <<= 1 )
{
p = s >> 1, wp = comp( cos( type * PI / p ), sin( type * PI / p ) );
for( int i = 0 ; i < len ; i += s )
{
w = comp( 1, 0 );
for( int j = 0 ; j < p ; j ++, w *= wp )
{
we = coe[i + j], wo = coe[i + j + p];
coe[i + j] = we + wo * w;
coe[i + j + p] = we - wo * w;
}
}
}
if( ~ type ) return ;
for( int i = 0 ; i <= len ; i ++ ) coe[i] /= len;
}
喜闻乐见的NTT，使用 998244353 作为模数，顺手预处理了单位根：
void init( const int len )
{
int lg2 = log2( len );
for( int i = 0 ; i < len ; i ++ )
for( int j = 0 ; j < lg2 ; j ++ )
rev[i] |= ( i >> j & 1 ) << ( lg2 - j - 1 );
wp[0] = 1, wp[1] = qkpow( g, phi / len );
wpinv[0] = 1, wpinv[1] = qkpow( g, phi - phi / len );
for( int i = 2 ; i < len ; i ++ )
wp[i] = 1ll * wp[i - 1] * wp[1] % mod,
wpinv[i] = 1ll * wpinv[i - 1] * wpinv[1] % mod;
}
void NTT( int *coe, const int len, const int t )
{
#define p ( s >> 1 )
for( int i = 0 ; i < len ; i ++ )
if( rev[i] < i )
swapp( coe[i], coe[rev[i]] );
int w, wo, we;
for( int s = 2 ; s <= len ; s <<= 1 )
for( int i = 0 ; i < len ; i += s )
for( int j = 0 ; j < s >> 1 ; j ++ )
{
w = t > 0 ? wp[len / s * j] : wpinv[len / s * j];
we = coe[i + j], wo = 1ll * coe[i + j + p] * w % mod;
coe[i + j] = ( we + wo ) % mod, coe[i + j + p] = ( we - wo + mod ) % mod;
}
#undef p
if( t > 0 ) return; int inver = inv( len );
for( int i = 0 ; i < len ; i ++ ) coe[i] = 1ll * coe[i] * inver % mod;
}
多项式牛顿迭代法
其实这个不应该叫做“牛顿迭代”，而应该叫做套路倍增。
一般情况下的牛顿迭代法可以求出(f(x)=0)的根，多项式牛顿迭代就可以求出(G(F(x))=0)的根。
方法是，使用倍增，考虑我们已经求出了：
[G(F_0(x))equiv 0pmod {x^{lceilfrac n 2
ceil}}
]
现在我们想要得到：
[G(F(x))equiv 0pmod {x^n}
]
那么，我们用泰勒公式把(G)展开：
[G(F(x))=sum_{k=0}^{+infty}frac{G^{(k)}(F_0(x))}{k!}(F(x)-F_0(x))^k
]
注意到，(F(x))在低次下也是方程的根：
[G(F(x))equiv 0pmod {x^{lceilfrac n 2
ceil}}
]
那么做差就可以得到：
[G(F(x))-G(F_0(x))equiv 0pmod {x^{lceilfrac n 2
ceil}}
]
(G)不应该是一个常函数，那么就应该有：
[F(x)-F_0(x)equiv 0pmod {x^{lceilfrac n 2
ceil}}
]
考虑左右平方。由于(F(x)-F_0(x))的后(lceilfrac n 2 ceil)项都是 0 ，那么平方后的式子的后(n)项，卷积得到结果也都应该是 0 （乘法时只要一个因数为 0 结果都是 0）。
于是就有：
[(F(x)-F_0(x))^2equiv 0pmod {x^n}
]
同样可知，更高次的幂在模意义下也应该是 0 。
然后泰勒公式就可以化简得到：
[G(F(x))equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))equiv 0pmod {x^n}
]
简单地整理一下得到：
[F(x)equiv F_0(x)-frac{G(F_0(x))}{G'(F_0(x))}pmod {x^n}
]
奇妙吧，我也觉得。
多项式乘法逆
问题是，给定一个(A(x))，求一个(A^{-1}(x))，使得(A(x)A^{-1}(x)equiv 1pmod {x^n})。
尝试倍增，假设我们已经有了：
[A(x)F_0(x)equiv 1pmod {x^{lceilfrac n 2
ceil}}
]
然后要求：
[A(x)F(x)equiv 1pmod {x^n}
]
然后通过牛顿迭代法部分的分析，我们知道应该有：
[(F(x)-F_0(x))^2equiv 0pmod {x^n}
]
展开，并且两边同时乘上(A(x))有：
[F(x)-2F_0(x)+A(x)F_0^2(x)equiv 0pmod {x^n}
]
简单变化得到：
[F(x)equiv F_0(x)(2-A(x)F_0(x))pmod {x^n}
]
时间复杂度应该是(T(n)=T(frac n 2)+O(nlog_2n)=O(nlog_2n))。
但是常数比较大，下面粘一份代码：
int P[MAXN], F0[MAXN], F[MAXN], B[MAXN];
void PolyInv( const int n )
{
if( n == 1 ) { F[0] = inv( P[0] ); return; }
int p = n + 1 >> 1, len = 1;
PolyInv( p );
while( len < n << 1 ) len <<= 1; init( len );
for( int i = 0 ; i < n ; i ++ ) B[i] = P[i];
for( int i = 0 ; i < len ; i ++ ) F0[i] = F[i], F[i] = 0;
NTT( B, len, 1 ), NTT( F0, len, 1 );
for( int i = 0 ; i < len ; i ++ ) F[i] = 1ll * F0[i] * ( 2 - 1ll * B[i] * F0[i] % mod + mod ) % mod;
NTT( F, len, -1 );
for( int i = n ; i < len ; i ++ ) F[i] = 0;
}
多项式带余除法
给定(n)次多项式(A(x))和(m)次多项式(B(x))，(m<n)，要求(n-m)次多项式(C(x))和低于(m)次的多项式(R(x))，使得以下等式成立：
[A(x)equiv B(x)C(x)+R(x)pmod {x^n}
]
这还不简单，大除法谁不会？
不难感受到，最恶心的部分是(R(x))，因而我们得想办法把它踢走。
观察一下，我们发现，(R(x))至少前(n-m+1)次的系数都是 0 ，如果我们把系数在(pmod {x^n})意义下翻转，(R(x))就自然而然地消失了。
翻转系数是比较简单的操作，大概就是：
[x^nA(frac 1 x)equiv x^mB(frac 1 x)cdot x^{n-m}C(frac 1 x)+x^nR(frac 1 x)pmod {x^n}
]
再见，(R(x))！
[x^nA(frac 1 x)equiv x^mB(frac 1 x)cdot x^{n-m}C(frac 1 x)pmod {x^n}
]
然后用 多项式求逆 + 多项式乘法 就可以得到(x^{n-m}C(frac 1 x))，翻转之后再做一下乘法就可以得到(R(x))。
时间复杂度是喜闻乐见的(O(nlog_2n))。
int G[MAXN], H[MAXN], A[MAXN], B[MAXN];
int N, M;
void PolyMul( int *ret, int *a, const int La, int *b, const int Lb )
{
int len = 1;
while( len < La + Lb ) len <<= 1;
init( len );
for( int i = 0 ; i < len ; i ++ ) ret[i] = A[i] = B[i] = 0;
for( int i = 0 ; i < La ; i ++ ) A[i] = a[i];
for( int i = 0 ; i < Lb ; i ++ ) B[i] = b[i];
NTT( A, len, 1 ), NTT( B, len, 1 );
for( int i = 0 ; i < len ; i ++ ) ret[i] = 1ll * A[i] * B[i] % mod;
NTT( ret, len, -1 );
for( int i = La + Lb ; i < len ; i ++ ) ret[i] = 0;
}
int main()
{
...
std :: reverse( G, G + N );
std :: reverse( H, H + M );
PolyInv :: PolyInv( HInv, H, N - M + 1 );
PolyMul :: PolyMul( C, G, N, HInv, N - M + 1 );
for( int i = N - M + 1 ; i < 2 * N - M + 1 ; i ++ ) C[i] = 0;
std :: reverse( C, C + N - M + 1 );
std :: reverse( H, H + M );
std :: reverse( G, G + N );
PolyMul :: PolyMul( R, H, M, C, N - M + 1 );
for( int i = 0 ; i < N ; i ++ ) R[i] = ( G[i] - R[i] + mod ) % mod;
...
}
多项式 ln
给定多项式(A(x))，求(B(x)equiv ln A(x)pmod {x^n})。
没想到多项式还有初等函数吧，我也没有想到。
这个还比较简单，众所周知：
[ln x=int frac{mathrm d x}{x}
]
因此，对原式子两边求导：
[B'(x)equiv frac{A'(x)}{A(x)}pmod {x^n}
]
其中多项式的积分和导数可以直接算，再套上一个求逆，我们就得到了(ln A(x))。时间复杂度(O(nlog_2n))。
设多项式(A(x)=sum_{i=0} a_ix^i)，它的导数是：
[A'(x)=sum_{i=1} ia_ix^{i-1}
]
积分是：
[int A(x) mathrm d x=sum_{i=0}frac {a_i} {i+1} x^{i+1} + C
]
模板题里面，保证了(a_0=1)，因此积分中(C=0)。不然我也不知道怎么计算
int inv[MAXN], B[MAXN], C[MAXN], D[MAXN];
void PolyDer( int *ret, int *A, const int len )
{
for( int i = 0 ; i < len - 1 ; i ++ )
ret[i] = 1ll * A[i + 1] * ( i + 1 ) % mod;
ret[len - 1] = 0;
}
void PolyInt( int *ret, int *A, const int len )
{
inv[1] = 1;
for( int i = 2 ; i <= len ; i ++ )
inv[i] = 1ll * inv[mod % i] * ( mod - mod / i ) % mod;
for( int i = 1 ; i <= len ; i ++ )
ret[i] = 1ll * A[i - 1] * inv[i] % mod;
ret[0] = 0;
}
void PolyLn( int *ret, int *A, const int n )
{
PolyDer :: PolyDer( B, A, n );
PolyInv :: PolyInv( C, A, n );
int len = 1;
while( len <= n << 1 ) len <<= 1;
init( len );
NTT( B, len, 1 ), NTT( C, len, 1 );
for( int i = 0 ; i < len ; i ++ ) D[i] = 1ll * B[i] * C[i] % mod;
NTT( D, len, -1 );
for( int i = n ; i < len ; i ++ ) D[i] = 0;
PolyInt :: PolyInt( ret, D, n - 1 );
for( int i = n ; i < len ; i ++ ) ret[i] = 0;
}
多项式 exp
给定多项式(A(x))，求(B(x)equiv e^{A(x)}pmod {x^n})。
两边求导：
[ln B(x)equiv A(x)pmod {x^n}
]
然后就可以尝试解方程：
[G(F(x))equiv ln F(x)-A(x)equiv 0pmod {x^n}
]
套用牛顿迭代法：
[egin{aligned}
&F(x)equiv F_0(x)-frac{ln F_0(x)-A(x)}{F_0^{-1}(x)}&pmod {x^n}\
Rightarrow & F(x)equiv F_0(x)(1-ln F_0(x)+A(x))&pmod {x^n}
end{aligned}
]
然后倍增计算，时间复杂度(T(n)=T(frac n 2)+O(nlog_2n)=O(nlog_2n))。
常数嘛......你懂的
int P[MAXN], F0[MAXN], F[MAXN], B[MAXN];
void PolyExp( const int n )
{
if( n == 1 ) { F[0] = 1; return ; }
int p = n + 1 >> 1, len = 1;
PolyExp( p );
while( len < ( n << 1 ) ) len <<= 1; init( len );
for( int i = 0 ; i < len ; i ++ ) F0[i] = F[i], F[i] = 0;
PolyLn :: PolyLn( B, F0, n );
for( int i = 0 ; i < n ; i ++ ) B[i] = ( P[i] - B[i] + mod ) % mod;
B[0] = ( B[0] + 1 ) % mod, NTT( B, len, 1 ), NTT( F0, len, 1 );
for( int i = 0 ; i < len ; i ++ ) F[i] = 1ll * F0[i] * B[i] % mod;
NTT( F, len, -1 );
for( int i = n ; i < len ; i ++ ) F[i] = 0;
}
多项式快速幂
给定(A(x))和正整数(k)，求(B(x)equiv (A(x))^kpmod {x^n})。
两边(ln)一下就有：
[ln B(x)equiv kA(x)pmod {x^n}
]
套板子计算一下就好，时间虽然是(O(nlog_2n))。
但常数就更大了
int A[MAXN], LnA[MAXN];
void PolyQkpow( int *ret, int *a, const int l1, const int indx )
{
for( int i = 0 ; i < l1 ; i ++ ) A[i] = a[i];
PolyLn :: PolyLn( LnA, A, l1 );
for( int i = 0 ; i < l1 ; i ++ ) LnA[i] = 1ll * LnA[i] * indx % mod;
PolyExp :: PolyExp( ret, LnA, l1 );
}
多项式开根
给定(A(x))，求(B(x)equiv sqrt{A(x)}pmod {x^n})。
问题等价于求下面函数的一个根：
[G(F(x))=(F(x))^2-A(x)
]
套用牛顿迭代公式：
[egin{aligned}
F(x)&equiv F_0(x)-frac{(F_0(x))^2-A(x)}{2F_0(x)}&pmod {x^n}\
Rightarrow F(x)&equiv frac{F_0(x)+frac{A(x)}{F_0(x)}}{2}&pmod {x^n}
end{aligned}
]
时间复杂度(O(nlog_2n))。
需要注意的是，如果(A)的常数项不是 1 ，那么我们就需要用二次剩余计算出(B)的常数项。代码中给出的是(A)的常数项是 1 的情况。
代码：
int P[MAXN], F0[MAXN], F[MAXN], B[MAXN], Finv[MAXN];
void PolySqrt( const int n )
{
if( n == 1 ) { F[0] = 1; return ; }
int half = n + 1 >> 1, len = 1;
PolySqrt( half );
while( len < ( n << 1 ) ) len <<= 1;
for( int i = 0 ; i < n ; i ++ ) B[i] = P[i];
for( int i = 0 ; i < len ; i ++ ) F0[i] = F[i], F[i] = 0;
PolyInv :: PolyInv( Finv, F0, n );
NTT( Finv, len, 1 ), NTT( B, len, 1 );
for( int i = 0 ; i < len ; i ++ ) F[i] = 1ll * Finv[i] * B[i] % mod;
NTT( F, len, -1 );
for( int i = n ; i < len ; i ++ ) F[i] = 0;
for( int i = 0 ; i < n ; i ++ ) F[i] = 1ll * ( F[i] + F0[i] ) % mod * inv2 % mod;
}
更多
咕咕咕
看！鸽子！