有关多项式处理的各种算法总结
前言
最近频繁接触了各种各样的生成函数,以及各种十(sang)分(xin)有(bing)趣(kuang)的多项式相关算法。
由于这部分内容太过于丧病外加内容较多,因此决定写篇总结来辅助复习。
阅读本文所需前置技能:快速数论变换
以下内容中代码片涉及的多项式卷积(NTT)及各种函数使用代码如下:
typedef long long ll;
const ll md=998244353;
inline int add(int a,int b){a+=b;if(a>=md)return a-md;return a;}
inline int mul(int a,int b){return (ll)a*b%md;}
inline void initrev(int n)
{
for(int i=0;i<n;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)*(n>>1));
}
inline void NTT(int *a,int n,bool f)
{
for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int h=2;h<=n;h<<=1)
{
int w=qpow(3,(md-1)/h);
if(f)w=qpow(w,md-2);
for(int j=0;j<n;j+=h)
{
int wn=1;
for(int k=j;k<j+(h>>1);k++)
{
int x=a[k],y=mul(wn,a[k+(h>>1)]);
a[k]=add(x,y);
a[k+(h>>1)]=add(x-y,md);
wn=mul(wn,w);
}
}
}
if(f)
for(int i=0,invn=inv[n];i<n;i++)
a[i]=mul(a[i],invn);
}
inline void mult(int *a,int n,int *b,int m,int *c)
{
static int x[N],y[N],l;
for(l=1;l<=(n+m-1);l<<=1);
initrev(l);
memcpy(x,a,sizeof(a[0])*n);
memcpy(y,b,sizeof(b[0])*m);
fill(x+n,x+l,0);fill(y+m,y+l,0);
NTT(x,l,0);NTT(y,l,0);
for(int i=0;i<l;i++)
c[i]=mul(x[i],y[i]);
NTT(c,l,1);
}
同时,绝大多数代码中传入的参数 n n 为多项式长度,要求为一个不小于所需答案长度的一个的次幂。
1.一些定义的简要介绍
非零多项式 G(x) G ( x ) 的次数为其最高项的次数,通常记为 deg(G(x)) d e g ( G ( x ) ) 。本文本着尽量直观的原则将避免使用此符号。
形如 [xn]G(x) [ x n ] G ( x ) 的符号代表多项式 G(x) G ( x ) 的 xn x n 一项的系数。
多项式之间的乘积叫做卷积,可以理解为初中学习的因式分解后展开的方法。
举个栗子就明白了:
多项式通常以生成函数的形式出现。这里简单介绍一下生成函数。
所谓生成函数,就是把一个数列用一个函数的各项系数来表示……
比如,这里以著名的斐波那契数列做例子:
假设
fib(1)=1,fib(2)=1,fib(n)=fib(n−1)+fib(n−2)
f
i
b
(
1
)
=
1
,
f
i
b
(
2
)
=
1
,
f
i
b
(
n
)
=
f
i
b
(
n
−
1
)
+
f
i
b
(
n
−
2
)
得到数列:
令其生成函数为 F(x) F ( x ) ,得到的函数如下:
相信大家这就大致能明白生成函数是个什么意思了~
对于一个无穷长的数列的生成函数,理论上会是个无穷长的多项式,而因此难以表示出来。
因此定义
mod xn
m
o
d
x
n
意义下的多项式为仅保留前
n
n
项的的原多项式。
以上关于生成函数的内容和剩下的内容貌似没有什么太大关系。
2.牛顿迭代
问题:
给出多项式,求解多项式
F(x)
F
(
x
)
,满足:
需要给出 F(x) mod xn F ( x ) m o d x n 意义下的答案,即前 n n 个系数。
解决方法(没看懂或不想看可直接到下面记忆结论,下面的内容同理):
假设咱们现在有
要求的是
也就是利用倍增的思想,将长度加倍。
注意如果求解的是特定长度的答案,可以考虑直接求出向上取最近的 2 2 的次幂长度的答案,并保留所需的长度即可。
对在
F0(x)
F
0
(
x
)
处来一发泰勒展开:
可以发现,在 mod x2t+1 m o d x 2 t + 1 意义下, F(x)−F0(x) F ( x ) − F 0 ( x ) 的最低非 0 0 次项大于,于是从第三项开始,后面的项贡献均为 0 0 。
那么得到一个简单的式子:
又由咱们的目标知 G(F(x))≡0 (mod x2t+1) G ( F ( x ) ) ≡ 0 ( m o d x 2 t + 1 )
那么最终移项得到:
这个形式正是咱们想要的~
3.多项式求逆
提出给定多项式
F(x)
F
(
x
)
,求多项式
G(x)
G
(
x
)
,使得下式成立:
解决方法:
考虑使用倍增,边界为
G0=Inv(F0)
G
0
=
I
n
v
(
F
0
)
,即直接求逆元。
考虑使用牛顿迭代公式,令
G(x)
G
(
x
)
为变量,
Gt(x)
G
t
(
x
)
代表
mod x2t
m
o
d
x
2
t
意义下的
G(x)
G
(
x
)
。
假设已经有了
Gt(x)
G
t
(
x
)
,要求
Gt+1(x)
G
t
+
1
(
x
)
,将
Gt(x)∗F(x)−1
G
t
(
x
)
∗
F
(
x
)
−
1
带入牛顿迭代公式,得:
最后的式子便是结果!
代码实现(个人习惯较为明显,仅供参考):
inline void cinv(int *a,int *b,int n)
{
static int C[N];
if(n==0)
{
b[0]=qpow(a[0],md-2);
return;
}
cinv(a,b,n>>1);
memcpy(C,a,sizeof(a[0])*n);
memset(C+n,0,sizeof(a[0])*n);
n<<=1;initrev(n);
NTT(C,n,0);NTT(b,n,0);
for(int i=0;i<n;i++)
b[i]=add(mul(2,b[i])-mul(C[i],mul(b[i],b[i])),md);
NTT(b,n,1);n>>=1;
memset(b+n,0,sizeof(b[0])*n);
}
4.多项式求 ln ln
问题:
给出多项式
F(x)
F
(
x
)
,求多项式
G(x)
G
(
x
)
,满足
解决方法:
考虑对两边同时求导:
那么有
考虑多项式如何求导……
每一项同时左移一位并乘上原来下标即可~
积分同理,右移一位并把下标乘回去即可~
配合多项式求逆食用,复杂度为 O(nlogn) O ( n l o g n ) ~
代码实现:
inline void cln(int *f,int *g,int n)
{
static int D[N],A[N],m;
memset(D,0,sizeof(D));
memset(A,0,sizeof(A));
for(int i=0;i<n;i++)
D[i]=mul(f[i+1],(i+1));
D[n]=0;
for(m=1;m<=n;m<<=1);
cinv(f,A,n);initrev(m);
NTT(A,m,0);NTT(D,m,0);
for(int i=0;i<m;i++)
A[i]=mul(A[i],D[i]);
NTT(A,m,1);
for(int i=1;i<=n;i++)
g[i]=mul(A[i-1],inv[i]);
g[0]=0;
}
5.多项式开方
问题:
给出多项式
F(x)
F
(
x
)
,求多项式
G(x)
G
(
x
)
,满足
解决方法:
考虑和多项式求逆一样使用倍增,初值
G0=F(0)‾‾‾‾√ (mod x)
G
0
=
F
(
0
)
(
m
o
d
x
)
。
同样将
G2(x)−F(x)
G
2
(
x
)
−
F
(
x
)
带入牛顿迭代公式:
配合多项式求逆食用,分析可知复杂度为 O(nlogn) O ( n l o g n )
代码实现:
void csqrt(int *a,int *b,int n)
{
static int C[N],D[N];
if(n==1)
{
b[0]=sqrt(a[0]);
return;
}
csqrt(a,b,n>>1);
memset(D,0,sizeof(D));
cinv(b,D,n);
memcpy(C,a,sizeof(a[0])*n);
memset(C+n,0,sizeof(C[0])*n);
n<<=1;NTT(C,n,0);NTT(b,n,0);NTT(D,n,0);
for(int i=0;i<n;i++)
b[i]=mul(inv2,add(b[i],mul(C[i],D[i])));
NTT(b,n,1);n>>=1;
memset(b+n,0,sizeof(b[0])*n);
}
6.多项式求指数函数(exp)
问题:
给出多项式
F(x)
F
(
x
)
,求多项式
G(x)
G
(
x
)
,满足
解决方法:
依旧考虑使用倍增和牛顿迭代。
然而发现一个问题:对于指数函数,无法知道常数项。
那么考虑通过同时取
ln
ln
来变形:
然后把这个式子带入牛顿迭代式:
配合多项式求逆和多项式求 ln ln 共同食用,复杂度仍为 O(nlogn) O ( n l o g n ) 。
代码实现:
inline void cexp(int *a,int *b,int n)
{
static int D[N];
if(n==1)
{
b[0]=1;
return;
}
cexp(a,b,n>>1);
memset(D,0,sizeof(D));
cln(b,D,n);
for(int i=0;i<n;i++)
D[i]=add(a[i]-D[i],md);
D[0]=add(1,D[0]);
n<<=1;initrev(n);
NTT(D,n,0);NTT(b,n,0);
for(int i=0;i<n;i++)
b[i]=mul(D[i],b[i]);
NTT(b,n,1);n>>=1;
memset(b+n,0,sizeof(b[0])*n);
}
7.多项式快速幂
问题:
给出多项式
F(x)
F
(
x
)
和正整数
k
k
,求多项式,满足
解决方法:
对等式求
ln
ln
得:
因此
配合多项式求 ln ln 和多项式求 exp e x p 共同食用,复杂度依旧为 O(nlogn) O ( n l o g n ) 。
代码实现:
inline void cpow(int *f,int *g,int n,int k)
{
static int C[N];
memset(C,0,sizeof(C));
cln(f,C,m);
for(int i=0;i<n;i++)
C[i]=mul(C[i],k);
cexp(C,g,m);
}
8.拉格朗日反演
问题:
给出多项式
F(x)
F
(
x
)
,求多项式
G(x)
G
(
x
)
的
n
n
次项系数,满足
解决方法:
采用拉格朗日反演:
本蒟蒻不会证明0_0
但是看完了上面的内容后显然还是会用的~
配合多项式求逆和多项式快速幂食用即可。
代码实现:
inline int lagrange(int *f,int n,int x)
{
static int C[N],D[N];
memset(C,0,sizeof(C));
memset(D,0,sizeof(D));
cinv(f,C,n);
cpow(C,D,n,x);
return mul(D[x-1],qpow(x,md-2));
}
9.多项式除法和取模
问题:
给出
n
n
次多项式和
m
m
次多项式,求多项式
G(x)
G
(
x
)
和
P(x)
P
(
x
)
,满足
且
解决方法:
定义运算
AR(x)
A
R
(
x
)
,满足对于
deg(A(x))=n
d
e
g
(
A
(
x
)
)
=
n
的多项式
A(x)
A
(
x
)
:
或者说
再举个例子
也就是系数反转的意思~
那么现在来考虑将原多项式反转一下:
注意此处默认 deg(G(x))=n−m d e g ( G ( x ) ) = n − m , deg(R(x))=m−1 d e g ( R ( x ) ) = m − 1 ,高次项补 0 0 。
此时,最低次项为
xn−m+1xm−1(1x)m−1=xn−m+1
x
n
−
m
+
1
x
m
−
1
(
1
x
)
m
−
1
=
x
n
−
m
+
1
,而
deg(G(x))=n−m
d
e
g
(
G
(
x
)
)
=
n
−
m
。
则放到
mod n−m+1
m
o
d
n
−
m
+
1
意义下,
R(x)
R
(
x
)
的影响被完全抵消,而不会影响
G(x)
G
(
x
)
的结果,更不会影响已知的
F(x)
F
(
x
)
和
D(x)
D
(
x
)
。
那么有
于是来一发多项式求逆+reverse函数即可求出 G(x) G ( x ) ,进而求出 R(x) R ( x ) 。复杂度 O(nlogn) O ( n l o g n ) 。
代码实现:
inline void cdiv(int *a,int n,int *b,int m,int *ret)
{
static int c[N],d[N],e[N],l;
if(n<m)return;
for(l=1;l<=n-m+1;l<<=1);
reverse_copy(a,a+n,c);
reverse_copy(b,b+m,d);
for(int i=m;i<n+m+1;i++)//不能用fill(),否则RE!!!
d[i]=0;
memset(e,0,sizeof(e[0])*(l*2));
cinv(d,e,l);
mult(e,n-m+1,c,n-m+1,d);
reverse_copy(d,d+n-m+1,ret);
}
inline void cmod(int *a,int n,int *b,int m,int *ret)
{
static int c[N];
cdiv(a,n,b,m,c);
mult(c,n-m+1,b,m,c);
for(int i=0;i<m-1;i++)
ret[i]=add(a[i]-c[i],md);
}
10.多项式的多点求值
问题:
给出多项式
F(x)
F
(
x
)
和点集
X={x0,x1,⋯,xn−1}
X
=
{
x
0
,
x
1
,
⋯
,
x
n
−
1
}
,求
解决办法:
考虑分治。
如何把问题变成两个规模减半的子问题?
首先把所求的值划成两半:
此时,答案集合分别为 Y0 Y 0 和 Y1 Y 1 。
然后系数也要划成两半。
考虑构造两个多项式:
可以发现,当带入的点 xi∈X0 x i ∈ X 0 ,那么 G1(xi)=0 G 1 ( x i ) = 0 。同理当 xi∈X1 x i ∈ X 1 ,那么 G0(xi)=0 G 0 ( x i ) = 0 。
令点集
X0
X
0
在知道
Y0
Y
0
的情况下可以解出的多项式为
F0(x)
F
0
(
x
)
,
X1
X
1
在知道
Y1
Y
1
的情况下可以解出的多项式为
F1(x)
F
1
(
x
)
。
可以发现,根据定义,
F0(x)
F
0
(
x
)
和
F1(x)
F
1
(
x
)
正是咱们需要的子问题系数,因为它们分别满足了带入
X0
X
0
和
X1
X
1
能解出正确的
Y0
Y
0
和
Y1
Y
1
。
现在可以用如下方式表示
F(x)
F
(
x
)
:
其中 D0(x) D 0 ( x ) 和 D1(x) D 1 ( x ) 为任意多项式。
可以发现在上式中,若带入的值 xi∈X0 x i ∈ X 0 ,则 D0(x)∗G0(x)=0 D 0 ( x ) ∗ G 0 ( x ) = 0 ,此时 F(x)=F0(x) F ( x ) = F 0 ( x ) 。
同理,若带入的值 xi∈X1 x i ∈ X 1 ,则 D1(x)∗G1(x)=0 D 1 ( x ) ∗ G 1 ( x ) = 0 ,此时 F(x)=F1(x) F ( x ) = F 1 ( x ) 。
也就是说,如下关系成立:
这样得到的子问题系数
F0(x)
F
0
(
x
)
和
F1(x)
F
1
(
x
)
满足次数不超过递归区间长度,同时可以求出正确的结果,因此这样做是正确的。
最后递归到底层时,余下来的系数便是答案!
配合多项式取模食用,复杂度
O(nlog2n)
O
(
n
l
o
g
2
n
)
。
注意,使用前需要先使用分治预处理出每一层分治区间的
G(x)
G
(
x
)
,预处理复杂度同样为
O(nlog2n)
O
(
n
l
o
g
2
n
)
。
代码实现:
inline void pre(int dep,int l,int r,int *a,int b[M][N])
{
if(l==r-1)
{
b[dep][l<<1]=add(md,-a[i]),b[dep][l<<1|1]=1;
return;
}
int mid=l+r>>1;
pre(dep+1,l,mid,a,b);
pre(dep+1,mid,r,a,b);
mult(b[dep+1]+(l<<1),mid-l+1,b[dep+1]+(mid<<1),r-mid+1,b[dep]+(l<<1));
}
inline void calc(int dep,int l,int r,int *a,int b[M][N])
{
static int c[N];
if(l==r-1)return;int mid=l+r>>1;
memcpy(c,a+l,sizeof(a[0])*(r-l));
cmod(c,r-l,b[dep+1]+(l<<1),mid-l+1,a+l);
cmod(c,r-l,b[dep+1]+(mid<<1),r-mid+1,a+mid);
calc(dep+1,l,mid,a,b);
calc(dep+1,mid,r,a,b);
}
不出意外还是未完待续…..
(因为还没学快速插值)
写到这里才惊讶地发现咱最近都学的是些什么毒瘤东西……
不过有用就是坠吼的~
一些比较好的练手题
BZOJ 4555 求和 (多项式求逆)
BZOJ 3456 城市规划 (多项式求
ln
ln
)
BZOJ 3625 小朋友和二叉树 (多项式开方)
BZOJ 3684 大朋友和多叉树 (拉格朗日反演)标程
COGS 2189 帕秋莉的超级多项式 (使用到的方法为第2-7条)标程