写法:
vector写有什么好处?
1.分治NTT的时候不用处理麻烦的下标。
2.传进函数的时候不用传指针而不会爆栈。
3.copy,swap,reverse,clear,resize等vector的函数可以帮助您快速地处理。
vector写的时候注意什么?
由于本人非常非常地懒,所以全部都是用vector写的。
根据测试,在开了 O2优化的情况下,是没有什么区别的。
但是不开O2,你懂的……
还是想用vector怎么办呢?
又经过一轮测试,只要在dft的时候把vector里的东西copy到一个数组,用这个数组做,再copy回去,是差不多的。
Dft:
不说了太累了。
我们设一次dft的常数为1,以便后面分析复杂度。
求逆:
设对多项式A求逆。
倍增法。
设有两个多项式 B 、 B 0 B、B0 B、B0
B
0
∗
A
=
1
(
m
o
d
x
n
/
2
)
B0 * A=1(mod~x^{n/2})
B0∗A=1(mod xn/2)
B
∗
A
=
1
(
m
o
d
x
n
)
B*A=1(mod~x^{n})
B∗A=1(mod xn)
其实可以用 B 0 B0 B0去推 B B B,推导如下:
首先B0的初值即n=1时 B 0 = a [ 0 ] − 1 B0=a[0]^{-1} B0=a[0]−1
B
0
∗
A
=
1
(
m
o
d
x
n
/
2
)
B0 * A=1(mod~x^{n/2})
B0∗A=1(mod xn/2)
B
0
∗
A
−
1
=
0
(
m
o
d
x
n
/
2
)
B0 * A-1=0(mod~x^{n/2})
B0∗A−1=0(mod xn/2)
等式两边同时平方:
(
B
0
∗
A
−
1
)
2
=
0
(
m
o
d
x
n
)
(B0 * A-1)^2=0(mod~x^{n})
(B0∗A−1)2=0(mod xn)
B
0
2
∗
A
2
−
2
B
0
∗
A
+
1
=
0
(
m
o
d
x
n
)
B0^2*A^2-2B0*A+1=0(mod~x^{n})
B02∗A2−2B0∗A+1=0(mod xn)
2
B
0
∗
A
−
B
0
2
∗
A
2
=
1
(
m
o
d
x
n
)
2B0*A-B0^2*A^2=1(mod~x^{n})
2B0∗A−B02∗A2=1(mod xn)
按A整理:
A
∗
(
2
B
0
−
B
0
2
∗
A
)
=
1
(
m
o
d
x
n
)
A*(2B0-B0^2*A)=1(mod~x^{n})
A∗(2B0−B02∗A)=1(mod xn)
所以
B
=
2
B
0
−
B
0
2
∗
A
(
m
o
d
x
n
)
B=2B0-B0^2*A(mod~x^{n})
B=2B0−B02∗A(mod xn)
我们来思考一下怎么做会快一点:
对B0求dft,对A求dft,然后算一算,再IDFT回去。
那么需要三次,但是每层的次数界必须开到4*n
复杂度为
T
(
n
)
=
n
l
o
g
n
+
T
(
n
/
2
)
≈
n
l
o
g
n
T(n)=nlogn+T(n/2)≈nlogn
T(n)=nlogn+T(n/2)≈nlogn
因为
n
+
n
/
2
+
n
/
4
+
…
<
2
n
n+n/2+n/4+…<2n
n+n/2+n/4+…<2n
所以复杂度大概为6次dft
开二次根:
方法类似于求逆,也是用倍增的思想。
已知 B 0 2 = A ( m o d x n / 2 ) B0^2=A(mod~x^{n/2}) B02=A(mod xn/2),推 B 2 = A ( m o d x n ) B^2=A(mod~x^{n}) B2=A(mod xn)
B0的初值可能需要二次剩余,一般多项式开根的题的常数项是个常数,可以直接预处理。
B
0
2
=
B
2
(
m
o
d
x
n
/
2
)
B0^2=B^2(mod~x^{n/2})
B02=B2(mod xn/2)
B
0
2
−
B
2
=
0
(
m
o
d
x
n
/
2
)
B0^2-B^2=0(mod~x^{n/2})
B02−B2=0(mod xn/2)
(
B
0
2
−
B
2
)
2
=
0
(
m
o
d
x
n
)
(B0^2-B^2)^2=0(mod~x^{n})
(B02−B2)2=0(mod xn)
B
0
4
−
2
B
0
2
∗
B
2
+
B
4
=
0
(
m
o
d
x
n
)
B0^4-2B0^2*B^2+B^4=0(mod~x^{n})
B04−2B02∗B2+B4=0(mod xn)
B
0
4
+
2
B
0
2
∗
B
2
+
B
4
=
4
B
0
2
∗
B
2
(
m
o
d
x
n
)
B0^4+2B0^2*B^2+B^4=4B0^2*B^2(mod~x^{n})
B04+2B02∗B2+B4=4B02∗B2(mod xn)
(
B
0
2
+
B
2
)
2
=
(
2
B
0
B
)
2
(
m
o
d
x
n
)
(B0^2+B^2)^2=(2B0B)^2(mod~x^{n})
(B02+B2)2=(2B0B)2(mod xn)
B
0
2
+
B
2
=
2
B
0
B
(
m
o
d
x
n
)
B0^2+B^2=2B0B(mod~x^{n})
B02+B2=2B0B(mod xn)
B
0
2
+
A
=
2
B
0
B
(
m
o
d
x
n
)
B0^2+A=2B0B(mod~x^{n})
B02+A=2B0B(mod xn)
B
=
B
0
/
2
+
A
/
(
2
B
0
)
(
m
o
d
x
n
)
B=B0/2+A/(2B0)(mod~x^{n})
B=B0/2+A/(2B0)(mod xn)
每层做的时候注意次数界是2n
每层要 d f t dft dft两次和求逆一次,复杂度大概是2*(3+6)=18次dft?
取模:
给出
A
(
x
)
,
B
(
x
)
,
d
e
g
(
A
)
>
=
d
e
g
(
B
)
A(x),B(x),deg(A)>=deg(B)
A(x),B(x),deg(A)>=deg(B)
求满足
A
(
x
)
=
B
(
x
)
∗
C
(
x
)
+
D
(
x
)
(
d
e
g
(
D
)
<
d
e
g
(
B
)
)
A(x)=B(x)*C(x)+D(x)(deg(D)<deg(B))
A(x)=B(x)∗C(x)+D(x)(deg(D)<deg(B))的
C
、
D
C、D
C、D
设
n
=
d
e
g
(
A
)
,
m
=
d
e
g
(
B
)
n=deg(A),m=deg(B)
n=deg(A),m=deg(B)
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
(
1
x
)
=
B
(
1
x
)
∗
C
(
1
x
)
+
D
(
1
x
)
A({1\over x})=B({1\over x})*C({1\over x})+D({1\over x})
A(x1)=B(x1)∗C(x1)+D(x1)
x
n
∗
A
(
1
x
)
=
x
m
∗
B
(
1
x
)
∗
x
n
−
m
∗
C
(
1
x
)
+
x
n
∗
D
(
1
x
)
x^n*A({1\over x})=x^m*B({1\over x})*x^{n-m}*C({1\over x})+x^n*D({1\over x})
xn∗A(x1)=xm∗B(x1)∗xn−m∗C(x1)+xn∗D(x1)
我们来理一下次数的关系:
x
n
∗
A
(
1
x
)
∈
[
0..
n
]
x^n*A({1\over x})∈[0..n]
xn∗A(x1)∈[0..n]
x
m
∗
B
(
1
x
)
∗
x
n
−
m
∗
C
(
1
x
)
∈
[
0..
n
]
x^m*B({1\over x})*x^{n-m}*C({1\over x})∈[0..n]
xm∗B(x1)∗xn−m∗C(x1)∈[0..n]
x
n
∗
D
(
1
x
)
∈
[
n
−
m
+
1..
n
]
x^n*D({1\over x})∈[n-m+1..n]
xn∗D(x1)∈[n−m+1..n]
对等式两边同时 m o d x n − m + 1 mod~x^{n-m+1} mod xn−m+1,则 x n ∗ D ( 1 x ) x^n*D({1\over x}) xn∗D(x1)没了。
x n ∗ A ( 1 x ) = x m ∗ B ( 1 x ) ∗ x n − m ∗ C ( 1 x ) ( m o d x n − m + 1 ) x^n*A({1\over x})=x^m*B({1\over x})*x^{n-m}*C({1\over x})(mod~x^{n-m+1}) xn∗A(x1)=xm∗B(x1)∗xn−m∗C(x1)(mod xn−m+1)
设 A T A^T AT就是把A反过来。
那么:
A
T
=
B
T
∗
C
T
(
m
o
d
x
n
−
m
+
1
)
A^T=B^T*C^T(mod~x^{n-m+1})
AT=BT∗CT(mod xn−m+1)
C
=
(
A
T
/
B
T
m
o
d
x
n
−
m
+
1
)
T
C=(A^T/{B^T}~mod~x^{n-m+1})^T
C=(AT/BT mod xn−m+1)T
坑点:必须要先 m o d x n − m + 1 mod~x^{n-m+1} mod xn−m+1再倒置。
求C要一次求逆和一次NTT,需要6+3=9次dft。
求了C后, D = A − B ∗ C D=A-B*C D=A−B∗C,需要一次NTT,即3次dft。
对数:
前置技能:
1.复合函数求导:
F ( G ( x ) ) ′ = F ′ ( G ( x ) ) ∗ G ′ ( x ) F(G(x))'=F'(G(x))*G'(x) F(G(x))′=F′(G(x))∗G′(x)
2.ln函数求导:
l n ( x ) ′ = 1 x ln(x)'={1\over x} ln(x)′=x1
l
n
(
A
)
ln(A)
ln(A)
=
∫
l
n
(
x
)
′
 
