多项式全家桶

每次复习完下一次都会忘，这次下定决心一定要记下来！！！

FFT 和 NTT

做法

太简单了，看原题吧，不会退役吧。

代码

int n, m, rev[maxn];
int qpow(int x, int k, int p) {
	int res = 1;
	while(k) {
		if(k & 1)
			res = res * x % p;
		k >>= 1, x = x * x % p;
	}
	return res;
}
void prepare(int len) {
	for (int i = 0; i < len; i++) {
		rev[i] = rev[i >> 1] >> 1;
		if(i & 1)
			rev[i] |= len >> 1;
	}
}
struct Poly {
	vector<int> a;
	int size() {
		return a.size();
	}
	int& operator[](int x) {
		return a[x];
	}
	void resize(int N) {
		a.resize(N);
	}
	void NTT(int f) {
		for (int i = 1; i < a.size(); i++)
			if(i < rev[i])
				swap(a[rev[i]], a[i]);
		for (int h = 2; h <= a.size(); h <<= 1) {
			int d = qpow((f == 1 ? gb : gi), (mod - 1) / h, mod);
			for (int i = 0; i < a.size(); i += h) {
				int nw = 1;
				for (int j = i; j < i + h / 2; j++) {
					int a0 = a[j], a1 = a[j + h / 2] * nw % mod;
					a[j] = (a0 + a1) % mod, a[j + h / 2] = (a0 - a1 + mod) % mod;
					nw = nw * d % mod;
				}
			}
		}
		if(f == -1) {
			int inv = qpow(a.size(), mod - 2, mod);
			for (int i = 0; i < a.size(); i++)
				a[i] = a[i] * inv % mod;
		}
	}
	void read() {
		for (int i = 0; i < a.size(); i++)
			cin >> a[i];
	}
	void print() {
		for (int i = 0; i < a.size(); i++)
			cout << a[i] << " ";
		cout << endl;
	}
	friend Poly operator*(Poly f, Poly g) {
		int len = 1, t = f.size() + g.size() - 1;
		while(len < t)
			len <<= 1;
		prepare(len);
		f.resize(len), g.resize(len);
		f.NTT(1), g.NTT(1);
		for (int i = 0; i < len; i++)
			f[i] = f[i] * g[i] % mod;
		f.NTT(-1); f.resize(t);
		return f;
	}
}

多项式求逆

做法

假设我们已经求出来了：

$$A\times B^{\prime }\equiv 1(\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })$$

那么我们要求：

$$A\times B\equiv 1(\mathrm{mod}{x}^n)$$

我们先做差，得到：

$$(B-B^{\prime })\equiv 0(\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })$$

然后平方，得到：

$$(B-B^{\prime }{)}^2\equiv 0(\mathrm{mod}{x}^n)$$

$$B^2-2B{B}^{\prime }+B^{\prime 2}\equiv 0(\mathrm{mod}{x}^n)$$

左右同时乘 $A$，得到：

$$B-2B^{\prime }+A{B}^{\prime 2}\equiv 0(\mathrm{mod}{x}^n)$$

$$B\equiv -A{B}^{\prime 2}+2B^{\prime }(\mathrm{mod}{x}^n)$$

按这个柿子计算即可，注意没有必要把 NTT 出来的结果直接 NTT 回来，可以直接乘起来加起来。

代码

Poly get_inv(Poly f, int lim) {
	if(lim == 1) {
		Poly ans; ans.resize(1);
		ans[0] = qpow(f[0], mod - 2, mod);
		return ans;
	}
	Poly ans = get_inv(f, lim + 1 >> 1);
	int len = 1;
	while(len < lim * 2)
		len <<= 1;
	prepare(len);
	f.resize(lim), f.resize(len), ans.resize(len);
	f.NTT(1), ans.NTT(1);
	for (int i = 0; i < len; i++)
		ans.a[i] = ans.a[i] * (2 - f.a[i] * ans.a[i] % mod + mod) % mod;
	ans.NTT(-1);
	ans.resize(lim);
	return ans;
}

注意，这里需要先 f.resize(lim)，因为我们需要排除多余的位。

多项式开根

做法

我们假设 $B^{\prime 2}(x)=F(x)(\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil }),B^2=F(x)(\mathrm{mod}{x}^n)$。

