多项式（Ⅱ）：进阶工业

书接上回多项式（Ⅰ）：基础工业。这部分主要写一下进阶的一些模板。

这部分你可能需要一点点导数和积分的知识。

多项式牛顿迭代

参考资料：

牛顿迭代法，就是用来找函数零点的一个理论复杂度同二分但是实际上快而 nb 多的东西-nyy

其实这个东西在高数上册微分在求方程近似解上的应用里面出现过，其中的切线法，就是牛顿迭代。

下面我们讨论一下多项式的牛顿迭代。

给定多项式 \(G(x)\)，已知有 \(F(x)\) 满足：

\[G(F(x)) \equiv 0(\bmod \ x^n) \]

求出模 \(x^n\) 意义下的 \(F(x)\)。

其实它跟我们下面要讲的多项式求逆的倍增法本质上是一样的，都是倍增的思想。

假设我们已知 \(F_0(x) \equiv F(x)(\bmod \ x^{\lceil \frac{n}{2} \rceil})\)，我们要求 \(F_1(x) \equiv F(x)(\bmod \ x^n)\)

下面我们先给出牛顿迭代的公式：

\[G(F(x)) = F_0(x) - \frac{G(F_0(x))}{G'(F_0(x))} \]

接下来我们简单证明一下：

我们在 \(F_0(x)\) 处用泰勒公式展开 \(G(F(x))\)

\[G(F(x)) = \sum_{i=0}^\infty \frac{G^{(i)}(F_0(x))}{i!}(F(x)-F_0(x))^i \]

根据给出的条件，我们可以知道 \(\displaystyle F(x) - F_0(x) \equiv 0(\bmod \ x^{\lceil \frac{n}{2} \rceil})\) 是成立的，所以说 \(\displaystyle (F(x)-F_0(x))^k \equiv 0(\bmod x^n)\) 在 \(k \ge 2\) 的时候成立。

所以说，我们泰勒展开后从二次项开始都是 \(0\)，有意义的只有前两项，因此：

\[G(F_1(x)) = G(F_0(x))+ G'(F_0(x))(F_1(x)-F_0(x))(\bmod \ x^n) \]

又因为 \(G(F_1(x)) \equiv 0(\bmod \ x^n)\)，我们把方程整理一下，便有：

\[F_1(x) = F_0(x) - \frac{G(F_0(x))}{G'(F_0(x))} \]

自此，公式得证。

可以看出，从 \(\lceil\frac{n}{2}\rceil\) 到 \(n\)，是一个 \(\log(n)\) 的过程，所以它的复杂度与二分和倍增是一样的，接下来的一些题目中，我们会运用到它来解题。

多项式求逆（Polyinv）

给定一个多项式 \(F(x)\)，求出一个多项式 \(G(x)\)，满足 \(F(x) * G(x) \equiv 1(\bmod \ x^n)\)。

多项式乘法逆

我们先讨论比较简单的模数是 \(998244353\) 这样类型的。

首先有个结论：若一个多项式有乘法逆，当且仅当它的常数项有逆元。

方法1：分治 FFT

目前还没学，咕掉先。时间复杂度是 \(O(n\log^2 n)\) 的。

方法2：倍增法

参考资料：多项式求逆-litble

如果多项式 \(F(x)\) 只有一项，那么显然 \(G_0\) 这个常数项就是 \(F_0\) 的逆元。

若有 \(n\) 项，考虑递归求解，像类似于归纳法一样往后推。

假设我们已经知道

\[F(x)H(x) \equiv 1(\bmod \ x^{\lceil \frac{n}{2}\rceil}) \]

因为 \(F(x) * G(x) \equiv 1(\bmod \ x^n)\)，这就说明对于第 \(1 \sim n-1\) 次方项，两个多项式乘起来后系数为 \(0\)，这就说明

\[F(x)G(x) \equiv 1(\bmod \ x^{\lceil \frac{n}{2}\rceil}) \]

两个式子相减

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

即

\[G(x) - H(x) \equiv 0(\bmod \ x^{\lceil \frac{n}{2}\rceil}) \]

我们对上面这个式子进行平方。由于 \(G(x) - H(x)\) 在模 \(x^{\lceil \frac{n}{2}\rceil}\)
意义下为 \(0\)，说明从 \(0\sim \lceil \frac{n}{2}\rceil-1\) 次项都为 \(0\)。设平方后的这个多项式为 \(P\)，则

\[P_i = \sum_{j=0}^i(G(x)-H(x))_j (G(x)-H(x))_{i-j} \]

那么对于 \(i \leq n\) 的项，\(j\) 和 \(i-j\) 至少有一项的次数 \(< \lceil \frac{n}{2 }\rceil\)，那么就是 \(0\)，所以

\[G(x)^2 - 2G(x)H(x) + H(x)^2 \equiv 0(\bmod \ x^n) \]

乘上一个 \(F(x)\)，根据 \(F(x)G(x) \equiv 1(\bmod x^n)\)，再移一下项，就得到

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

转化到了这个地方便可以用 NTT 求解了。

复杂度用主定理分析一下：

\[T(n) = T(n/2) + O(n\log n) = O(n \log n) \]

此外还有一些常用的递归复杂度：

\[T(n) = T(n/2) + O(n) = O(n) \]

\[T(n) = 2T(n/2) + O(n) = O(n\log n) \]

\[T(n) = 2T(n/2) + O(n\log n) = O(n\log^2 n) \]

可以用主定理分析一下，感觉基本的这么多就可以了。

代码实现

有一些小细节需要注意：

代码中 (n+1)>>1 就是向上取整，至于为什么是向上取整，这样就能保证 \(\lceil \frac{n}{2} \rceil\) 的平方一定是能包含到 \(n\) 的。
为什么 lim < (n<<1) ，因为这样相乘以后就能保证两个 \(n\) 项的多项式相乘后的多项式能存的开。
为什么要把 \(n \sim \text{lim}-1\) 设成 \(0\)，因为模 \(x^n\) 意义下这些地方是同余于 \(0\) 的，为了避免干扰计算，就把系数也设成 \(0\) 了。后面的那个 \(b_i\) 设 \(0\) 也是同理的。

il void Polyinv(int n,ll a[],ll b[])
	{
		if(n == 1) { b[0] = inv(a[0]); return ; }
		Polyinv((n+1)>>1,a,b);
		static ll c[N];
		lim = 1;
		while(lim < (n<<1)) lim <<= 1;
		for(re int i=0;i<lim;i++) R[i] = (R[i>>1]>>1) + ((i&1) ? lim/2 : 0);
		for(re int i=0;i<n;i++) c[i] = a[i];
		for(re int i=n;i<lim;i++) c[i] = 0;
		NTT(c,lim,1) , NTT(b,lim,1);
		for(re int i=0;i<lim;i++) b[i] = 1ll * ((2 * b[i] % mod) - (c[i] * b[i] % mod * b[i] % mod) + mod) % mod;
		NTT(b,lim,-1);
		for(re int i=n;i<lim;i++) b[i] = 0;
		return ;
	}

方法3：牛顿迭代法

其实牛顿迭代法推出来的式子是和倍增法一样的，不过牛顿迭代法用途更为广泛。

我们先把 \(F(x) * G(x) \equiv 1(\bmod \ x^n)\) 来转化一下，把它们拆成相加的形式，那么就有

\[F(x) - \frac{1}{G(x)} \equiv 0(\bmod \ x^n) \]

你会发现这符合牛顿迭代法的样子，我们像牛顿迭代法那样的构造一个外层函数 \(H(t)\)，使它
能把上面这个式子表达出来，因为我们是要求 \(G(x)\) 的，所以 \(G(x)\) 应该放在 \(H(t)\) 里，那么就构造

\[H(t) = F(x) - \frac{1}{t} \]

令 \(t = G(x)\)，我们便得到了我们想要的式子：\(H(G(x)) \equiv 0(\bmod \ x^n)\)

类似的，假设我们已知 \(G_0(x) \equiv G(x)(\bmod \ x^{\lceil \frac{n}{2} \rceil})\)，我们要求 \(G_1(x) \equiv G(x)(\bmod \ x^n)\)

我们将它套入牛顿迭代公式，则有

\[G_1(x) = G_0(x) - \frac{H(G_0(x))}{H'(G_0(x))} \]

把 \(H\) 的解析式套进去，然后求导的时候要注意，我们是对 \(G(x)\) 求导，所以这里的 \(F(x)\) 是一个常数，常数求导后等于 \(0\)，那么有

\[G_1(x) = G_0(x) - \frac{F(x) - \frac{1}{G_0(x)}}{\frac{1}{G_0(x)^2}} \]

化简一下，则有

\[G_1(x) = G_0(x)(2-G_0(x)F(x)) \]

你会发现这和我们倍增得到的结果是一样的，这也说明它们本质上其实是相同的。

任意模数多项式乘法逆

把 NTT 换成 MTT ，然后再注意亿点细节即可。

感谢万能的 UOJ 群友解答。

三模 NTT 版本

三模 NTT 版本常数巨大无比，由于相乘会爆 long long，所以在一次递归中，三模 NTT 需要调用 \(18\) 次 DFT。常数可想而知，又加上题解区三模 NTT 版本的少之又少，导致我调了半上午 + 一下午才有了一个具体的理解。下面进入正文。

首先，因为 CRT 的合并会导致数很大，所以我们分成两次计算。

首先，我们计算出三个模数下的 \(F(x) \times H(x)\)，然后我们直接把它用 exCRT 合并起来，目的是直接把这个转化成 \(\bmod \ 10^9+7\) 的意义下(模板题给定的模数是 \(10^9+7\))进行，这样就不会出锅了。

之后，我们把它取反，变成 \(-F(x)H(x)\)，然后再将其变成 \(2 - F(x)H(x)\)，也就是在这个多项式的 \(x^0\) 项加上一个 \(2\)。这么简单的问题我竟然当时没绕过弯来，其实就是这个 \(2\) 其实是 \(x^0\) 次方的系数，你不能在所有的 \(x^k\) 的系数前都加上 \(2\)，这样显然不对。

然后我们在让这个多项式与 \(H(x)\) 相乘，这样便得到了 \(2H(x) - F(x)H(x)^2\)，也就是我们的 \(G(x)\)。

还是，注意取模。

代码

只截取了主要部分，其中 inv 函数求逆元，是用快速幂实现的。

#include<bits/stdc++.h>
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e5 + 5;
const ll mod1 = 469762049;
const ll mod2 = 998244353;
const ll mod3 = 1004535809;
const ll Mod = 1e9 + 7;
const ll g = 3;
using namespace std;

int n,m,lim;
ll invg,invn,inv1,inv12;
ll a[N],b1[N],b2[N],b3[N],R[N],B[N],aa[N],bb[N];

namespace Poly
{
	il void NTT(ll A[],int n,int type,ll mod)
	{
		//do something
	}
	il ll crt(ll a,ll b,ll c)//采用 exCRT 的合并方法就不会爆 long long了
	{
		ll k = ((b-a) % mod2 + mod2) % mod2 * inv1 % mod2;
		ll x = k * mod1 + a;
		k = ((c-x) % mod3 + mod3) % mod3 * inv12 % mod3;
		x = (x + k * mod1 % Mod * mod2 % Mod) % Mod;
		return x;
	}
	il void Mul(ll *A,ll *B,ll *C,int n,ll modx)
	{
		for(re int i=0;i<n;i++) aa[i] = A[i] , bb[i] = B[i];
		for(re int i=n;i<lim;i++) aa[i] = bb[i] = 0;
		NTT(aa,lim,1,modx) , NTT(bb,lim,1,modx);
		for(re int i=0;i<lim;i++) C[i] = aa[i] * bb[i] % modx;
		NTT(C,lim,-1,modx);
		for(re int i=n;i<lim;i++) C[i] = 0;
		return ;
	}
	il void Polyinv(int n)
	{
		if(n == 1) { B[0] = inv(a[0],Mod); return ; }
		Polyinv((n+1)>>1);
		lim = 1;
		while(lim < (n<<1)) lim <<= 1;
		for(re int i=0;i<lim;i++) R[i] = (R[i>>1]>>1) + ((i&1) ? lim/2 : 0);
		Mul(a,B,b1,n,mod1) , Mul(a,B,b2,n,mod2) , Mul(a,B,b3,n,mod3);//求 F(x)H(x)	
		for(re int i=0;i<n;i++) b1[i] = b2[i] = b3[i] = (Mod - crt(b1[i],b2[i],b3[i])) % Mod;//取反
		b1[0] += 2 , b2[0] += 2 , b3[0] += 2;//x^0项加2
		Mul(b1,B,b1,n,mod1) , Mul(b2,B,b2,n,mod2) , Mul(b3,B,b3,n,mod3);//求 2H(x) - F(x)H(x)^2
		for(re int i=0;i<n;i++) B[i] = crt(b1[i],b2[i],b3[i]);//合并之后就是答案
		return ;
	}
}
using namespace Poly;

signed main()
{
	n = read();
	for(re int i=0;i<n;i++) a[i] = read() % Mod;
	inv1 = inv(mod1 % mod2,mod2) , inv12 = inv(mod1 * mod2 % mod3,mod3);
	Poly::Polyinv(n);
	for(re int i=0;i<n;i++) cout << B[i] << " ";
	return 0;
}

拆系数 FFT 版本

发现拆系数 FFT 还是比三模 NTT 好写+好理解，相比于正常的拆系数 FFT 没有什么大变动。

仍然只保留主要部分。我采用的是 \(4\) 次 FFT 版本，相比于三模 NTT 来说，它每轮只需要 \(8\) 次调用 DFT，所以最后的速度比 NTT 快了一倍左右，FFT 是 10s，NTT 接近 20s。

#include<bits/stdc++.h>
#define int long long
#define ll long long
#define complex comple
#define double long double
#define re register
#define il inline
const int N = 4e5 + 5;
const int mod = 1e9 + 7;
const double Pi = acos(-1.0);
const ll bas = (1<<15) , bas2 = bas * bas;
using namespace std;
int max(int x,int y){return x > y ? x : y;}
int min(int x,int y){return x < y ? x : y;}

int n,m,lim;
ll a[N],b[N],c[N],R[N];
ll x,AC,AD,BC,BD;

namespace Poly
{
	struct complex
	{
		double x,y;
		complex (double xx = 0 , double yy = 0) { x = xx; y = yy; }
	}f[N],g[N],h[N];
    //运算符重载
    
	il ll ksm(ll a,ll b)
	{
		//do something
	} ll inv(ll x) { return ksm(x,mod-2); }
	il void FFT(complex A[],int n,int type)
	{
		// do something
	}
	il void Mul(ll a[],ll b[],ll c[],int n)
	{
		for(re int i=0;i<lim;i++) f[i] = g[i] = h[i] = {0,0};
		for(re int i=0;i<n;i++)
		{
			a[i] %= mod , b[i] %= mod;
			f[i].x = (a[i]>>15) , f[i].y = a[i] & 32767;
			h[i].x = (b[i]>>15) , h[i].y = b[i] & 32767;
		}//优化 4 次 FFT
		FFT(f,lim,1) , FFT(h,lim,1);
		g[0] = {f[0].x,-f[0].y};
		for(re int i=1;i<lim;i++) g[i] = {f[lim-i].x,-f[lim-i].y};
		for(re int i=0;i<lim;i++) h[i].x /= lim , h[i].y /= lim , f[i] = f[i] * h[i] , g[i] = g[i] * h[i];
		FFT(f,lim,-1) , FFT(g,lim,-1);
		for(re int i=0;i<n;i++)
		{
			AC = ((ll)((f[i].x+g[i].x)/2+0.5)) % mod;
			AD = ((ll)((f[i].y+g[i].y)/2+0.5)) % mod;
			BD = ((ll)(g[i].x-AC+0.5) % mod + mod) % mod;
			BC = ((ll)(f[i].y-AD+0.5) % mod + mod) % mod;
			AC = bas2 % mod * AC % mod , AD = bas % mod * (AD+BC) % mod;
			c[i] = (AC + AD + BD) % mod;
		}
	}
	il void Polyinv(int n)
	{
		if(n == 1) { b[0] = inv(a[0]); return ; }
		Polyinv((n+1)>>1);
		lim = 1;
		while(lim < (n<<1)) lim <<= 1;
		for(re int i=0;i<lim;i++) R[i] = (R[i>>1]>>1) + ((i&1) ? lim/2 : 0) , c[i] = 0;
		Mul(a,b,c,n);//算第一次
		for(re int i=0;i<n;i++) c[i] = (mod - c[i]) % mod;
		c[0] = (c[0] + 2) % mod;//取反求第二次
		Mul(b,c,b,n);
		for(re int i=n;i<lim;i++) b[i] = 0;
		return ;
	}
}
using namespace Poly;

signed main()
{
	n = read();
	for(re int i=0;i<n;i++) a[i] = read();
	Poly::Polyinv(n);
	for(re int i=0;i<n;i++) cout << b[i] << " ";
	return 0;
}

多项式对数函数（Polyln）

给出 \(n-1\) 次多项式 \(A(x)\)，求一个 \(\bmod{\:x^n}\) 下的多项式 \(B(x)\)，满足 \(B(x) \equiv \ln A(x)\).

在 \(\text{mod } 998244353\) 下进行，且 \(a_i \in [0, 998244353) \cap \mathbb{Z}\)

首先若 \(\ln f(x)\) 存在，则必有 \([x^0]f(x) = 1\)，因为 \(e\) 在模意义下没有对应的值，只有 \(\ln 1 = 0\) 是确定的。

在这种情况下，我们考虑给 \(f\) 数组进行变换。

因为一个多项式先求导再积分就是原函数加个常数，我们令这个常数为 \(0\)，然后就有了下面的变换：

\[\begin{aligned} \ln f(x) &= \int \text{d}(\ln f(x))\\ &= \int (\ln f(x))'\text{d}x \end{aligned} \]

这是一个复合函数的形式，我们先给它导一下

\[\ln f(x)= \int \frac{f'(x)}{f(x)} \]

到了这步其实就很简单了。

运用导数的知识，我们可以把分子求出来，多项式的导数还是很简单的，就是每项分别求导再相加就行了，运用 \((x^a)' = ax^{a-1}\)，然后相加就行了，然后注意第 \(n-1\) 项是 \(0\) 就行。

然后下面的这个式子我们可以用多项式求逆整出来，这个我们也会。

然后就是把这个求出来的多项式积分一下，这个也好整，运用 \(\displaystyle \int kx^a = \frac{k}{a+1}x^{a+1} + C\)（\(C\) 是常数）这个式子往里代就行，同样的也是对每项积分再相加，这样就得到了最终的式子。

时间复杂度是多项式求逆和 NTT 的 \(O(n\log n)\)

代码

#include<bits/stdc++.h>
//#define int long long
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e5 + 5;
const int mod = 998244353;
const ll g = 3;
using namespace std;

int n,m,lim,invg,invn;
int R[N];
ll a[N],b[N];

namespace Poly
{
	il ll ksm(ll a,ll b)
	{
		//快速幂
	} ll inv(ll x) { return ksm(x,mod-2); }
	il void NTT(ll A[],int n,int type)
	{
		// do something
	}
	il void Polymul(int n,ll a[],ll b[])
	{
		//多项式相乘
	}
	il void Polyinv(int n,ll a[],ll b[])
	{
		//多项式求逆
	}
	il void Polyderiv(int n,ll a[],ll b[])
	{
		for(re int i=0;i<n-1;i++)
			b[i] = a[i+1] * (i+1) % mod;
		b[n-1] = 0;
	}
	il void Polyinte(int n,ll a[],ll b[])
	{
		b[0] = 0;
		for(re int i=1;i<=n;i++) b[i] = a[i-1] * inv(i) % mod;
	}
	il void Polyln(int n,ll a[],ll b[])
	{
		static ll F[N],G[N];
		memset(F , 0 , sizeof F);
		memset(G , 0 , sizeof G);
		Polyderiv(n,a,F);
		Polyinv(n,a,G);
		Polymul(n,F,G);
		Polyinte(n,F,b);
	}
}
using namespace Poly;

signed main()
{
	n = read();
	for(re int i=0;i<n;i++) a[i] = read();
	Poly::Polyln(n,a,b);
	for(re int i=0;i<n;i++) cout << b[i] << " ";
	return 0;
}

多项式指数函数（Polyexp）

给出 \(n-1\) 次多项式 \(F(x)\)，求一个 \(\bmod{\:x^n}\) 下的多项式 \(G(x)\)，满足 \(G(x) \equiv \text e^{F(x)}\)。系数对 \(998244353\) 取模。

如果你细心的话，你会发现洛谷上对数函数和指数函数的样例刚好是反过来的，这也恰恰说明了这俩是逆运算（其实这是废话）。

和对数函数相同的，我们要保证 \([x^0]A(x)= 0\)，这样才能使得 \(B(x)\) 的常数项在模意义下有意义。

在这里我们会用到牛顿迭代法来解决。

首先，我们可以给左右两边取个 \(\ln\)，那么就有：

\[\ln(G(x)) \equiv F(x)(\bmod \ x^n) \]

移项，

\[\ln(G(x)) - F(x)\equiv 0(\bmod \ x^n) \]

我们发现这又变成了牛顿迭代的形式了，我们再令 \(G(x)\) 为自变量，给它套一个外层函数 \(H(t)\)，那么就是

\[H(t) = \ln t - F(x) \]

将 \(G(x)\) 代入

\[H(G(x)) = \ln(G(x)) - F(x) \]

还是一样的，假设已知 \(G_0(x) \equiv F(x)(\bmod \ x^{\lceil \frac{n}{2} \rceil})\)，我们要求 \(G_1(x) \equiv F(x)(\bmod \ x^n)\)

代入牛顿迭代公式，就有

\[G_1(x) = G_0(x) - \frac{H(G_0(x))}{H'(G_0(x))} \]

把 \(H\) 的解析式套进去

\[G_1(x) = G_0(x) - \frac{\ln(G_0(x)) - F(x)}{\frac{1}{G_0(x)}}(\bmod \ x^n) \]

化简一下

\[G_1(x) = G_0(x)(1-\ln G_0(x)+F(x)) \]

其中，\(G_0(x)\) 我们是已知的，这也就是说，我们可以采取类似于多项式求逆的递归，从小到大的把 \(G(x)\) 求出来。

时间复杂度 \(T(n) = T(n/2) + O(n\log n) = O(n\log n)\)

代码

#include<bits/stdc++.h>
//#define int long long
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e5 + 5;
const int mod = 998244353;
const ll g = 3;
using namespace std;

int n,m,lim,invg,invn;
int R[N];
ll a[N],b[N];

il int read()
{
	int f=0,s=0;
	char ch=getchar();
	for(;!isdigit(ch);ch=getchar()) f |= (ch=='-');
	for(; isdigit(ch);ch=getchar()) s = (s<<1) + (s<<3) + (ch^48);
	return f ? -s : s;
}

namespace Poly
{
	il ll ksm(ll a,ll b)
	{
		//do something
	} ll inv(ll x) { return ksm(x,mod-2); }
	il void NTT(ll A[],int n,int type)
	il void Polymul(int n,ll a[],ll b[])
	il void Polyinv(int n,ll a[],ll b[])
	il void Polyderiv(int n,ll a[],ll b[])
	il void Polyinte(int n,ll a[],ll b[])
	
   il void Polyln(int n,ll a[],ll b[])
	{
		static ll F[N],G[N];
		memset(F , 0 , sizeof F);
		memset(G , 0 , sizeof G);
		Polyderiv(n,a,F);
		Polyinv(n,a,G);
		Polymul(n,F,G);
		Polyinte(n,F,b);
	}
	il void Polyexp(int n,ll a[],ll b[])
	{
		if(n == 1) { b[0] = 1; return ; }
		Polyexp((n+1)>>1,a,b);
		static ll lnb[N];
		Polyln(n,b,lnb);
		lim = 1;
		while(lim < (n<<1)) lim <<= 1;
		for(re int i=0;i<n;i++) lnb[i] = ((a[i] - lnb[i]) % mod + mod) % mod;
		for(re int i=n;i<lim;i++) lnb[i] = b[i] = 0;
		lnb[0]++;
		Polymul(n,b,lnb);
		for(re int i=n;i<lim;i++) b[i] = 0;
	}
}
using namespace Poly;

signed main()
{
	n = read();
	for(re int i=0;i<n;i++) a[i] = read();
   invg = inv(g);
	Poly::Polyexp(n,a,b);
	for(re int i=0;i<n;i++) cout << b[i] << " ";
	return 0;
}

多项式幂函数（Polypow）

因为第一遍跑代码的时候 invg 不会被预处理导致调了 3h+。

普通版

给定一个 \(n-1\) 次多项式 \(F(x)\)，求一个在 \(\bmod\ x^n\) 意义下的多项式 \(G(x)\)，使得 \(G(x) \equiv (F(x))^k \ (\bmod\ x^n)\)。

多项式的系数在 \(\bmod\ 998244353\) 的意义下进行运算。

保证 \([x^0]F(x)=0\)。

首先朴素暴力快速幂肯定是 \(O(n\log n \log k)\) 的，但实际上我们可以继续优化它，这需要用到对数函数的性质：\(\ln x^a = a\ln x\)。

我们先把 \(F(x)^k\) 取个 \(\ln\)，变成 \(\ln(F(x)^k)\)，然后再应用上面的这个性质，就变成了 \(k\ln(F(x))\)，然后再对它求个 exp，于是这题就做完了，没错，就这么简单，并且它保证了 \([x^0]F(x) = 0\)，意味着可以取 \(\ln\)。

然后就有我上面说的，因为第一遍跑代码的时候 invg 不会被预处理导致调了 3h+，警钟敲烂。

还有一个问题是 \(k\) 的大小，在原题中 \(k \leq 10^{10^5}\)，很大。然后 EI 太阳神好像证明了一下当 \(n<p\) 时，一个多项式的幂次 \(k\) 直接对 \(p\) 取模是对的，具体可以看题解里的证明。

代码

#include<bits/stdc++.h>
#define int long long
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e5 + 5;
const int mod = 998244353;
const ll g = 3;
using namespace std;
int max(int x,int y){return x > y ? x : y;}
int min(int x,int y){return x < y ? x : y;}

int n,m,lim,invg,invn;
int R[N];
ll k,a[N],b[N];

namespace Poly
{
	il ll ksm(ll a,ll b)
	il ll inv(ll x) { return ksm(x,mod-2); }
	il void NTT(ll A[],int n,int type)
	il void Polymul(int n,ll a[],ll b[])
	il void Polyinv(int n,ll a[],ll b[])
	il void Polyderiv(int n,ll a[],ll b[])
	il void Polyinte(int n,ll a[],ll b[])
	il void Polyln(int n,ll a[],ll b[])
	il void Polyexp(int n,ll a[],ll b[])
    
	il void Polypow(int n,ll k,ll a[],ll b[])
	{
		static ll lna[N];
		Polyln(n,a,lna);
		for(re int i=0;i<n;i++) lna[i] = lna[i] * k % mod;
		Polyexp(n,lna,b);
		memset(lna , 0 , sizeof lna);
	}
}
using namespace Poly;

signed main()
{
	n = read() , k = read();
	invg = inv(g);//because of you
	for(re int i=0;i<n;i++) a[i] = read();
	Poly::Polypow(n,k,a,b);
	for(re int i=0;i<n;i++) cout << b[i] << " ";
	return 0;
}

加强版

加强版和普通版就差了一个地方，并不保证 \(a_0 = 1\) 。

为此，我们应该把这个多项式转化成 \(a_0=1\) 的样子。

具体的，我们从 \(F(x)\) 中提取一个公因式 \(b \times x^k\)，其中 \(k\) 满足 \(a_k \not = 0\)，且 \(\forall 0 \leq i <k,a_i =0\)，\(b\) 就是 \(x^k\) 这一项的系数，这样就可以保证最低位是 \(1\) 了。

举个例子，我们现在是 \(0 + 2x + 2x^2\)，我们需要提取出一个 \(2x\)，式子就变成了 \(2x(1+x)\)，然后对右边的多项式进行快速幂，然后再乘上左边的系数的幂次，再右移即可。

注意一种特殊情况：\(k \times m \ge n\)（这里的 \(k\) 是上文的 \(k\)，\(m\) 就是题目中给定的 \(k\)，也就是幂次），这也就是说前 \(k\times m\) 项都是 \(0\) 了，我们也就没必要往后做了，输出 \(n\) 个 \(0\) 即可。

还有这个题的取模问题，我们提取出来的公因式的次幂，其实是对 \(\varphi(p)\) 取模的，因为它是一个单项式，应该用扩展欧拉定理来降幂，也就是模 \(\varphi(p)\)，右边的那个多项式还是同理的对 \(p\) 取模。

代码

根据上文，你需要存三个 \(k\)：

\(k\)：对 \(p\) 取模。

\(k_1\)：对 \(\varphi(p)\) 取模。

\(k_2\)：用来求在不模 \(p\) 的意义下，每个项会往右平移多少，这个不好说清楚，它的作用就是和 \(m\) 相乘判断前 \(n\) 位是不是 \(0\) 的，也就是看前 \(m\) 个系数为 \(0\) 的项在 \(k_2\) 次方后会产生前多少项系数是 \(0\)。

#include<bits/stdc++.h>
#define int long long
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e5 + 5;
const int mod = 998244353;
const ll g = 3;
using namespace std;

int n,m,lim,pos,val;
int invg,invn,invval;
int R[N];
ll k,k1,k2,a[N],b[N];

namespace Poly
{
	il void Getk()
	{
		char ch=getchar();
		for(;!isdigit(ch);ch=getchar());
		for(; isdigit(ch);ch=getchar())
		{
			k = ((k<<1) + (k<<3) + (ch^48)) % mod;
			k1 = ((k1<<1) + (k1<<3) + (ch^48)) % (mod-1);
			if(10*k2+(ch^48) <= mod) k2 = 10*k2+(ch^48);
		}
	}
	il ll ksm(ll a,ll b)
	ll inv(ll x) { return ksm(x,mod-2); }
	il void NTT(ll A[],int n,int type)
	il void Polymul(int n,ll a[],ll b[])
	il void Polyinv(int n,ll a[],ll b[])
	il void Polyderiv(int n,ll a[],ll b[])
	il void Polyinte(int n,ll a[],ll b[])
	il void Polyln(int n,ll a[],ll b[])
	il void Polyexp(int n,ll a[],ll b[])
	il void Polypow(int n,ll k,ll a[],ll b[])
    
	il void Polyexpow(int n,ll k,ll k1,ll k2,ll a[],ll b[])
	{
		pos = n;
		for(re int i=0;i<n;i++)
			if(a[i]) { pos = i , val = a[i]; break; }
		if(pos * k2 >= n)//超过n次项了
		{
			for(re int i=0;i<n;i++) cout << "0 ";
			exit(0);
		}
		invval = inv(val);
		for(re int i=0;i<n-pos;i++) a[i] = a[i+pos] * invval % mod;
		for(re int i=n-pos;i<n;i++) a[i] = 0;
		//for(re int i=0;i<n-pos;i++) cout << a[i] << " ";
		//cout << endl;
		Polypow(n-pos,k,a,b);
		//for(re int i=0;i<n-pos;i++) cout << b[i] << " ";
		//cout << endl;
		val = ksm(val,k1) , pos = min(k2*pos,n);//左边的单项式的k1次幂
		for(re int i=n-1;i>=pos;i--) b[i] = b[i-pos] * val % mod;//看看向右平移后有多少不是 0 的项
		for(re int i=pos-1;i>=0;i--) b[i] = 0;//前k2*pos项都是0
		return ;
	}
}
using namespace Poly;

signed main()
{
	n = read();
	Getk();
	invg = inv(g);
	for(re int i=0;i<n;i++) a[i] = read();
	Poly::Polyexpow(n,k,k1,k2,a,b);
	for(re int i=0;i<n;i++) cout << b[i] << " ";
	return 0;
}

多项式开根（Polysqrt）

前置知识：多项式快速幂。

给定一个 \(n-1\) 次多项式 \(F(x)\)，求一个在 \({} \bmod x^n\) 意义下的多项式 \(G(x)\)，使得 \(G^2(x) \equiv F(x) \pmod{x^n}\)。若有多解，请取零次项系数较小的作为答案。

多项式的系数在 \({}\bmod 998244353\) 的意义下进行运算。

普通版

保证 \(a_0 =1\)。

你会发现，这其实就是 \(G(x) \equiv F(x)^{\frac{1}{2}}\)，这个 \(\displaystyle \frac{1}{2}\) 其实就是 \(2\) 的逆元，那也就是求 \(G(x) \equiv F(x)^{\text{inv}(2)}\)，并且保证 \(a_0 =1\)，所以多项式快速幂是可行的。

实际上根据这个，你开几次根都是可行的，只需要跑它的逆元次幂即可。

代码

开 \(k\) 次根，代码就这么简单。

il void Polysqrt(int n,ll k,ll a[],ll b[])
{
	Polypow(n,k,a,b);
}

signed main()
{
	n = read() , k = read();
	invg = inv(g);
	for(re int i=0;i<n;i++) a[i] = read();
	Polysqrt(n,inv(k),a,b);
	for(re int i=0;i<n;i++) cout << b[i] << " ";
	return 0;
}

加强版

不保证 \(a_0 = 1\)。

参考博客：

题解 P5277 多项式开根（加强版）- Great_Influence

以开平方根为例。

我们先同理的跟幂函数加强版的思路相同的提取出一个公因数 \(a \times x^k\)，那么最后式子转化为

\[\large F(x)^{\frac{1}{2}} = a^{\frac{1}{2}}x^{k \cdot \frac{1}{2}}e^{\frac{1}{2}\ln(\frac{F(x)}{b\times x^k})} \]

右边的 \(e\) 的那一坨其实就用逆元跑快速幂即可。但是左边 \(b^{\frac{1}{2}}\) 要用到的其实是二次剩余。

这里我们不考虑二次剩余，我们给出一个 \(\mathcal O(\sqrt{\bmod})\) 的 BSGS 解法。

因为 \(\displaystyle \sqrt[k]{a} = a^{\frac{1}{k}}\)，所以我们实际上求的就是 \(\frac{1}{k}\) 次方根，而在模意义下，这个数等价于 \(a^{k^{-1}}\)。

而根据扩展欧拉定理，\(k^{-1}\) 应该是在模 \(\bmod-1\) 意义下的，也就是说 \(k^{-1}\) 有可能没有逆元。

接下来我们给出一个求 \(k\) 次剩余的方法。

因为我们做的是多项式题，所以一般模数就是类似于 \(998244353\) 这样的 NTT 模数，所以一般是有原根的。

那么 \(a = 3^t\) 是一定存在的，我们要求的就是 \(3^{\frac{t}{k}}\)。

如果 \(\frac{t}{k}\) 在模 \(\bmod -1\) 意义下不存在，那么说明 \(a\) 不存在 \(k\) 次剩余。

实际实现就是在 BSGS 后得到 \(t\)，然后判断一下是否 \(k \mid t\)，如果是的话，用 exgcd 求出 \(k\) 在模 \(998244352\) 意义下的逆元，再求出 \(3^{xk^{-1}}\) 作为 \(k\) 次剩余。

因为模板题给的是二次剩余，所以不用 exgcd 也行了。

其实我也不是很懂，明天去 U 群问一下，得到答案再补充。

代码

#include<bits/stdc++.h>
#define int long long
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e5 + 5;
const int mod = 998244353;
const ll g = 3;
using namespace std;
int max(int x,int y){return x > y ? x : y;}
int min(int x,int y){return x < y ? x : y;}

int n,m,lim,pos,val;
int invg,invn,invval;
int R[N];
ll k,k1,k2,a[N],b[N];

map <int,int> Hash;

il ll BSGS(ll x)
{
	m = ceil(sqrt(mod));
	Hash[x] = 0 ;ll t = x;
	for(re int i=1;i<m;i++)
	{
		t = t * 3 % mod;
		Hash[t] = i;
	}
	ll Val = ksm(3,m); t = 1;
	for(re int i=1;;i++)
	{
		t  = t * Val % mod;
		if(Hash.count(t)) {  return ksm(3,(i*m-Hash[t])/2); }
	}
}
	il void Polyexsqrt(int n,ll k,ll a[],ll b[])
{
	pos = n;
	for(re int i=0;i<n;i++)
		if(a[i]) { pos = i , val = a[i]; break; }
	if(pos == n)
	{
		for(re int i=0;i<n;i++) cout << "0 ";
		exit(0);
	}
	invval = inv(val);
	for(re int i=0;i<n-pos;i++) a[i] = a[i+pos] * invval % mod;
	for(re int i=n-pos;i<n;i++) a[i] = 0;
	Polypow(n,inv(k),a,b);
	val = BSGS(val) , pos = pos / 2;
	if(val > (mod+1)/2) val = mod - val;
	for(re int i=n-1;i>=pos;i--) b[i] = b[i-pos] * val % mod;
   for(re int i=pos-1;i>=0;i--) b[i] = 0;
	}
}

signed main()
{
	n = read() , k = 2;
	invg = inv(g);
	for(re int i=0;i<n;i++) a[i] = read();
	Polyexsqrt(n,k,a,b);
	for(re int i=0;i<n;i++) cout << b[i] << " ";
	return 0;
}

总结

以上就是多项式计算的一些进阶计算了，还有一些更进阶的计算等以后有时间提升一下数学基础并且有时间了再继续学吧。

posted @ 2023-06-26 20:42 Bloodstalk 阅读(32) 评论(0) 编辑收藏举报

刷新页面返回顶部

bloodstalk