多项式
由若干个单项式相加组成的代数式叫做多项式
形如:
f(x)=∑ni=0aixi
f
(
x
)
=
∑
i
=
0
n
a
i
x
i
,
deg f(x)
deg
f
(
x
)
称为
f
f
的度,是最高次项的次数。
生成函数
形如
∑∞i=0aixi
∑
i
=
0
∞
a
i
x
i
生成函数又称母函数,往往和多项式算法联系起来起到优化转移的作用
多项式算法
加减法
多项式加减法比较简单
若
f(x)=∑ni=0aixi
f
(
x
)
=
∑
i
=
0
n
a
i
x
i
,
g(x)=∑ni=0bixi
g
(
x
)
=
∑
i
=
0
n
b
i
x
i
则 (f±g)(x)=∑ni=0(ai±bi)xi ( f ± g ) ( x ) = ∑ i = 0 n ( a i ± b i ) x i
代码实现也比较简单
乘法
多项式乘法是所有多项式算法的基础,也是生成函数的卷积运算
若 f(x)=∑ni=0aixi f ( x ) = ∑ i = 0 n a i x i , g(x)=∑ni=0bixi g ( x ) = ∑ i = 0 n b i x i
则 (f×g)(x)=∑2ni=0∑ij=0aj×bi−jxi ( f × g ) ( x ) = ∑ i = 0 2 n ∑ j = 0 i a j × b i − j x i
朴素的多项式乘法是 O(n2) O ( n 2 ) 的,但FFT(快速傅里叶变换)运用复数根的特性可以在 O(nlogn) O ( n log n ) 的复杂度内实现多项式的系数表达和点值表达的转化,运用点值表达实现 O(n) O ( n ) 的相乘,总复杂度就是 O(nlogn) O ( n log n ) ,但是因为常数过大,往往小范围会使用暴力。
inline void FFT(E *a,int r){
for(int i=0;i<n;i++) if(rev[i]>i) swap(a[rev[i]],a[i]);
for(int i=1;i<n;i<<=1){
E wn(cos(M_PI/i),r*sin(M_PI/i));
for(int j=0;j<n;j+=(i<<1)){
E w(1,0);
for(int k=0;k<i;k++,w=w*wn){
E x=a[j+k],y=w*a[j+k+i];
a[j+k]=x+y; a[j+k+i]=x-y;
}
}
}
if(r==-1) for(int i=0;i<n;i++) a[i].real/=n;
}
运用生成函数的计数问题往往会对一个质数取模,如果这个质数可以表示成 k×2n+1 k × 2 n + 1 ,其中 k k 为奇数大于多项式的度数,那么就可以用这个质数的原根代替复数根,实现多项式的系数对一个质数取模,这个算法就是NTT(快速数论变换),这种模数称为NTT模数。
inline void Pre(int n){
num=n;
int g=Pow(3,(P-1)/num);
w[0][0]=w[1][0]=1; for(int i=1;i<num;i++) w[0][i]=1LL*w[0][i-1]*g%P;
for(int i=1;i<num;i++) w[1][i]=w[0][num-i];
}
inline void NTT(int *a,int n,int r){
for(int i=1;i<n;i++) if(rev[i]>i) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1)
for(int j=0;j<n;j+=i<<1)
for(int k=0;k<i;k++){
int x=a[j+k],y=1LL*a[j+k+i]*w[r][num/(i<<1)*k]%P;
a[j+k]=(x+y)%P; a[j+k+i]=(x+P-y)%P;
}
if(!r) for(int i=0,inv=Pow(n,P-2);i<n;i++) a[i]=1LL*a[i]*inv%P;
}
不过有些丧心病狂的题的模数不是NTT模数,这个时候需要CRT合并。
多项式求逆
比如BZOJ4555
最后会得到 f(x)=f(x)g(x)+1 f ( x ) = f ( x ) g ( x ) + 1 ,这是一个非齐次线性递推,因为求的第n项不大,只要直接求出 f(x) f ( x ) 就行了
很容易可以得到
f(x)=11−g(x)
f
(
x
)
=
1
1
−
g
(
x
)
,如果令
h(x)=1−g(x)
h
(
x
)
=
1
−
g
(
x
)
,那么
f(x)=(h(x))−1
f
(
x
)
=
(
h
(
x
)
)
−
1
,所以就要求
h(x)
h
(
x
)
的逆元
求逆元的话网上讲解也有很多,鏼爷集训队论文里也有,大概是用了倍增的做法,摘一下
求
A(x)
A
(
x
)
的逆元,即是求
B(x)
B
(
x
)
,使
A(x)B(x)≡1(modx2k)
A
(
x
)
B
(
x
)
≡
1
(
mod
x
2
k
)
如果我们已经求出
G(x)
G
(
x
)
,
A(x)G(x)≡1(modxn)
A
(
x
)
G
(
x
)
≡
1
(
mod
x
n
)
那么
n=0 n = 0 时,求一下 A(x) A ( x ) 常数项的乘法逆元就可以了
那么用
O(nlogn)
O
(
n
log
n
)
的复杂度把
B(x)
B
(
x
)
赋成
2G(x)−A(x)G2(x)
2
G
(
x
)
−
A
(
x
)
G
2
(
x
)
进入下一层运算就可以了
总复杂度
T(n)=T(n/2)+O(nlogn)=O(nlogn)
T
(
n
)
=
T
(
n
/
2
)
+
O
(
n
log
n
)
=
O
(
n
log
n
)
void Inv(int *a,int *b,int n){
if(n==1){
b[0]=Pow(a[0],P-2); return ;
}
Inv(a,b,n>>1);
for(int i=0;i<n;i++) tmp[i]=a[i],tmp[i+n]=0;
int L=0; while(!(n>>L&1)) L++;
for(int i=1;i<(n<<1);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
FFT(tmp,n<<1,1); FFT(b,n<<1,1);
for(int i=0;i<(n<<1);i++)
tmp[i]=(1LL*b[i]*2+P-1LL*tmp[i]*b[i]%P*b[i]%P)%P;
FFT(tmp,n<<1,0);
for(int i=0;i<n;i++) b[i]=tmp[i],b[i+n]=0;
}
应用的话,就是一些解一些现行递推式,或者求伯努利数
预处理伯努利数:http://blog.youkuaiyun.com/coldef/article/details/72876762
学习自:http://blog.miskcoo.com/2015/05/polynomial-inverse
多项式开根
就是对
A(x)
A
(
x
)
,求
B(x)
B
(
x
)
,使
B2(x)≡A(x)(modxn)
B
2
(
x
)
≡
A
(
x
)
(
mod
x
n
)
同样是用倍增的方法。
已求得
G(x)
G
(
x
)
,使
G2(x)≡A(x)(modxn)
G
2
(
x
)
≡
A
(
x
)
(
mod
x
n
)
那么
n=0 n = 0 的时候就是 B(x) B ( x ) 的常数项就是1(我也不知道为什么,Manchery说的)
时间复杂度 T(n)=T(n/2)+多项式求逆的复杂度+O(nlogn)=O(nlogn) T ( n ) = T ( n / 2 ) + 多 项 式 求 逆 的 复 杂 度 + O ( n log n ) = O ( n log n )
void Sqrt(int *a,int *b,int n){
if(n==1) return void(b[0]=1);
Sqrt(a,b,n>>1);
memset(invb,0,sizeof(int)*n); Inv(b,invb,n);
int L=0; while(!(n>>L&1)) L++;
for(int i=1;i<(n<<1);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
for(int i=0;i<n;i++) tmp[i]=a[i],tmp[i+n]=0;
FFT(invb,n<<1,1); FFT(tmp,n<<1,1);
for(int i=0,inv2=P+1>>1;i<(n<<1);i++) tmp[i]=1LL*tmp[i]*inv2%P*invb[i]%P;
FFT(tmp,n<<1,0);
for(int i=0,inv2=P+1>>1;i<n;i++) b[i]=(1LL*b[i]*inv2+tmp[i])%P;
}
学习自: http://www.cnblogs.com/acha/p/6472798.html
多项式取模
此取模不是指系数取模,而是多项式模多项式
UPD:
暴力的多项式取模比较简单,
设模多项式为
M(x)=∑ni=0aixi
M
(
x
)
=
∑
i
=
0
n
a
i
x
i
,化一下可以得到
xn=∑n−1i=0aianxi
x
n
=
∑
i
=
0
n
−
1
a
i
a
n
x
i
,将所有大于
n
n
的项带入就行
对于,
B(x)
B
(
x
)
,设
degA(x)=n,degB(x)=m(m≤n)
deg
A
(
x
)
=
n
,
deg
B
(
x
)
=
m
(
m
≤
n
)
求
D(x)
D
(
x
)
,
R(x)
R
(
x
)
,满足
A(x)=D(x)B(x)+R(x)
A
(
x
)
=
D
(
x
)
B
(
x
)
+
R
(
x
)
只要求出
D(x)
D
(
x
)
,通过多项式加减乘就可以得到
R(x)
R
(
x
)
把
1x
1
x
带入多项式
两边乘 xn x n
设 fR(x)=xdegf(x)f(1x) f R ( x ) = x deg f ( x ) f ( 1 x )
则
发现 degR(x)<m deg R ( x ) < m
那么
这样多项式求个逆就可以求出
D(x)
D
(
x
)
啦
带回去加加减减就可以得到
R(x)
R
(
x
)
inline void Div(int *a,int n,int *b,int m){
static int tmp[N],A[N],B[N];
for(int i=0;i<m;i++) t[i]=b[m-i-1];
for(int i=0;i<n;i++) A[i]=a[n-i-1];
int nn=1,d=n-m+1; for(;nn<d<<1;nn<<=1);
for(int i=n;i<nn;i++) A[i]=0;
for(int i=d;i<nn;i++) t[i]=0;
for(int i=0;i<nn;i++) B[i]=0;
Inv(t,B,nn);
for(int i=d;i<nn;i++) B[i]=0;
Mul(A,B,max(n,d));
for(int i=d;i<=n<<1;i++) A[i]=0;
for(int i=0;i<d;i++) if(i>d-i-1) swap(A[i],A[d-i-1]);
for(int i=0;i<m;i++) t[i]=b[i];
Mul(t,A,max(d,m));
for(int i=0;i<n;i++) a[i]=(a[i]+P-t[i])%P;
}
学习自:http://blog.miskcoo.com/2015/05/polynomial-division
大致就是这些算法
还有多项式求exp什么的…还不会…先挖个坑
多项式exp
可以参考金策的2015年国家集训队论文
两边同求导
对于每一项的系数
也就是说
分治FFT
牛顿迭代不会~~
void exp(int *a,int *b,int l,int r){
if(l==r){
add(b[l],1LL*inv[l]*a[l]%P); return ;
}
int mid=l+r>>1;
exp(a,b,l,mid);
static int tmpa[N],tmpb[N];
for(int i=l;i<=mid;i++) tmpa[i-l]=b[i];
for(int i=1;i<=r-l;i++) tmpb[i]=a[i];
int m,L=0;
for(m=1;m<r-l+1;m<<=1,L++); m<<=1;
for(int i=1;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
NTT(tmpa,m,1); NTT(tmpb,m,1);
for(int i=0;i<m;i++) tmpa[i]=1LL*tmpa[i]*tmpb[i]%P;
NTT(tmpa,m,0);
for(int i=mid+1;i<=r;i++) add(b[i],1LL*inv[i]*tmpa[i-l]%P);
for(int i=0;i<m;i++) tmpa[i]=tmpb[i]=0;
exp(a,b,mid+1,r);
}
应用
根据生成函数的卷积运算优化转移什么的都很常见
还有就是多项式的多点求值,这玩意看http://blog.miskcoo.com/2015/05/polynomial-multipoint-eval-and-interpolation
我实现了一下…然而自带大常数,跑的不如暴力…
NOI2017D1T3,就运用了多项式取模来解决线性递推数列的第n项(n很大),也就是 f(x)=f(x)g(x)+1 f ( x ) = f ( x ) g ( x ) + 1 的第n项,
具体看叉姐论文和Picks的博客
http://www.docin.com/p-724323397.html
http://picks.logdown.com/posts/197262-polynomial-division
codechef rng是模板题,可以在上面看大神代码
其他的用处…惯例先挖个坑
本人因能力限制,代码仅供参考,如有疑问,欢迎与本人讨论