那么做差再平方得到：

$$(B^{\prime }(x)-B(x){)}^2=0(\mathrm{mod}{x}^n)$$

$$B^{\prime }(x{)}^2-2B^{\prime }(x)B(x)+A(x)=0(\mathrm{mod}{x}^n)$$

$$B(x)=\frac{B^{\prime }(x{)}^2+A(x)}{2B^{\prime }(x)}(\mathrm{mod}{x}^n)$$

直接做就可以。

代码

Poly get_sqrt(Poly f, int lim) {
	if(lim == 1) {
		Poly ans; ans.resize(1);
		ans.a[0] = 1;
		return ans;
	}
	Poly ans = get_sqrt(f, lim + 1 >> 1), inv = get_inv(ans, lim);
	int len = 1;
	while(len < lim * 2)
		len <<= 1;
	f.resize(lim), f.resize(len), ans.resize(len), inv.resize(len);
	prepare(len);
//	cout << lim << endl;
	inv.NTT(1), f.NTT(1);
	for (int i = 0; i < len; i++)
		inv.a[i] = inv.a[i] * f.a[i] % mod;
	inv.NTT(-1);
	for (int i = 0; i < len; i++)
		ans.a[i] = (mod + 1) / 2 * (inv.a[i] + ans.a[i]) % mod;
	ans.resize(lim);
	return ans;
}

一点记忆的小技巧：求逆和开根都是单式做的，所以要用倍增的想法去构造。

多项式除法

做法

我们记 $F_R(x)=x^{n}F(\frac{1}{x})$，那么就有：

$$F(x)=G(x)Q(x)+R(x)$$

$$F(\frac{1}{x})=G(\frac{1}{x})Q(\frac{1}{x})+R(\frac{1}{x})$$

$$x^{n}F(\frac{1}{x})=x^{m}G(\frac{1}{x})x^{n-m}Q(\frac{1}{x})+x^{n-m+1}\times x^{m-1}R(\frac{1}{x})$$

$$F_R(\frac{1}{x})=G_R(\frac{1}{x})Q_R(\frac{1}{x})+x^{n-m+1}{R}_R(\frac{1}{x})$$

$$F_R(\frac{1}{x})=G_R(\frac{1}{x})Q_R(\frac{1}{x})(\mathrm{mod}{x}^{n-m+1})$$

$$Q_R(x)=G_R^{-1}(x)F_R(x)(\mathrm{mod}{x}^{n-m+1})$$

然后做就可以了。

代码

这里只给出核心代码，重载的运算符就没有了。

pair<Poly, Poly> get_div(Poly f, Poly g) {
	int n = f.size(), m = g.size();
	if(f.size() < g.size()) {
		Poly q, r;
		q.resize(n - m + 1); r = g;
		return make_pair(q, r);
	}
	Poly q, r;
	q = get_rev(f) * (get_inv(get_rev(g), n - m + 1));
	q.resize(n - m + 1); q = get_rev(q);
	r = f - g * q;
	r.resize(m - 1);
	return make_pair(q, r);
}

多项式 ln

做法

我们要求 $G(x)=\ln F(x)(\mathrm{mod}{x}^n)$。

你可能会疑问为什么多项式会有 ln，我也不知道，但是就是对 EGF，OGF 的优化非常有用。

如何求呢？我们对两边同时求导，得到：

$$G^{\prime }(x)=\frac{F^{\prime }(x)}{F(x)}$$

我们直接对 $F$ 求逆然后乘上 $F^{\prime }$ 就可以了。

代码

int inv[maxn];
void prework_inv() {
	inv[1] = 1;
	for (int i = 2; i <= n; i++)
		inv[i] = (mod - mod / i) * inv[mod % i] % mod;
}
Poly get_ji(Poly f) {
	Poly g; g.resize(f.size() + 1);
	g.a[0] = 0;
	for (int i = 1; i < g.size(); i++)
		g[i] = f[i - 1] * inv[i] % mod;
	return g;
}
Poly get_dao(Poly f) {
	Poly g; g.resize(f.size() - 1);
	for (int i = 0; i < g.size(); i++)
		g[i] = f[i + 1] * (i + 1) % mod;
	return g;
}
Poly get_ln(Poly f) {
	prework_inv();
	Poly ans = get_ji(get_dao(f) * get_inv(f, f.size())); ans.resize(f.size());
	return ans;
}

