多项式专题
除法、开根、exp、多点插值、快速求值、拉格朗日插值的板子还没有整理,重心拉格朗日插值法还不会,正在补坑中。
\(\color{red}{\text{约定:}}\)
\(1.F(x)\)表示一个普通的项数为\(2\)的幂次多项式,\(F_D(x)\)是他的点值表示。
\(2.w\)代表单位根,\(w_m\)表示\(m\)次单位根。
\(3.A\)代表一个数列。
\(4.g\)表示原根。
\(5.F^{(z)}(x)\)表示\(F(x)\)的\(z\)次导数。
\(\color{red}{\text{多项式系列之零 底层知识:}}\)
多项式的表示:
多项式可以通过系数数列\(A\)表示,\(a_i\)是\(x_i\)的系数。
多项式可以通过点值表示,对于一个\(n\)次多项式,取\(n\)种不同的\(x\)取值带入\(F(x)\),得到\(n\)个值,在取相同这\(n\)个数的意义下,可以唯一的表示这个多项式。
多项式乘法:
定义\(F(x)\oplus G(x)=\sum_{i=0}^n\sum_{j=0}^i f_ix^i\times g_{j-i}x^{j-i}\),在系数表示之下相乘复杂度\(\Theta(n^2)\),在点值表示之下\(F(x)\oplus G(x)=A_f\times A_g=\sum_{i=1}^n a_{fi}\times a_{gi}\),复杂度\(\Theta(n)\)。
复数:
复数一般情况下可以表示成\(a+bi\)的形式,\(a,b\)是实数,\(i=\sqrt{-1}\)。
复数的幅角:平面直角坐标系上点\((a,b)\)所在的任意角。
复数的模长:\(\sqrt{a^2+b^2}\)
两个复数相乘:\((a+bi)\times(c+di)=ac+adi+bci-bd=(ac-bd)+(ad+bc)i\),复数相乘之后,模长等于原来两个复数的模长的乘积,幅角的角度等于原来两个幅角的和。
复数可以加减乘除,可以和实数一样的带入\(F(x)\)。
单位根:
在单位圆上从\(w_m^0=(1,0)\)开始平均取\(m\)个点,从\(0\)开始编号,分别是\(w_m^0,w_m^1,w_m^2,w_m^3\cdots w_m^{m-1}\)。
画图观察可得:
\(1.w_m^k\;=(cos(\frac k m 2\pi),sin(\frac k m 2\pi))\)所代表的复数
\(2.w_m^{-k}=(cos(\frac {-k} n 2\pi),sin(\frac {-k} n 2\pi))\)所代表的复数
\(3.w_m^m\;=w_m^0=(1,0)\)
\(4.w_m^m\;=-w_m^{\frac m 2}\)
\(5.w_{2m}^{2k}\;=w_m^k\)
\(6.w_m^{k+\frac m 2}\!=-w_m^k\)
DFT&IDFT:
在科学的数学函数意义上DFT是讲一个函数转化成三角函数的加减乘除的形式,三角函数的系数是原函数系数与点值之间的变换规律。IDFT是DFT的逆变换。
原根:
(NTT用)
\(1.\)什么是\(g\):在\(mod~p\)意义下\(g^i(i\in[0,p-1])\)互不相同,即\(g\)可以张成整个\(mod~p\)下的域。
\(2.g\)存在的条件:\(p=2,4,q^a,2q^a\),\(q\)是奇素数。
\(3.\)如何求\(g\):把\(\phi(p)\)进行质因数分解\(\phi(p)=\prod p_i^{a_i}\),如果对于任意的\(p_i\),总有\(g^{\frac {\phi(p)} {p_i}}\neq 1(mod~p)\),暴力枚举即可。
CRT合并:
(NTT用)
求解\({\begin{cases}x\equiv a_1 (mod~p_1)\\x\equiv a_2 (mod~p_2)\end{cases}},\gcd(p_1,p_2)=1\)
由\(x\equiv a_1 (mod~p_1)\),得
带入二式,得
若\(\gcd(p1,p2)==1\),用逆元直接除便可;否则通过\(exgcd\)可求得\(y\),若无解则方程组无解。
最后\(x=p_1y+a_1(mod~p_1p_2)\)。
泰勒展开:
(牛顿迭代用)
简单来说,\(F(x)\)在\(x_0\)处展开就是
牛顿迭代:
(开根、exp用)
牛顿迭代是用来求解零点问题的一种方法。
已知\(G(z)\) (这是一个函数,如\(G(z)=z^2-H(x)\)那么我们就是对\(H(x)\)开根) 求\(G(F(x))\equiv 0(mod\;x^n)\)的\(F(x)\)。我们尝试递归求解。
\(1.\)结束位置:
当\(n==1\)时,\(G(F(x))\equiv0(mod\;x)\)直接求解即可,与\(G(x)\)的具体形式有关。
\(2.\)转移:
假设我们已经得知\(G(J(x))\equiv0(mod\;x^{\lceil\frac n 2\rceil})\):
我们运用泰勒展开把\(G(F(x))\)在\(J(x)\)处展开,得
由于\(F(x)\)和\(J(x)\)的前\(\lceil\frac n 2\rceil\)相同,所以\(F(x)-J(x)\)的最低非零项大于\(\lceil\frac n 2\rceil\),所以\((F(x)-J(x))^2\)的最低非零项大于\(2\lceil\frac n 2\rceil\)取模后为\(0\)。
由于\(G(F(x))\equiv 0\;(mod\;x^n):\)
拉格朗日插值法:
(多项式快速插值用)
\(0.\)问题描述:
给定点集\(X=\{(x_0,y_0),(x_1,y_1)\cdots\}\),求一个函数\(F(x)\),使得\(\forall x_i,F(x_i)=y_i\)
\(1.\)朴素的拉格朗日插值法:
构造函数\(L_i(x)=\prod_{i\ne j}\frac{x-x_j}{x_i-x_j}\)。
当\(x=x_i\)时,\(L_i(x)=1\)。
当\(x\in X,x\ne x_i\)时,\(L_i(x)=0\)
这样,\(F(x)=\sum_{i}L_i(x)\times y_i\)
复杂度\(\Theta(n^2)\)
\(2.\)重心拉格朗日插值法:
在朴素拉格朗日插值法中,如果我们新插入一个节点,修改的复杂度是\(\Theta(n^2)\)的,因为每个函数都要修改,而函数的修改要先减去它再修改再加回来,加减的复杂度是\(\Theta(n)\)的。
我们推一下式子:
令\(w_i=\prod_{i\ne j}\frac{1}{x_i-x_j},l(x)=\prod_jx-x_j\)
修改时,我们只需要
\(\color{red}{\text{多项式全集之一 FFT:}}\)
什么是FFT:
FFT是利用DFT的特殊性质,把\(w\)带入\(x\)从而\(\Theta(nlogn)\)求一个系数多项式的点值表示,所以叫FDFT。
\(w\)的具体应用:
\(1.\)可以方便的IDFT:
设\(F(x)\)的系数是\(A\),在\(w_m\)的DFT下点值是\(B\),\(G(x)\)的系数是\(B\),在\(w_m^{-1}\)的DFT下点值是\(C\)。
当\(j-k=0\)时\(\sum_{i=0}^{m-1}w_m^{(j-k)i}=m\),否则根据等比数列求和公式得
由此可得:\(c_k=m\times a_k\),\(a_k=\frac{c_k}{m}\)
综上所述,对于点值取的\(w\)相反数做DFT再除以\(m\)可得到系数。
\(2.\)可以快速的DFT:
直接将\(w\)带入多项式做DFT需要复杂度\(\Theta(n^2)\),我们利用\(w\)的性质优化:
把\(F(x)\)按照奇偶分裂,\(F(x)=(a_0x^0+a_2x^2+a_4x^4+a_6x^6\cdots)+(a_1x^1+a_3x^3+a_5x^5+a_7x^7\cdots)\)
我们令
我们可以发现\(F(x)=F_0(x^2)+xF_1(x^2)\)
现在我们把\(w_m\)带入,令\(k<\frac m 2\):
我们知道取\(w_{\frac m 2}\)时,\(A_0,A_1\)的取值,就可以算出\(A(w_m)\),而\(A_0,A_1\)的长度都为\(A\)的一半,所以可以递归计算。
非递归优化FFT:
\(1.\)优化原理:
画图可知,递归版FFT最底层结束状态第\(i\)个位置的项是\(i\)二进制翻转后的结果。我们可以\(\Theta(n)\)的得到最底层的结果,然后向上模拟回溯合并即可。
\(2.\)蝴蝶变换:
由上述式子:
可得在迭代时\(w_m^k\)与\(w_m^{k+\frac m 2}\)都只与\(w_{\frac m 2}^k,w_{\frac m 2}^{k+\frac m 2}\)有关,所以我们可以用临时变量记录下一层的两个信息向上迭代。
共轭复数优化FFT:
\(1.\)优化原理:
(在DFT时)
令
那么
\(2.\)证明:
\(step~1\):
\(step~2\):为方便起见,我们用\(X\)代替\(\frac{2\pi j k} {m}\)
而在IDFT时,我们需要
数论优化FFT(NTT):
\(g^{\frac {p-1} m}\)与\(w_m\)的共性:
\(1.(g^{\frac {p-1} m})^k\)和\(w_m^k\)都互不相同\((k\in[0,m-1])\)。
\(2.(g^{\frac {p-1} m})^m=g^{p-1}=1\),\(w_m^m=1\)。
\(3.(g^{\frac {p-1} m})^{\frac m 2}=g^{\frac{p-1} 2}=\sqrt{g^{p-1}}\),由于原根的互不相同,\((g^{\frac {p-1} m})^{\frac m 2}=-1=-g^{p-1}=-(g^{\frac {p-1} m})^m\),\(w_m^m=-w_m^{\frac m 2}\)
\(4.(g^{\frac {p-1} {2m}})^{2k}=(g^{\frac {p-1} {m}})^{k}\),\(w_{2m}^{2k}=w_m^k\)
\(5.(g^{\frac {p-1} {m}})^{k+\frac m 2}=(g^{\frac {p-1} {m}})^{k}\times (g^{\frac {p-1} {m}})^{\frac m 2}=-(g^{\frac {p-1} {m}})^{k}\),\(w_m^{k+\frac m 2}=-w_m^k\)
因为有这些共性,所以\(g^{\frac {p-1} m}\)可以代替\(w_m\)。
喜闻乐见的模板:
FFT模板(共轭优化)
namespace FFT{
const double pi = acos(-1);
struct cp{
double x, y;
cp() {x = y = 0;}
cp(double X,double Y) {x = X; y = Y; }
cp conj() {return (cp) {x, -y};}
}a[3000005], b[3000005], c[3000005], I(0, 1), d[3000005];
cp operator+ (const cp &a, const cp &b) {return (cp){a.x + b.x, a.y + b.y}; }
cp operator- (const cp &a, const cp &b) {return (cp){a.x - b.x, a.y - b.y}; }
cp operator* (const cp &a, const cp &b) {return (cp){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}
cp operator* (const cp &a, double b) {return (cp){a.x * b, a.y * b};}
cp operator/ (const cp &a, double b) {return (cp){a.x / b, a.y / b};}
struct p_l_e{
int wz[3000005];
void FFT(cp *a, int N, int op) {
for(int i = 0; i < N; i++)
if (i<wz[i]) swap(a[i],a[wz[i]]);
for(int le = 2; le <= N; le <<= 1) {
int mid = le >> 1;
for(int i = 0; i < N; i += le) {
cp x, y, w = (cp) {1, 0};
cp wn = (cp){cos(op * pi / mid), sin(op * pi / mid)};
for(int j = 0 ; j < mid; j++) {
x = a[i + j]; y = a[i + j + mid] * w;
a[i + j] = x + y;
a[i + j + mid] = x - y;
w = w * wn;
}
}
}
}
void D_FFT(cp *a, cp *b, int N, int op){
for(int i = 0; i < N; i++) d[i] = a[i] + I * b[i];
FFT(d, N, op);
d[N] = d[0];
for(int i = 0; i < N; i++){
a[i] = (d[i] + d[N - i].conj()) / 2;
b[i] = I * (-1) * (d[i] - d[N - i].conj()) / 2;
}
d[N] = cp(0, 0);
}
void mult(cp *a, cp *b, cp *c, int M){
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
D_FFT(a, b, N, 1);
for(int i = 0; i < N; i++) c[i] = a[i] * b[i];
FFT(c, N, -1);
for(int i = 0; i < N; i++) c[i].x = c[i].x / N;
}
}PLE;
int n, m;
void main() {
scanf("%d%d", &n, &m); n++; m++;
for(int i = 0; i < n; i++) scanf("%lf", &a[i].x);
for(int i = 0; i < m; i++) scanf("%lf", &b[i].x);
PLE.mult(a, b, c, n + m - 1);
for(int i = 0; i < n + m - 1; i++) printf("%d ", (int)round(c[i].x));
return;
}
}
NTT模板:
namespace NTT{
typedef long long LL;
const int mod = 998244353;
const int g = 3;
LL a[3000005], b[3000005], c[3000005];
int n, m;
LL qpow(LL a, LL b){
LL ans = 1;
while(b){
if(b & 1) ans = ans * a % mod;
a = a * a % mod;
b >>= 1;
}
return ans;
}
struct p_l_e{
int wz[3000005];
void NTT(LL *a, int N, int op) {
for(int i = 0; i < N; i++)
if(i < wz[i]) swap(a[i], a[wz[i]]);
for(int le = 2; le <= N; le <<= 1) {
int mid = le >> 1;
LL wn = qpow(g, (mod - 1) / le);
if(op == -1) wn = qpow(wn, mod - 2);
for(int i = 0; i < N; i += le) {
int w = 1, x, y;
for(int j = 0; j < mid; j++) {
x = a[i + j];
y = a[i + j + mid] * w % mod;
a[i + j] = (x + y) % mod;
a[i + j + mid] = (x - y + mod) % mod;
w = w * wn % mod;
}
}
}
}
void mult(LL *a, LL *b, LL *c, int M) {
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
NTT(a, N, 1); NTT(b, N, 1);
for(int i = 0; i < N; i++) c[i] = a[i] * b[i] % mod;
NTT(c, N, -1);
LL t = qpow(N, mod - 2);
for(int i = 0; i < N; i++) c[i] = c[i] * t % mod;
}
}PLE;
void main() {
scanf("%d%d", &n, &m); n++; m++;
for(int i = 0; i < n; i++) scanf("%lld", &a[i]);
for(int i = 0; i < m; i++) scanf("%lld", &b[i]);
PLE.mult(a, b, c, n + m - 1);
for(int i = 0; i < n + m - 1; i++) printf("%lld ", c[i]);
}
}
\(\color{red}{\text{多项式全集之二 任长任模的FFT:}}\)
三模NTT现任意模数FFT(MTT):
\(1.\)为什么要用MTT:当\(p\)不是NTT模数或者多项式长度大于模数限制时,就要使用MTT。
\(2.\)MTT的使用原理:我们对初始多项式取模,那么如果在不取模卷积情况下,答案\(x\)不会超过\(N\times p^2\)。我们取三个NTT模数\(p_1,p_2,p_3\),分别做多项式乘法,得到\(x\)分别\(mod~p_1,p_2,p_3\)的答案,通过CRT合并可以得到\(x~mod~p_1p_2p_3\)的答案,如果\(x<p_1p_2p_3\)那么就可以得到准确的答案,再对\(p\)取模即可。
\(3.\)CRT合并的小优化:
\(step~0:\)初始式子
\(step~1:\)把一式二式合并(LL范围内)。
\(step~2:\)再次合并(不需要\(long~double\) 快速乘)。
\(4.\)常用NTT模数:
以下模数的共同\(g=3189\)
\(p=r\times 2^k+1\) | \(k\) | \(g\) |
---|---|---|
\(104857601\) | \(22\) | \(3\) |
\(167772161\) | \(25\) | \(3\) |
\(469762049\) | \(26\) | \(3\) |
\(950009857\) | \(21\) | \(7\) |
\(998244353\) | \(23\) | \(3\) |
\(1004535809\) | \(21\) | \(3\) |
\(2013265921\) | \(27\) | \(31\) |
\(2281701377\) | \(27\) | \(3\) |
\(3221225473\) | \(30\) | \(5\) |
拆系数FFT(CFFT)实现任模FFT:
\(1.\)实现原理:运用实数FFT不取模做乘法,然后取模回归到整数。但是由于误差较大(值域是\(10^{23}\)),我们令\(t=\sqrt{m}\)把系数\(a_i=k_it+b_i\),对\(k_i,t_i\)交叉做四遍卷积,求出答案按系数贡献取模加入。
\(2.\)可按合并DFT的方法优化DFT次数。
\(bluestein\)算法实现任意长长度FFT:
当\(m\)不是\(2\)的幂次的时候,我们从式子入手:
令\(X_i=a_iw_m^{\frac {i^2} 2},Y_i=w_m^{\frac{-i^2}2}\)
喜闻乐见的模板:
三模NTT模板(注意:不可以MTT回来,因为系数会取模)
namespace MTT{
typedef long long LL;
int n, m;
LL p, mod;
const LL p1 = 998244353;
const LL p2 = 1004535809;
const LL p3 = 104857601;
const int g = 3189;
LL a[300005], b[300005], c[300005], cpa[300005], cpb[300005];
LL c3[300005], c1[300005], c2[300005];
LL qpow(LL a, LL b, LL mod) {
LL ans = 1;
while(b) {
if(b & 1) ans = ans * a % mod;
a = a * a % mod;
b >>= 1;
}
return ans;
}
const LL inv12 = qpow(p1, p2 - 2, p2);
const LL inv123 = qpow(p1 * p2 % p3, p3 - 2, p3);
struct p_l_e{
int wz[300005];
void MTT(LL *a, int N, int op) {
for(int i = 0; i < N; i++)
if(i < wz[i]) swap(a[i], a[wz[i]]);
for(int le = 2; le <= N; le <<= 1) {
int mid = le >> 1;
LL wn = qpow(g, (mod - 1) / le, mod);
if(op == -1) wn = qpow(wn, mod - 2, mod);
for(int i = 0; i < N ;i += le) {
LL w = 1, x, y;
for(int j = 0; j < mid; j++) {
x = a[i + j];
y = a[i + j + mid] * w % mod;
a[i + j] = (x + y) % mod;
a[i + j + mid] = (x - y + mod) % mod;
w = w * wn % mod;
}
}
}
}
void mult(LL *a, LL *b, LL *c, int M) {
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
MTT(a, N, 1); MTT(b, N, 1);
for(int i = 0; i < N; i++) c[i] = a[i] * b[i] % mod;
MTT(c, N, -1);
LL t = qpow(N, mod - 2, mod);
for(int i = 0; i < N; i++) c[i] = c[i] * t % mod;
}
}PLE;
LL CRT(LL c1, LL c2, LL c3) {
LL x = (c1 + p1 * ((c2 - c1 + p2) % p2 * inv12 % p2));
LL y = (x % p + p1 * p2 % p * ((c3 - x % p3 + p3) % p3 * inv123 % p3) % p) % p;
return y;
}
void merge(LL *c1, LL *c2, LL *c3, LL *c, int N) {
for(int i = 0; i < N; i++)
c[i] = CRT(c1[i], c2[i], c3[i]);
return;
}
void main() {
scanf("%d%d%lld", &n, &m, &p); n++; m++;
for(int i = 0; i < n; i++) scanf("%lld", &a[i]);
for(int i = 0; i < m; i++) scanf("%lld", &b[i]);
mod = p1; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c1, n + m - 1);
mod = p2; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c2, n + m - 1);
mod = p3; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c3, n + m - 1);
merge(c1, c2, c3, c, n + m - 1);
for(int i = 0; i < n + m - 1; i++) printf("%lld ", (c[i] % p + p) % p);
return;
}
}
拆系数FFT模板(注意:相同系数的两项可以合并一起IDFT。采用共轭优化法,只进行四次DFT)
namespace CFFT{
typedef long long LL;
int n, m, p ,sqrp;
int a[300005], b[300005];
const long double pi = acos(-1);
struct cp{
long double x, y;
cp() {x = y = 0;}
cp(long double X,long double Y) {x = X; y = Y; }
cp conj() {return (cp) {x, -y};}
}ka[300005], kb[300005], ta[300005], tb[300005], kk[300005], kt[300005], tt[300005], c[300005], I(0, 1), d[300005];
cp operator+ (const cp &a, const cp &b) {return (cp){a.x + b.x, a.y + b.y}; }
cp operator- (const cp &a, const cp &b) {return (cp){a.x - b.x, a.y - b.y}; }
cp operator* (const cp &a, const cp &b) {return (cp){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}
cp operator* (const cp &a, long double b) {return (cp){a.x * b, a.y * b};}
cp operator/ (const cp &a, long double b) {return (cp){a.x / b, a.y / b};}
struct p_l_e{
int wz[300005];
void FFT(cp *a, int N, int op){
for(int i = 0; i < N; i++)
if(i < wz[i]) swap(a[i], a[wz[i]]);
for(int le = 2; le <= N; le <<= 1){
int mid = le >> 1;
cp x, y, w, wn = (cp){cos(op * 2 * pi / le), sin(op * 2 * pi / le)};
for(int i = 0; i < N; i += le){
w = (cp){1, 0};
for(int j = 0; j < mid; j++){
x = a[i + j];
y = a[i + j + mid] * w;
a[i + j] = x + y;
a[i + j + mid] = x - y;
w = w * wn;
}
}
}
}
void D_FFT(cp *a, cp *b, int N, int op){
for(int i = 0; i < N; i++) d[i] = a[i] + I * b[i];
FFT(d, N, op);
d[N] = d[0];
if(op == 1){
for(int i = 0; i < N; i++){
a[i] = (d[i] + d[N - i].conj()) / 2;
b[i] = I * (-1) * (d[i] - d[N - i].conj()) / 2;
}
} else {
for(int i = 0; i < N; i++){
a[i] = cp(d[i].x, 0);
b[i] = cp(d[i].y, 0);
}
}
d[N] = cp(0, 0);
}
void mult(int *a, int *b, int M){
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
for(int i = 0; i < N; i++){
ka[i].x = a[i] >> 15;
kb[i].x = b[i] >> 15;
ta[i].x = a[i] & 32767;
tb[i].x = b[i] & 32767;
}
D_FFT(ta, ka, N, 1); D_FFT(tb, kb, N, 1);
for(int i = 0; i < N; i++){
kk[i] = ka[i] * kb[i];
kt[i] = ka[i] * tb[i] + ta[i] * kb[i];
tt[i] = ta[i] * tb[i];
}
D_FFT(tt, kk, N, -1); FFT(kt, N, -1);
for(int i = 0; i < N; i++){
tt[i] = tt[i] / N;
kt[i] = kt[i] / N;
kk[i] = kk[i] / N;
}
}
}PLE;
void main() {
scanf("%d%d%d", &n, &m, &p); n++; m++;
for(int i = 0; i < n; i++) scanf("%d", &a[i]),a[i] = a[i] % p;
for(int i = 0; i < m; i++) scanf("%d", &b[i]),b[i] = b[i] % p;
PLE.mult(a, b, n + m - 1);
for(int i = 0; i < n + m - 1; i++)
printf("%lld ",(((((LL)round(kk[i].x)) % p) << 30) + ((((LL)round(kt[i].x)) % p) << 15) + ((LL)round(tt[i].x)) % p) % p);
}
}
\(blue\_stein\)模板:
struct polynie {
CP getw(int m, int k, int op) {
return CP(cos(2 * pi * k / m), op * sin(2 * pi * k / m));
}
int wz[MAXN];
CP A[MAXN], B[MAXN], C[MAXN];
void FFT(CP *a, int N, int op) {
rop(i, 0, N) if(i < wz[i]) swap(a[i], a[wz[i]]);
for(int l = 2; l <= N; l <<= 1) {
int mid = l >> 1;
CP x, y, w, wn = CP(cos(pi / mid), sin(op * pi / mid));
for(int i = 0; i < N; i += l) {
w = CP(1, 0);
rop(j, 0, mid) {
x = a[i + j];
y = w * a[i + j + mid];
a[i + j] = x + y;
a[i + j + mid] = x - y;
w = w * wn;
}
}
}
}
void mult(CP *a, CP *b, CP *c, int M) {
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
FFT(a, N, 1); FFT(b, N, 1);
rop(i, 0, N) c[i] = a[i] * b[i];
FFT(c, N, -1);
rop(i, 0, N) c[i].x = c[i].x / N, c[i].y = c[i].y / N;
}
void blue_stein(CP *a, int M, int op) {
int M2 = M << 1;
memset(A, 0, sizeof(A));
memset(B, 0, sizeof(B));
memset(C, 0, sizeof(C));
rop(i, 0, M) A[i] = a[i] * getw(M2, 1ll * i * i % M2, op);
rop(i, 0, M2) B[i] = getw(M2, 1ll * (i - M) * (i - M) % M2, -op);
mult(A, B, C, M2 + M - 1);
rop(i, 0, M) a[i] = C[i + M] * getw(M2, 1ll * i * i % M2, op);
if(op == -1) rop(i, 0, M) a[i].x = a[i].x / M, a[i].y = a[i].y / M;
}
}PLE;
\(\color{red}{\text{多项式全集之三 多项式求逆与除法:}}\)
多项式求逆:
\(1.\)问题描述:
已知\(A(x)\),且\(F(x)A(x)\equiv 1 (mod~x^n)\),求\(F(x)\)
\(2.\)推导过程:
两边平方,得:
由于\([G(x)-B(x)]^2\)的第\(k<n\)项为
\(i,k-i\)一定有一项\(<\frac n 2\),所以
两边同乘\(A(x)\),得:
多项式除法:
\(1.\)问题描述:
已知一个\(n\)次多项式\(A(x)\),一个\(m\)次多项式\(B(x)\),且\(A(x)=B(x)C(x)+D(x)\),求\(n-m\)次多项式\(C(x)\),\(<m\)次多项式\(D(x)\)。
\(2.\)推导过程:
由\(A(x) = B(x)C(x)+D(x)\)得:
求逆可得\(B_r(x)\),再反转得\(B(x)\),然后乘\(C(x)\)去减\(A(x)\)得\(D(x)\).
喜闻乐见的模板:
namespace INV{
typedef long long LL;
int n, a[300005], b[300005];
const int mod = 998244353;
const int g = 3189;
int qpow(int a, int b){
int ans = 1;
while(b){
if(b & 1) ans = 1ll * ans * a % mod;
a = 1ll * a * a % mod;
b >>= 1;
}
return ans;
}
struct p_l_e{
int wz[300005], i_c[300005];
void NTT(int *a, int N, int op){
for(int i = 0; i < N; i++)
if(i < wz[i]) swap(a[i], a[wz[i]]);
for(int le = 2; le <= N; le <<= 1){
int mid = le >> 1, wn = qpow(g, (mod - 1) / le);
if(op == -1) wn = qpow(wn, mod - 2);
for(int i = 0; i < N; i += le){
LL w = 1; int x, y;
for(int j = 0; j < mid; j++){
x = a[i + j];
y = w * a[i + j + mid] % mod;
a[i + j] = (x + y) % mod;
a[i + j + mid] = (x - y + mod) % mod;
w = w * wn % mod;
}
}
}
}
int init(int M){
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
return N;
}
void INV(int *a, int *b, int deg){
if(deg == 1){b[0] = qpow(a[0], mod - 2); return;}
INV(a, b, (deg + 1) >> 1);
int N = init(deg + deg - 1);
for(int i = 0; i < deg; i++) i_c[i] = a[i];
for(int i = deg; i < N; i++) i_c[i] = 0;
NTT(b, N, 1);NTT(i_c, N, 1);
for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * (2 - 1ll * b[i] * i_c[i] % mod + mod) % mod;
NTT(b, N, -1);
int t = qpow(N, mod - 2);
for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * t % mod;
for(int i = deg; i < N; i++) b[i] = 0;
}
}PLE;
void main(){
scanf("%d", &n);
for(int i = 0; i < n; i++) scanf("%d", &a[i]);
PLE.INV(a, b, n);
for(int i = 0; i < n; i++) printf("%d ",b[i]);
}
}
\(\color{red}{\text{多项式全集之四 多项式ln与exp:}}\)
多项式ln:
\(1.\)做法:
设
两边求导得
积分回去即可。
\(2.\)应用:
这个的组合意义是:无序组合。
设\(F(x)\),\(f_i\)表示一些东西,那么这些东西有序组合的方案数为
而无序组成的方案数为:
如果无序组合方案数好求,那么求\(\ln\)就能得到\(F(x)\)。
多项式\(\exp\):
我们就是求\(F(x)=e^{A(x)}\),变形一下得\(\ln(F(x))-A(x)=0\),就是求函数\(G(z)=\ln(z)-A(x)\)的零点,带入牛顿迭代的公式得:
由于\(G(F(x))\equiv 0 (mod\;x^n)\)
把\(F(x)\)看成自变量\(x,A(x)\)看成常数来求导
对于终止情况:
\(\ln F(x)-A(x)=0(mod\;x)\)
因为都只剩下了常数项,
\(F(x)=e^{A(x)}\)
喜闻乐见的代码:
多项式\(\ln\):
namespace PLE_ln{
struct polyme {
int li[SZ], wz[SZ];
void NTT(int *a, int N, int op) {
rop(i, 0, N) if(i < wz[i]) swap(a[i], a[wz[i]]);
for(int l = 2; l <= N; l <<= 1) {
int mid = l >> 1;
int x, y, w, wn = qpow(g, (mod - 1) / l);
if(op) wn = qpow(wn, mod - 2);
for(int i = 0; i < N; i += l) {
w = 1;
for(int j = 0; j < mid; ++j) {
x = a[i + j]; y = 1ll * w * a[i + j + mid] % mod;
a[i + j] = (x + y) % mod;
a[i + j + mid] = (x - y + mod) % mod;
w = 1ll * w * wn % mod;
}
}
}
}
void qd(int *a, int *b, int n) {
rop(i, 0, n) b[i] = 1ll * a[i + 1] * (i + 1) % mod;
}
void jf(int *a, int *b, int n) {
rop(i, 1, n) b[i] = 1ll * a[i - 1] * qpow(i, mod - 2) % mod;
}
void mult(int *a, int *b, int *c, int M) {
int N = 1, len = 0;
while(N < M) N <<= 1, len ++;
rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
NTT(a, N, 0); NTT(b, N, 0);
rop(i, 0, N) c[i] = 1ll * a[i] * b[i] % mod;
NTT(c, N, 1);
int t = qpow(N, mod - 2);
rop(i, 0, N) c[i] = 1ll * c[i] * t % mod;
}
void inv(int *a, int *b, int deg) {
if(deg == 1) {b[0] = qpow(a[0], mod - 2) % mod; return;}
inv(a, b, (deg + 1) >> 1);
rop(i, 0, deg) li[i] = a[i];
int N = 1, len = 0;
while(N < deg + deg - 1) N <<= 1, len ++;
rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
rop(i, deg, N) li[i] = 0;
NTT(li, N, 0); NTT(b, N, 0);
rop(i, 0, N) b[i] = 1ll * b[i] * (2 - 1ll * li[i] * b[i] % mod + mod) % mod;
NTT(b, N, 1);
int t = qpow(N, mod - 2);
for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * t % mod;
rop(i, deg, N) b[i] = 0;
}
}PLE;
int a[SZ], da[SZ], ia[SZ], dla[SZ], la[SZ], n;
void main() {
scanf("%d", &n);
rop(i, 0, n) scanf("%d", &a[i]);
PLE.qd(a, da, n);
PLE.inv(a, ia, n);
PLE.mult(ia, da, dla, n + n - 1);
PLE.jf(dla, la, n);
rop(i, 0, n) printf("%d ", la[i]);
}
}
\(\color{red}{\text{多项式全集之五 多项式快速幂与开根:}}\)
多项式快速幂方法一:
直接套用快速幂模板,重载*运算符,\(O(n\log n\log W)\)。
多项式快速幂方法二:
\(F(x)^k=e^{\ln F^k(x)}=e^{k\ln F(x)}\)多项式\(ln\)再\(exp\)即可。
多项式开根:
已知\(F^2(x)\equiv A(x)\;\;\;(mod\;x^n)\),求\(F(x)\)。
移项得:\(F^2(x)-A(x)\equiv 0\;\;\;(mod\;x^n)\)
问题转化成求函数\(G(z)=z^2-A(x)\equiv 0\;\;\;(mod\;x^n)\)的零点,带入牛顿迭代得:
多项式求逆即可。
结束状态:\(F(x)=\sqrt{A_0}\)二次剩余一下。
\(\color{red}{\text{多项式全集之六 多项式多点求值与快速插值:}}\)
多项式多点求值 :
给定多项式\(F(x)\),\(X=\{{x_0, x_1, x_2\cdots x_m}\}\),求\(F(x_0),F(x_1)\cdots\)。
设:
显然当\(x\in X_0\)时,\(P_0(x)=0\),\(P_1\)同理。
我们设
那么
当\(x\in X_0\)时,\(F(x)=A_0(x)\),我们求出\(A_0(x)\)的值就可以得到\(F(x)\)的值了,而且\(A_0(x)\)长度小于\(P_0(x)\)长度(\(\lfloor\frac n 2\rfloor\))。\(P_1,A_1\)同理,这样我们把\(P\)预处理一下,就可以递归解决了。
多项式快速插值 :
\(\color{red}{\text{多项式全集之八 本文提到知识的部分扩展:}}\)
麦克劳林级数:
泰勒公式中,取\(x_0=0\)得到的式子叫麦克劳林级数: