文章目录
所谓多项式,就是很多个项组成的式子
(逃)
前言
多项式可以优化卷积、快速的计算式子,以及做很多奇怪的事。
因此,了解多项式知识是必要的。(牵强附会)
过程中会涉及一些导数和积分内容,但都会给出较为详细的讲解或证明,想我一样的导数了解较少的同学请放心使用。
多项式
如果一个函数可以写成: f ( x ) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 + ⋯ + a n x n f(x)=a_0+a_1x+a_2x^2+a_3x^3+ \dots +a_{n}x^{n} f(x)=a0+a1x+a2x2+a3x3+⋯+anxn,那么就称这个函数是一个多项式,也可写作 f ( x ) = ∑ i n a i x i f(x)=\sum_i^na_ix^i f(x)=∑inaixi。
注意!多项式的项数必须是有限的!比如,正弦函数 sin ( x ) \sin(x) sin(x) 的项数是无穷的,所以不是一个多项式。
快速傅立叶变换(FFT)
FFT,用于在
n
l
o
g
n
nlogn
nlogn的复杂度内求两个多项式相乘(卷积)。
算是多项式的入门吧。
很妙的算法。
代码写成迭代,精炼之后很优美。
系数表示法和点值表示法
系数表示法就是我们常用的函数的表示形式。
对于一个n次的函数,点值表示法则是给出了n+1个互异的点,通过这n+1个互异的点就可以唯一确定的描述出这个函数(可以想想高斯消元)。
这东西有啥用?
考虑如果两个多项式都表示成了点值形式,且选取的点的横坐标相同,我把它们的纵坐标对应乘起来,就能得到它们乘积的多项式的点值表示
乘起来的复杂度降到了
O
(
n
)
O(n)
O(n)
但是这个点值表示法求起来暴力似乎还是n方啊…
所以快速傅立叶变换其实解决的就是系数表示法和点值表示法快速互相转换的问题
单位根
规定:
若复数 ω \omega ω满足 ω n = 1 \omega^n=1 ωn=1,就称 ω \omega ω是n次单位根
对于一个n项的多项式,我们让它的点值表示取的横坐标为
ω
n
k
(
0
≤
k
<
n
)
\omega_n^k(0\leq k <n)
ωnk(0≤k<n)
然后就会有很多很妙的性质:
- ω n m k m = ω n k \omega_{nm}^{km}=\omega_n^k ωnmkm=ωnk
- ω n k + n = ω n k \omega_{n}^{k+n}=\omega_n^k ωnk+n=ωnk
- ω n k + n 2 = − w n k \omega_n^{k+\frac{n}{2}}=-w_n^k ωnk+2n=−wnk
都可以结合几何意义较为显然的证明。
离散傅立叶变换(DFT)
那么回到问题。
首先,我们要把多项式填零直到
n
=
2
k
n=2^k
n=2k 的长度。
接下来,我们要求:
f
(
x
)
=
a
0
+
a
1
x
+
a
2
x
2
+
a
3
x
3
+
.
.
.
+
a
n
−
1
x
n
−
1
f(x)=a_0+a_1x+a_2x^2+a_3x^3+...+a_{n-1}x^{n-1}
f(x)=a0+a1x+a2x2+a3x3+...+an−1xn−1
让我们把这个函数分奇偶拆成两个函数:
A
(
x
)
=
a
0
+
a
2
x
+
a
4
x
2
+
a
6
x
3
+
.
.
.
+
a
n
−
1
x
n
2
A(x)=a_0+a_2x+a_4x^2+a_6x^3+...+a_{n-1}x^{\frac{n}{2}}
A(x)=a0+a2x+a4x2+a6x3+...+an−1x2n
B
(
x
)
=
a
1
+
a
3
x
+
a
5
x
2
+
a
7
x
3
+
.
.
.
+
a
n
−
2
x
n
2
B(x)=a_1+a_3x+a_5x^2+a_7x^3+...+a_{n-2}x^{\frac{n}{2}}
B(x)=a1+a3x+a5x2+a7x3+...+an−2x2n
那么就有:
f
(
x
)
=
A
(
x
2
)
+
x
B
(
x
2
)
f(x)=A(x^2)+xB(x^2)
f(x)=A(x2)+xB(x2)
我们把
x
=
ω
n
k
x=\omega_n^k
x=ωnk 带入。
分情况讨论:
(设:
0
≤
k
<
n
2
0\le k<\dfrac{n}{2}
0≤k<2n)
当指数
<
n
2
< \dfrac{n}{2}
<2n 时:
f
(
ω
n
k
)
=
A
(
(
ω
n
k
)
2
)
+
ω
n
k
B
(
(
ω
n
k
)
2
)
f(\omega_n^k)=A((\omega_n^k)^2)+\omega_n^kB((\omega_n^k)^2)
f(ωnk)=A((ωnk)2)+ωnkB((ωnk)2)
=
A
(
ω
n
2
k
)
+
ω
n
k
B
(
ω
n
2
k
)
=A(\omega_n^{2k})+\omega_n^kB(\omega_n^{2k})
=A(ωn2k)+ωnkB(ωn2k)
=
A
(
ω
n
2
k
)
+
ω
n
k
B
(
ω
n
2
k
)
=A(\omega_{\frac{n}{2}}^{k})+\omega_n^kB(\omega_{\frac{n}{2}}^{k})
=A(ω2nk)+ωnkB(ω2nk)
类似的,当指数
≥
n
2
\ge\dfrac{n}{2}
≥2n 时:
f
(
ω
n
k
+
n
2
)
=
A
(
ω
n
2
k
+
n
)
+
ω
n
k
+
n
2
B
(
ω
n
2
k
+
n
)
f(\omega_n^{k+\frac{n}{2}})=A(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}B(\omega_n^{2k+n})
f(ωnk+2n)=A(ωn2k+n)+ωnk+2nB(ωn2k+n)
=
A
(
ω
n
2
k
)
−
ω
n
k
B
(
ω
n
2
k
)
=A(\omega_n^{2k})-\omega_n^{k}B(\omega_n^{2k})
=A(ωn2k)−ωnkB(ωn2k)
=
A
(
ω
n
2
k
)
−
ω
n
k
B
(
ω
n
2
k
)
=A(\omega_{\frac{n}{2}}^{k})-\omega_n^{k}B(\omega_{\frac{n}{2}}^{k})
=A(ω2nk)−ωnkB(ω2nk)
这样,我们就可以在
O
(
l
e
n
)
O(len)
O(len) 的复杂度内合并两个长度为
l
e
n
len
len 的多项式,就可以利用分治把复杂度降到
O
(
n
log
n
)
O(n\log n)
O(nlogn)。
位逆序置换(蝴蝶变换)
这个也很妙
有点线性求逆元那味了
通过观察发现,系数 i 调换位置之后会去往它的二进制表示反过来的位置(记为
r
i
r_i
ri)
而上面那个东西可以通过:
r[i]=r[i>>1]>>1+[i&1]*(n>>1)
来线性的求
这个还是挺好理解的画一画就行了
离散傅立叶逆变换(IDFT)
我们变完之后当然还要变回来。
我会
n
3
n^3
n3 高斯消元!
设之前DFT得到的点值表示
b
k
=
f
(
ω
n
k
)
=
∑
i
=
0
n
−
1
a
i
(
ω
n
k
)
i
b_k=f(\omega_n^k)=\sum_{i=0}^{n-1}a_i(\omega_{n}^k)^i
bk=f(ωnk)=∑i=0n−1ai(ωnk)i
构造一个函数:
g
(
x
)
=
b
0
+
b
1
x
+
b
2
x
2
+
b
3
x
3
+
.
.
.
+
b
n
−
1
x
n
−
1
g(x)=b_0+b_1x+b_2x^2+b_3x^3+...+b_{n-1}x_{n-1}
g(x)=b0+b1x+b2x2+b3x3+...+bn−1xn−1
然后我们用
ω
n
0
,
ω
n
−
1
,
ω
n
−
2
,
.
.
.
,
ω
n
−
n
\omega_n^{0},\omega_n^{-1},\omega_n^{-2},...,\omega_n^{-n}
ωn0,ωn−1,ωn−2,...,ωn−n 作为横坐标带入。(不难发现刚才 DFT 需要的性质依然成立,所以可以用同样的方法求解)
那么我们在第
k
k
k 位得到的新的点值就是:
g
(
ω
n
−
k
)
=
∑
i
=
0
n
−
1
b
i
(
ω
n
−
k
)
i
g(\omega_n^{-k})=\sum_{i=0}^{n-1}b_i(\omega_n^{-k})^i
g(ωn−k)=i=0∑n−1bi(ωn−k)i
=
∑
i
=
0
n
−
1
∑
j
=
0
n
−
1
a
j
(
ω
n
i
)
j
(
ω
n
−
k
)
i
=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_{n}^i)^j(\omega_n^{-k})^i
=i=0∑n−1j=0∑n−1aj(ωni)j(ωn−k)i
=
∑
j
=
0
n
−
1
a
j
∑
i
=
0
n
−
1
(
ω
n
j
−
k
)
i
=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i
=j=0∑n−1aji=0∑n−1(ωnj−k)i
考虑
∑
i
=
0
n
−
1
(
ω
n
j
−
k
)
i
\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i
∑i=0n−1(ωnj−k)i,这一项,当
j
=
k
j=k
j=k 时,显然等于
n
n
n;当
j
≠
k
j\ne k
j=k 时,根据等比公式,有:
∑
i
=
0
n
−
1
(
ω
n
j
−
k
)
i
=
(
w
n
j
−
k
)
n
−
1
w
n
j
−
k
−
1
\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i=\frac{(w_n^{j-k})^n-1}{w_n^{j-k}-1}
i=0∑n−1(ωnj−k)i=wnj−k−1(wnj−k)n−1
=
(
w
n
n
)
j
−
k
−
1
w
n
j
−
k
−
1
=
1
j
−
k
−
1
w
n
j
−
k
−
1
=
0
=\frac{(w_n^{n})^{j-k}-1}{w_n^{j-k}-1}=\frac{1^{j-k}-1}{w_n^{j-k}-1}=0
=wnj−k−1(wnn)j−k−1=wnj−k−11j−k−1=0
所以原来的式子只有当
j
=
k
j=k
j=k 的时候有值,等于
n
a
k
na_k
nak。
所以我们直接 FFT 完除
n
n
n 即可。
struct node{
double x,y;
node(double a=0,double b=0){
x=a;y=b;
}
}A[N],B[N];
il node operator * (node a,node b){
return (node){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
}
il node operator + (node a,node b){
return (node){a.x+b.x,a.y+b.y};
}
il node operator - (node a,node b){
return (node){a.x-b.x,a.y-b.y};
}
int r[N];
il void fft(node *x,int flag){
for(int i=1;i<=len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));
for(int i=0;i<=len;i++){
if(i<r[i]) swap(x[i],x[r[i]]);
}
for(int l=1;l<len;l<<=1){
node o(cos(pi/l),flag*sin(pi/l));
for(int st=0;st<len;st+=l<<1){
node t(1,0);
for(int j=0;j<l;j++,t=t*o){
node u=x[st+j],v=t*x[st+j+l];
x[st+j]=u+v;
x[st+j+l]=u-v;
}
}
}
return;
}
快速数论变换
学完这个就可以把 FFT 忘掉了。
然而,由于 NTT 的思想是建立在 FFT 的基础之上,学习 FFT 还是有必要的。
原根
NTT说简单也简单说难也难。
简单是因为板子几乎和FFT一模一样。
难是因为那个群的概念完全看不懂qwq。
然后我就选择调过了证明部分。
背下来背下来!——宋队
总的来说,写一写对背板子有帮助的我可以理解的。
对于一个质数P,若其存在原根(基本全是3)且可以表示为:
P
=
k
∗
2
n
+
1
P=k*2^n+1
P=k∗2n+1
它就是一个ntt模数。
g
g
g为什么叫原根?
因为它有一些关键性质:
1. g k ∗ n = 1 ( m o d P ) 1. g^{k*n}=1 \pmod P 1.gk∗n=1(modP)
2. g i m o d P g^i \bmod P gimodP 在 0 ≤ i ≤ p 0\leq i\leq p 0≤i≤p 的范围内取值两两不同。
设:
g
k
g^k
gk 为单位根,记为
g
n
g_n
gn。
那么它就有很多熟悉的性质:
- g n n ≡ 1 ( m o d P ) g_n^n \equiv 1 \pmod P gnn≡1(modP)
- g 2 k 2 k = g k k g_{2k}^{2k}=g_k^k g2k2k=gkk
- g n k + n 2 = − g n k g_n^{k+ \frac n 2 }=-g_n^k gnk+2n=−gnk
没错,推FFT时
ω
\omega
ω 需要有的性质它全有!
所以就可以直接上FFT的板子啦!
然后烦人的复数、浮点运算和精度问题全都拜拜了。
针不戳。
宋队诚不欺我。
逆变换还是倒数,从FFT的求共轭改成求逆元。
il void ntt(ll *x,int flag){
for(int i=0;i<lim;i++){
if(i<r[i]) swap(x[i],x[r[i]]);
}
for(register int l=1;l<lim;l<<=1){
register ll g=ksm(3,(mod-1)/(l<<1));if(flag==-1) g=ksm(g,mod-2);
for(register int st=0;st<lim;st+=l<<1){
register ll now=1;
for(register int i=0;i<l;i++,now=now*g%mod){
ll u=x[st+i],v=x[st+i+l]*now%mod;
x[st+i]=(u+v)%mod;x[st+i+l]=(u-v+mod)%mod;
}
}
}
if(flag==-1){
register ll ni=ksm(lim,mod-2);
for(register int i=0;i<lim;i++) x[i]=x[i]*ni%mod;
}
}
快速莫比乌斯变换 & 快速沃尔什变换(FMT&FWT)
正常卷积的定义:
c
k
=
∑
i
+
j
=
k
a
i
b
j
c_k=\sum_{i+j=k}a_ib_j
ck=∑i+j=kaibj。
可以用FFT或者NTT在
O
(
n
log
n
)
O(n\log n)
O(nlogn) 的复杂度内解决。
然而,有些时候,我们计算卷积的时候下标关系并不是喜闻乐见的加法,而是形如
c
k
=
∑
i
⊕
j
=
k
a
i
b
j
c_k=\sum_{i\oplus j=k}a_ib_j
ck=∑i⊕j=kaibj,其中的
⊕
\oplus
⊕ 是 and or xor
中的一种位运算。
这个时候就需要快速莫比乌斯变换(and or
)和快速沃尔什变换(xor
),同样可以把复杂度降到
O
(
n
log
n
)
O(n\log n)
O(nlogn) 级别。
思想
由于三个运算思想比较类似,因此不对两个算法进行特别的区分,统称为FWT。
采用类比的思想,FFT先把多项式转化为点值表示,乘到一起,最后再反演回去。
类似的,FWT也是先把多项式
A
,
B
A,B
A,B 转化为某种变换
FWT
(
A
)
,
FWT
(
B
)
\operatorname{FWT}(A),\operatorname{FWT}(B)
FWT(A),FWT(B),然后乘起来得到
FWT
(
C
)
\operatorname{FWT}(C)
FWT(C),最后再反演回去得到
C
C
C。
根据不同位运算的性质,对应的变换有所不同。
OR
定义
求:
c
k
=
∑
i
∣
j
=
k
a
i
b
j
c_k=\sum_{i\operatorname{|}j=k}a_ib_j
ck=∑i∣j=kaibj
考虑或运算有如下性质:若
i
∣
k
=
k
,
j
∣
k
=
k
i|k=k,j|k=k
i∣k=k,j∣k=k,那么
(
i
∣
j
)
∣
k
=
k
(i|j)|k=k
(i∣j)∣k=k,反之亦然。
那么我们就按照这个充要条件设计变换:
FWT
o
r
(
A
)
k
=
∑
i
∣
k
=
k
A
i
\operatorname{FWT}_{or}(A)_k=\sum_{i|k=k}A_i
FWTor(A)k=∑i∣k=kAi。
自然语言:
i
i
i 必须是
k
k
k 的子集才能转移。
考虑两个变换后的多项式逐项系数相乘:
FWT
o
r
(
A
)
∗
FWT
o
r
(
B
)
=
∑
k
(
∑
i
∣
k
=
k
a
i
)
×
(
∑
j
∣
k
=
k
b
j
)
\operatorname{FWT}_{or}(A)*\operatorname{FWT}_{or}(B)=\sum_{k}(\sum_{i|k=k}a_i)\times (\sum_{j|k=k}b_j)
FWTor(A)∗FWTor(B)=k∑(i∣k=k∑ai)×(j∣k=k∑bj)
=
∑
k
∑
i
∣
k
=
k
∑
j
∣
k
=
k
a
i
b
j
=
∑
k
∑
(
i
∣
j
)
∣
k
=
k
a
i
b
j
=
FWT
o
r
(
C
)
=\sum_{k}\sum_{i|k=k}\sum_{j|k=k}a_ib_j=\sum_{k}\sum_{(i|j)|k=k}a_ib_j=\operatorname{FWT}_{or}(C)
=k∑i∣k=k∑j∣k=k∑aibj=k∑(i∣j)∣k=k∑aibj=FWTor(C)
所以我们就证出了系数直接相乘是对的,现在只需要一种快速的变换得到
FWT
o
r
\operatorname{FWT}_{or}
FWTor 以及将其反演的方法。
变换:
还是和FFT类似的,考虑分治,假设我们要合并两个长度为
2
n
2^n
2n 的序列,设它们不同的最高位为1的序列为
A
1
A_1
A1,为0的为
A
0
A_0
A0。
考虑 or
的性质,
A
0
A_0
A0 可以转移到
A
1
A_1
A1,而
A
1
A_1
A1 无法转移到
A
0
A_0
A0,因此有:
FWT
o
r
(
A
)
=
(
FWT
o
r
(
A
0
)
,
FWT
o
r
(
A
0
)
+
FWT
o
r
(
A
1
)
)
\operatorname{FWT}_{or}(A)=(\operatorname{FWT}_{or}(A_0),\operatorname{FWT}_{or}(A_0)+\operatorname{FWT}_{or}(A_1))
FWTor(A)=(FWTor(A0),FWTor(A0)+FWTor(A1))
其中
+
+
+ 表示系数逐项相加,
(
f
(
x
)
,
g
(
x
)
)
(f(x),g(x))
(f(x),g(x)) 表示把
f
(
x
)
f(x)
f(x) 和
g
(
x
)
g(x)
g(x) 顺次写出。
逆变换
那么对应的,还原成原序列只需要逆向操作即可:
UFWT
o
r
(
A
)
=
(
UFWT
o
r
(
A
0
)
,
UFWT
o
r
(
A
0
)
−
UFWT
o
r
(
A
1
)
)
\operatorname{UFWT}_{or}(A)=(\operatorname{UFWT}_{or}(A_0),\operatorname{UFWT}_{or}(A_0)-\operatorname{UFWT}_{or}(A_1))
UFWTor(A)=(UFWTor(A0),UFWTor(A0)−UFWTor(A1))
代码写起来非常简洁:
void Or(ll *x,int n,int op){
for(int l=1;l<n;l<<=1){
for(int st=0;st<n;st+=l*2){
for(int i=0;i<l;i++){
ll u=x[st+i],v=x[st+l+i];
if(op==1) x[st+i]=u,x[st+l+i]=(v+u)%mod;
else x[st+i]=u,x[st+l+i]=(v+mod-u)%mod;
}
}
}
}
AND
和 or
的整个定义、证明几乎都是一样的。
求:
c
k
=
∑
i
&
j
=
k
a
i
b
j
c_k=\sum_{i\operatorname{\&}j=k}a_ib_j
ck=∑i&j=kaibj
与运算有如下性质:若
i
&
k
=
k
,
j
&
k
=
k
i\&k=k,j\&k=k
i&k=k,j&k=k,那么
(
i
&
j
)
&
k
=
k
(i\&j)\&k=k
(i&j)&k=k,反之亦然。
按照这个充要条件设计变换:
FWT
o
r
(
A
)
k
=
∑
i
∣
k
=
k
A
i
\operatorname{FWT}_{or}(A)_k=\sum_{i|k=k}A_i
FWTor(A)k=∑i∣k=kAi。
自然语言:
i
i
i 必须是
k
k
k 的超集才能转移。
后面的证明一模一样,变换也可以同理的得到:
FWT
a
n
d
(
A
)
=
(
FWT
a
n
d
(
A
0
)
+
FWT
a
n
d
(
A
1
)
,
FWT
a
n
d
(
A
1
)
)
\operatorname{FWT}_{and}(A)=(\operatorname{FWT}_{and}(A_0)+\operatorname{FWT}_{and}(A_1),\operatorname{FWT}_{and}(A_1))
FWTand(A)=(FWTand(A0)+FWTand(A1),FWTand(A1))
UFWT
a
n
d
(
A
)
=
(
UFWT
a
n
d
(
A
0
)
−
UFWT
a
n
d
(
A
1
)
,
UFWT
a
n
d
(
A
1
)
)
\operatorname{UFWT}_{and}(A)=(\operatorname{UFWT}_{and}(A_0)-\operatorname{UFWT}_{and}(A_1),\operatorname{UFWT}_{and}(A_1))
UFWTand(A)=(UFWTand(A0)−UFWTand(A1),UFWTand(A1))
void And(ll *x,int n,int op){
for(int l=1;l<n;l<<=1){
for(int st=0;st<n;st+=l*2){
for(int i=0;i<l;i++){
ll u=x[st+i],v=x[st+l+i];
if(op==1) x[st+i]=(u+v)%mod,x[st+l+i]=v;
else x[st+i]=(u+mod-v)%mod,x[st+l+i]=v;
}
}
}
}
XOR
这个和之前的有所不同。
定义
求:
c
k
=
∑
i
⊕
j
=
k
a
i
b
j
c_k=\sum_{i\operatorname{\oplus}j=k}a_ib_j
ck=∑i⊕j=kaibj
设
d
(
x
)
d(x)
d(x) 为二进制下
x
x
x 的1的个数的奇偶性(模2的结果)。
有性质:
d
(
i
&
k
)
⊕
d
(
j
&
k
)
=
d
(
(
i
⊕
j
)
&
k
)
d(i\& k)\oplus d(j\& k)=d((i\oplus j)\& k)
d(i&k)⊕d(j&k)=d((i⊕j)&k)。
证明:
由于与运算,只考虑
k
k
k 为1的那些位。如果某一位
i
,
j
i,j
i,j 都是1或0,那么等号两边的奇偶性都不改变;如果某一位
i
,
j
i,j
i,j 一个是1一个是0,那么等号两边的奇偶性都改变;所以等号两遍的奇偶性始终同步改变,最后就一定是相等的。
设计:
FWT
x
o
r
(
A
)
k
=
∑
i
(
−
1
)
d
(
i
&
k
)
a
i
\operatorname{FWT}_{xor}(A)_k=\sum_{i} (-1)^{d(i\& k)}a_i
FWTxor(A)k=i∑(−1)d(i&k)ai
那么还是尝试逐项相乘:
FWT
x
o
r
(
A
)
∗
FWT
x
o
r
(
B
)
=
∑
k
(
∑
i
(
−
1
)
d
(
i
&
k
)
a
i
)
×
(
∑
j
(
−
1
)
d
(
j
&
k
)
b
j
)
\operatorname{FWT}_{xor}(A)*\operatorname{FWT}_{xor}(B)=\sum_k(\sum_{i} (-1)^{d(i\& k)}a_i)\times (\sum_{j} (-1)^{d(j\& k)}b_j)
FWTxor(A)∗FWTxor(B)=k∑(i∑(−1)d(i&k)ai)×(j∑(−1)d(j&k)bj)
=
∑
k
∑
i
,
j
(
−
1
)
d
(
i
&
k
)
+
d
(
j
&
k
)
a
i
b
j
=
∑
k
∑
i
,
j
(
−
1
)
d
(
i
&
k
)
⊕
d
(
j
&
k
)
a
i
b
j
=\sum_k\sum_{i,j} (-1)^{d(i\& k)+d(j\& k)}a_ib_j=\sum_k\sum_{i,j} (-1)^{d(i\& k)\oplus d(j\& k)}a_ib_j
=k∑i,j∑(−1)d(i&k)+d(j&k)aibj=k∑i,j∑(−1)d(i&k)⊕d(j&k)aibj
=
∑
k
∑
i
,
j
(
−
1
)
d
(
(
i
⊕
j
)
&
k
)
a
i
b
j
=
FWT
x
o
r
(
C
)
=\sum_k\sum_{i,j} (-1)^{d((i\oplus j)\& k)}a_ib_j=\operatorname{FWT}_{xor}(C)
=k∑i,j∑(−1)d((i⊕j)&k)aibj=FWTxor(C)
得证。
变换
还是分治的思想,考虑到增加一位后,只有下标最高位带1的数贡献到了下标最高位带1的位置时,最高位与运算结果为1,1的个数增加一个,奇偶性改变,所有的符号都变为相反。其它的转移都不变。
也就是:
FWT
x
o
r
(
A
)
=
(
FWT
x
o
r
(
A
0
)
+
FWT
x
o
r
(
A
1
)
,
FWT
x
o
r
(
A
0
)
−
FWT
x
o
r
(
A
1
)
)
\operatorname{FWT}_{xor}(A)=(\operatorname{FWT}_{xor}(A_0)+\operatorname{FWT}_{xor}(A_1),\operatorname{FWT}_{xor}(A_0)-\operatorname{FWT}_{xor}(A_1))
FWTxor(A)=(FWTxor(A0)+FWTxor(A1),FWTxor(A0)−FWTxor(A1))
逆变换
反过来即可。
FWT
x
o
r
(
A
)
=
(
FWT
x
o
r
(
A
0
)
+
FWT
x
o
r
(
A
1
)
2
,
FWT
x
o
r
(
A
0
)
−
FWT
x
o
r
(
A
1
)
2
)
\operatorname{FWT}_{xor}(A)=(\frac{\operatorname{FWT}_{xor}(A_0)+\operatorname{FWT}_{xor}(A_1)}{2},\frac{\operatorname{FWT}_{xor}(A_0)-\operatorname{FWT}_{xor}(A_1)}{2})
FWTxor(A)=(2FWTxor(A0)+FWTxor(A1),2FWTxor(A0)−FWTxor(A1))
void Xor(ll *x,int n,int op){
for(int l=1;l<n;l<<=1){
for(int st=0;st<n;st+=l*2){
for(int i=0;i<l;i++){
ll u=x[st+i],v=x[st+l+i];
if(op==1) x[st+i]=(u+v)%mod,x[st+l+i]=(u+mod-v)%mod;
else x[st+i]=(u+v)%mod*499122177%mod,x[st+l+i]=(u+mod-v)%mod*499122177%mod;
}
}
}
}
多项式求导
对于一个多项式
f
(
x
)
=
∑
i
=
0
n
a
i
x
i
f(x)=\sum_{i=0}^na_ix^i
f(x)=∑i=0naixi,有:
f
′
(
x
)
=
∑
i
=
0
n
−
1
(
i
+
1
)
a
i
+
1
x
i
f'(x)=\sum_{i=0}^{n-1}(i+1)a_{i+1}x^i
f′(x)=i=0∑n−1(i+1)ai+1xi
通俗的说就是:每一项的系数乘上自己的下标后往下一位传。
证明:
我们可以把多项式的导数写成每一个单项的导数相加,因此我们只需要证出每个单项的
f
(
x
)
=
a
k
x
k
(
k
>
0
)
f(x)=a_kx^k(k>0)
f(x)=akxk(k>0) 的导数是
f
′
(
x
)
=
k
a
k
x
k
−
1
f'(x)=k a_{k}x^{k-1}
f′(x)=kakxk−1 即可。
直接上导数定义:
f
′
(
x
)
=
lim
d
x
→
0
f
(
x
+
d
x
)
−
f
(
x
)
d
x
=
lim
d
x
→
0
a
k
(
x
+
d
x
)
k
−
a
k
x
k
d
x
f'(x)=\lim_{\operatorname d \! x\to0}\frac{f(x+dx)-f(x)}{dx}=\lim_{dx\to0}\frac{a_k(x+dx)^k-a_kx^k}{dx}
f′(x)=dx→0limdxf(x+dx)−f(x)=dx→0limdxak(x+dx)k−akxk
=
a
k
×
lim
d
x
→
0
(
x
+
d
x
)
k
−
x
k
d
x
=
a
k
×
lim
d
x
→
0
∑
i
=
0
k
(
k
i
)
x
i
d
x
k
−
i
−
x
k
d
x
=a_k\times\lim_{\operatorname d \! x\to0}\frac{(x+dx)^k-x^k}{dx}=a_k\times\lim_{\operatorname d \! x\to0}\frac{\sum_{i=0}^k\binom{k}{i}x^idx^{k-i}-x^k}{dx}
=ak×dx→0limdx(x+dx)k−xk=ak×dx→0limdx∑i=0k(ik)xidxk−i−xk
=
a
k
×
lim
d
x
→
0
∑
i
=
0
k
(
k
i
)
x
i
d
x
k
−
i
−
x
k
d
x
=
a
k
×
lim
d
x
→
0
∑
i
=
0
k
−
1
(
k
i
)
x
i
d
x
k
−
i
d
x
=a_k\times\lim_{dx\to0}\frac{\sum_{i=0}^k\binom{k}{i}x^idx^{k-i}-x^k}{dx}=a_k\times\lim_{dx\to0}\frac{\sum_{i=0}^{k-1}\binom{k}{i}x^id x^{k-i}}{dx}
=ak×dx→0limdx∑i=0k(ik)xidxk−i−xk=ak×dx→0limdx∑i=0k−1(ik)xidxk−i
=
a
k
×
lim
d
x
→
0
∑
i
=
0
k
−
1
(
k
i
)
x
i
d
x
k
−
i
−
1
=a_k\times\lim_{dx\to0}\sum_{i=0}^{k-1}\binom{k}{i}x^idx^{k-i-1}
=ak×dx→0limi=0∑k−1(ik)xidxk−i−1
由于
d
x
→
0
\operatorname d \! x \to0
dx→0,
d
x
dx
dx 次数大于0的项都可以忽略,因此这个式子就只剩下:
a
k
(
k
k
−
1
)
x
k
−
1
=
a
k
k
x
k
−
1
a_k\binom{k}{k-1}x^{k-1}=a_kkx^{k-1}
ak(k−1k)xk−1=akkxk−1
证毕。
我们顺便也可以写出多项式求导的代码了。
多项式积分就是求导的逆运算。
积分也可以线性求逆元后做到线性。
(不过后面需要用积分函数的地方都不会使快速幂成为算法复杂度瓶颈,这里就偷懒了)
void dao(ll *a,ll *b,int n){
static ll tmp[N];
memcpy(tmp,a,sizeof(ll)*n);
for(int i=0;i<n-1;i++) b[i]=tmp[i+1]*(i+1)%mod;
fill(b+n,b+lim,0);
return;
}
void jifen(ll *a,ll *b,int n){
static ll tmp[N];
memcpy(tmp,a,sizeof(ll)*n);
for(int i=1;i<n;i++) b[i]=tmp[i-1]*ksm(i,mod-2)%mod;
fill(b+n,b+lim,0);
b[0]=0;
return;
}
泰勒展开
为了后面的式子,我们有必要补充一些前置知识。
思想
泰勒展开的思想:
如果存在两个函数
f
(
x
)
,
g
(
x
)
f(x),g(x)
f(x),g(x),在某个位置
x
0
x0
x0 满足:
f
(
x
0
)
=
g
(
x
0
)
,
f
(
1
)
(
x
0
)
=
g
(
1
)
(
x
0
)
,
f
(
2
)
(
x
0
)
=
g
(
2
)
(
x
0
)
.
.
.
f
(
n
)
(
x
0
)
=
g
(
n
)
(
x
0
)
f(x_0)=g(x_0),f^{(1)}(x_0)=g^{(1)}(x_0),f^{(2)}(x_0)=g^{(2)}(x_0)...f^{(n)}(x_0)=g^{(n)}(x_0)
f(x0)=g(x0),f(1)(x0)=g(1)(x0),f(2)(x0)=g(2)(x0)...f(n)(x0)=g(n)(x0)
其中
f
(
k
)
(
x
0
)
f^{(k)}(x_0)
f(k)(x0) 表示
f
(
x
)
f(x)
f(x) 在
x
0
x_0
x0 处的
k
k
k 阶导数,上面的式子也就是指
f
(
x
)
,
g
(
x
)
f(x),g(x)
f(x),g(x) 在
x
0
x_0
x0 处的
0
−
n
0-n
0−n 阶导分别相等。
那么随着
n
n
n 的增大,
g
(
x
)
g(x)
g(x) 就会越来越与
f
(
x
)
f(x)
f(x) 相似,当且仅当其每一阶导数都相同时(即
n
→
∞
n\to \infty
n→∞),两个函数会完全相等。
利用这个思想,我们如果想要仿制某个函数
f
(
x
)
f(x)
f(x) 时,可以通过使仿制函数
g
(
x
)
g(x)
g(x) 的尽量多阶导数都与
f
(
x
)
f(x)
f(x) 相同。
泰勒展开式
利用泰勒展开的思想,我们如果想要仿制某个函数
f
(
x
)
f(x)
f(x) 时,可以通过使仿制函数
g
(
x
)
g(x)
g(x) 的尽量多阶导数都与
f
(
x
)
f(x)
f(x) 相同。
具体如何求呢?
设
g
(
x
)
=
∑
i
n
a
i
x
i
g(x)=\sum_i^na_ix_i
g(x)=∑inaixi。
先考虑
x
0
=
0
x_0=0
x0=0 的情况。
考虑如何求一个特定的系数
a
k
a_k
ak。
注意到,对于多项式
g
(
x
)
g(x)
g(x),在0点时,
g
(
k
)
(
0
)
=
k
!
a
k
g^{(k)}(0)=k!a_k
g(k)(0)=k!ak,因为
k
k
k 阶导0点的函数值只有从
k
k
k 下标一直移到常数项的
a
k
a_k
ak,往下移的过程中一直乘下标,因此总共乘了
k
!
k!
k!。
那么我们就有:
g
(
k
)
(
0
)
=
k
!
a
k
=
f
(
k
)
(
0
)
g^{(k)}(0)=k!a_k=f^{(k)}(0)
g(k)(0)=k!ak=f(k)(0)
也就是:
a
k
=
f
(
k
)
(
0
)
k
!
a_k=\frac{f^{(k)}(0)}{k!}
ak=k!f(k)(0)
我们可以用这个式子求出所有的系数,因此就有
g
(
x
)
=
∑
i
=
0
n
f
(
i
)
(
0
)
i
!
x
i
g(x)=\sum_{i=0}^n\frac{f^{(i)}(0)}{i!}x^i
g(x)=i=0∑ni!f(i)(0)xi
n
n
n 越大,
g
(
x
)
g(x)
g(x) 就越接近
f
(
x
)
f(x)
f(x),当
n
→
∞
n\to\infty
n→∞ 时,
g
(
x
)
g(x)
g(x) 就与
f
(
x
)
f(x)
f(x) 完全相同了,因此任何一个函数
f
(
x
)
f(x)
f(x) 都可以写成:
f
(
x
)
=
∑
i
=
0
∞
f
(
i
)
(
0
)
i
!
x
i
f(x)=\sum_{i=0}^{\infty}\frac{f^{(i)}(0)}{i!}x^i
f(x)=i=0∑∞i!f(i)(0)xi
这就是
f
(
x
)
f(x)
f(x) 在0点展开的泰勒展开式。
如果在一个比较一般的位置
x
0
x_0
x0 展开,可以理解为函数平移,那么式子就变成:
f
(
x
)
=
∑
i
=
0
∞
f
(
i
)
(
x
0
)
i
!
(
x
−
x
0
)
i
f(x)=\sum_{i=0}^{\infty}\frac{f^{(i)}(x_0)}{i!}(x-x_0)^i
f(x)=i=0∑∞i!f(i)(x0)(x−x0)i
这就是泰勒展开的一般形式。
值得注意的是,对于一个 n n n 次多项式,其 f ( i ) ( x ) ( i > n ) f^{(i)}(x)(i>n) f(i)(x)(i>n) 是等于 0 的常函数,所以对于一个 n n n 次多项式,我们只需要展开到 n n n 阶导即可。
牛顿迭代
比较虚空的强大科技…
如果某个地方你没有特别理解(尤其是
g
(
f
(
x
)
)
g(f(x))
g(f(x))),可以结合后面的具体例子,会更好理解一些。
考虑递归倍增法,假设对于给定函数
h
(
x
)
h(x)
h(x) 我们已知了在
(
m
o
d
x
⌈
n
2
⌉
)
\pmod {x^{\lceil\frac{n}{2}\rceil}}
(modx⌈2n⌉) 下的要求函数(比如逆、开根等)
f
0
(
x
)
f_0(x)
f0(x),现在要求在
(
m
o
d
x
n
)
\pmod {x^n}
(modxn) 下的要求函数。
设计一个复合函数
g
(
f
(
x
)
)
g(f(x))
g(f(x)),使其零点位置得到的
f
(
x
)
f(x)
f(x) 的解就是要求的函数。
举个例子:在求平方根的时候,就可以设计 g ( f ( x ) ) = f 2 ( x ) − h ( x ) g(f(x))=f^2(x)-h(x) g(f(x))=f2(x)−h(x)。
(这里的自变量是 f ( x ) f(x) f(x),因此请把 h ( x ) h(x) h(x) 当成常数(或者说:常量)来看)
将
g
(
f
(
x
)
)
g(f(x))
g(f(x)) 在
f
0
(
x
)
f_0(x)
f0(x) 的位置泰勒展开:
g
(
f
(
x
)
)
≡
∑
i
=
0
∞
g
(
i
)
(
f
0
(
x
)
)
i
!
(
f
(
x
)
−
f
0
(
x
)
)
i
(
m
o
d
x
n
)
g(f(x))\equiv\sum_{i=0}^{\infty}\frac{g^{(i)}(f_0(x))}{i!}(f(x)-f_0(x))^i\pmod {x^n}
g(f(x))≡i=0∑∞i!g(i)(f0(x))(f(x)−f0(x))i(modxn)
这个玩意看起来根本不可求,但是我们有一个很好的性质。
考虑
(
f
(
x
)
−
f
0
(
x
)
)
i
(f(x)-f_0(x))^i
(f(x)−f0(x))i 这一项,由于
f
0
(
x
)
≡
f
(
x
)
(
m
o
d
x
⌈
n
2
⌉
)
f_0(x)\equiv f(x) \pmod {x^{\lceil\frac{n}{2}\rceil}}
f0(x)≡f(x)(modx⌈2n⌉),所以这一项的
0
−
⌊
n
2
⌋
0-\lfloor\frac{n}{2}\rfloor
0−⌊2n⌋ 项的系数全是0。
那么这一项平方或者更高次方后在
[
0
,
n
−
1
]
[0,n-1]
[0,n−1] 项之间就没有系数了!
所以我们只需要统计上面那个式子
i
=
0
/
1
i=0/1
i=0/1 的贡献即可。
那么就是:
f
(
x
)
=
g
(
f
0
(
x
)
)
+
g
′
(
f
0
(
x
)
)
(
f
(
x
)
−
f
0
(
x
)
)
≡
0
(
m
o
d
x
n
)
f(x)=g(f_0(x))+g'(f_0(x))(f(x)-f_0(x))\equiv 0\pmod{x^n}
f(x)=g(f0(x))+g′(f0(x))(f(x)−f0(x))≡0(modxn)
移项即得:
f
(
x
)
=
f
0
(
x
)
−
g
(
f
0
(
x
)
)
g
′
(
f
0
(
x
)
)
f(x)=f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}
f(x)=f0(x)−g′(f0(x))g(f0(x))
这个式子是重中之重!
我们把对应的
g
(
f
0
(
x
)
)
g(f_0(x))
g(f0(x)) 带入即可求出我们想要的求特定多项式函数的式子了。
你要是看不明白这个推导背下来上面这个式子也行。
多项式求逆
牛迭的无脑开始体现出来了。
设计 :
g
(
f
(
x
)
)
=
1
f
(
x
)
−
h
(
x
)
g(f(x))=\frac{1}{f(x)}-h(x)
g(f(x))=f(x)1−h(x)
想要带入牛迭的式子,我们需要
f
(
x
)
=
1
x
f(x)=\frac{1}{x}
f(x)=x1 的导数。
直接上定义即可:
f
′
(
x
)
=
lim
d
x
→
0
(
1
x
+
d
x
−
1
x
)
d
x
=
lim
d
x
→
0
x
−
(
x
+
d
x
)
x
2
+
x
d
x
d
x
=
f'(x)=\lim_{dx\to0}\frac{(\dfrac{1}{x+dx}-\dfrac{1}{x})}{dx}=\lim_{dx\to0}\frac{\dfrac{x-(x+dx)}{x^2+xdx}}{dx}=
f′(x)=dx→0limdx(x+dx1−x1)=dx→0limdxx2+xdxx−(x+dx)=
lim
d
x
→
0
−
1
x
2
+
x
d
x
=
−
1
x
2
\lim_{dx\to0}\frac{-1}{x^2+xdx}=-\frac{1} {x^2}
dx→0limx2+xdx−1=−x21
(其实可以发现,如果把这个函数理解成 f ( x ) = x − 1 f(x)=x^{-1} f(x)=x−1,那么我们之前推导的多项式的导数的结论在这里还是适用的。)
那么我们现在带回牛迭:
f
(
x
)
=
f
0
(
x
)
−
g
(
f
0
(
x
)
)
g
′
(
f
0
(
x
)
)
=
f
0
(
x
)
−
1
f
0
(
x
)
−
h
(
x
)
−
1
f
0
2
(
x
)
f(x)= f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}=f_0(x)-\frac{\dfrac{1}{f_0(x)}-h(x)}{-\dfrac{1}{f_0^2(x)}}
f(x)=f0(x)−g′(f0(x))g(f0(x))=f0(x)−−f02(x)1f0(x)1−h(x)
=
2
f
0
(
x
)
−
h
(
x
)
f
0
2
(
x
)
=2f_0(x)-h(x)f_0^2(x)
=2f0(x)−h(x)f02(x)
需要卷积。
时间复杂度:
T
(
n
)
=
T
(
n
2
)
+
O
(
n
log
n
)
=
O
(
n
log
n
)
T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n)
T(n)=T(2n)+O(nlogn)=O(nlogn)
void inv(ll *h,ll *f,int n){
static ll tmp[N];
if(n==1){
f[0]=ksm(h[0],mod-2);return;
}
inv(h,f,(n+1)>>1);
init(n<<1);
fill(f+n,f+lim,0);
for(int i=0;i<n;i++) tmp[i]=h[i];
fill(tmp+n,tmp+lim,0);
NTT(tmp,lim,1);NTT(f,lim,1);
for(int i=0;i<lim;i++) f[i]=(2*f[i]+mod-tmp[i]*f[i]%mod*f[i]%mod)%mod;
NTT(f,lim,-1);
fill(f+n,f+lim,0);
}
多项式开根
还是一样的套路,设计:
g
(
f
(
x
)
)
=
f
0
2
(
x
)
−
h
(
x
)
g(f(x))=f_0^2(x)-h(x)
g(f(x))=f02(x)−h(x)
这个函数就是个一般的多项式,它的导数就是:
g
′
(
f
(
x
)
)
=
2
f
0
(
x
)
g'(f(x))=2f_0(x)
g′(f(x))=2f0(x)
带入牛迭:
f
(
x
)
=
f
0
(
x
)
−
g
(
f
0
(
x
)
)
g
′
(
f
0
(
x
)
)
=
f
0
(
x
)
−
f
0
2
(
x
)
−
h
(
x
)
2
f
0
(
x
)
f(x)= f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}=f_0(x)-\frac{f_0^2(x)-h(x)}{2f_0(x)}
f(x)=f0(x)−g′(f0(x))g(f0(x))=f0(x)−2f0(x)f02(x)−h(x)
=
f
0
2
(
x
)
+
h
(
x
)
2
f
0
(
x
)
=\frac{f_0^2(x)+h(x)}{2f_0(x)}
=2f0(x)f02(x)+h(x)
需要卷积和求逆。
时间复杂度:
T
(
n
)
=
T
(
n
2
)
+
O
(
n
log
n
)
=
O
(
n
log
n
)
T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n)
T(n)=T(2n)+O(nlogn)=O(nlogn)
void Sqrt(ll *h,ll *f,int n){
static ll tmp[N];
if(n==1){f[0]=1;return;}
Sqrt(h,f,(n+1)>>1);
memset(tmp,0,sizeof(ll)*lim);
inv(f,tmp,n);mul(tmp,h,tmp,n,n);
init(n<<1);
NTT(f,lim,1);NTT(tmp,lim,1);
for(int i=0;i<lim;i++){
f[i]=(tmp[i]+f[i])*499122177%mod;
}
NTT(f,lim,-1);
fill(f+n,f+lim,0);
return;
}
多项式求 ln
读者尝试一下就会发现,如果想要用类似的牛迭的方法分别求
ln
\ln
ln 和
exp
\exp
exp,那么就会出现求
ln
\ln
ln 需要
exp
\exp
exp,求
exp
\exp
exp 需要
ln
\ln
ln 的闭环窘境. . . . . .
所以至少有一个我们需要别的方法来求。
考虑用别的方法求
ln
\ln
ln。
首先我们求一下
ln
x
\ln x
lnx 的导数。
我们需要一个
e
e
e 的定义式:
lim
n
→
∞
(
1
+
1
n
)
n
=
e
\lim_{n\to \infty}(1+\frac{1}{n})^n=e
limn→∞(1+n1)n=e
然后还是直接用导数的定义式:
f
′
(
x
)
=
lim
d
x
→
0
ln
(
x
+
d
x
)
−
ln
(
x
)
d
x
=
lim
d
x
→
0
ln
x
+
d
x
x
d
x
f'(x)=\lim_{dx\to0}\frac{\ln(x+dx)-\ln(x)}{dx}=\lim_{dx\to0}\frac{\ln\dfrac{x+dx}{x}}{dx}
f′(x)=dx→0limdxln(x+dx)−ln(x)=dx→0limdxlnxx+dx
=
1
x
lim
d
x
→
0
x
d
x
ln
x
+
d
x
x
=
1
x
lim
d
x
→
0
ln
(
x
+
d
x
x
)
x
d
x
=\frac{1}{x}\lim_{dx\to0}\frac{x}{dx}\ln\dfrac{x+dx}{x}=\frac{1}{x}\lim_{dx\to0}\ln(\dfrac{x+dx}{x})^{\dfrac{x}{dx}}
=x1dx→0limdxxlnxx+dx=x1dx→0limln(xx+dx)dxx
=
1
x
lim
x
d
x
→
∞
ln
(
x
+
d
x
x
)
x
d
x
=
1
x
lim
x
d
x
→
∞
ln
e
=
1
x
=\frac{1}{x}\lim_{\frac{x}{dx}\to\infty}\ln(\dfrac{x+dx}{x})^{\dfrac{x}{dx}}=\frac{1}{x}\lim_{\frac{x}{dx}\to\infty}\ln e=\frac{1}{x}
=x1dxx→∞limln(xx+dx)dxx=x1dxx→∞limlne=x1
证毕。
根据复合函数求导法则:
f
′
(
g
(
x
)
)
=
f
′
(
g
(
x
)
)
g
′
(
x
)
f'(g(x))=f'(g(x))g'(x)
f′(g(x))=f′(g(x))g′(x)
就有:
ln
′
(
f
(
x
)
)
=
f
′
(
x
)
f
(
x
)
=
d
ln
x
d
x
\ln'(f(x))=\frac{f'(x)}{f(x)}=\frac{d\ln x}{dx}
ln′(f(x))=f(x)f′(x)=dxdlnx
d
ln
x
=
f
′
(
x
)
f
(
x
)
d
x
d\ln x=\frac{f'(x)}{f(x)}dx
dlnx=f(x)f′(x)dx
还有一个恒等变换:
f
(
x
)
=
∫
d
f
(
x
)
f(x)=\int df(x)
f(x)=∫df(x)
大概就是函数值等于它所有的变化量加在一起。
因此有:
ln
f
(
x
)
=
∫
d
ln
f
(
x
)
\ln f(x)=\int d\ln f(x)
lnf(x)=∫dlnf(x)
再将刚才求的
d
ln
x
=
f
′
(
x
)
f
(
x
)
d
x
d\ln x=\dfrac{f'(x)}{f(x)}dx
dlnx=f(x)f′(x)dx 带入:
ln
x
=
∫
d
ln
x
=
∫
f
′
(
x
)
f
(
x
)
d
x
\ln x=\int d\ln x=\int \dfrac{f'(x)}{f(x)}dx
lnx=∫dlnx=∫f(x)f′(x)dx
因此,我们求出
f
(
x
)
f(x)
f(x) 的导数和逆,卷积后积分即可。
需要求导、积分、卷积、求逆。
时间复杂度:
O
(
n
log
n
)
O(n\log n)
O(nlogn)。
void Ln(ll *a,ll *b,int n){
static ll tmp[N],u[N],v[N];
init(n<<1);
memcpy(tmp,a,sizeof(ll)*n);
fill(tmp+n,tmp+lim,0);
memset(u,0,sizeof(ll)*lim);
memset(v,0,sizeof(ll)*lim);
dao(tmp,u,n);
inv(tmp,v,n);
fill(v+n,v+lim,0);
mul(u,v,tmp,n,lim);
jifen(tmp,tmp,n);
memcpy(b,tmp,sizeof(ll)*n);
fill(b+n,b+lim,0);
return;
}
多项式exp
有了多项式求
l
n
ln
ln,我们就可以继续用牛迭耍赖了。
设计 :
g
(
f
(
x
)
)
=
ln
(
f
0
(
x
)
)
−
h
(
x
)
g(f(x))=\ln (f_0(x))-h(x)
g(f(x))=ln(f0(x))−h(x)
ln
x
\ln x
lnx 的导数我们刚才已经求完了,是
1
x
\dfrac{1}{x}
x1
带入牛迭:
f
(
x
)
=
f
0
(
x
)
−
g
(
f
0
(
x
)
)
g
′
(
f
0
(
x
)
)
=
f
0
(
x
)
−
ln
(
f
0
(
x
)
)
−
h
(
x
)
1
f
0
(
x
)
f(x)= f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}=f_0(x)-\frac{\ln (f_0(x))-h(x)}{\dfrac{1}{f_0(x)}}
f(x)=f0(x)−g′(f0(x))g(f0(x))=f0(x)−f0(x)1ln(f0(x))−h(x)
=
f
0
(
x
)
(
1
+
h
(
x
)
−
ln
(
f
0
(
x
)
)
)
=f_0(x)(1+h(x)-\ln(f_0(x))\space)
=f0(x)(1+h(x)−ln(f0(x)) )
需要卷积和求
ln
\ln
ln
时间复杂度:
T
(
n
)
=
T
(
n
2
)
+
O
(
n
log
n
)
=
O
(
n
log
n
)
T(n)=T(\frac{n}{2})+O(n\log n)=O(n\log n)
T(n)=T(2n)+O(nlogn)=O(nlogn)
void Exp(ll *h,ll *f,int n){
static ll tmp[N];
if(n==1){f[0]=1;return;}
Exp(h,f,(n+1)>>1);
init(n<<1);
//printf(" tmp: ");for(int i=0;i<n;i++){printf("%lld ",f[i]);} putchar('\n');
Ln(f,tmp,n);
//printf(" Ln: ");for(int i=0;i<n;i++){printf("%lld ",tmp6[i]);} putchar('\n');
for(int i=0;i<n;i++) tmp[i]=(h[i]+(i==0)+mod-tmp[i])%mod;
fill(tmp+n,tmp+lim,0);
fill(f+n,f+lim,0);
//printf(" g: ");for(int i=0;i<n;i++){printf("%lld ",tmp6[i]);} putchar('\n');
NTT(f,lim,1);NTT(tmp,lim,1);
for(int i=0;i<lim;i++){
f[i]=f[i]*tmp[i]%mod;
}
NTT(f,lim,-1);
fill(f+n,f+n+lim,0);
//for(int i=0;i<n;i++){printf("%lld ",f[i]);} putchar('\n');putchar('\n');
return;
}
分治FFT
当卷积中包含自身时的操作技巧。
思想类似于CDQ分治,比较容易理解。
给出情境:
给出一个多项式 g ( x ) = ∑ i = 1 n − 1 a i x i g(x)=\sum_{i=1}^{n-1}a_ix^i g(x)=∑i=1n−1aixi,要求求出 n − 1 n-1 n−1 次多项式 f ( x ) = ∑ i = 0 n − 1 b i x i f(x)=\sum_{i=0}^{n-1}b_ix^i f(x)=∑i=0n−1bixi,使得 b 0 = 1 b_0=1 b0=1,对于 i > 0 i>0 i>0,有 b i = ∑ j = 1 i a j b i − j b_i=\sum_{j=1}^ia_jb_{i-j} bi=∑j=1iajbi−j。
考虑分治求解。
类似于CDQ分治,对于区间
(
l
,
r
)
(l,r)
(l,r) 先递归求出左侧答案,然后计算左侧对右侧的贡献,最后递归计算右侧答案。
左侧对于右侧的贡献仍然是一个卷积:
f
l
+
i
×
g
j
→
f
l
+
i
+
j
(
l
+
i
≤
m
i
d
,
m
i
d
<
l
+
i
+
j
≤
r
)
f_{l+i}\times g_j\to f_{l+i+j}(l+i\le mid,mid<l+i+j\le r)
fl+i×gj→fl+i+j(l+i≤mid,mid<l+i+j≤r)
注意一些加一减一的边界条件即可。
void solve(ll *g,ll *f,int l,int r){
static ll a[N],b[N],c[N];
if(l==r) return;
int mid=(l+r)>>1;
solve(g,f,l,mid);
int n=r-l+1,m=mid-l+1;
memcpy(a,g,sizeof(ll)*n);
memcpy(b,f+l,sizeof(ll)*m);
init(n+m);
fill(a+n,a+lim,0);
fill(b+m,b+lim,0);
fill(c,c+lim,0);
NTT(a,lim,1);NTT(b,lim,1);
for(int i=0;i<lim;i++) c[i]=a[i]*b[i]%mod;
NTT(c,lim,-1);
for(int i=1;i<=m;i++) (f[mid+i]+=c[m+i-1])%=mod;
solve(g,f,mid+1,r);
return;
}
source
完整代码
int r[N];
ll a[N],b[N],c[N];
void NTT(ll *x,int lim,int op){
for(int i=0;i<lim;i++) if(i<r[i]) swap(x[i],x[r[i]]);
for(int l=1;l<lim;l<<=1){
ll w=ksm(3,(mod-1)/(l<<1));if(op==-1) w=ksm(w,mod-2);
for(int st=0;st<lim;st+=(l<<1)){
for(ll i=0,now=1;i<l;i++,now=now*w%mod){
ll u=x[st+i],v=now*x[st+i+l]%mod;
x[st+i]=(u+v)%mod;x[st+i+l]=(u+mod-v)%mod;
}
}
}
if(op==-1){
ll ni=ksm(lim,mod-2);
for(int i=0;i<lim;i++) x[i]=x[i]*ni%mod;
}
return;
}
void Or(ll *x,int n,int op){
for(int l=1;l<n;l<<=1){
for(int st=0;st<n;st+=l*2){
for(int i=0;i<l;i++){
ll u=x[st+i],v=x[st+l+i];
if(op==1) x[st+i]=u,x[st+l+i]=(v+u)%mod;
else x[st+i]=u,x[st+l+i]=(v+mod-u)%mod;
}
}
}
}
void And(ll *x,int n,int op){
for(int l=1;l<n;l<<=1){
for(int st=0;st<n;st+=l*2){
for(int i=0;i<l;i++){
ll u=x[st+i],v=x[st+l+i];
if(op==1) x[st+i]=(u+v)%mod,x[st+l+i]=v;
else x[st+i]=(u+mod-v)%mod,x[st+l+i]=v;
}
}
}
}
//499122177
void Xor(ll *x,int n,int op){
for(int l=1;l<n;l<<=1){
for(int st=0;st<n;st+=l*2){
for(int i=0;i<l;i++){
ll u=x[st+i],v=x[st+l+i];
if(op==1) x[st+i]=(u+v)%mod,x[st+l+i]=(u+mod-v)%mod;
else x[st+i]=(u+v)%mod*499122177%mod,x[st+l+i]=(u+mod-v)%mod*499122177%mod;
}
}
}
}
void init(int n){
lim=1,L=0;
while(lim<n) lim<<=1,L++;
for(int i=1;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(L-1));
return;
}
//499122177
void inv(ll *h,ll *f,int n){
static ll tmp[N];
if(n==1){
f[0]=ksm(h[0],mod-2);return;
}
inv(h,f,(n+1)>>1);
init(n<<1);
fill(f+n,f+lim,0);
for(int i=0;i<n;i++) tmp[i]=h[i];
fill(tmp+n,tmp+lim,0);
NTT(tmp,lim,1);NTT(f,lim,1);
for(int i=0;i<lim;i++) f[i]=(2*f[i]+mod-tmp[i]*f[i]%mod*f[i]%mod)%mod;
NTT(f,lim,-1);
fill(f+n,f+lim,0);
}
void mul(ll *a,ll *b,ll *c,int n,int m){
static ll u[N],v[N];
init(n+m+1);
memcpy(u,a,sizeof(ll)*n);
memcpy(v,b,sizeof(ll)*m);
fill(u+n,u+lim,0);
fill(v+m,v+lim,0);
NTT(u,lim,1);NTT(v,lim,1);
for(int i=0;i<lim;i++) c[i]=u[i]*v[i]%mod;
NTT(c,lim,-1);
for(int i=n+m+1;i<lim;i++) c[i]=0;
return;
}
void Sqrt(ll *h,ll *f,int n){
static ll tmp[N];
if(n==1){f[0]=1;return;}
Sqrt(h,f,(n+1)>>1);
memset(tmp,0,sizeof(ll)*lim);
inv(f,tmp,n);mul(tmp,h,tmp,n,n);
init(n<<1);
NTT(f,lim,1);NTT(tmp,lim,1);
for(int i=0;i<lim;i++){
f[i]=(tmp[i]+f[i])*499122177%mod;
}
NTT(f,lim,-1);
fill(f+n,f+lim,0);
return;
}
void dao(ll *a,ll *b,int n){
static ll tmp[N];
memcpy(tmp,a,sizeof(ll)*n);
for(int i=0;i<n-1;i++) b[i]=tmp[i+1]*(i+1)%mod;
fill(b+n,b+lim,0);
return;
}
void jifen(ll *a,ll *b,int n){
static ll tmp[N];
memcpy(tmp,a,sizeof(ll)*n);
for(int i=1;i<n;i++) b[i]=tmp[i-1]*ksm(i,mod-2)%mod;
b[0]=0;
return;
}
void Ln(ll *a,ll *b,int n){
static ll tmp[N],u[N],v[N];
init(n<<1);
memcpy(tmp,a,sizeof(ll)*n);
fill(tmp+n,tmp+lim,0);
memset(u,0,sizeof(ll)*lim);
memset(v,0,sizeof(ll)*lim);
dao(tmp,u,n);
inv(tmp,v,n);
fill(v+n,v+lim,0);
mul(u,v,tmp,n,lim);
jifen(tmp,tmp,n);
memcpy(b,tmp,sizeof(ll)*n);
fill(b+n,b+lim,0);
return;
}
void Exp(ll *h,ll *f,int n){
static ll tmp[N];
if(n==1){f[0]=1;return;}
Exp(h,f,(n+1)>>1);
init(n<<1);
Ln(f,tmp,n);
for(int i=0;i<n;i++) tmp[i]=(h[i]+(i==0)+mod-tmp[i])%mod;
fill(tmp+n,tmp+lim,0);
fill(f+n,f+lim,0);
NTT(f,lim,1);NTT(tmp,lim,1);
for(int i=0;i<lim;i++){
f[i]=f[i]*tmp[i]%mod;
}
NTT(f,lim,-1);
fill(f+n,f+n+lim,0);
return;
}
void solve(ll *g,ll *f,int l,int r){
static ll a[N],b[N],c[N];
if(l==r) return;
int mid=(l+r)>>1;
solve(g,f,l,mid);
int n=r-l+1,m=mid-l+1;
memcpy(a,g,sizeof(ll)*n);
memcpy(b,f+l,sizeof(ll)*m);
init(n+m);
fill(a+n,a+lim,0);
fill(b+m,b+lim,0);
fill(c,c+lim,0);
NTT(a,lim,1);NTT(b,lim,1);
for(int i=0;i<lim;i++) c[i]=a[i]*b[i]%mod;
NTT(c,lim,-1);
for(int i=1;i<=m;i++) (f[mid+i]+=c[m+i-1])%=mod;
solve(g,f,mid+1,r);
return;
}
int main(){
return 0;
}