「学习笔记」多项式科技
做多项式题就像嗑药，出多项式题就像贩毒。 ——某著名 OIer
多项式乘法
快速傅里叶变换 (FFT)
直接上链接（
总的来说就是先 DFT 从系数表示法到点值表示法，再 IDFT 从点值表示法到系数表示法。
简单说一下不太理解的，在 DFT 中 \(\omega_n^k = -\omega_n^{k+\frac{n}{2}}\)，其实就是在单位圆上旋转了 \(180°\)。
感性理解 DFT 和 IDFT 是逆运算，所以 \(\omega_n\) 就是相反数。
Code
#include <bits/stdc++.h>
#define ll long long
#define db double
#define gc getchar
#define pc putchar
using namespace std;
namespace IO
{
template <typename T>
void read(T &x)
{
x = 0; bool f = 0; char c = gc();
while(!isdigit(c)) f |= c == '-', c = gc();
while(isdigit(c)) x = x * 10 + c - '0', c = gc();
if(f) x = -x;
}
template <typename T>
void write(T x)
{
if(x < 0) pc('-'), x = -x;
if(x > 9) write(x / 10);
pc('0' + x % 10);
}
}
using namespace IO;
struct Complex
{
db x, y;
Complex(db _x = 0.0, db _y = 0.0)
{
x = _x, y = _y;
}
Complex operator + (const Complex b) const
{
return Complex(x + b.x, y + b.y);
}
Complex operator - (const Complex b) const
{
return Complex(x - b.x, y - b.y);
}
Complex operator * (const Complex b) const
{
return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
}
};
const int N = 1e6 + 5;
const db pi = acos(-1.0);
int n, m;
Complex f[N << 2], g[N << 2];
int len = 1, bit, rev[N << 2];
void FFT(Complex a[], int type)
{
for(int i = 0; i < len; i++)
if(i < rev[i]) swap(a[i], a[rev[i]]);
for(int mid = 1; mid < len; mid <<= 1)
{
Complex wn(cos(pi / mid), type * sin(pi / mid));
for(int i = 0; i < len; i += (mid << 1))
{
Complex w(1, 0);
for(int j = 0; j < mid; j++, w = w * wn)
{
Complex x = a[i + j], y = w * a[i + mid + j];
a[i + j] = x + y;
a[i + mid + j] = x - y;
}
}
}
if(type == -1)
for(int i = 0; i < len; i++)
a[i].x /= len;
return;
}
int main()
{
read(n), read(m);
for(int i = 0; i <= n; i++) read(f[i].x);
for(int i = 0; i <= m; i++) read(g[i].x);
while(len <= n + m) len <<= 1, bit++;
for(int i = 0; i < len; i++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
FFT(f, 1);
FFT(g, 1);
for(int i = 0; i < len; i++) f[i] = f[i] * g[i];
FFT(f, -1);
for(int i = 0; i <= n + m; i++)
printf("%lld ", (ll)(f[i].x + 0.5));
pc('\n');
return 0;
}
// A.S.
快速数论变换 (NTT)
继续上链接（
就是把 FFT 中的 \(\omega_n\) 改为了 \(g^{\frac{p-1}{n}}\)，其中 \(g\) 为模数的原根。
证明就算了
然后在 IDFT 时就感性理解地取一下逆元就行了（
Code
#include <bits/stdc++.h>
#define ll long long
#define db double
#define gc getchar
#define pc putchar
#define swap(a, b) a ^= b ^= a ^= b
using namespace std;
namespace IO
{
template <typename T>
void read(T &x)
{
x = 0; bool f = 0; char c = gc();
while(!isdigit(c)) f |= c == '-', c = gc();
while(isdigit(c)) x = x * 10 + c - '0', c = gc();
if(f) x = -x;
}
template <typename T>
void write(T x)
{
if(x < 0) pc('-'), x = -x;
if(x > 9) write(x / 10);
pc('0' + x % 10);
}
}
using namespace IO;
const int N = 1e6 + 5;
const int p = 998244353;
const int G = 3;
const int Gi = 332748118;
int n, m, len = 1, bit;
ll f[N << 2], g[N << 2];
int rev[N << 2];
ll qpow(ll a, int b)
{
ll res = 1;
while(b)
{
if(b & 1) res = res * a % p;
a = a * a % p, b >>= 1;
}
return res % p;
}
void NTT(ll a[], int type)
{
for(int i = 0; i < len; i++)
if(i < rev[i]) swap(a[i], a[rev[i]]);
for(int mid = 1; mid < len; mid <<= 1)
{
ll wn = qpow(type == 1 ? G : Gi, (p - 1) / (mid << 1));
for(int i = 0; i < len; i += (mid << 1))
{
ll w = 1;
for(int j = 0; j < mid; j++, w = w * wn % p)
{
ll x = a[i + j], y = w * a[i + mid + j] % p;
a[i + j] = (x + y) % p;
a[i + mid + j] = (x - y + p) % p;
}
}
}
ll inv = qpow(len, p - 2);
if(type == -1)
for(int i = 0; i < len; i++)
f[i] = f[i] * inv % p;
return;
}
int main()
{
read(n), read(m);
for(int i = 0; i <= n; i++) read(f[i]);
for(int i = 0; i <= m; i++) read(g[i]);
while(len <= n + m) len <<= 1, bit++;
for(int i = 0; i < len; i++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
NTT(f, 1);
NTT(g, 1);
for(int i = 0; i < len; i++) f[i] = f[i] * g[i] % p;
NTT(f, -1);
for(int i = 0; i <= n + m; i++)
write(f[i]), pc(' ');
pc('\n');
return 0;
}
// A.S.
多项式求逆
对于一个多项式 \(F(x)\) 求 一个多项式 \(G(x)\)，满足 \(F(x)*G(x)\equiv 1\ (\bmod x^n)\)，系数对 \(998244353\) 取模。
就是多项式的逆元。
推一下式子
设
那么
我们要把模数改为 \(x^n\)，只需要平方一下
将等式左边乘上 \(A(x)\) 不影响等号
因为
所以可以将 \(A(x)*B(x)\) 都去掉
我们只需要求出 \(C(x)\)，而 \(C(x)\) 与 \(B(x)\) 的形式是一样的，只是模数不一样，所以可以递归求解，当然也可以递推。
多项式 ln
给定一个多项式 \(A(x)\)，求一个多项式 \(B(x)\)，满足 \(B(x)\equiv \ln A(x)\quad (\bmod x^n)\)。
设 \(f(x) = \ln x\)
\(\ln\) 不好处理，但是对 \(\ln\) 求导后就很好算了，\(f'(x)=\dfrac{1}{x}\)
所以将同余号两边同时求导，\(B'(x)\equiv f'(A(x))*A'(x)\) （复合函数求导）
因为 \(f'(A(x))=\dfrac{1}{A(x)}\)
所以 \(B'(x)\equiv \dfrac{A'(x)}{A(x)}\)
有了 \(B'(x)\) 后再积分求出 \(B(x)\)
求导公式：\((x^a)'=ax^{a-1}\)
积分公式：\(\int x^adx=\dfrac{1}{a+1}x^{a+1}\)
积分就是求导的逆运算，你会发现把 \(\dfrac{1}{a+1}x^{a+1}\) 求导后就是 \(x^a\)。
多项式开根
给定一个多项式 \(A(x)\)，求一个多项式 \(B(x)\)，满足 \(B^2(x)\equiv A(x)\quad(\bmod x^n)\)。
设
那么
然后就可以递归求解了，只需要多项式求逆。
多项式除法
给定 \(n\) 次多项式 \(F(x)\) 和 \(m\) 次多项式 \(G(x)\)，求多项式 \(Q(x)\) 和 \(R(x)\)，满足 ：
\(Q(x)\) 次数为 \(n - m\) ，\(R(x)\) 次数小于 \(m\)
\(F(x)= Q(x)*G(x)+R(x)\)
大力推式子
我们发现 \(x^nF(\frac{1}{x})\) 其实就是 \(F(x)\) 翻转一下，将其设为 \(F_R(x)\)，其余同理。
多项式求逆即可
有了 \(Q_R(x)\) 也就有了 \(Q(x)\)，然后
多项式 exp
给定一个多项式 \(A(x)\)，求多项式 \(B(x)\)，满足 \(B(x)\equiv e^{A(x)}\quad (\bmod x^n)\)。系数对 \(998244353\) 取模。
牛顿迭代
对于形如 \(F(G(x))\) 的复合函数，求满足 \(F(G(x))=0\) 的 \(G(x)\)。
通过这个式子可以迅速逼近真实值。
在本题中
令
只需要使
其中 \(G_0(x)\) 满足
因为对于 \(F(G(x))\) 来说 \(G(x)\) 为自变量，\(A(x)\) 相当于一个常量，所以
那么
然后就可以递归处理了，需要多项式 ln。
多项式快速幂
给定一个多项式 \(A(x)\)，求多项式 \(B(x)\)，满足 \(B(x)\equiv A^k(x)\quad (\bmod x^n)\)，系数对 \(998244353\) 取模。
普通版
保证 \(a_0=1\)。
推一下式子
只需要多项式 exp 和多项式 ln。
\(k\) 能取模就感性理解一下吧
加强版
不保证 \(a_0=1\)
但是我们需要让 \(a_0=1\) 才能进行多项式的计算，所以我们可以把多项式都除以 \(a_0\)，最后再乘回来。
这样在 \(a_0\neq 0\) 时是可以的，但是如果 \(a_0=0\) 呢，其实只需要把前缀的 \(0\) 都删掉，相当于降次，记录一下最终结果前缀会有多少个 \(0\)，向右平移即可。
代码细节比较多（
还有就是指数的取模，对于多项式的指数可以直接 \(\bmod p\)，但是对于一开始除掉的 \(a_0\)，指数需要 \(\bmod \varphi(p)\)
Code
ll lna[N << 2], lowa[N << 2];
void Qpow(ll *a, int k, ll *b, int n, int k2) // k = k % p, k2 = k % phi(p)
{
Clear(b, 0, n);
int shift = 0;
while(!a[shift]) shift++;
if((ll)shift * k >= n) return;
n -= shift;
for(int i = 0; i < n; i++) lowa[i] = a[i + shift];
int low0 = lowa[0], inv0 = inv(low0);
for(int i = 0; i < n; i++) lowa[i] = lowa[i] * inv0 % p;
Ln(lowa, lna, n);
for(int i = 0; i < n; i++) lna[i] = lna[i] * k % p;
Exp(lna, b, n);
shift *= k;
for(int i = n + shift - 1; i >= shift; i--)
b[i] = b[i - shift] * qpow(low0, k2) % p;
Clear(b, 0, shift);
return;
}
