文章目录
头函数
#define poly vector<int>
#define bg begin
#define pb push_back
const int mod=998244353,g=3;
inline int add(int a,int b){return a+b>=mod?a+b-mod:a+b;}
inline void Add(int &a,int b){a=add(a,b);}
inline int dec(int a,int b){return a>=b?a-b:a-b+mod;}
inline void Dec(int &a,int b){a=dec(a,b);}
inline int mul(int a,int b){return 1ll*a*b>=mod?1ll*a*b%mod:a*b;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));return res;}
多项式加减点乘点除
幼儿园小朋友应该都会了吧
inline poly operator +(poly a,poly b){
poly c;int lim=max(a.size(),b.size());c.resize(lim);
for(int i=0;i<lim;i++)c[i]=add(a[i],b[i]);return c;
}
inline poly operator -(poly a,poly b){
poly c;int lim=max(a.size(),b.size());c.resize(lim);
for(int i=0;i<lim;i++)c[i]=dec(a[i],b[i]);return c;
}
inline poly operator *(poly a,int b){
for(int i=0;i<a.size();i++)Mul(a[i],b);return a;
}
inline poly operator /(poly a,int b){
for(int i=0,inv=ksm(b,mod-2);i<a.size();i++)Mul(a[i],inv);
return a;
}
多项式乘法
FFT:
前置
多项式的点值和系数表示法:
对于一个
n
n
n次多项式
f
(
x
)
f(x)
f(x)
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……a_nx^n
f(x)=a0+a1x+a2x2+a3x3……anxn被称作该多项式的系数表示
我们可以通过带任意一个
x
x
x都可以的到唯一的一个
f
(
x
)
f(x)
f(x)
但
o
i
oi
oi中一般
x
x
x一般都只是一个不定元,不会带入特定值计算,比如用作表示下标之类的
而如果我们把不同的
n
+
1
n+1
n+1带入进去得到的
n
+
1
n+1
n+1个点值叫做点值表示法
由
n
+
1
n+1
n+1个点也可以还原出唯一一个
n
n
n次多项式
虚数
即
i
,
−
1
i,\sqrt{-1}
i,−1
考虑在平面直角坐标系内
将
y
y
y轴用
i
i
i来表示
复数更准确的定义是
e i k = cos ( k ) + i sin ( k ) e^{ik}=\cos(k)+i\sin(k) eik=cos(k)+isin(k)
这样平面上一个点就是
(
a
,
b
i
)
(a,bi)
(a,bi)的形式
而2个虚数相乘,就对应平面直角坐标系上2个向量模长相乘,极角相加
由于
C
+
+
C++
C++自带的复数很慢
所以我们一般手写复数结构体
const double pi=acos(-1);
struct complex{
double x,y;
complex(double _x=0,double _y=0):x(_x),y(_y){}
friend inline complex operator +(const complex &a,const complex &b){
return complex(a.x+b.x,a.y+b.y);
}
friend inline complex operator -(const complex &a,const complex &b){
return complex(a.x-b.x,a.y-b.y);
}
friend inline complex operator /(const complex &a,const double &b){
return complex(a.x/b,a.y/b);
}
friend inline complex operator *(const complex &a,const complex &b){
return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
};
单位根
f
f
t
fft
fft的根本
n
n
n次单位根指的是满足
w
n
=
1
w^n=1
wn=1的复数
n
n
n次单位根有
n
n
n个,分别表示为
w
n
k
,
k
=
0
,
1
,
2
,
…
…
n
−
1
w_{n}^{k},k=0,1,2,……n-1
wnk,k=0,1,2,……n−1
实际上对应的将平面单位圆周上的
n
n
n个点
这些点将单位圆周均分成
n
n
n块,且构成一个正
n
n
n边形
更准确的表示为
w
n
k
=
e
2
π
i
k
/
n
w_n^k=e^{2\pi ik/n}
wnk=e2πik/n
即
w
n
k
=
cos
(
2
π
k
/
n
)
+
i
sin
(
2
π
k
/
n
)
w_n^k=\cos(2\pi k/n)+i\sin(2\pi k/n)
wnk=cos(2πk/n)+isin(2πk/n)表示
即单位圆上的
n
n
n点
单位根几个重要的性质:
1、消去引理
对于任何整数
n
≥
0
,
k
≥
0
,
d
≥
0
n\geq0,k\geq 0,d\geq 0
n≥0,k≥0,d≥0
有
w
d
n
d
k
=
w
n
k
w_{dn}^{dk}=w_{n}^{k}
wdndk=wnk
证明: w d n d k = e 2 π i d k / d n = e 2 π i k / n = w n k w_{dn}^{dk}=e^{2\pi idk/dn}=e^{2\pi ik/n}=w_{n}^{k} wdndk=e2πidk/dn=e2πik/n=wnk
2、折半引理
如果
n
≥
0
n\geq 0
n≥0且
n
n
n为偶数,那么
对所有
n
n
n次单位根平方,得到的集合是
n
/
2
n/2
n/2个
n
/
2
n/2
n/2次单位根的集合
说白了就是
w
n
k
+
n
2
=
−
w
n
k
w_{n}^{k+\frac n 2}=-w_n^{k}
wnk+2n=−wnk或者就是
w
n
n
2
=
−
1
w_n^{\frac n 2}=-1
wn2n=−1
证明:画个单位圆,
n
2
\frac n 2
2n就是旋转了
180
°
180°
180°,自然取反
3、求和引理
对于任意整数
n
≥
1
n\geq 1
n≥1和
n
̸
∣
k
n\not | k
n̸∣k
有
∑
j
=
0
n
−
1
(
w
n
k
)
j
=
0
\sum_{j=0}^{n-1}(w_{n}^k)^j=0
∑j=0n−1(wnk)j=0
考虑这其实是一个等比数列,
w
n
k
̸
=
1
w_{n}^{k}\not= 1
wnk̸=1时
原式
=
(
w
n
k
)
n
−
1
w
n
k
−
1
=
(
w
n
n
)
k
−
1
w
n
k
−
1
=
0
=\frac{(w_{n}^{k})^{n}-1}{w_{n}^{k}-1}=\frac{(w_{n}^n)^k-1}{w_{n}^{k}-1}=0
=wnk−1(wnk)n−1=wnk−1(wnn)k−1=0
而当 n ∣ k n|k n∣k时,原式显然为 n n n
这个性质在后面很重要
算法
考虑对于两个多项式
A
(
x
)
=
a
0
+
a
1
x
+
a
2
x
2
+
a
3
x
3
…
…
a
n
x
n
A(x)=a_0+a_1x+a_2x^2+a_3x^3……a_nx^n
A(x)=a0+a1x+a2x2+a3x3……anxn
B
(
x
)
=
b
0
+
b
1
x
+
b
2
x
2
+
b
3
x
3
…
…
b
m
x
m
B(x)=b_0+b_1x+b_2x^2+b_3x^3……b_mx^m
B(x)=b0+b1x+b2x2+b3x3……bmxm
我们要求一个
n
+
m
n+m
n+m次多项式
C
(
x
)
=
A
(
x
)
∗
B
(
x
)
C(x)=A(x)*B(x)
C(x)=A(x)∗B(x)
更具体的
C
(
x
)
C(x)
C(x)满足
c
i
=
∑
j
=
0
i
a
j
∗
b
i
−
j
c_i=\sum_{j=0}^{i}a_j*b_{i-j}
ci=∑j=0iaj∗bi−j
也就是2个数列倒着乘的和,所谓的卷积
考虑如果我们直接拆开暴力做是
O
(
n
2
)
O(n^2)
O(n2)的
当然也有一种分治乘法可以做到
n
l
o
g
2
3
n^{log_23}
nlog23(大概就是大常数
n
n
n\sqrt n
nn)
而 F F T FFT FFT可以做到 O ( n l o g n ) O(nlogn) O(nlogn)求出 C C C
下面为了方便假设 n = m n=m n=m, n ̸ = m n\not =m n̸=m也没有区别
考虑如果我们有
n
+
1
n+1
n+1个值
x
1
,
…
…
,
x
n
+
1
x_1,……,x_{n+1}
x1,……,xn+1
如果已经求出
A
(
x
1
)
,
…
…
A
(
x
n
+
1
)
A(x_1),……A(x_{n+1})
A(x1),……A(xn+1)
和
B
(
x
1
)
,
…
…
B
(
x
n
+
1
)
B(x_1),……B(x_{n+1})
B(x1),……B(xn+1)
也就是分别求出
x
1
,
…
…
,
x
n
+
1
x_1,……,x_n+1
x1,……,xn+1的点值
那么我们可以直接在
O
(
n
)
O(n)
O(n)的时间内求出
C
(
x
1
)
=
A
(
x
1
)
∗
B
(
x
1
)
,
…
…
,
C
(
x
n
+
1
)
=
A
(
x
n
+
1
)
∗
B
(
x
n
+
1
)
C(x_1)=A(x_1)*B(x_1),……,C(x_{n+1})=A(x_{n+1})*B(x_{n+1})
C(x1)=A(x1)∗B(x1),……,C(xn+1)=A(xn+1)∗B(xn+1)
现在我们考虑这样
先对
A
(
x
)
A(x)
A(x),
B
(
x
)
B(x)
B(x)求出
n
+
1
n+1
n+1个点的点值,乘起来得到
C
(
x
)
C(x)
C(x)的点值
又对于一个
n
n
n次多项式,如果我们知道其
n
+
1
n+1
n+1个不同点值
就一定可以还原出一个唯一的多项式
于是最后再由点值表示还原出原来的
C
(
x
)
C(x)
C(x)
注意由于实际乘出来
C
(
x
)
C(x)
C(x)最高系数是
2
n
−
1
2n-1
2n−1
所以我们需要带入
2
n
2n
2n个点求值
于是我们会将
A
(
x
)
,
B
(
x
)
A(x),B(x)
A(x),B(x)补0到
2
n
2n
2n次项
实际上由于
f
f
t
fft
fft的特殊性
我们会将项数补充到2的整数次幂次
当然非二的整数次幂项的多项式乘法也是可以做的(见下面补充的混合基 f f t fft fft和 B l u e s t e i n Bluestein Bluestein)
第一步将系数转成点值是正变换,称作
D
F
T
DFT
DFT,单次复杂度
O
(
n
l
o
g
n
)
O(nlogn)
O(nlogn)
由点值还原系数为逆变换,称作
I
D
F
T
IDFT
IDFT,单次复杂度
O
(
n
l
o
g
n
)
O(nlogn)
O(nlogn)
于是总复杂度就是
O
(
n
l
o
g
n
)
O(nlogn)
O(nlogn)
整个操作被称为快速傅里叶变换
(
F
F
T
)
(FFT)
(FFT)
DFT:
考虑我们带入
n
n
n次单位根
w
n
w_n
wn
A
(
w
n
k
)
=
∑
i
=
0
n
a
i
(
w
n
k
)
i
A(w_n^k)=\sum_{i=0}^na_i(w_n^k)^i
A(wnk)=∑i=0nai(wnk)i
考虑将下标按照奇偶分类
A
(
w
n
k
)
=
∑
i
=
0
n
2
−
1
a
2
i
(
w
n
k
)
2
i
+
∑
i
=
0
n
2
−
1
a
2
i
+
1
(
w
n
k
)
2
i
+
1
A(w_n^k)=\sum_{i=0}^{\frac n 2-1}a_{2i}(w_n^{k})^{2i}+\sum_{i=0}^{\frac n2 -1}a_{2i+1}(w_{n}^k)^{2i+1}
A(wnk)=∑i=02n−1a2i(wnk)2i+∑i=02n−1a2i+1(wnk)2i+1
= ∑ i = 0 n 2 − 1 a 2 i ( w n 2 k i ) + w n k ∑ i = 0 n 2 − 1 a 2 i + 1 ( w n 2 k i ) =\sum_{i=0}^{\frac n 2-1}a_{2i}(w_n^{2ki})+w_{n}^k\sum_{i=0}^{\frac n2 -1}a_{2i+1}(w_{n}^{2ki}) =∑i=02n−1a2i(wn2ki)+wnk∑i=02n−1a2i+1(wn2ki)
= ∑ i = 0 n 2 − 1 a 2 i ( w n 2 k ) + w n k ∑ i = 0 n 2 − 1 a 2 i + 1 ( w n 2 k ) =\sum_{i=0}^{\frac n 2-1}a_{2i}(w_{\frac n 2}^{k})+w_n^k\sum_{i=0}^{\frac n2 -1}a_{2i+1}(w_{\frac n 2}^k) =∑i=02n−1a2i(w2nk)+wnk∑i=02n−1a2i+1(w2nk)
= A o ( w n 2 k ) + w n k A 1 ( w n 2 k ) =A_o(w_{\frac n 2}^k)+{w_{n}^k}A_1(w_{\frac n2 }^k) =Ao(w2nk)+wnkA1(w2nk)
这样
A
0
A_0
A0和
A
1
A_1
A1都只有
n
2
\frac n2
2n项了
可以继续递归去做
尽管现在复杂度并没有变化
考虑对于
A
(
w
n
k
+
n
2
)
A(w_n^{k+\frac n 2})
A(wnk+2n)
A
(
w
n
k
+
n
2
)
=
∑
i
=
0
n
2
−
1
a
2
i
(
w
n
k
+
n
2
)
2
i
+
w
n
k
+
n
2
∑
i
=
0
n
2
−
1
a
2
i
+
1
(
w
n
k
+
n
2
)
2
i
A(w_n^{k+\frac n 2})=\sum_{i=0}^{\frac n 2-1}a_{2i}(w_n^{k+\frac n 2})^{2i}+w_{n}^{k+\frac n 2}\sum_{i=0}^{\frac n2 -1}a_{2i+1}(w_{n}^{k+\frac n2 })^{2i}
A(wnk+2n)=∑i=02n−1a2i(wnk+2n)2i+wnk+2n∑i=02n−1a2i+1(wnk+2n)2i
考虑单位根的消去引理
w n k + n 2 = − w n k w_{n}^{k+\frac n 2}=-w_n^k wnk+2n=−wnk
则 A ( w n k + n 2 ) = ∑ i = 0 n 2 − 1 a 2 i ( − w n k ) 2 i + ( − w n k ) ∑ i = 0 n 2 − 1 a 2 i + 1 ( − w n k ) 2 i A(w_n^{k+\frac n 2})=\sum_{i=0}^{\frac n 2-1}a_{2i}(-w_n^k)^{2i}+(-w_n^k)\sum_{i=0}^{\frac n2 -1}a_{2i+1}(-w_n^k)^{2i} A(wnk+2n)=∑i=02n−1a2i(−wnk)2i+(−wnk)∑i=02n−1a2i+1(−wnk)2i
= ∑ i = 0 n 2 − 1 a 2 i ( w n k ) 2 i − w n k ∑ i = 0 n 2 − 1 a 2 i + 1 ( w n k ) 2 i =\sum_{i=0}^{\frac n 2-1}a_{2i}(w_n^k)^{2i}-w_n^k\sum_{i=0}^{\frac n2 -1}a_{2i+1}(w_n^k)^{2i} =∑i=02n−1a2i(wnk)2i−wnk∑i=02n−1a2i+1(wnk)2i
= ∑ i = 0 n 2 − 1 a 2 i ( w n 2 k ) − w n k ∑ i = 0 n 2 − 1 a 2 i + 1 ( w n 2 k ) =\sum_{i=0}^{\frac n 2-1}a_{2i}(w_{\frac n 2}^{k})-w_n^k\sum_{i=0}^{\frac n2 -1}a_{2i+1}(w_{\frac n 2}^k) =∑i=02n−1a2i(w2nk)−wnk∑i=02n−1a2i+1(w2nk)
= A o ( w n 2 k ) − w n k A 1 ( w n 2 k ) =A_o(w_{\frac n 2}^k)-{w_{n}^k}A_1(w_{\frac n2 }^k) =Ao(w2nk)−wnkA1(w2nk)
我们发现唯一不同的就是第二项的符号
也就是说如果我们知道
A
o
(
w
n
2
k
)
A_o(w_{\frac n 2}^k)
Ao(w2nk)和
A
1
(
w
n
2
k
)
A_1(w_{\frac n2 }^k)
A1(w2nk)
我们就可以同时知道
A
(
w
n
k
+
n
2
)
A(w_n^{k+\frac n 2})
A(wnk+2n)和
A
(
w
n
k
)
A(w_n^k)
A(wnk)
考虑对于
k
∈
[
0
,
n
−
1
]
,
k\in[0,n-1],
k∈[0,n−1],我们求
A
(
w
n
k
)
A(w_n^k)
A(wnk)
就只需要知道
k
∈
[
0
,
n
2
−
1
]
,
A
0
(
w
n
2
k
)
k\in[0,\frac n2-1 ],A_0(w_{\frac n2 }^k)
k∈[0,2n−1],A0(w2nk)和
A
1
(
w
n
2
k
)
A_1(w_{\frac n 2}^k)
A1(w2nk)
就可以在
O
(
n
)
O(n)
O(n)的时间内得到
A
A
A
而
A
0
A_0
A0,
A
1
A_1
A1系数都只有
n
2
\frac n 2
2n个,所以规模只有原来的一般
递归求解即可
时间复杂度 T ( n ) = 2 ∗ T ( n 2 ) + O ( n ) = O ( n l o g n ) T(n)=2*T(\frac n2 )+O(n)=O(nlogn) T(n)=2∗T(2n)+O(n)=O(nlogn)
这里也是之所以要把项数补充到
2
k
2^k
2k
因为每一次都把
n
n
n项分成
n
2
\frac n 2
2n项
如果
n
n
n是奇数,那就没法分开了
由于递归常数比较大,而一般
f
f
t
fft
fft有关的题时间瓶颈就在
D
F
T
DFT
DFT上面
D
F
T
DFT
DFT不知道为什么 常数也很大
写的差的
f
f
t
fft
fft甚至可以跑
1
e
6
1e6
1e6的数据
T
T
T掉
所以我们考虑迭代做
由于每个数最终在的位置和原来不一样
所以我们要预处理出最终的位置上
据说找规律得到了
O
(
n
)
O(n)
O(n)预处理的方法
代码如下:
没看懂,选择全文背诵
当然也可以
n
l
o
g
n
nlogn
nlogn模拟最终位置
int rev[N<<2];
inline void init(int lim){
for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
}
我们先把每个数放到最终应该在的地方然后一步步迭代回去就是了
inline void DFT(complex f[],int lim){
for(int i=0;i<lim;i++)if(rev[i]>i)swap(f[i],f[rev[i]]);
for(int mid=1;mid<lim;mid<<=1){
complex now=plx(cos(pi/mid),sin(pi/mid));
for(int i=0;i<lim;i+=(mid<<1)){
complex w=plx(1,0);
for(int j=0;j<mid;j++,w=w*now){
plx p0=f[i+j],p1=w*f[i+j+mid];
f[i+j]=p0+p1,f[i+j+mid]=p0-p1;
}
}
}
}
IDFT:
以下摘抄自 m i s k c o o miskcoo miskcoo神犇( m a r k d o w n markdown markdown太难写了QAQ)
代码实现
const double pi=acos(-1);
struct plx{
double x,y;
plx(double _x=0,double _y=0):x(_x),y(_y){}
friend inline plx operator +(const plx &a,const plx &b){
return plx(a.x+b.x,a.y+b.y);
}
friend inline plx operator -(const plx &a,const plx &b){
return plx(a.x-b.x,a.y-b.y);
}
friend inline plx operator /(const plx &a,const double &b){
return plx(a.x/b,a.y/b);
}
friend inline plx operator *(const plx &a,const plx &b){
return plx(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
};
inline void fft(plx f[],int lim,int kd){//kd表示在做正变换还是逆变换
for(int i=0;i<lim;i++)if(rev[i]>i)swap(f[i],f[rev[i]]);
for(int mid=1;mid<lim;mid<<=1){
plx now=plx(cos(pi/mid),kd*sin(pi/mid));
for(int i=0;i<lim;i+=(mid<<1)){
plx w=plx(1,0);
for(int j=0;j<mid;j++,w=w*now){
plx p0=f[i+j],p1=w*f[i+j+mid];
f[i+j]=p0+p1,f[i+j+mid]=p0-p1;
}
}
}
if(kd==-1)for(int i=0;i<lim;i++)f[i]=f[i]/lim;
}
NTT:
由于
f
f
t
fft
fft涉及复数和实数运算,实际会出现精度误差
在整数运算时难免会出锅
所以我们考虑一个能在模意义下的变换
这就是
n
t
t
ntt
ntt
首先引入原根的概念
对于一个素数
p
p
p,其原根
g
g
g定义为满足
g
0
,
g
1
,
…
…
,
g
p
−
2
g^0,g^1,……,g^{p-2}
g0,g1,……,gp−2互不相同的数
又由于费马定理,对一个素数
p
p
p,有
a
p
−
1
≡
1
(
m
o
d
p
)
a^{p-1}\equiv 1(mod\ p)
ap−1≡1(mod p)
这个和单位根很相似
我们考虑单位根之所以能够做
f
f
t
fft
fft,是因为
w
w
w的三个性质
考虑如果我们对于一个素数
p
=
2
n
∗
k
+
1
p=2^n*k+1
p=2n∗k+1,我们令
g
n
=
g
k
g_n=g^k
gn=gk
这样我们就可以满足
g
n
0
,
g
n
1
,
g
n
2
…
…
g
n
n
−
1
g_n^0,g_n^1,g_n^2……g_n^{n-1}
gn0,gn1,gn2……gnn−1互不相同且满足这些性质(不想证了)
但这也就限制了
n
t
t
ntt
ntt的模数必须是
k
∗
2
n
+
1
k*2^n+1
k∗2n+1的形式
否则就要做任意模数
n
t
t
ntt
ntt
做法就是选出几个模数分别做之后
C
r
t
Crt
Crt合并答案(并不会)
代码和
f
f
t
fft
fft类似
由于
w
n
−
k
=
w
n
n
−
k
w_n^{-k}=w_n^{n-k}
wn−k=wnn−k
所以我们可以先做的时候不管正逆变换
最后把
1
1
1~
n
−
1
n-1
n−1反转一下就可以了
代码实现
inline void ntt(poly &f,int lim,int kd){
for(int i=0;i<lim;i++)if(i<rev[i])swap(f[i],f[rev[i]]);
for(int mid=1;mid<lim;mid<<=1){
int now=ksm(g,(mod-1)/(mid<<1));
for(int i=0;i<lim;i+=(mid<<1)){
int w=1;
for(int j=0;j<mid;j++,w=mul(w,now)){
int a0=f[i+j],a1=mul(w,f[i+j+mid]);
f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
}
}
}
if(kd==-1&&(reverse(f.bg()+1,f.bg()+lim),1))for(int i=0,inv=ksm(lim,mod-2);i<lim;i++)f[i]=mul(f[i],inv);
}
由于
n
t
t
ntt
ntt中每次乘起来很耗费常数
我也不知道为什么,但就是很耗时间
于是我们可以预处理出原根优化常数
预处理原根
const int N=1000005,C=21;
int *w[22];
inline void init_w(){
for(int i=1;i<=C;i++)
w[i]=new int[1<<(i-1)];
int wn=ksm(g,(mod-1)/(1<<C));
w[C][0]=1;
for(int i=1;i<(1<<(C-1));i++)w[C][i]=mul(w[C][i-1],wn);
for(int i=C-1;i;i--)
for(int j=0;j<(1<<(i-1));j++)
w[i][j]=w[i+1][j<<1];
}
速度比不预处理快了差不多 25 % 25\% 25%到 45 % 45\% 45%不等
其实
n
t
t
ntt
ntt本身常数不算很大,运算常数大概也只有5、6左右
不过下标不连续可能会导致慢一些
inline void ntt(poly &f,int lim,int kd){
for(int i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
for(int mid=1,l=1;mid<lim;mid<<=1,l++)
for(int i=0;i<lim;i+=(mid<<1))
for(int j=0,a0,a1;j<mid;j++)
a0=f[i+j],a1=mul(f[i+j+mid],w[l][j]),
f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
if(kd==-1&&(reverse(f.begin()+1,f.begin()+lim),1))
for(int inv=ksm(lim,mod-2),i=0;i<lim;i++)Mul(f[i],inv);
}
乘法
可以在比较小的时候暴力加循环展开 优化常数
inline poly operator *(poly a,poly b){
int deg=a.size()+b.size()-1,lim=1;
if(deg<=128){
poly c(deg,0);
for(int i=0;i<a.size();i++)
for(int j=0;j<b.size();j++)
Add(c[i+j],mul(a[i],b[j]));
return c;
}
while(lim<deg)lim<<=1;init(lim);
a.resize(lim),ntt(a,lim,1);
b.resize(lim),ntt(b,lim,1);
for(int i=0;i<lim;i++)Mul(a[i],b[i]);
ntt(a,lim,-1),a.resize(deg);
return a;
}
多项式求逆:
已知一个 n − 1 n-1 n−1次多项式 A ( x ) A(x) A(x),求多项式 B ( x ) B(x) B(x)满足:
A ( x ) B ( x ) ≡ 1 ( m o d x n ) A(x)B(x)\equiv 1(mod\ x^n) A(x)B(x)≡1(mod xn)
求解过程
倍增:
若已知 A ( x ) B ( x ) ′ ≡ 1 ( m o d x ⌈ n 2 ⌉ ) A(x)B(x)'\equiv 1(mod \ x^{\lceil \frac n 2\rceil}) A(x)B(x)′≡1(mod x⌈2n⌉)
A ( x ) B ( x ) ′ − A ( x ) B ( x ) ≡ 0 ( m o d x ⌈ n 2 ⌉ ) A(x)B(x)'-A(x)B(x)\equiv 0(mod \ x^{\lceil \frac n 2\rceil}) A(x)B(x)′−A(x)B(x)≡0(mod x⌈2n⌉)
B ( x ) ′ − B ( x ) ≡ 0 ( m o d x ⌈ n 2 ⌉ ) B(x)'-B(x)\equiv 0(mod \ x^{\lceil \frac n 2\rceil}) B(x)′−B(x)≡0(mod x⌈2n⌉)
平方:
B
(
x
)
′
2
+
B
(
x
)
2
−
2
B
(
x
)
′
B
(
x
)
≡
0
(
m
o
d
x
n
)
B(x)'^2+B(x)^2-2B(x)'B(x)\equiv 0(mod \ x^{ n})
B(x)′2+B(x)2−2B(x)′B(x)≡0(mod xn)
A
(
x
)
(
B
(
x
)
′
2
+
B
(
x
)
2
−
2
B
(
x
)
′
B
(
x
)
)
≡
0
(
m
o
d
x
n
)
A(x)(B(x)'^2+B(x)^2-2B(x)'B(x))\equiv 0(mod \ x^{ n})
A(x)(B(x)′2+B(x)2−2B(x)′B(x))≡0(mod xn)
B
(
x
)
≡
2
B
(
x
)
′
−
A
(
x
)
B
(
x
)
′
2
(
m
o
d
x
n
)
B(x)\equiv 2B(x)'-A(x)B(x)'^2(mod\ x^n)
B(x)≡2B(x)′−A(x)B(x)′2(mod xn)
复杂度 T ( n ) = T ( n 2 ) + O ( n l o g n ) = O ( n l o g n ) T(n)=T(\frac n 2)+O(nlogn)=O(nlogn) T(n)=T(2n)+O(nlogn)=O(nlogn)
注意次数,最高到3倍,开的4倍
代码实现
inline poly Inv(poly a,int deg){
poly c,b(1,ksm(a[0],mod-2));
for(int lim=4;lim<(deg<<2);lim<<=1){
init(lim);
c=a,c.resize(lim>>1);
c.resize(lim),ntt(c,lim,1);
b.resize(lim),ntt(b,lim,1);
for(int i=0;i<lim;i++)Mul(b[i],dec(2,mul(b[i],c[i])));
ntt(b,lim,-1),b.resize(lim>>1);
}b.resize(deg);return b;
}
多项式开方:
已知一个 n − 1 n-1 n−1次多项式 A ( x ) A(x) A(x),求一个 m o d x n mod\ x^n mod xn意义下的多项式 B ( x ) B(x) B(x)满足
B ( x ) 2 ≡ A ( x ) ( m o d x n ) B(x)^2\equiv A(x)(mod\ x^n) B(x)2≡A(x)(mod xn)
满足 A [ 0 ] = 1 A[0]=1 A[0]=1
求解过程
倍增
首先当
n
=
1
n=1
n=1时要满足
A
[
0
]
=
1
A[0]=1
A[0]=1(否则要二次剩余解,老子不会)
假设已知
B
(
x
)
′
2
≡
A
(
x
)
(
m
o
d
x
⌈
n
2
⌉
)
B(x)'^2\equiv A(x)(mod\ x^{\lceil \frac n 2\rceil})
B(x)′2≡A(x)(mod x⌈2n⌉)
由于
⌈
n
2
⌉
\lceil\frac n 2\rceil
⌈2n⌉以上的项是不会有影响的
所以
B
(
x
)
′
≡
B
(
x
)
(
m
o
d
x
⌈
n
2
⌉
)
B(x)'\equiv B(x) (mod\ x^{\lceil \frac n 2\rceil})
B(x)′≡B(x)(mod x⌈2n⌉)
移项平方得:
B
(
x
)
2
+
B
(
x
)
′
2
−
2
B
(
x
)
B
(
x
)
′
≡
0
(
m
o
d
x
n
)
B(x)^2+B(x)'^2-2B(x)B(x)'\equiv 0(mod \ x^n)
B(x)2+B(x)′2−2B(x)B(x)′≡0(mod xn)
A
(
x
)
+
B
(
x
)
′
2
≡
2
B
(
x
)
B
(
x
)
′
(
m
o
d
x
n
)
A(x)+B(x)'^2\equiv 2B(x)B(x)'(mod \ x^n)
A(x)+B(x)′2≡2B(x)B(x)′(mod xn)
B
(
x
)
≡
A
(
x
)
2
B
(
x
)
′
+
B
(
x
)
2
(
m
o
d
x
n
)
B(x)\equiv \frac{A(x)}{2B(x)'}+\frac{B(x)}{2}(mod\ x^n)
B(x)≡2B(x)′A(x)+2B(x)(mod xn)
复杂度 T ( n ) = T ( n 2 ) + O ( n l o g n ) = O ( n l o g n ) T(n)=T(\frac n 2)+O(nlogn)=O(nlogn) T(n)=T(2n)+O(nlogn)=O(nlogn)
代码实现
inline poly Sqrt(poly a,int deg){
poly b(1,1),c,d;
for(int lim=4;lim<(deg<<2);lim<<=1){
c=a,c.resize(lim>>1);
init(lim),d=Inv(b,lim>>1),
c.resize(lim),ntt(c,lim,1);
d.resize(lim),ntt(d,lim,1);
for(int i=0;i<lim;i++)Mul(c[i],d[i]);
ntt(c,lim,-1),b.resize(lim>>1);
for(int i=0;i<(lim>>1);i++)b[i]=mul(inv2,add(b[i],c[i]));
}b.resize(deg);return b;
}
多项式除法和取模:
给定一个
n
n
n次多项式
A
(
x
)
A(x)
A(x)和一个
m
m
m次多项式
B
(
x
)
B(x)
B(x)
求一个
n
−
m
n-m
n−m次多项式
C
(
x
)
C(x)
C(x)和
m
−
1
m-1
m−1次多项式
D
(
x
)
D(x)
D(x)满足:
A
(
x
)
=
B
(
x
)
C
(
x
)
+
D
(
x
)
A(x)=B(x)C(x)+D(x)
A(x)=B(x)C(x)+D(x)
求解过程
考虑对一个最高次数为
n
n
n的多项式
B
(
x
)
B(x)
B(x)操作
B
R
(
x
)
=
x
n
B
(
1
x
)
B^R(x)=x^nB(\frac 1 x)
BR(x)=xnB(x1)
会发现
B
R
(
x
)
B^R(x)
BR(x)只是
B
(
x
)
B(x)
B(x)的系数反转的柿子
考虑
A
(
x
)
=
B
(
x
)
C
(
x
)
+
D
(
x
)
A(x)= B(x)C(x)+D(x)
A(x)=B(x)C(x)+D(x)
A
(
x
)
A(x)
A(x)最高项为
n
n
n,
B
(
x
)
B(x)
B(x)最高项为
m
m
m,则
C
(
x
)
C(x)
C(x)最高项为
n
−
m
n-m
n−m,
D
(
x
)
D(x)
D(x)为
m
−
1
m-1
m−1
两边同时乘一个
x
n
x^n
xn,并带入
1
x
\frac 1 x
x1
x
n
A
(
x
)
=
x
m
B
(
x
)
x
n
−
m
C
(
x
)
+
x
n
−
m
+
1
∗
x
m
−
1
D
(
x
)
x^nA(x)= x^mB(x)x^{n-m}C(x)+x^{n-m+1}*x^{m-1}D(x)
xnA(x)=xmB(x)xn−mC(x)+xn−m+1∗xm−1D(x)
A
R
(
x
)
=
B
R
(
x
)
C
R
(
x
)
+
x
n
−
m
+
1
D
R
(
x
)
A^R(x)=B^R(x)C^R(x)+x^{n-m+1}D^R(x)
AR(x)=BR(x)CR(x)+xn−m+1DR(x)
考虑在
m
o
d
x
n
−
m
+
1
mod\ x^{n-m+1}
mod xn−m+1意义下,
A
,
B
A,B
A,B已知,
C
C
C最高为
n
−
m
n-m
n−m不影响,而
D
D
D被消去
则
A
R
(
x
)
≡
B
R
(
x
)
C
R
(
x
)
(
m
o
d
x
n
−
m
+
1
)
A^R(x)\equiv B^R(x)C^R(x)(mod\ x^{n-m+1})
AR(x)≡BR(x)CR(x)(mod xn−m+1)
多项式求逆就可以了
复杂度 O ( n l o g n ) O(nlogn) O(nlogn)
代码实现
inline poly operator /(poly a,poly b){
int lim=1,deg=a.size()-b.size()+1;
reverse(a.bg(),a.end());
reverse(b.bg(),b.end());
while(lim<=deg)lim<<=1;
b=Inv(b,lim),b.resize(deg);
a=a*b,a.resize(deg);
reverse(a.bg(),a.end());
return a;
}
inline poly operator %(poly a,poly b){
poly c=a-(a/b)*b;
c.resize(b.size()-1);
return c;
}
多项式求导与积分
我怕是自学的是一个假的微积分
假装自己会的差不多了
复杂度 O ( n ) O(n) O(n)
代码实现
inline poly deriv(poly a){
for(int i=0;i<a.size()-1;i++)a[i]=mul(a[i+1],i+1);
a.pob();return a;
}
inline poly integ(poly a){
a.pb(0);
for(int i=a.size()-1;i;i--)a[i]=mul(a[i-1],inv[i]);
a[0]=0;
return a;
}
多项式Ln
已知一个
n
−
1
n-1
n−1次多项式
g
(
x
)
g(x)
g(x),求一个
m
o
d
x
n
mod\ x^n
mod xn意义下的多项式
f
(
x
)
f(x)
f(x),满足:
f
(
x
)
≡
L
n
(
g
(
x
)
)
(
m
o
d
x
n
)
f(x)\equiv Ln(g(x))(mod\ x^n)
f(x)≡Ln(g(x))(mod xn)
求解过程
由于有公式
若
f
(
x
)
=
l
n
(
g
(
x
)
)
f(x)=ln(g(x))
f(x)=ln(g(x))
则
g
′
(
x
)
=
f
′
(
x
)
g
(
x
)
g'(x)=f'(x)g(x)
g′(x)=f′(x)g(x)
f
′
(
x
)
=
g
′
(
x
)
g
(
x
)
f'(x)=\frac{g'(x)}{g(x)}
f′(x)=g(x)g′(x)
且要满足
g
[
0
]
=
1
g[0]=1
g[0]=1否则老子不会
求导求逆最后再积分就可以了
复杂度 O ( n l o g n ) O(nlogn) O(nlogn)
代码实现
/*
if f(x)=Ln(g(x))
then g'(x)=f'(x)g(x)
*/
inline poly Ln(poly a,int lim){
a=integ(deriv(a)*Inv(a,lim)),a.resize(lim);
return a;
}
多项式Exp
前置知识
泰勒展开
f ( x ) = f ( x 0 ) + f 1 ( x 0 ) 1 ! ( x − x 0 ) + f 2 ( x 0 ) 2 ! ( x − x 0 ) 2 + … … + f n ( x 0 ) n ! ( x − x 0 ) n + … … f(x)=f(x_0)+\frac{f^1{(x_0)}}{1!}(x-x_0)+\frac{f^2(x_0)}{2!}(x-x_0)^2+……+\frac{f^n(x_0)}{n!}(x-x_0)^n+…… f(x)=f(x0)+1!f1(x0)(x−x0)+2!f2(x0)(x−x0)2+……+n!fn(x0)(x−x0)n+……
考虑我们要构造一个函数
g
(
x
)
g(x)
g(x)使
g
(
x
)
g(x)
g(x)完全拟合
f
(
x
)
f(x)
f(x)
那首先初始点
x
0
x_0
x0的值要和
f
(
x
0
)
f(x_0)
f(x0)一样
在此基础上只需要保证一阶导数,二阶导数……都完全相同即可
即
∀
i
,
f
i
(
x
)
=
g
i
(
x
)
\forall i,f^i(x)=g^i(x)
∀i,fi(x)=gi(x)
由于求
g
g
g第
i
i
i阶导数时为
i
!
a
i
=
f
i
(
x
)
i!a_i=f^i(x)
i!ai=fi(x)
即
a
i
=
f
i
(
x
)
i
!
a_i=\frac{f^i(x)}{i!}
ai=i!fi(x)
所以得证
牛顿迭代
在多项式中一般用来解这类问题:
假设有函数
f
f
f和一个多项式
g
(
x
)
m
o
d
x
n
g(x)mod\ x^n
g(x)mod xn
满足
f
(
g
(
x
)
)
≡
0
(
m
o
d
x
n
)
f(g(x))\equiv 0 ( mod\ x^n)
f(g(x))≡0(mod xn)
已知
f
f
f,求
g
g
g
说白了就是用来解 f ( g ( x ) ) ≡ 0 ( m o d x n ) f(g(x))\equiv 0(mod\ x^n) f(g(x))≡0(mod xn)之类的方程
首先在 n = 1 n=1 n=1的时候即常数时单独求
假设已经知道
f
(
g
′
(
x
)
)
≡
0
(
m
o
d
x
⌈
n
2
⌉
)
f(g'(x))\equiv 0(mod\ x^{\lceil \frac n2 \rceil})
f(g′(x))≡0(mod x⌈2n⌉)
要求
f
(
g
(
x
)
)
≡
0
(
m
o
d
x
n
)
f(g(x))\equiv0(mod \ x^n)
f(g(x))≡0(mod xn)
考虑将
f
(
g
(
x
)
)
f(g(x))
f(g(x))在
g
′
(
x
)
g'(x)
g′(x)处泰勒展开
f
(
g
(
x
)
)
=
f
(
g
′
(
x
)
)
+
f
1
(
g
′
(
x
)
)
1
!
(
g
(
x
)
−
g
′
(
x
)
)
+
f
2
(
g
′
(
x
)
)
2
!
(
g
(
x
)
−
g
(
x
)
′
)
2
+
…
…
f(g(x))=f(g'(x))+\frac{f^1(g'(x))}{1!}(g(x)-g'(x))+\frac{f^2(g'(x))}{2!}(g(x)-g(x)')^2+……
f(g(x))=f(g′(x))+1!f1(g′(x))(g(x)−g′(x))+2!f2(g′(x))(g(x)−g(x)′)2+……
首先显然有
g
(
x
)
≡
g
′
(
x
)
(
m
o
d
x
⌈
n
2
⌉
)
g(x)\equiv g'(x)(mod \ x^{\lceil \frac n 2 \rceil})
g(x)≡g′(x)(mod x⌈2n⌉)
所以
g
(
x
)
−
g
′
(
x
)
g(x)-g'(x)
g(x)−g′(x)最低项次数一定大于
⌈
n
2
⌉
\lceil \frac n 2\rceil
⌈2n⌉
则在
m
o
d
x
n
mod \ x^n
mod xn意义下,整个式子从
(
g
(
x
)
−
g
(
x
)
′
)
2
(g(x)-g(x)')^2
(g(x)−g(x)′)2开始都为
0
0
0了
则
f
(
g
(
x
)
)
≡
f
(
g
′
(
x
)
)
+
f
1
(
g
′
(
x
)
)
(
g
(
x
)
−
g
′
(
x
)
)
(
m
o
d
x
n
)
f(g(x))\equiv f(g'(x))+{f^1(g'(x))}(g(x)-g'(x))(mod\ x^n)
f(g(x))≡f(g′(x))+f1(g′(x))(g(x)−g′(x))(mod xn)
又 f ( g ( x ) ) ≡ 0 ( m o d x n ) f(g(x))\equiv 0 (mod\ x^n) f(g(x))≡0(mod xn)
g ( x ) ≡ g ′ ( x ) − f ( g ′ ( x ) ) f 1 ( g ′ ( x ) ) ( m o d x n ) g(x)\equiv g'(x)-\frac{f(g'(x))}{f^1(g'(x))}(mod\ x^n) g(x)≡g′(x)−f1(g′(x))f(g′(x))(mod xn)
这就大功告成了
例:
比如多项式开根
就是解
B
(
x
)
2
−
A
(
x
)
≡
0
(
m
o
d
x
n
)
B(x)^2-A(x)\equiv 0(mod\ x^n)
B(x)2−A(x)≡0(mod xn)
假设已知
B
′
(
x
)
2
−
A
(
x
)
≡
0
(
m
o
d
x
⌈
n
2
⌉
)
B'(x)^2-A(x)\equiv 0(mod\ x^{\lceil \frac n 2\rceil})
B′(x)2−A(x)≡0(mod x⌈2n⌉)
这时候
f
(
B
(
x
)
)
=
B
(
x
)
2
−
A
(
x
)
,
f(B(x))=B(x)^2-A(x),
f(B(x))=B(x)2−A(x),则
f
1
(
B
(
x
)
)
=
2
B
(
x
)
f^1(B(x))=2B(x)
f1(B(x))=2B(x)
带入
B
(
x
)
≡
B
′
(
x
)
−
B
′
(
x
)
2
−
A
(
x
)
2
B
′
(
x
)
≡
B
′
(
x
)
2
+
A
(
x
)
2
B
′
(
x
)
B(x)\equiv B'(x)-\frac{B'(x)^2-A(x)}{2B'(x)}\equiv \frac{B'(x)^2+A(x)}{2B'(x)}
B(x)≡B′(x)−2B′(x)B′(x)2−A(x)≡2B′(x)B′(x)2+A(x)
就是我们推出来的式子
已知一个
n
−
1
n-1
n−1次多项式
g
(
x
)
g(x)
g(x),求一个
m
o
d
x
n
mod\ x^n
mod xn意义下的多项式
f
(
x
)
f(x)
f(x),满足:
f
(
x
)
≡
e
g
(
x
)
(
m
o
d
x
n
)
f(x)\equiv e^{g(x)}(mod\ x^n)
f(x)≡eg(x)(mod xn)
也就是 L n ( f ( x ) ) ≡ g ( x ) ( m o d x n ) Ln(f(x))\equiv g(x)(mod\ x^n) Ln(f(x))≡g(x)(mod xn)
求解过程
倍增:
考虑原问题,即求解方程 L n ( f ( x ) ) − g ( x ) ≡ 0 ( m o d x n ) Ln(f(x))-g(x)\equiv 0(mod\ x^n) Ln(f(x))−g(x)≡0(mod xn)
同样假设已经知道 L n ( f ′ ( x ) ) − g ( x ) ≡ 0 ( m o d x ⌈ n 2 ⌉ ) Ln(f'(x))-g(x)\equiv 0(mod\ x^{\lceil \frac n 2\rceil}) Ln(f′(x))−g(x)≡0(mod x⌈2n⌉)
令 P ( f ( x ) ) = L n ( f ( x ) ) − g ( x ) P(f(x))=Ln(f(x))-g(x) P(f(x))=Ln(f(x))−g(x),则 P 1 ( f ( x ) ) = 1 f ( x ) P^1(f(x))=\frac 1 {f(x)} P1(f(x))=f(x)1
则 f ( x ) ≡ f ′ ( x ) − L n ( f ′ ( x ) ) − g ( x ) 1 f ′ ( x ) ≡ f ′ ( x ) − f ′ ( x ) L n ( f ′ ( x ) ) + f ′ ( x ) g ( x ) f(x)\equiv f'(x)-\frac{Ln(f'(x))-g(x)}{\frac{1}{f'(x)}}\equiv f'(x)-f'(x)Ln(f'(x))+f'(x)g(x) f(x)≡f′(x)−f′(x)1Ln(f′(x))−g(x)≡f′(x)−f′(x)Ln(f′(x))+f′(x)g(x)
≡ f ′ ( x ) ( 1 − L n ( f ′ ( x ) ) + g ( x ) ) ( m o d x n ) \equiv f'(x)(1-Ln(f'(x))+g(x)) (mod\ x^n) ≡f′(x)(1−Ln(f′(x))+g(x))(mod xn)
复杂度 T ( n ) = T ( n 2 ) + O ( n l o g n ) = O ( n l o g n ) T(n)=T(\frac n 2)+O(nlogn)=O(nlogn) T(n)=T(2n)+O(nlogn)=O(nlogn)
代码实现
inline poly exp(poly a,int deg){
poly b(1,1),c;a.resize(deg<<1);
for(int lim=2;lim<(deg<<1);lim<<=1){
c=Ln(b,lim);
for(int i=0;i<lim;i++)c[i]=dec(a[i],c[i]);
Add(c[0],1),b=b*c;
b.resize(lim);
}b.resize(deg);return b;
}
多项式多点求值
已知一个 n n n次多项式 f ( x ) f(x) f(x),求 f ( a 1 ) … … f ( a m ) f(a_1)……f(a_m) f(a1)……f(am)
求解过程
考虑构造函数
P
(
x
)
=
∏
i
=
1
m
(
x
−
a
i
)
P(x)=\prod_{i=1}^{m}(x-a_i)
P(x)=∏i=1m(x−ai)
显然
∀
i
∈
[
1
,
m
]
,
P
(
a
i
)
=
0
\forall i\in[1,m],P(a_i)=0
∀i∈[1,m],P(ai)=0
假设
f
(
x
)
=
P
(
x
)
g
(
x
)
+
A
(
x
)
f(x)=P(x)g(x)+A(x)
f(x)=P(x)g(x)+A(x)
A
(
x
)
=
f
(
x
)
%
P
(
x
)
A(x)=f(x)\%P(x)
A(x)=f(x)%P(x)
那显然
f
(
a
i
)
=
A
(
a
i
)
f(a_i)=A(a_i)
f(ai)=A(ai)
但由于
P
(
x
)
P(x)
P(x)是
m
m
m次的,没有起到优化的作用
而考虑对于
k
,
P
(
x
)
=
(
x
−
a
k
)
k,P(x)=(x-a_k)
k,P(x)=(x−ak)
则
f
(
x
)
%
P
(
x
)
f(x)\% P(x)
f(x)%P(x)后就只剩下一个常数项,即
f
(
a
i
)
f(a_i)
f(ai)的值了
但是这样一次就
n
l
o
g
n
nlogn
nlogn了
考虑分治优化
P
1
(
x
)
=
∏
i
=
l
m
i
d
(
x
−
a
i
)
,
P
2
(
x
)
=
∏
i
=
m
i
d
+
1
r
(
x
−
a
i
)
P_1(x)=\prod_{i=l}^{mid}(x-a_i),P_2(x)=\prod_{i=mid+1}^{r}(x-a_i)
P1(x)=∏i=lmid(x−ai),P2(x)=∏i=mid+1r(x−ai)
∀
i
∈
[
l
,
m
i
d
]
,
f
(
x
)
%
P
1
(
x
)
=
f
(
a
i
)
\forall i\in[l,mid] , f(x)\%P_1(x)=f(a_i)
∀i∈[l,mid],f(x)%P1(x)=f(ai)
取模之后
f
(
x
)
f(x)
f(x)的次数减少了一半,继续递归求解即可
P
(
x
)
P(x)
P(x)可以先分治
n
t
t
ntt
ntt预处理出来
复杂度
T
(
n
)
=
2
∗
T
(
n
2
)
+
O
(
n
l
o
g
n
)
=
O
(
n
l
o
g
2
n
)
T(n)=2*T(\frac n 2)+O(nlogn)=O(nlog^2n)
T(n)=2∗T(2n)+O(nlogn)=O(nlog2n)
代码实现
第一次写
T
T
T掉了,预处理了一波单位根就跑过去了
也可以在比较小的时候暴力秦九韶展开,然并没写
poly f[N<<2];
int a[N];
int n,m;
#define lc (u<<1)
#define rc ((u<<1)|1)
#define mid ((l+r)>>1)
void build(int u,int l,int r,int *v){
if(l==r){f[u].pb(dec(0,v[l])),f[u].pb(1);return;}
build(lc,l,mid,v),build(rc,mid+1,r,v);
f[u]=f[lc]*f[rc];
}
void calc(int u,int l,int r,poly now,int *v){
if(l==r){v[l]=now[0];return;}
calc(lc,l,mid,now%f[lc],v),calc(rc,mid+1,r,now%f[rc],v);
}
#undef lc
#undef rc
#undef mid
多项式快速插值
考虑传统的拉格朗日插值法插多项式是
O
(
n
2
)
O(n^2)
O(n2)的
即构造函数
f
(
x
)
=
∑
i
=
1
n
y
i
∏
j
=
̸
i
(
x
−
x
j
)
x
i
−
x
j
f(x)=\sum_{i=1}^{n}y_i\prod_{j=\not i}\frac{(x-x_j)}{x_i-x_j}
f(x)=i=1∑nyij≠i∏xi−xj(x−xj)
化一下
f ( x ) = ∑ i = 1 n y i ∏ j = ̸ i ( x i − x j ) ∏ j = ̸ i ( x − x j ) f(x)=\sum_{i=1}^{n}\frac{y_i}{\prod_{j=\not i}(x_i-x_j)}\prod_{j=\not i}(x-x_j) f(x)=i=1∑n∏j≠i(xi−xj)yij≠i∏(x−xj)
考虑如何求出
G
=
∏
j
=
̸
i
(
x
i
−
x
j
)
G=\prod_{j=\not i}(x_i-x_j)
G=∏j≠i(xi−xj)
设
g
(
x
)
=
∏
j
(
x
−
x
j
)
g(x)=\prod_j(x-x_j)
g(x)=∏j(x−xj)
j
=
̸
i
j=\not i
j≠i就相当于除以了
x
−
x
i
x-x_i
x−xi
那就变成了
g
(
x
i
)
(
x
−
x
i
)
\frac{g(x_i)}{(x-x_i)}
(x−xi)g(xi)
但是这个分子分母就都是0了,没法直接求
根据洛必达法则:
若
lim
x
→
a
f
(
x
)
=
0
,
lim
x
→
a
g
(
x
)
=
0
\lim_{x\rightarrow a}f(x)=0,\lim_{x\rightarrow a}g(x)=0
x→alimf(x)=0,x→alimg(x)=0
有
lim
x
→
a
f
(
x
)
g
(
x
)
=
lim
x
→
a
f
′
(
x
)
g
′
(
x
)
\lim_{x\rightarrow a}\frac{f(x)}{g(x)}=\lim_{x\rightarrow a}\frac{f'(x)}{g'(x)}
x→alimg(x)f(x)=x→alimg′(x)f′(x)
同时取导得到
G
=
g
′
(
x
i
)
G=g'(x_i)
G=g′(xi)
接下来考虑对整个式子分治
设
f
l
,
r
f_{l,r}
fl,r表示分治
[
l
,
r
]
[l,r]
[l,r]得到的答案
f l , r = ∑ i = l r y i g ′ ( x i ) ∏ j = ̸ i ( x − x j ) f_{l,r}=\sum_{i=l}^{r}\frac{y_i}{g'(x_i)}\prod_{j=\not i}(x-x_j) fl,r=i=l∑rg′(xi)yij≠i∏(x−xj)
= ∏ k = m i d + 1 r ( x − x k ) ∑ i = l m i d y i g ′ ( x i ) ∏ j = ̸ i [ l , m i d ] ( x − x j ) + ∏ k = l m i d ( x − x k ) ∑ i = m i d + 1 r y i g ′ ( x i ) ∏ j = ̸ i [ m i d + 1 , r ] ( x − x j ) =\prod_{k=mid+1}^{r}(x-x_k)\sum_{i=l}^{mid}\frac{y_i}{g'(x_i)}\prod_{j=\not i}^{[l,mid]}(x-x_j)+\prod_{k=l}^{mid}(x-x_k)\sum_{i=mid+1}^{r}\frac{y_i}{g'(x_i)}\prod_{j=\not i}^{[mid+1,r]}(x-x_j) =k=mid+1∏r(x−xk)i=l∑midg′(xi)yij≠i∏[l,mid](x−xj)+k=l∏mid(x−xk)i=mid+1∑rg′(xi)yij≠i∏[mid+1,r](x−xj)
= ∏ i = m i d + 1 r ( x − x i ) f l , m i d + ∏ i = l m i d ( x − x i ) f m i d + 1 , r =\prod_{i=mid+1}^r(x-x_i)f_{l,mid}+\prod_{i=l}^{mid}(x-x_i)f_{mid+1,r} =i=mid+1∏r(x−xi)fl,mid+i=l∏mid(x−xi)fmid+1,r
先分治 n t t ntt ntt求出 g g g,多点求值把 g ′ ( x i ) g'(x_i) g′(xi)求出来再分治一波就完了
复杂度 O ( n l o g 2 n ) O(nlog^2n) O(nlog2n)
下降幂多项式乘法
首先考虑对于
x
n
‾
x^{\underline n}
xn构建
E
G
F
EGF
EGF
x
n
‾
=
∑
i
=
n
∞
i
n
‾
i
!
x
i
=
∑
i
=
n
∞
1
(
i
−
n
)
!
x
i
=
x
n
e
x
x^{\underline n}=\sum_{i=n}^{\infty}\frac{i^{\underline n}}{i!}x^i=\sum_{i=n}^{\infty}\frac{1}{(i-n)!}x^i=x^ne^x
xn=∑i=n∞i!inxi=∑i=n∞(i−n)!1xi=xnex
考虑对于 f ( x ) = ∑ i = 0 ∞ a i x i ‾ f(x)=\sum_{i=0}^{\infty}a_ix^{\underline i} f(x)=∑i=0∞aixi的点值构造 E G F EGF EGF
g ( x ) = ∑ i = 0 ∞ f ( i ) i ! x i = ∑ i = 0 ∞ x i i ! ∑ j = 0 ∞ a j i j ‾ g(x)=\sum_{i=0}^{\infty}\frac{f(i)}{i!}x^i=\sum_{i=0}^{\infty}\frac{x^i}{i!}\sum_{j=0}^{\infty}a_ji^{\underline j} g(x)=∑i=0∞i!f(i)xi=∑i=0∞i!xi∑j=0∞ajij
= ∑ j = 0 ∞ a j ∑ i = 0 ∞ i j ‾ i ! x i = ∑ j = 0 ∞ a j x j e x =\sum_{j=0}^{\infty}a_j\sum_{i=0}^{\infty}\frac{i^{\underline j}}{i!}x^i=\sum_{j=0}^{\infty}a_jx^je^x =∑j=0∞aj∑i=0∞i!ijxi=∑j=0∞ajxjex
= e x ∑ i = 0 ∞ a i x i =e^x\sum_{i=0}^{\infty}a_ix^i =ex∑i=0∞aixi
所以只需要用普通多项式的系数乘个
e
x
e^x
ex就得到了点值的
E
G
F
EGF
EGF
点值还原原多项式只需要乘一个
e
−
x
e^{-x}
e−x即可
其他技巧
多项式快速幂
直接快速幂要多个
l
o
g
log
log而且常数大(虽然
l
n
ln
ln和
e
x
p
exp
exp常数一样大死个仙人)
当
A
[
0
]
=
1
A[0]=1
A[0]=1时(
L
n
Ln
Ln要求保证这个),
A
(
x
)
k
=
E
x
p
(
l
n
(
A
(
x
)
)
∗
k
)
A(x)^k=Exp(ln(A(x))*k)
A(x)k=Exp(ln(A(x))∗k)
洛谷板子读入时取模
inline poly ksm(poly a,int deg,int k){
a=exp(Ln(a,deg)*k,deg),a.resize(deg);
return a;
}
分治NTT
Bluestein
混合基NTT
MTT
模板合集
const int mod=998244353,g=3;
inline int add(int a,int b){return a+b>=mod?a+b-mod:a+b;}
inline void Add(int &a,int b){a=add(a,b);}
inline int dec(int a,int b){return a>=b?a-b:a-b+mod;}
inline void Dec(int &a,int b){a=dec(a,b);}
inline int mul(int a,int b){return 1ll*a*b>=mod?1ll*a*b%mod:a*b;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));return res;}
const int N=100005,C=17;
int *w[18];
int rev[N<<2];
#define poly vector<int>
#define pb push_back
#define bg begin
inline void init_rev(int lim){
for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
}
inline void init_w(){
for(int i=1;i<=C;i++)
w[i]=new int[1<<(i-1)];
int wn=ksm(g,(mod-1)/(1<<C));
w[C][0]=1;
for(int i=1;i<(1<<(C-1));i++)w[C][i]=mul(w[C][i-1],wn);
for(int i=C-1;i;i--)
for(int j=0;j<(1<<(i-1));j++)
w[i][j]=w[i+1][j<<1];
}
inline void ntt(poly &f,int lim,int kd){
for(int i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
for(int mid=1,l=1;mid<lim;mid<<=1,l++)
for(int i=0;i<lim;i+=(mid<<1))
for(int j=0,a0,a1;j<mid;j++)
a0=f[i+j],a1=mul(f[i+j+mid],w[l][j]),
f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
if(kd==-1&&(reverse(f.bg()+1,f.bg()+lim),1))
for(int inv=ksm(lim,mod-2),i=0;i<lim;i++)Mul(f[i],inv);
}
inline poly operator +(poly a,poly b){
poly c;int lim=max(a.size(),b.size());c.resize(lim);
for(int i=0;i<lim;i++)c[i]=add(a[i],b[i]);return c;
}
inline poly operator -(poly a,poly b){
poly c;int lim=max(a.size(),b.size());c.resize(lim);
for(int i=0;i<lim;i++)c[i]=dec(a[i],b[i]);return c;
}
inline poly operator *(poly a,int b){
for(int i=0;i<a.size();i++)Mul(a[i],b);return a;
}
inline poly operator /(poly a,int b){
for(int i=0,inv=ksm(b,mod-2);i<a.size();i++)Mul(a[i],inv);
return a;
}
inline poly operator *(poly a,poly b){
int deg=a.size()+b.size()-1,lim=1;
if(deg<=128){
poly c(deg,0);
for(int i=0;i<a.size();i++)
for(int j=0;j<b.size();j++)
Add(c[i+j],mul(a[i],b[j]));
return c;
}
while(lim<deg)lim<<=1;init(lim);
a.resize(lim),ntt(a,lim,1);
b.resize(lim),ntt(b,lim,1);
for(int i=0;i<lim;i++)Mul(a[i],b[i]);
ntt(a,lim,-1),a.resize(deg);
return a;
}
inline poly Inv(poly a,int deg){
poly c,b(1,ksm(a[0],mod-2));
for(int lim=4;lim<(deg<<2);lim<<=1){
init(lim);
c=a,c.resize(lim>>1);
c.resize(lim),ntt(c,lim,1);
b.resize(lim),ntt(b,lim,1);
for(int i=0;i<lim;i++)Mul(b[i],dec(2,mul(b[i],c[i])));
ntt(b,lim,-1),b.resize(lim>>1);
}b.resize(deg);return b;
}
inline poly Sqrt(poly a,int deg){
poly b(1,1),c,d;
for(int lim=4;lim<(deg<<2);lim<<=1){
c=a,c.resize(lim>>1);
init(lim),d=Inv(b,lim>>1),
c.resize(lim),ntt(c,lim,1);
d.resize(lim),ntt(d,lim,1);
for(int i=0;i<lim;i++)Mul(c[i],d[i]);
ntt(c,lim,-1),b.resize(lim>>1);
for(int i=0;i<(lim>>1);i++)b[i]=mul(inv[2],add(b[i],c[i]));
}b.resize(deg);return b;
}
inline poly operator /(poly a,poly b){
int lim=1,deg=a.size()-b.size()+1;
reverse(a.bg(),a.end());
reverse(b.bg(),b.end());
while(lim<=deg)lim<<=1;
b=Inv(b,lim),b.resize(deg);
a=a*b,a.resize(deg);
reverse(a.bg(),a.end());
return a;
}
inline poly operator %(poly a,poly b){
poly c=a-(a/b)*b;
c.resize(b.size()-1);
return c;
}
inline poly deriv(poly a){
for(int i=0;i<a.size()-1;i++)a[i]=mul(a[i+1],i+1);
a.pob();return a;
}
inline poly integ(poly a){
a.pb(0);
for(int i=a.size()-1;i;i--)a[i]=mul(a[i-1],inv[i]);
a[0]=0;
return a;
}
inline poly Ln(poly a,int lim){
a=integ(deriv(a)*Inv(a,lim)),a.resize(lim);
return a;
}
inline poly exp(poly a,int deg){
poly b(1,1),c;int n=a.size();
for(int lim=2;lim<(deg<<1);lim<<=1){
c=Ln(b,lim);
for(int i=0;i<lim;i++)c[i]=dec(i<n?a[i]:0,c[i]);
Add(c[0],1),b=b*c;
b.resize(lim);
}b.resize(deg);return b;
}
inline poly ksm(poly a,int deg,int k){
a=exp(Ln(a,deg)*k,deg),a.resize(deg);
return a;
}
poly f[N<<2];
#define mid ((l+r)>>1)
#define lc (u<<1)
#define rc ((u<<1)|1)
inline void build(int u,int l,int r,int *a){
if(l==r){
f[u].clear();
f[u].pb(dec(0,a[l]));
f[u].pb(1);return;
}build(lc,l,mid,a),build(rc,mid+1,r,a);
f[u]=f[lc]*f[rc];
}
inline void calc(int u,int l,int r,poly g,int *a){
if(l==r){a[l]=g[0];return;}
calc(lc,l,mid,g%f[lc],a);
calc(rc,mid+1,r,g%f[rc],a);
}
inline void getans(poly a,int *b,int num){
build(1,1,num,b),calc(1,1,num,a,b);
}
#undef mid
#undef lc
#undef rc