多项式 exp

牛顿迭代

我们先来看一个经典的问题：给定一个数 $a$，求 $b^2=a$，精确到小数点后若干位，并且手算。

那么你会说我会二分啊，但是比如 $a={10}^9$ 的时候二分其实挺难做的。

那么我们就可以请出牛顿迭代了。

我们可以把上面这个问题抽象为求 $f(x)=x^2-a$ 的零点。

牛顿迭代的做法是我们先选取 $t=a$ 作为我们的初始点，然后每次做一条 $f(x)$ 在 $(t,f(t))$ 处的切线，找其与 $x$ 轴的交点 $(t^{\prime },0)$，然后令 $t=t^{\prime }$。

根据很简单的知识我们知道这条直线为 $y=f^{\prime }(x)(x-t)+f(t)$。那么零点为 $x=t-\frac{f(t)}{f^{\prime }(t)}$。

我们发现，这样会一步步逼近我们的答案，并且次数非常小。

做法

我们对于多项式同样做牛顿迭代。

我们假设我们要求 $ln(B(x))=A(x)(\mathrm{mod}{x}^n)$，那么我们记 $F(G(x))=\ln (G(x))-A(x)$。

那么我们如果求出了 $ln(B_0(x))=A(x)(\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })$，那么我们就可以去做一次牛顿迭代就行。

那么此时有：

$$B(x)=B_0(x)-\frac{F(B_0(x))}{F^{\prime }(B_0(x))}$$

$$B(x)=B_0(x)-(\ln B_0(x)-A(x))B_0(x)$$

$$B(x)=B_0(x)(1-ln{B}_0(x)+A(x))$$

注意这里是求 $F$ 的导，其实并不是一个复合函数，你可以把 $G(x)$ 视为一个变量，所以不需要对里面求导。

按这个柿子做即可。

代码

Poly get_exp(Poly f, int lim) {
	if(lim == 1) {
		Poly ans; ans.resize(1); ans[0] = 1;
		return ans;
	}
	Poly ans = get_exp(f, lim + 1 >> 1);
	f.resize(lim), ans.resize(lim);
	ans = ans * (1 - get_ln(ans) + f);
	ans.resize(lim);
	return ans;
}

多项式快速幂

做法

你也许会说直接类似于正常快速幂，但是复杂度 $O(n{\log }^{2}n)$。

其实我们可以变成 $e^{k\ln F(x)}$ 这么求值就可以了，是单 $\log$ 的。

代码

Poly qpow(Poly f, int k) {
	return get_exp(get_ln(f) * k, f.size());
}

多项式多点求值

做法

我们考虑其实 $f(a_i)=f(x)(\mathrm{mod}x-a_i)$。

那么我们就可以对这个东西分治，就是我们记 $g(l,r)=\prod _{i=l}^r(x-a_i)$，然后每次对左儿子和右儿子分别取模，一步步下去，复杂度 $O(n{\log }^{2}n)$ 但是常数极大，需要一点逆天的优化（比如区间长度不大于一定值就暴力（））。

据说有常数更小的厉害做法，但是我看不懂（），所以给出简单的这个做法，卡卡能过。

代码

int a[maxn];
Poly f[maxn], g;
void build(int l, int r, int t) {
	if(l == r) {
		f[t].resize(2), f[t][1] = 1, f[t][0] = mod - a[l];
		return ;
	}
	int mid = l + r >> 1;
	build(l, mid, t << 1), build(mid + 1, r, t << 1 | 1);
	f[t] = f[t << 1] * f[t << 1 | 1];
}
void solve(int l, int r, Poly &g, int t) {
	if(l + 2000 >= r) {
		for (int i = l; i <= r; i++) {
			int res = 0;
			for (int j = g.size() - 1; j >= 0; j--) {
				res = 1ll * res * a[i] % mod + g[j]; res %= mod;
			}
			cout << res << endl;
		}
		return ;
	}
	int mid = l + r >> 1;
	Poly l1 = get_div(g, f[t << 1]).second;
	solve(l, mid, l1, t << 1); 
	l1 = get_div(g, f[t << 1 | 1]).second;
	solve(mid + 1, r, l1, t << 1 | 1);
}