多项式膜法（咕咕咕？）

前言

本篇文章为基础的多项式计算的集合。

文章中涉及到的多项式计算，通常会在生成函数中使用到。

另外，模板题可以去洛谷上面找，一搜就有。

代码有一点需要注意：如果你看到上下文中代码的调用方式不一样，不要惊慌。

这是因为，我的写法是，给每个运算开辟一个 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 0\pmod {x^{\lceil\frac n 2\rceil}} \]

现在我们想要得到：

\[G(F(x))\equiv 0\pmod {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 0\pmod {x^{\lceil\frac n 2\rceil}} \]

那么做差就可以得到：

\[G(F(x))-G(F_0(x))\equiv 0\pmod {x^{\lceil\frac n 2\rceil}} \]

\(G\)不应该是一个常函数，那么就应该有：

\[F(x)-F_0(x)\equiv 0\pmod {x^{\lceil\frac n 2\rceil}} \]

考虑左右平方。由于\(F(x)-F_0(x)\)的后\(\lceil\frac n 2\rceil\)项都是 0 ，那么平方后的式子的后\(n\)项，卷积得到结果也都应该是 0 （乘法时只要一个因数为 0 结果都是 0）。

于是就有：

\[(F(x)-F_0(x))^2\equiv 0\pmod {x^n} \]

同样可知，更高次的幂在模意义下也应该是 0 。

然后泰勒公式就可以化简得到：

\[G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\equiv 0\pmod {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 1\pmod {x^n}\)。

尝试倍增，假设我们已经有了：

\[A(x)F_0(x)\equiv 1\pmod {x^{\lceil\frac n 2\rceil}} \]

然后要求：

\[A(x)F(x)\equiv 1\pmod {x^n} \]

然后通过牛顿迭代法部分的分析，我们知道应该有：

\[(F(x)-F_0(x))^2\equiv 0\pmod {x^n} \]

展开，并且两边同时乘上\(A(x)\)有：

\[F(x)-2F_0(x)+A(x)F_0^2(x)\equiv 0\pmod {x^n} \]

简单变化得到：

\[F(x)\equiv F_0(x)(2-A(x)F_0(x))\pmod {x^n} \]

时间复杂度应该是\(T(n)=T(\frac n 2)+O(n\log_2n)=O(n\log_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(n\log_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(n\log_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 0\pmod {x^n} \]

套用牛顿迭代法：

\[\begin{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(n\log_2n)=O(n\log_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))^k\pmod {x^n}\)。

两边\(\ln\)一下就有：

\[\ln B(x)\equiv kA(x)\pmod {x^n} \]

套板子计算一下就好，时间虽然是\(O(n\log_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) \]

套用牛顿迭代公式：

\[\begin{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(n\log_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;
}

更多

咕咕咕

看！鸽子！