怎么感觉讲的这么快,我好菜啊
终于差不多把多项式水过去了……来波总结多项式真好玩
1.FFT
这里主要搞一个裸的FFT,原理已经讲过了
放个模板跑:
# include <algorithm>
# include <complex>
# include <cstdio>
# include <cmath>
# define N 262144
using namespace std;
const double PI=acos(-1);
typedef complex<double> E;
int rev[N];
int n,m,L;
E a[N+N],b[N+N],c[N+N];
void FFT(E f[],int s){
for (int i=0;i<n;++i) if (i>rev[i]) swap(f[i],f[rev[i]]);
for (int i=1;i<n;i<<=1){
E wn(cos(PI/i),sin(PI/i)); wn*=s;
for (int j=0;j<n;j+=(i<<1)){
E w(1,0);
for (int k=0;k<i;++k,w*=wn){
E x=f[j+k],y=w*f[j+k+i];
f[j+k]=x+y,f[j+k+i]=x-y;
}
}
}
if (s == -1) for (int i=0;i<n;++i) f[i] /= n;
}
int main(){
scanf("%d",&n); --n;
m=2*n; for (n=1;n<=m;n<<=1) ++L;
for (int i=0;i<n;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
fft(a,1),fft(b,1);
for (int i=0;i<n;++i) c[i]=a[i]*b[i];
fft(c,-1);
for (int i=0;i<(m>>1);++i) printf("%.8lf ",c[i].real());
return 0;
}
这个写法是经典的蝴蝶操作的写法,而且稍微改进过了一下,常数变得更小(交换了循环,于是访问更加连续)
2.NTT
特殊
好像比较慢?反正原理也讲过了(在原根意义下可以搞出相应的单位根,所以把复数部分搞一搞就可以了),就是由于
gP−1≡1(modP)
,故
e−2πiN
相当于
gP−1N(modP)
进一步地,变换公式里面写的那点东西也相应乘的是
gnk
和
g−nk
,就搞完了
这里的质数
P
一般需要满足
其他模数
但是怎样做对于其他的模数的东西呢?
首先合数就不管了……CRT直接做啊……所以只考虑质数
做法1
我们可以利用中国剩余定理合并:
首先选取三个模数,一般为
998244353
,
1004535809
,
46976204
,或
9985661441=235∗222+1
,其中
3
都是上面模数的原根
然后对三个模数做DFT,得出结果之后乘起来再取模就可以了,只要三个模数乘积大于
可是要
9
次点值和插值……好慢啊
做法2
把两个多项式
令
然后考虑乘起来之后怎么样?
对于 Ai 和 Bj (注意不要看错了,这里下标不是一样的符号)有
于是我们只要先搞四次点值,求出 a′ib′j , a′′ib′j , a′ib′′j , a′′ib′′j 就可以做了
然后插值回去只要三次就可以了
这样少了两次但是还是很慢
做法3
这种做法好像是myy论文的方法……表示非常地强啊
我们考虑FFT的时候一开始虚部是 0 ,这是一种浪费,我们考虑如果我们是对一个复数序列DFT会出现些什么奇怪的东西来
先考虑某个其中第
其中 θ 的意义就是那个单位根的转角
然后再考虑对于另一边的 A(ω−wn) 的影响
然后就很明显了?我们如果分别设上面两个结果的实部和虚部为 A(ωwn)=x1+iy1 和 A(ω−wn)=x2+iy2
我们可以解出
然后我们就可以将两个DFT合并成一个了!跑两遍DFT就可以求出四个点值序列!
然后插值回来有没有什么更简单的办法呢?考虑IDFT之后,一个复数变成了实数,而如果我们把这个复数乘上一个 i ,就会变成纯虚数,那么我们只要把多项式的点值表达
因此其实我们又可以将两个插值合成一个,然后就只需要一半的次数就可以搞完了!
这种方法是比较普适的,所以属于一个很好的优化!
3.多项式的各种运算
我们对数列操作的时候经常会搞到生成函数,然后就经常需要对于由题目给出数据所决定的多项式进行奇怪的运算,例如求逆,除,模,开方,取对数,指数,牛顿迭代等等的操作,其中因为一般题目求的是像“第
n
项的值”这样的东西,那么一般就是在模
a.多项式求逆
就是对于一个给定的
A(x)
,要找出另一个
B(x)
,使得
A(x)B(x)≡1(modxn)
,那么这里
B(x)
记为
A−1(x)
利用倍增算法解决,
n=1
的时候可以直接求常数项的逆元
n>1
的时候,只需要考虑怎么将模
xn2
意义下的逆元
B(x)
转成这里我们需要的模
xn
意义下的逆元
G(x)
有关系式
A(x)B(x)≡1(modxn2)
就是
A(x)B(x)−1≡0(modxn2)
两边平方,我们发现
modxn2
意义下就变成了
modxn
意义下的答案(次数变成两倍),那么这时就可以化了
得到
已经大致有眉目了,我们稍微化一下式子就得到 A(x)[2B(x)−A(x)B2(x)]≡1(modxn)
于是就有 G(x)≡2B(x)−A(x)B2(x)(modxn)
倍增即可,复杂度由主定理还是 O(nlogn)
b.多项式开根
现在我们继续求另一个
B2(x)≡A(x)(modxn)
类似尺规作图中考虑如何开方的思路,最终可以转成求逆的问题来解决
仍然考虑倍增,
n=1
的时候考虑Cipolla算法即可
然后
n>1
的时候考虑怎样从
B(x)modxn2
推
G(x)modxn
类似的套路,但是我们先要考虑到因为是平方,所以
G(x)
对于
modxn2
的情况也一样是成立的
写出关系式
B2(x)≡A(x)(modxn2)
,
G2(x)≡A(x)(modxn2)
两式相减,然后再同样的套路两边平方得到
G4(x)+B4(x)−2G2(x)B2(x)≡0(modxn)
实际上这时我们就可以考虑把这个多项式配起来,两边加上
4G2(x)B2(x)
得到
[G2(x)+B2(x)]2≡(2G(x)B(x))2(modxn)
这时可以拆开平方(因为里面两个多项式的次数仍然是达到
xn
的)
得到
G2(x)+B2(x)≡2G(x)B(x)(modxn)
利用定义式
G2(x)≡A(x)(modxn)
就基本上可以得到答案了
可以发现
G(x)≡B(x)2+A(x)2B(x)(modxn)
每一次求个逆就可以了,时间复杂度还是
O(nlogn)
的
但是常数巨大,一般不会出什么题目(除非是一些生成函数的方程是二次方程需要求解的情况)
另外需要注意的是根据类似二次剩余的思路,多项式开方里面也是有一半不可做的
c.多项式除法&多项式取模
做法
我们做出多项式除法就可以求取模了,所以两个东西是一样的
顺便澄清一下,多项式除法不能用求逆做,因为求逆相当于没有余式
考虑就是对于某个
A(x)
,要求分解出
A(x)=B(x)C(x)+D(x)
,其中
degD<degB
不妨先设
degA=n
,
degB=m<degD
,那么
degC>n−m+1
将系数向量反过来搞(代入
1x
)并乘回一个
xn
就得到
这时候考虑下各个项的指数的范围, xnA(1x) 的在 [0,n] ; xmB(1x) 的在 [0,m] ; xn−mC(1x) 的在 [0,n−m] ; xnD(1x) 的在 [n−m+1,m]
这时就有点思路了,如果我们要求 C ,那直接在
然后就有 xn−mC(1x)≡xnA(1x)xmB(1x) ,求逆即可求得
等等好像还没做完啊这时候左边不是还没消掉吗?那些 1x 怎么求啊?
实际上我们会发现,我们代进去 1x ,但是外面也乘了相应的次数上去,实际上相当于把系数向量翻转然后移了下位而已,那就很可以做了
时间复杂度还是 O(nlogn)
(非常不常用的)应用
- 多项式欧几里得,这时和原来几乎是一样的,复杂度 O(n2logn) (基本没什么题目用)
- 多项式扩展欧几里得,也是一样的,可以用来求多项式模另一个多项式意义下的逆元(有情况下有用,可用来加速常系数线性递推)
然后下面有两个常用的应用:多点求值和多点插值
多点求值
给出多项式
F(x)
,求出
F(x)
在
x1,x2,x3,⋯
这些点上的值
考虑分治,设
Gl,r=∏i=lr(x−xi)
,
Ql,r=F(x)modGl,r(x)
很明显
F(xi)=F(x)modGi,i(x)
(由余数定理可知)
然后利用分治每次暴力往下取模就可以了
因为式子的规模是逐渐减小的,每次取模
O(nlogn)
,总的时间就是
O(nlog2n)
常数巨大,所以感觉可能暴力秦九韶求值都和这差不多
多点插值
主要就是搞拉格朗日插值,在求和号里面把那个多项式分子分母分开来算
F(x)=∑i=1m[∏1≤j≤m,j≠i(x−xj)][yi∏1≤j≤m,j≠i1x−xj]
设右边括号里的是
Gi
,考虑怎么求
我们考虑设
M(x)=∏i=1m(x−xi)
有
考虑事实上有 M(xi)=0 ,那么可以放上去:
由于这里考虑的是 x→xi 的时候,所以说有 Gi=yiM′(xi)
对 M′(x) 做多点求值就可以求出 Gi 了
剩下的就是 F(x)=∑i=1mGi[∏1≤j≤m,j≠i(x−xi)]
这下就可以直接分治FFT解决了,时间复杂度显然是 O(nlog2n)
d.多项式求对数、牛顿迭代、多项式求指数
实际上下面的问题都是在模某个
xn
意义下才能进行的
求对数就很简单了,考虑
ln′(A(x))=A′(x)A(x)
,然后再积分回来就可了,时间复杂度
O(nlogn)
牛顿迭代的原理大致就是利用Taylor展开,可以证明这在模多项式的意义下也可以做
现在要搞得就是给出一个
g(x)
,求一个
f(x)
的前
n
项,使得
其中
f
是一个
那么 有
然后从第三项起最低次项次数都至少是 x2n ,就会被整个模掉,只需要展开到一阶就可以了
又 g(f)=0 ,故
然后移项得到 f≡f0−g(f0)g′(f0)(modx2n)
就可以牛顿迭代辣!我们就是用 f0→f 不断迭代就可以了
顺便说一句,这里的求导是相当于将 f 看做变量,其余
但是一般这样的过程很慢……一般就是要用题目的特殊性来求解(picks表示查到的最好时间复杂度为 O((nlogn)1.5) )
那么这时候考虑回怎么求指数,实际上只要解决 B(x)≡exp(A(x))(modxn) 的情况就可以了
考虑设 g(x)=ln(x)−A ,最终要使 g(B)=0
稍微化一下式子, B≡B0−g(B0)g′(B0)(modxn)
而又有 g′(B0)=1B0 ,那么代回去只需要求一下对数就可以了,时间复杂度 O(nlogn)
顺便提一下,牛顿迭代还可以用来回顾我们之前写的几个式子
多项式求逆:对 f(x)p(x)−1≡0(modxn) ,我们将 f(x) 视作变量就有
是不是很眼熟?事实上,这就是我们之前得到的式子,因为牛顿迭代法本身也相当于一种类似分治的东西
更多应用参见picks的博客: http://picks.logdown.com/posts/209226-newtons-method-of-polynomial
e.FWT
最后突然发现这个毒瘤东西没讲……
据zfr大佬说这实际上只是个小tricks?
主要就是求
∑i⊕jAiBj
其中
⊕
表示一些二进制运算,例如xor/and/or
然后和FFT差不多,定义变换和逆变换,搞到另一个域去直积再变换回来
不加证明地给出几个变换……
然后代码就都差不多了……这么对称的东西随便记一记就可以了吧