d
x
=\int {ln(x)'} \,{\rm d}x
=∫ln(x)′dx
=
∫
l
n
′
(
x
)
∗
x
′
 
d
x
=\int {ln'(x)*x'} \,{\rm d}x
=∫ln′(x)∗x′dx
=
∫
x
′
x
 
d
x
=\int {{x' \over x}} \,{\rm d}x
=∫xx′dx
只需要求逆一次和NTT一次,9次dft
指数:
前置技能:
牛顿迭代:
https://blog.csdn.net/ccnt_2012/article/details/81837154
可以发现这就是个求函数零点的玩意儿,虽然有些局限性,但是面对exp函数这么优美的图像肯定是没有问题的。
现在来思考多项式牛顿迭代,即求 f ( B ) = 0 f(B)=0 f(B)=0,其中 B B B是一个多项式。
思路是倍增,即知道了前 n / 2 n/2 n/2项的答案,去推前 n n n项的答案。
我们已知前 n / 2 n/2 n/2项为多项式B0,想要求出前 n n n项B。
可以直接由牛顿迭代得出:
B = B 0 − f ( B 0 ) f ′ ( B 0 ) ( m o d x n ) B=B0-{f(B0)\over f'(B0)}(mod~x^{n}) B=B0−f′(B0)f(B0)(mod xn)
下面介绍一下用泰勒展开推的方法:
f ( B ) = ∑ i = 0 ∞ f ( i ) ( B 0 ) ∗ ( B − B 0 ) i / i ! f(B)=\sum_{i=0}^{∞}f^{(i)}(B0)*(B-B0)^i/i! f(B)=∑i=0∞f(i)(B0)∗(B−B0)i/i!
注意当i>1时,每一项最低次项都是 x n x^{n} xn
即:
f
(
B
)
=
f
(
B
0
)
+
f
′
(
B
0
)
∗
(
B
−
B
0
)
(
m
o
d
x
n
)
f(B)=f(B0)+f'(B0)*(B-B0)(mod~x^n)
f(B)=f(B0)+f′(B0)∗(B−B0)(mod xn)
∵
f
(
B
)
=
0
(
m
o
d
x
n
)
∵f(B)=0(mod~x^n)
∵f(B)=0(mod xn)
∴
f
(
B
0
)
+
f
′
(
B
0
)
∗
(
B
−
B
0
)
=
0
(
m
o
d
x
n
)
∴f(B0)+f'(B0)*(B-B0)=0(mod~x^n)
∴f(B0)+f′(B0)∗(B−B0)=0(mod xn)
B
=
B
0
−
f
(
B
0
)
f
′
(
B
0
)
(
m
o
d
x
n
)
B=B0-{f(B0)\over f'(B0)}(mod~x^n)
B=B0−f′(B0)f(B0)(mod xn)
会了牛顿迭代之后思考如何把
e
x
p
(
A
)
exp(A)
exp(A)换出牛顿迭代的形式:
设
e
x
p
(
A
)
=
B
exp(A)=B
exp(A)=B
等式两边同时取
l
n
ln
ln,得
A
=
l
n
(
B
)
A=ln(B)
A=ln(B)
那么
f
(
B
)
=
l
n
(
B
)
−
A
f(B)=ln(B)-A
f(B)=ln(B)−A,现在就是求
f
f
f的零点。
套多项式牛顿迭代式得:
f
(
B
)
=
B
0
−
f
(
B
0
)
f
′
(
B
0
)
(
m
o
d
x
n
)
f(B)=B0-{f(B0)\over f'(B0)}(mod~x^n)
f(B)=B0−f′(B0)f(B0)(mod xn)
把
f
(
B
)
=
l
n
(
B
)
−
A
f(B)=ln(B)-A
f(B)=ln(B)−A代入,得:
f
(
B
)
=
B
0
−
(
l
n
(
B
0
)
−
A
)
/
(
l
n
(
B
0
)
−
A
)
′
(
m
o
d
x
n
)
f(B)=B0-(ln(B0)-A)/(ln(B0)-A)'(mod~x^n)
f(B)=B0−(ln(B0)−A)/(ln(B0)−A)′(mod xn)
f
(
B
)
=
B
0
−
(
l
n
(
B
0
)
−
A
)
/
(
l
n
(
B
0
)
′
−
A
′
)
(
m
o
d
x
n
)
f(B)=B0-(ln(B0)-A)/(ln(B0)'-A')(mod~x^n)
f(B)=B0−(ln(B0)−A)/(ln(B0)′−A′)(mod xn)
注意这里的函数的自变量是一个多项式,因变量也是。
A
A
A的话其实算作常数,求导之后是没有的。
因为自变量是多项式,所以
l
n
(
B
0
)
′
ln(B0)'
ln(B0)′不算复合函数求导,只是对
l
n
ln
ln求导,
l
n
(
B
0
)
′
=
1
/
B
0
ln(B0)'=1/B0
ln(B0)′=1/B0。
f
(
B
)
=
B
0
−
(
l
n
(
B
0
)
−
A
)
∗
B
0
(
m
o
d
x
n
)
f(B)=B0-(ln(B0)-A)*B0(mod~x^n)
f(B)=B0−(ln(B0)−A)∗B0(mod xn)
f
(
B
)
=
B
0
∗
(
1
−
l
n
(
B
0
)
+
A
)
(
m
o
d
x
n
)
f(B)=B0*(1-ln(B0)+A)(mod~x^n)
f(B)=B0∗(1−ln(B0)+A)(mod xn)
每层需要用到 l n ln ln和一次乘法,复杂度为 ( 9 + 3 ) ∗ 2 = 24 (9+3)*2=24 (9+3)∗2=24次dft
模板:
#include<ctime>
#include<vector>
#include<cstdio>
#include<algorithm>
#define ll long long
#define pp printf
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i < B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define pb push_back
#define si size()
using namespace std;
typedef vector<ll> V;
const int N = (1 << 19) + 5;
const int mo = 998244353;
const int ni2 = (mo + 1) / 2;
ll ksm(ll x, ll y) {
ll s = 1;
for(; y; y /= 2, x = x * x % mo)
if(y & 1) s = s * x % mo;
return s;
}
int r[N]; ll aa[N];
void dft(V &b, int f) {
int n = b.si;
ff(i, 0, n) aa[i] = b[i];
ff(i, 0, n) {
r[i] = r[i / 2] / 2 + (i & 1) * (n / 2);
if(i < r[i]) swap(aa[i], aa[r[i]]);
}
for(int h = 1; h < n; h *= 2) {
ll wn = ksm(ksm(3, (mo - 1) / 2 / h), f == -1 ? mo - 2 : 1);
for(int j = 0; j < n; j += 2 * h) {
ll A, W = 1, *l = aa + j, *r = aa + j + h;
ff(i, 0, h) {
A = *r * W, *r = (*l - A) % mo, *l = (*l + A) % mo;
W = W * wn % mo; l ++; r ++;
}
}
}
if(f == -1) {
ll v = ksm(n, mo - 2);
ff(i, 0, n) aa[i] = (aa[i] + mo) * v % mo;
}
ff(i, 0, n) b[i] = aa[i];
}
V operator * (V a, V b) {
int z = a.si + b.si - 1;
int n = 1; while(n < z) n *= 2;
a.resize(n); b.resize(n);
dft(a, 1); dft(b, 1);
ff(i, 0, n) a[i] = a[i] * b[i] % mo;
dft(a, -1); a.resize(z); return a;
}
V operator + (V a, V b) {
int sz = max(a.si, b.si); a.resize(sz); b.resize(sz);
ff(i, 0, sz) a[i] = (a[i] + b[i] + mo) % mo;
return a;
}
V operator - (V a, V b) {
int sz = max(a.si, b.si); a.resize(sz); b.resize(sz);
ff(i, 0, sz) a[i] = (a[i] - b[i] + mo) % mo;
return a;
}
V qni(V a) {
int n0 = 1; while(n0 < a.si) n0 *= 2;
V b; b.clear(); b.pb(ksm(a[0], mo - 2));
for(int n = 2; n <= n0; n *= 2) {
b.resize(n * 2); dft(b, 1);
V c = a; c.resize(n); c.resize(n * 2); dft(c, 1);
ff(i, 0, n * 2) b[i] = (b[i] * 2 - b[i] * b[i] % mo * c[i]) % mo;
dft(b, -1); b.resize(n);
}
b.resize(a.si); return b;
}
void fan(V &a) {
reverse(a.begin(), a.end());
}
V div(V a, V b) {
int n = a.size() - b.size() + 1;
fan(a); fan(b);
b.resize(a.size()); b = qni(b);
a = a * b; a.resize(n);
fan(a);
return a;
}
V qmo(V a, V b) {
return a - (b * (div(a, b)));
}
ll w;
struct P {
ll x, y;
P(){}
P(ll _x, ll _y) {x = _x, y = _y;}
};
P operator *(P a, P b) { return P((a.x * b.x + a.y * b.y % mo * w) % mo, (a.x * b.y + a.y * b.x) % mo);}
int pd(ll x) {
return ksm(x, (mo - 1) / 2) == 1;
}
P ksm(P x, ll y) {
P s = P(1, 0);
for(; y; y /= 2, x = x * x)
if(y & 1) s = s * x;
return s;
}
ll Sqrt(ll n) {
if(!n) return 0;
while(1) {
ll a = rand() % mo;
w = (a * a - n + mo) % mo;
if(pd(w)) continue;
return ksm(P(a, 1), (mo + 1) / 2).x;
}
}
V Sqrt(V a) {
int n0 = 1; while(n0 < a.si) n0 *= 2;
V b; b.clear(); b.push_back(Sqrt(a[0])); b[0] = 1;
for(int n = 2; n <= n0; n *= 2) {
V c = a; c.resize(n);
V d = b; d.resize(n);
ff(i, 0, d.si) d[i] = d[i] * 2 % mo;
ff(i, 0, b.si) b[i] = b[i] * ni2 % mo;
c = c * qni(d);
b = b + c;
}
b.resize(n0); return b;
}
V qd(V a) {
a[0] = 0;
ff(i, 1, a.si) a[i - 1] = a[i] * i % mo;
return a;
}
V jf(V a) {
a.pb(0);
fd(i, a.si, 1) a[i] = a[i - 1] * ksm(i, mo - 2) % mo;
a[0] = 0;
return a;
}
V ln(V a) {
int sa = a.si;
V b = a; b = qni(b); a = qd(a);
a = a * b; a = jf(a); a.resize(sa);
return a;
}
V exp(V a) {
int n0 = 1; while(n0 < a.si) n0 *= 2;
V b; b.clear(); b.pb(1);
for(int n = 2; n <= n0; n *= 2) {
V c = b; c.resize(n); c = ln(c);
V d = a; d.resize(n);
c = c - d;
c.resize(2 * n); dft(c, 1);
b.resize(2 * n); dft(b, 1);
ff(i, 0, 2 * n) b[i] = (b[i] - c[i] * b[i]) % mo;
dft(b, -1); b.resize(n);
}
b.resize(a.si);
return b;
}
int n, m;
V a, b;
int main() {
srand(time (0));
}