施工现场
代码就不放了太难写了。。这里都是一些工具,说白了还是要靠脑子的QUQ
多项式
在数学中,由若干个单项式相加组成的代数式叫做多项式(若有减法:减一个数等于加上它的相反数)。多项式中的每个单项式叫做多项式的项,这些单项式中的最高项次数,就是这个多项式的次数。其中多项式中不含字母的项叫做常数项。
为了方便接下来瞎逼逼,我们定义n次多项式 A ( x ) = ∑ i = 1 n a i ⋅ x i A(x)=\sum\limits_{i=1}^n a_i\cdot x^i A(x)=i=1∑nai⋅xi,记 [ x n ] F ( x ) [x^n]F(x) [xn]F(x)表示多项式 F ( x ) F(x) F(x)的n次项系数。在这里x并不是很重要,我们更关注的是系数a[]
多项式加法&减法
我是来搞笑的。。
对应系数相加即可,取界次较高的多项式为和多项式的界次
多项式乘法
多项式的乘法,是指把一个多项式中的每个单项式与另一个多项式中的每个单项式相乘之后合并同类项。
实现就是FFT和NTT,这个随便一搜都是嘛
任意膜数的多项式乘法
众所周知NTT需要一些特殊条件才能写,当膜数比较一般的时候就莫得办法了
一个比较暴力的做法就是我们取两个膜数分别做6次dft,然后把它们用CRT合并起来。这样子显然是不够优秀的
于是这里有一个被广大群众称为MTT的做法。具体推导可以看这里
这里只放一个板子
由于一些奇怪的原因我们得用longdouble来保证精度,包括定义的常量
π
\pi
π
Code
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <math.h>
#define copy(x,t) memcpy(x,t,sizeof(x))
typedef long long LL;
typedef long double ld;
const ld pi=acos(-1);
const int M=32767;
const int N=524299;
struct com {ld r,i;} ;
com operator +(const com &a,const com &b) {return (com) {a.r+b.r,a.i+b.i};}
com operator -(const com &a,const com &b) {return (com) {a.r-b.r,a.i-b.i};}
com operator *(const com &a,const com &b) {return (com) {a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r};}
com operator /(const com &a,const ld &b) {return (com) {a.r/b,a.i/b};}
com conj(com a) {return (com) {a.r,-a.i};}
int rv[N],a[N],b[N],c[N],MOD;
void FFT(com *a,int n) {
for (int i=0;i<n;++i) if (i<rv[i]) std:: swap(a[i],a[rv[i]]);
for (int i=1;i<n;i<<=1) {
com wn=(com) {cos(pi/i),sin(pi/i)};
for (int j=0;j<n;j+=(i<<1)) {
com w=(com) {1,0};
for (int k=0;k<i;++k) {
com u=a[j+k],v=a[j+k+i]*w;
a[j+k]=u+v; a[j+k+i]=u-v;
w=w*wn;
}
}
}
}
void mul(int *A,int *B,int *C,int n) {
static com a[N],b[N],Da[N],Db[N],Dc[N],Dd[N];
for (int i=0;i<n;++i) A[i]=(A[i]+MOD)%MOD;
for (int i=0;i<n;++i) B[i]=(B[i]+MOD)%MOD;
for (int i=0;i<n;++i) a[i]=(com) {(ld)(A[i]&M),(ld)(A[i]>>15)};
for (int i=0;i<n;++i) b[i]=(com) {(ld)(B[i]&M),(ld)(B[i]>>15)};
FFT(a,n),FFT(b,n);
for (int i=0;i<n;++i) {
com da,db,dc,dd;
int j=(n-1)&(n-i);
da=(a[i]+conj(a[j]))*(com) {0.5,0};
db=(a[i]-conj(a[j]))*(com) {0,-0.5};
dc=(b[i]+conj(b[j]))*(com) {0.5,0};
dd=(b[i]-conj(b[j]))*(com) {0,-0.5};
Da[j]=da*dc,Db[j]=da*dd,Dc[j]=db*dc,Dd[j]=db*dd;
}
for (int i=0;i<n;++i) a[i]=Da[i]+Db[i]*(com) {0,1};
for (int i=0;i<n;++i) b[i]=Dc[i]+Dd[i]*(com) {0,1};
FFT(a,n),FFT(b,n);
for (int i=0;i<n;++i) {
LL da=(LL)(a[i].r/n+0.5)%MOD;
LL db=(LL)(a[i].i/n+0.5)%MOD;
LL dc=(LL)(b[i].r/n+0.5)%MOD;
LL dd=(LL)(b[i].i/n+0.5)%MOD;
C[i]=((LL)da+((LL)(db+dc)<<15)+((LL)dd<<30))%MOD;
}
}
int main(void) {
freopen("data.in","r",stdin);
int n,m; scanf("%d%d%d",&n,&m,&MOD);
for (int i=0;i<=n;++i) scanf("%d",&a[i]);
for (int i=0;i<=m;++i) scanf("%d",&b[i]);
int len,lg;
for (len=1,lg=0;len<=n+m;) len<<=1,lg++;
for (int i=0;i<len;++i) rv[i]=(rv[i>>1]>>1)|((i&1)<<(lg-1));
mul(a,b,c,len);
for (int i=0;i<=n+m;++i) printf("%d ", (c[i]+MOD)%MOD);
return 0;
}
多项式积分&求导
记
∫
A
(
x
)
=
B
(
x
)
\int A(x)=B(x)
∫A(x)=B(x),那么
B
(
x
)
=
∑
i
=
1
n
a
i
−
1
i
⋅
x
i
B(x)=\sum\limits_{i=1}^n \frac{a_{i-1}}{i}\cdot x^i
B(x)=i=1∑niai−1⋅xi
记
A
′
(
x
)
=
C
(
x
)
A'(x)=C(x)
A′(x)=C(x),那么
C
(
x
)
=
∑
i
=
1
n
(
i
+
1
)
a
i
+
1
⋅
x
i
C(x)=\sum\limits_{i=1}^n (i+1)a_{i+1}\cdot x^i
C(x)=i=1∑n(i+1)ai+1⋅xi
具体推导可以问高中数学老师。。。我是认真的
实现起来需要注意for的顺序,以及不存在的位为0。这样做是O(n)的。
多项式求逆
忘的差不多了。。
记
A
(
x
)
B
(
x
)
≡
1
(
m
o
d
x
n
)
A(x)B(x)\equiv1\pmod {x^n}
A(x)B(x)≡1(modxn),且
A
(
x
)
G
(
x
)
≡
1
(
m
o
d
x
n
2
)
A(x)G(x)\equiv1\pmod{x^{\frac{n}{2}}}
A(x)G(x)≡1(modx2n)
显然 A ( x ) B ( x ) ≡ 1 ( m o d x n 2 ) A(x)B(x)\equiv1\pmod{x^{\frac{n}{2}}} A(x)B(x)≡1(modx2n),于是两柿子可以相减得到
G ( x ) − B ( x ) ≡ 0 ( m o d x n 2 ) G(x)-B(x)\equiv0\pmod{x^{\frac{n}{2}}} G(x)−B(x)≡0(modx2n)
平方得到 G 2 ( x ) + B 2 ( x ) − 2 B ( x ) G ( x ) ≡ 0 ( m o d x n ) G^2(x)+B^2(x)-2B(x)G(x)\equiv0\pmod{x^n} G2(x)+B2(x)−2B(x)G(x)≡0(modxn),注意这里模数变了
然后乘上 A ( x ) A(x) A(x)得到 G 2 ( x ) A ( x ) + B ( x ) − 2 G ( x ) ≡ 0 ( m o d x n ) G^2(x)A(x)+B(x)-2G(x)\equiv0\pmod{x^n} G2(x)A(x)+B(x)−2G(x)≡0(modxn)
移项就是 B ( x ) ≡ 2 G ( x ) − G 2 ( x ) A ( x ) ( m o d x n ) B(x)\equiv 2G(x)-G^2(x)A(x)\pmod{x^n} B(x)≡2G(x)−G2(x)A(x)(modxn)
于是倍增然后NTT或FFT就行了,要把多项式补齐为2的幂次项。这样做是O(nlogn)的。
多项式开根
这个要用到求逆
记
B
2
(
x
)
≡
A
(
x
)
(
m
o
d
x
n
2
)
B^2(x)\equiv A(x)\pmod{x^{\frac{n}{2}}}
B2(x)≡A(x)(modx2n),且
G
2
(
x
)
≡
A
(
x
)
(
m
o
d
x
n
)
G^2(x)\equiv A(x)\pmod{x^{n}}
G2(x)≡A(x)(modxn)
显然 G 2 ( x ) ≡ A ( x ) ( m o d x n 2 ) G^2(x)\equiv A(x)\pmod{x^{\frac{n}{2}}} G2(x)≡A(x)(modx2n),于是两柿子可以相减得到
[ G ( x ) − B ( x ) ] [ G ( x ) + B ( x ) ] ≡ 0 ( m o d x n 2 ) \left[G(x)-B(x)\right]\left[G(x)+B(x)\right]\equiv0\pmod{x^{\frac{n}{2}}} [G(x)−B(x)][G(x)+B(x)]≡0(modx2n)
显然后面不可能为0,我们把前面单独考虑就是
G ( x ) − B ( x ) ≡ 0 ( m o d x n 2 ) G(x)-B(x)\equiv0\pmod{x^{\frac{n}{2}}} G(x)−B(x)≡0(modx2n)
平方得到 G 2 ( x ) + B 2 ( x ) − 2 B ( x ) G ( x ) ≡ 0 ( m o d x n ) G^2(x)+B^2(x)-2B(x)G(x)\equiv0\pmod{x^n} G2(x)+B2(x)−2B(x)G(x)≡0(modxn),同样注意这里模数变了
根据 G ( x ) G(x) G(x)的定义得到 A ( x ) + B 2 ( x ) − 2 B ( x ) G ( x ) ≡ 0 ( m o d x n ) A(x)+B^2(x)-2B(x)G(x)\equiv0\pmod{x^n} A(x)+B2(x)−2B(x)G(x)≡0(modxn)
随便移项就有了 G ( x ) ≡ B 2 ( x ) + A ( x ) 2 B ( x ) ( m o d x n ) G(x)\equiv\frac{B^2(x)+A(x)}{2B(x)}\pmod{x^n} G(x)≡2B(x)B2(x)+A(x)(modxn),这个同样要倍增补齐,还要套一个求逆
多项式ln
常规的完了接下来就是一些比较不那么常规的了。。
记 l n ( A ( x ) ) = B ( x ) ln(A(x))=B(x) ln(A(x))=B(x)
两边求导再积分得到 B ( x ) = ∫ l n ( A ( x ) ) ′ = ∫ A ′ ( x ) A ( x ) B(x)=\int ln(A(x))'=\int \frac{A'(x)}{A(x)} B(x)=∫ln(A(x))′=∫A(x)A′(x)
然后参照上面的逆元和积分就可以了,这个也是O(nlogn)的
多项式exp
讲这个之前要科普一下牛顿迭代
牛顿迭代
将这个之前要科普一下泰勒展开
泰勒展开
数学中,泰勒公式是一个用函数在某点的信息描述其附近取值的公式。如果函数足够平滑的话,在已知函数在某一点的各阶导数值的情况之下,泰勒公式可以用这些导数值做系数构建一个多项式来近似函数在这一点的邻域中的值。泰勒公式还给出了这个多项式和实际的函数值之间的偏差。
大概长这个样子
f
(
x
)
=
f
(
x
0
)
0
!
(
x
−
x
0
)
0
+
f
′
(
x
0
)
1
!
(
x
−
x
0
)
1
+
f
′
′
(
x
0
)
2
!
(
x
−
x
0
)
2
+
.
.
+
f
(
n
)
(
x
0
)
n
!
(
x
−
x
0
)
n
+
R
n
(
x
)
f(x)=\frac{f(x_0)}{0!}{(x-x_0)}^0+\frac{f'(x_0)}{1!}{(x-x_0)}^1+\frac{f''(x_0)}{2!}{(x-x_0)}^2+..+\frac{f^{(n)}(x_0)}{n!}{(x-x_0)}^n+R_n(x)
f(x)=0!f(x0)(x−x0)0+1!f′(x0)(x−x0)1+2!f′′(x0)(x−x0)2+..+n!f(n)(x0)(x−x0)n+Rn(x)
其中
f
′
(
x
)
f'(x)
f′(x)是
f
(
x
)
f(x)
f(x)的一阶导函数,
f
′
′
(
x
)
f''(x)
f′′(x)是二阶导,也就是求导再求导。
f
(
n
)
(
x
)
f^{(n)}(x)
f(n)(x)则是n阶导。
R
n
(
x
)
R_n(x)
Rn(x)是泰勒公式的余项,展开很多次的情况下这东西接近高阶无穷小。如果不好理解暂时忽视它好像也没什么问题
牛顿迭代
牛顿真厉害。。
现在有一个函数
f
(
x
)
f(x)
f(x),我们要求它的零点
首先随便取一个位置
x
0
x_0
x0,然后在
x
0
x_0
x0处泰勒展开。为了方便(祖师爷规定)我们只搞出前两项线性项来
也就是
f
(
x
0
)
+
f
′
(
x
0
)
(
x
−
x
0
)
=
f
(
x
)
=
0
f(x_0)+f'(x_0)(x-x_0)=f(x)=0
f(x0)+f′(x0)(x−x0)=f(x)=0,化一下柿子就是
x
=
x
0
−
f
(
x
0
)
f
′
(
x
0
)
x=x_0-\frac{f(x_0)}{f'(x_0)}
x=x0−f′(x0)f(x0)
注意到这样得出来的并不是真正的零点,那么把x作为新的x0再做直到符合精度即可
假设我们现在有一个多项式函数
G
(
F
(
x
)
)
G(F(x))
G(F(x)),其中
F
(
x
)
F(x)
F(x)也是一个多项式,我们要求G的零点
考虑记
G
(
A
(
x
)
)
≡
0
(
m
o
d
x
n
2
)
G(A(x))\equiv0\pmod{x^{\frac{n}{2}}}
G(A(x))≡0(modx2n),
G
(
B
(
x
)
)
≡
0
(
m
o
d
x
n
)
G(B(x))\equiv0\pmod {x^n}
G(B(x))≡0(modxn)
由于 G ( B ( x ) ) ≡ 0 ( m o d x n ) G(B(x))\equiv0\pmod {x^n} G(B(x))≡0(modxn),所以 G ( B ( x ) ) ≡ 0 ( m o d x n 2 ) G(B(x))\equiv0\pmod {x^{\frac{n}{2}}} G(B(x))≡0(modx2n)
又因为 G ( x ) G(x) G(x)是一个多项式,所以 B ( x ) B(x) B(x)中次数较小的 n 2 \frac{n}{2} 2n项贡献为0,并且 B ( x ) B(x) B(x)的前 n 2 \frac{n}{2} 2n项恰好为 A ( x ) A(x) A(x)的前 n 2 \frac{n}{2} 2n项
即 B ( x ) − A ( x ) B(x)-A(x) B(x)−A(x)的前 n 2 \frac{n}{2} 2n项为0,那么 ( B ( x ) − A ( x ) ) 2 ≡ 0 ( m o d x n ) (B(x)-A(x))^2\equiv0 \pmod {x^n} (B(x)−A(x))2≡0(modxn)
讲了这么多有啥用呢?我们用
A
(
x
)
A(x)
A(x)做泰勒展开,根据上面可知只有前两项不为0,因此由牛顿迭代可由
A
(
x
)
A(x)
A(x)得到
B
(
x
)
B(x)
B(x)
B
(
x
)
≡
A
(
x
)
−
G
(
A
(
x
)
)
G
′
(
A
(
x
)
)
(
m
o
d
x
n
)
B(x) \equiv A(x)-\frac{G(A(x))}{G'(A(x))} \pmod{x^n}
B(x)≡A(x)−G′(A(x))G(A(x))(modxn)
于是我们倍增+求导+求逆就可以了。根据某主定理可知复杂度O(nlogn),但是常数超级大
多项式exp
容我去吃个宵夜
考虑设
e
F
(
x
)
=
A
(
x
)
e^{F(x)}=A(x)
eF(x)=A(x),设
G
(
A
(
x
)
)
=
F
(
x
)
−
ln
(
A
(
x
)
)
G(A(x))=F(x)-\ln(A(x))
G(A(x))=F(x)−ln(A(x)),我们只需要求这个函数的零点就可以了
根据上面的牛顿迭代可以知道
B
(
x
)
≡
A
(
x
)
(
1
−
l
n
(
A
(
x
)
)
+
B
(
x
)
)
(
m
o
d
x
n
)
B(x)\equiv A(x)\left(1-ln(A(x))+B(x)\right) \pmod{x^n}
B(x)≡A(x)(1−ln(A(x))+B(x))(modxn)
然后就没了
多项式除法&取模
考虑一个n次被除多项式
A
(
x
)
A(x)
A(x)和一个m次多项式
B
(
x
)
B(x)
B(x),我们要求一个n-m次的商多项式
C
(
x
)
C(x)
C(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),这里D(x)是一个m-1次的余多项式。除法和取模的不同在于我们的目的多项式不同(obivous
注意到求出了
C
(
x
)
C(x)
C(x)就可以一次NTT+一次减法求出
D
(
x
)
D(x)
D(x),因此我们想想要怎么弄出这个商
由那个除法的等柿可以知道
x
n
A
(
1
x
)
=
x
m
B
(
1
x
)
x
n
−
m
C
(
1
x
)
+
x
n
−
m
+
1
x
m
−
1
D
(
1
x
)
x^nA(\frac{1}{x})=x^mB(\frac{1}{x})x^{n-m}C(\frac{1}{x})+x^{n-m+1}x^{m-1}D(\frac{1}{x})
xnA(x1)=xmB(x1)xn−mC(x1)+xn−m+1xm−1D(x1)
观察这个
x
n
A
(
1
x
)
x^nA(\frac{1}{x})
xnA(x1),我们发现实际上就是把n次多项式
A
(
x
)
A(x)
A(x)的系数倒过来而已,记翻转的结果为
A
′
(
x
)
A'(x)
A′(x)
而通过这个巧妙的系数翻转操作我们成功地让
x
n
−
m
+
1
D
′
(
x
)
x^{n-m+1}D'(x)
xn−m+1D′(x)较低的
n
−
m
+
1
n-m+1
n−m+1项为0,于是
A
′
(
x
)
=
B
′
(
x
)
C
′
(
x
)
(
m
o
d
x
n
−
m
+
1
)
A'(x)=B'(x)C'(x)\pmod{x^{n-m+1}}
A′(x)=B′(x)C′(x)(modxn−m+1)
那么我们就可以先翻转系数,取较低的n-m+1位,然后求逆乘上再翻转回来即可
取模的话,商和除数都出来了你说要怎么办啊(滑稽
多点求值
感觉越来越毒瘤了。。这个好像适用面不是那么广啊
给定一个n次多项式
F
(
x
)
F(x)
F(x),和m个
x
x
x的取值,求
F
(
x
)
F(x)
F(x)在这m个x处的值
考虑构造一个多项式 G ( x ) = ∏ i = 1 ⌊ n 2 ⌋ ( x − x i ) G(x)=\prod\limits_{i=1}^{\lfloor\frac{n}{2}\rfloor}\left(x-x_i\right) G(x)=i=1∏⌊2n⌋(x−xi)很容易发现对于 x 0 ∈ [ 1 , ⌊ n 2 ⌋ ] x_0\in[1,\lfloor\frac{n}{2}\rfloor] x0∈[1,⌊2n⌋]有 G ( x 0 ) = 0 G(x_0)=0 G(x0)=0
设 F ( x ) = G ( x ) C ( x ) + D ( x ) F(x)=G(x)C(x)+D(x) F(x)=G(x)C(x)+D(x),根据上面的结论有 ∀ x 0 ∈ [ 1 , ⌊ n 2 ⌋ ] ,    D ( x 0 ) = F ( x 0 ) \forall x_0\in[1,\lfloor\frac{n}{2}\rfloor],\;D(x_0)=F(x_0) ∀x0∈[1,⌊2n⌋],D(x0)=F(x0),这一步比较显然
然后这个 D ( x ) D(x) D(x)正好就是余数,并且它是 ⌊ n 2 ⌋ \lfloor\frac{n}{2}\rfloor ⌊2n⌋次的,我们就成功地把规模缩小到了原本的一半,这样递归做就可以了
复杂度的话就是
T
(
n
)
=
T
(
n
/
2
)
+
O
(
n
log
n
)
=
O
(
n
log
2
n
)
T(n)=T(n/2)+O(n\log n)=O(n\log^2n)
T(n)=T(n/2)+O(nlogn)=O(nlog2n)的,常数很大,1e5大概跑2s的样子
luogu上有人用nm暴力+各种底层优化艹过去了,这也太会乱搞了。。
多点插值
先放着
拉格朗日反演
这个东西的证明有一点鬼畜。。我还不会证。。
对于两个多项式
F
(
x
)
F(x)
F(x)和
G
(
x
)
G(x)
G(x),满足它们没有常数项,且一次项系数互逆。若
G
(
F
(
x
)
)
=
x
G(F(x))=x
G(F(x))=x,且我们已知
G
(
x
)
G(x)
G(x),则有
[
x
n
]
F
(
x
)
=
1
n
[
x
n
−
1
]
(
x
G
(
x
)
)
n
[x^n]F(x)=\frac{1}{n}[x^{n-1}]{\left(\frac{x}{G(x)}\right)}^n
[xn]F(x)=n1[xn−1](G(x)x)n
这是一个在oi中比较常用的形式,我们把
F
(
x
)
F(x)
F(x)称为
G
(
x
)
G(x)
G(x)的复合逆
有一个结论就是 F ( G ( x ) ) = G ( F ( x ) ) F(G(x))=G(F(x)) F(G(x))=G(F(x)),因为复合运算构成了一个群
这个东西好像很稀有,目前见过的题只有大朋友和多叉树啥的。。