2022年7月7日和8日,day04
背景
- 计数问题利用生成函数的方法转化为形式幂级数的代数计算问题。
- 假设 C ( x ) = A ( x ) + B ( x ) C(x)=A(x)+B(x) C(x)=A(x)+B(x) 或 A ( x ) B ( x ) A(x)B(x) A(x)B(x) 或 A ( x ) B ( x ) \frac{A(x)}{B(x)} B(x)A(x) 或 A ( x ) \sqrt{A(x)} A(x) 等形式, C ( x ) C(x) C(x) 的前 n + 1 n+1 n+1 项的系数 c 0 , c 1 , … , c n c_0, c_1, \dots, c_n c0,c1,…,cn 仅仅需要 A ( x ) , B ( x ) A(x), B(x) A(x),B(x) 的前 n + 1 n+1 n+1 项系数就可以计算。
- 多项式 A ( x ) , B ( x ) m o d x n + 1 A(x), B(x) \mod x^{n+1} A(x),B(x)modxn+1 代替形式幂级数 A ( x ) , B ( x ) A(x), B(x) A(x),B(x)。
- 朴素多项式乘法、求逆、开根号等的复杂度为 O ( n 2 ) O(n^2) O(n2)。
- F F T / N T T FFT/NTT FFT/NTT 可以 O ( n l o g ( n ) ) O(nlog(n)) O(nlog(n)) 内实现两个 n n n 次多项式的乘积。
多项式的表示形式
系数表示与点值表示
假设 f ( x ) f(x) f(x) 是一个 n n n 次多项式,则 f ( x ) f(x) f(x) 的系数表示为 f ( x ) = a n x n + a n − 1 x n − 1 + ⋯ + a 0 f(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_0 f(x)=anxn+an−1xn−1+⋯+a0。 f ( x ) f(x) f(x) 的点值表示为 ( x 0 , f ( x 0 ) ) , ( x 1 , f ( x 1 ) ) , … , ( x n , f ( x n ) ) (x_0, f(x_0)), (x_1, f(x_1)), \dots, (x_n, f(x_n)) (x0,f(x0)),(x1,f(x1)),…,(xn,f(xn))。
- n + 1 n+1 n+1 个点值可以表示一个 n n n 次多项式。
- 在点值表示下
n
n
n 次多项式的乘法复杂度为
O
(
n
)
O(n)
O(n)。
h ( x ) = f ( x ) × g ( x ) , h ( x ) 是 n 次 多 项 式 ( x 0 , f ( x 0 ) ) , ( x 1 , f ( x 1 ) ) , … , ( x n , f ( x n ) ) ( x 0 , g ( x 0 ) ) , ( x 1 , g ( x 1 ) ) , … , ( x n , g ( x n ) ) ⇒ ( x 0 , h ( x 0 ) ) , ( x 1 , h ( x 1 ) ) , … , ( x n , h ( x n ) ) h ( x i ) = f ( x i ) × g ( x i ) h(x)=f(x)\times g(x),\ h(x)是\ n\ 次多项式\\ \begin{aligned} &(x_0, f(x_0)), (x_1, f(x_1)), \dots, (x_n, f(x_n))\\ &(x_0, g(x_0)), (x_1, g(x_1)), \dots, (x_n, g(x_n))\\ \Rightarrow &(x_0, h(x_0)), (x_1, h(x_1)), \dots, (x_n, h(x_n))\\ \end{aligned}\\ h(x_i)=f(x_i)\times g(x_i) h(x)=f(x)×g(x), h(x)是 n 次多项式⇒(x0,f(x0)),(x1,f(x1)),…,(xn,f(xn))(x0,g(x0)),(x1,g(x1)),…,(xn,g(xn))(x0,h(x0)),(x1,h(x1)),…,(xn,h(xn))h(xi)=f(xi)×g(xi)
复数与单位根
复数的指数形式: a + b i = r e i θ a+bi=re^{i\theta} a+bi=reiθ,其中 r = a 2 + b 2 , t a n ( θ ) = b a r=\sqrt{a^2+b^2},\ tan(\theta)=\frac{b}{a} r=a2+b2, tan(θ)=ab。
单位根: x n = 1 x^n=1 xn=1 在复数域上的根称为 n n n 次单位根。 n n n 次单位根有 n n n 个,形式为 w n k = e i 2 k π n w_n^k=e^{i\frac{2k\pi}{n}} wnk=ein2kπ。
欧拉公式:
e
i
x
=
c
o
s
x
+
i
s
i
n
x
e^{ix}=cosx+isinx
eix=cosx+isinx。
w
n
=
e
i
2
π
n
⇒
(
w
n
)
k
=
e
i
2
k
π
n
⇒
(
w
n
k
)
n
=
e
i
2
k
n
π
n
=
1
w_n=e^{i\frac{2\pi}{n}}\Rightarrow (w_n)^k=e^{i\frac{2k\pi}{n}}\Rightarrow (w_n^k)^n=e^{i\frac{2kn\pi}{n}}=1
wn=ein2π⇒(wn)k=ein2kπ⇒(wnk)n=ein2knπ=1
w
n
k
w_n^k
wnk 可以理解为
x
n
=
1
x^n=1
xn=1 在复数域上的第
k
+
1
k+1
k+1 个根。
单位根的性质:
- w n k = w 2 n 2 k w_n^k=w_{2n}^{2k} wnk=w2n2k
- w 2 n k + n = − w 2 n k w_{2n}^{k+n}=-w_{2n}^k w2nk+n=−w2nk
FFT(快速傅里叶变换,Fast Fourier Transform)
DFT(离散傅里叶变换,discrete Fourier Transform)
将多项式 A ( x ) = a 0 + a 1 x + ⋯ + a n − 1 x n − 1 A(x)=a_0+a_1x+\dots+a_{n-1}x^{n-1} A(x)=a0+a1x+⋯+an−1xn−1 转化为其点值表达式 ( w n k , A ( w n k ) ) (w_n^k,\ A(w_n^k)) (wnk, A(wnk)), k = 0 , 1 , … , n − 1 k=0,1,\dots,n-1 k=0,1,…,n−1。
不妨设
n
=
2
l
n=2^l
n=2l
A
(
x
)
=
a
0
+
a
1
x
+
⋯
+
a
n
−
1
x
n
−
1
A
(
x
)
=
(
a
0
+
a
2
x
2
+
⋯
+
a
n
−
2
x
n
−
2
)
+
(
a
1
x
+
a
3
x
3
+
⋯
+
a
n
−
1
x
n
−
1
)
\begin{aligned} A(x)=&a_0+a_1x+\cdots+a_{n-1}x^{n-1}\\ A(x)=&(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})\\ &+(a_1x+a_3x^3+\cdots+a_{n-1}x^{n-1}) \end{aligned}
A(x)=A(x)=a0+a1x+⋯+an−1xn−1(a0+a2x2+⋯+an−2xn−2)+(a1x+a3x3+⋯+an−1xn−1)
令
B
(
x
)
,
C
(
x
)
B(x), C(x)
B(x),C(x) 为
B
(
x
)
=
a
0
+
a
2
x
+
⋯
+
a
n
−
2
x
n
2
−
1
C
(
x
)
=
a
1
+
a
3
x
+
⋯
+
a
n
−
1
x
n
2
−
1
B(x)=a_0+a_2x+\cdots+a_{n-2}x^{\frac{n}{2}-1}\\ C(x)=a_1+a_3x+\cdots+a_{n-1}x^{\frac{n}{2}-1}\\
B(x)=a0+a2x+⋯+an−2x2n−1C(x)=a1+a3x+⋯+an−1x2n−1
那么
A
(
x
)
A(x)
A(x) 重写为:
A
(
x
)
=
B
(
x
2
)
+
x
C
(
x
2
)
A(x)=B(x^2)+xC(x^2)
A(x)=B(x2)+xC(x2)
分两种情况:
-
0
≤
k
<
n
2
0\le k <\frac{n}{2}
0≤k<2n
A ( w n k ) = B ( w n 2 k ) + w n k C ( w n 2 k ) = B ( w n 2 k ) + w n k C ( w n 2 k ) A(w_n^k)=B(w_n^{2k}) + w_n^kC(w_n^{2k})=B(w_{\frac{n}{2}}^k)+w_n^kC(w_{\frac{n}{2}}^k) A(wnk)=B(wn2k)+wnkC(wn2k)=B(w2nk)+wnkC(w2nk) -
n
2
≤
k
+
n
2
<
n
\frac{n}{2}\le k+\frac{n}{2}<n
2n≤k+2n<n
A ( w n k + n 2 ) = B ( w n 2 k ) + w n k + n 2 C ( w n 2 k ) = B ( w n 2 k ) − w n k C ( w n 2 k ) A(w_n^{k+\frac{n}{2}})=B(w_n^{2k}) + w_n^{k+\frac{n}{2}}C(w_n^{2k})=B(w_{\frac{n}{2}}^k)-w_n^kC(w_{\frac{n}{2}}^k) A(wnk+2n)=B(wn2k)+wnk+2nC(wn2k)=B(w2nk)−wnkC(w2nk)
时间复杂度: T ( n ) = 2 T ( n 2 ) + O ( n ) ⇒ T ( n ) = O ( n l o g n ) T(n)=2T(\frac{n}{2})+O(n)\Rightarrow T(n)=O(nlogn) T(n)=2T(2n)+O(n)⇒T(n)=O(nlogn)。
IDFT(逆离散傅里叶变换,inverse DFT)
将多项式的点值表示 ( w n k , b k ) (w_n^k,\ b_k) (wnk, bk), ( k = 0 , 1 , … , n − 1 ) (k=0,1,\dots,n-1) (k=0,1,…,n−1) 转化为其系数表示 A ( x ) = a 0 + a 1 x + ⋯ + a n − 1 x n − 1 A(x)=a_0+a_1x+\dots+a_{n-1}x^{n-1} A(x)=a0+a1x+⋯+an−1xn−1。
[ w n 0 ⋅ 0 w n 0 ⋅ 1 ⋯ w n 0 ⋅ ( n − 1 ) w n 1 ⋅ 0 w n 1 ⋅ 1 ⋯ w n 1 ⋅ ( n − 1 ) ⋮ ⋮ ⋮ w n ( n − 1 ) ⋅ 0 w n ( n − 1 ) ⋅ 1 ⋯ w n ( n − 1 ) ⋅ ( n − 1 ) ] [ a 0 a 1 ⋮ a n − 1 ] = [ b 0 b 1 ⋮ b n − 1 ] \begin{bmatrix} w_n^{0\cdot 0} & w_n^{0\cdot 1} & \cdots & w_n^{0\cdot (n-1)} \\ w_n^{1\cdot 0} & w_n^{1\cdot 1} & \cdots & w_n^{1\cdot (n-1)} \\ \vdots & \vdots & & \vdots\\ w_n^{(n-1)\cdot 0} & w_n^{(n-1)\cdot 1} & \cdots & w_n^{(n-1)\cdot (n-1)} \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{bmatrix} = \begin{bmatrix} b_0\\ b_1\\ \vdots\\ b_{n-1} \end{bmatrix} ⎣⎢⎢⎢⎢⎡wn0⋅0wn1⋅0⋮wn(n−1)⋅0wn0⋅1wn1⋅1⋮wn(n−1)⋅1⋯⋯⋯wn0⋅(n−1)wn1⋅(n−1)⋮wn(n−1)⋅(n−1)⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎡a0a1⋮an−1⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡b0b1⋮bn−1⎦⎥⎥⎥⎤
设
n
×
n
n\times n
n×n 的矩阵
Ω
\Omega
Ω,其中
Ω
i
,
j
=
w
n
i
j
\Omega_{i,j}=w_n^{ij}
Ωi,j=wnij,
(
0
≤
i
,
j
<
n
)
(0\le i,j<n)
(0≤i,j<n),设向量
a
=
(
a
0
,
a
1
,
⋯
,
a
n
−
1
)
\boldsymbol{a}=(a_0,\ a_1,\ \cdots,\ a_{n-1})
a=(a0, a1, ⋯, an−1),
b
=
(
b
0
,
b
1
,
⋯
,
b
n
−
1
)
\boldsymbol{b}=(b_0,\ b_1,\ \cdots,\ b_{n-1})
b=(b0, b1, ⋯, bn−1),则 IDFT 相当于求解方程
Ω
a
=
b
\Omega \boldsymbol{a} = \boldsymbol{b}
Ωa=b
记
Ω
ˉ
i
,
j
=
w
n
−
i
j
\bar\Omega_{i,j}=w_n^{-ij}
Ωˉi,j=wn−ij,那么
M
=
Ω
ˉ
⋅
Ω
M=\bar\Omega\cdot\Omega
M=Ωˉ⋅Ω 为
[
w
n
−
0
⋅
0
w
n
−
0
⋅
1
⋯
w
n
−
0
⋅
(
n
−
1
)
w
n
−
1
⋅
0
w
n
−
1
⋅
1
⋯
w
n
−
1
⋅
(
n
−
1
)
⋮
⋮
⋮
w
n
−
(
n
−
1
)
⋅
0
w
n
−
(
n
−
1
)
⋅
1
⋯
w
n
−
(
n
−
1
)
⋅
(
n
−
1
)
]
[
w
n
0
⋅
0
w
n
0
⋅
1
⋯
w
n
0
⋅
(
n
−
1
)
w
n
1
⋅
0
w
n
1
⋅
1
⋯
w
n
1
⋅
(
n
−
1
)
⋮
⋮
⋮
w
n
(
n
−
1
)
⋅
0
w
n
(
n
−
1
)
⋅
1
⋯
w
n
(
n
−
1
)
⋅
(
n
−
1
)
]
=
[
n
0
⋯
0
0
n
⋯
0
⋮
⋮
⋮
0
0
⋯
n
]
\begin{bmatrix} w_n^{-0\cdot 0} & w_n^{-0\cdot 1} & \cdots & w_n^{-0\cdot (n-1)} \\ w_n^{-1\cdot 0} & w_n^{-1\cdot 1} & \cdots & w_n^{-1\cdot (n-1)} \\ \vdots & \vdots & & \vdots\\ w_n^{-(n-1)\cdot 0} & w_n^{-(n-1)\cdot 1} & \cdots & w_n^{-(n-1)\cdot (n-1)} \end{bmatrix} \begin{bmatrix} w_n^{0\cdot 0} & w_n^{0\cdot 1} & \cdots & w_n^{0\cdot (n-1)} \\ w_n^{1\cdot 0} & w_n^{1\cdot 1} & \cdots & w_n^{1\cdot (n-1)} \\ \vdots & \vdots & & \vdots\\ w_n^{(n-1)\cdot 0} & w_n^{(n-1)\cdot 1} & \cdots & w_n^{(n-1)\cdot (n-1)} \end{bmatrix}= \begin{bmatrix} n & 0 & \cdots & 0 \\ 0 & n & \cdots & 0 \\ \vdots & \vdots & & \vdots\\ 0 & 0 & \cdots & n \end{bmatrix}
⎣⎢⎢⎢⎢⎡wn−0⋅0wn−1⋅0⋮wn−(n−1)⋅0wn−0⋅1wn−1⋅1⋮wn−(n−1)⋅1⋯⋯⋯wn−0⋅(n−1)wn−1⋅(n−1)⋮wn−(n−1)⋅(n−1)⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡wn0⋅0wn1⋅0⋮wn(n−1)⋅0wn0⋅1wn1⋅1⋮wn(n−1)⋅1⋯⋯⋯wn0⋅(n−1)wn1⋅(n−1)⋮wn(n−1)⋅(n−1)⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎡n0⋮00n⋮0⋯⋯⋯00⋮n⎦⎥⎥⎥⎤
由于
M
i
j
=
w
n
−
i
⋅
0
w
n
0
⋅
j
+
w
n
−
i
⋅
1
w
n
1
⋅
j
+
⋯
+
w
n
−
i
⋅
(
n
−
1
)
w
n
(
n
−
1
)
⋅
j
=
w
n
(
j
−
i
)
⋅
0
+
w
n
(
j
−
i
)
⋅
1
+
⋯
+
w
n
(
j
−
i
)
⋅
(
n
−
1
)
=
{
n
if
j
=
i
1
−
(
w
n
n
)
(
j
−
i
)
1
−
w
n
j
−
i
=
0
if
j
≠
i
\begin{aligned} M_{ij}&=w_n^{-i\cdot 0}w_n^{0\cdot j}+w_n^{-i\cdot 1}w_n^{1\cdot j}+\cdots+w_n^{-i\cdot (n-1)}w_n^{(n-1)\cdot j}\\ &=w_n^{(j-i)\cdot 0}+w_n^{(j-i)\cdot 1}+\cdots+w_n^{(j-i)\cdot (n-1)}\\ &= \begin{cases} n &\text{if} &j=i\\ \frac{1-(w_n^n)^{(j-i)}}{1-w_{n}^{j-i}}=0 &\text{if} &j\not= i\\ \end{cases} \end{aligned}
Mij=wn−i⋅0wn0⋅j+wn−i⋅1wn1⋅j+⋯+wn−i⋅(n−1)wn(n−1)⋅j=wn(j−i)⋅0+wn(j−i)⋅1+⋯+wn(j−i)⋅(n−1)={n1−wnj−i1−(wnn)(j−i)=0ififj=ij=i
因此有
Ω
ˉ
Ω
=
n
I
\bar\Omega\Omega=nI
ΩˉΩ=nI,那么
a
=
1
n
Ω
ˉ
b
\boldsymbol{a} = \frac{1}{n}\bar\Omega\boldsymbol{b}
a=n1Ωˉb
相当于给定
B
(
x
)
=
b
0
+
b
1
x
+
⋯
+
b
n
−
1
x
n
−
1
B(x)=b_0+b_1x+\cdots+b_{n-1}x^{n-1}
B(x)=b0+b1x+⋯+bn−1xn−1,求点值
B
(
w
n
−
k
)
B(w_n^{-k})
B(wn−k),
(
0
≤
k
<
n
)
(0\le k< n)
(0≤k<n),那么
a
k
=
1
n
B
(
w
n
−
k
)
a_k=\frac{1}{n}B(w_n^{-k})
ak=n1B(wn−k)
例题1:https://ac.nowcoder.com/acm/contest/26012/A
记
A
A
A 和
B
B
B 分别为
A
=
a
n
⋅
1
0
n
+
⋯
+
a
i
⋅
1
0
i
+
⋯
+
a
1
⋅
10
+
a
0
⇒
A
(
x
)
=
a
n
⋅
x
n
+
⋯
+
a
i
⋅
x
i
+
⋯
+
a
1
⋅
x
+
a
0
B
=
b
m
⋅
1
0
m
+
⋯
+
b
j
⋅
1
0
j
+
⋯
+
b
1
⋅
10
+
b
0
⇒
B
(
x
)
=
b
m
⋅
x
m
+
⋯
+
b
j
⋅
x
j
+
⋯
+
b
1
⋅
x
+
b
0
A=a_n\cdot 10^n+\cdots+a_i\cdot 10^i+\cdots+a_1\cdot 10+a_0\\ \Rightarrow A(x)=a_n\cdot x^n+\cdots+a_i\cdot x^i+\cdots+a_1\cdot x+a_0\\ B=b_m\cdot 10^m+\cdots+b_j\cdot 10^j+\cdots+b_1\cdot 10+b_0\\ \Rightarrow B(x)=b_m\cdot x^m+\cdots+b_j\cdot x^j+\cdots+b_1\cdot x+b_0\\
A=an⋅10n+⋯+ai⋅10i+⋯+a1⋅10+a0⇒A(x)=an⋅xn+⋯+ai⋅xi+⋯+a1⋅x+a0B=bm⋅10m+⋯+bj⋅10j+⋯+b1⋅10+b0⇒B(x)=bm⋅xm+⋯+bj⋅xj+⋯+b1⋅x+b0
那么
C
C
C 为
C
(
x
)
=
A
(
x
)
×
B
(
x
)
⇒
C
(
x
)
=
c
n
+
m
⋅
x
n
+
m
+
⋯
+
c
i
⋅
x
i
+
⋯
+
c
1
⋅
x
+
c
0
C(x)=A(x)\times B(x)\\ \Rightarrow C(x)=c_{n+m}\cdot x^{n+m}+\cdots+c_i\cdot x^i+\cdots+c_1\cdot x+c_0\\
C(x)=A(x)×B(x)⇒C(x)=cn+m⋅xn+m+⋯+ci⋅xi+⋯+c1⋅x+c0
#include <bits/stdc++.h>
using namespace std;
const double Pi = acos(-1.0);
const int N = 2e5 + 50;
struct com {
double x, y;
com operator + (com r) { return {x + r.x, y + r.y}; }
com operator - (com r) { return {x - r.x, y - r.y}; }
com operator * (com r) { return {x * r.x - y * r.y, x * r.y + y * r.x}; }
}f[N], g[N];
int n, m, ans[N];
int rev[N], lim, bit;
int read (com* c) {
static char s[N]; cin >> s;
int len = strlen (s);
for (int i = len - 1; i >= 0; i --)
c[i].x = s[len - i - 1] - '0';
return len;
}
void fft (com* h, int inv) {
for (int i = 0; i < lim; i ++)
if (rev[i] < i) swap (h[i], h[rev[i]]);
for (int k = 1; k < lim; k <<= 1) {
com wn = {cos(Pi/k), inv * sin(Pi/k)};
for (int i = 0; i < lim; i += k << 1) {
com w = {1, 0};
for (int j = 0; j < k; w = w * wn, j ++) {
com x = h[i + j], y = w * h[i + j + k];
h[i + j] = x + y;
h[i + j + k] = x - y;
}
}
}
if (inv == 1) return ;
for (int i = 0; i < lim; i ++) h[i].x /= lim;
}
int main () {
n = read (f), m = read (g);
for (lim = 1; lim <= n + m; lim <<= 1) bit ++;
for (int i = 0; i < lim; i ++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
fft (f, 1), fft (g, 1);
for (int i = 0; i < lim; i ++) f[i] = f[i] * g[i];
fft (f, -1);
for (int i = 0; i < lim; i ++) {
ans[i] += f[i].x + 0.5;
if (ans[i] > 9) {
ans[i + 1] += ans[i] / 10;
ans[i] %= 10;
}
}
while (lim > 0 && ! ans[lim]) lim --;
while (~ lim) cout << ans[lim --];
return 0;
}
例题2:https://ac.nowcoder.com/acm/contest/26012/B
要注意的是三元组
(
i
,
j
,
k
)
(i,j,k)
(i,j,k) 是有序的,因此直接进行FFT会导致方案数错误,考虑容斥
U
=
{
(
i
,
j
,
k
)
∣
1
≤
i
,
j
,
k
≤
n
}
A
=
{
(
i
,
j
,
k
)
∣
i
=
j
}
B
=
{
(
i
,
j
,
k
)
∣
j
=
k
}
C
=
{
(
i
,
j
,
k
)
∣
i
=
k
}
U=\{(i,j,k)\ |\ 1\le i,j,k \le n\}\\ A=\{(i,j,k)\ |\ i=j\}\\ B=\{(i,j,k)\ |\ j=k\}\\ C=\{(i,j,k)\ |\ i=k\}
U={(i,j,k) ∣ 1≤i,j,k≤n}A={(i,j,k) ∣ i=j}B={(i,j,k) ∣ j=k}C={(i,j,k) ∣ i=k}
那么答案为
1
3
!
f
(
U
−
A
∪
B
∪
C
)
=
1
3
!
(
f
(
U
)
−
f
(
A
)
−
f
(
B
)
−
f
(
C
)
+
f
(
A
∩
B
)
+
f
(
B
∩
C
)
+
f
(
A
∩
C
)
−
f
(
A
∩
B
∩
C
)
)
=
1
3
!
(
f
(
U
)
−
3
f
(
A
)
+
2
f
(
A
∩
B
)
)
\begin{aligned} &\frac{1}{3!}f(U-A\cup B\cup C)\\ =&\frac{1}{3!}\Big(f(U)-f(A)-f(B)-f(C)+f(A\cap B)+f(B\cap C)+f(A\cap C) -f(A\cap B\cap C)\Big)\\ =&\frac{1}{3!}\Big(f(U)-3f(A)+2f(A\cap B)\Big) \end{aligned}
==3!1f(U−A∪B∪C)3!1(f(U)−f(A)−f(B)−f(C)+f(A∩B)+f(B∩C)+f(A∩C)−f(A∩B∩C))3!1(f(U)−3f(A)+2f(A∩B))
记
S
(
x
)
=
∑
i
x
s
i
,
S
2
(
x
)
=
∑
i
x
2
s
i
,
S
3
(
x
)
=
∑
i
x
3
s
i
S(x)=\sum_ix^{s_i},S_2(x)=\sum_ix^{2s_i},S_3(x)=\sum_ix^{3s_i}
S(x)=∑ixsi,S2(x)=∑ix2si,S3(x)=∑ix3si ,那么
1
3
!
(
f
(
U
)
−
3
f
(
A
)
+
2
f
(
A
∩
B
)
)
=
1
6
(
S
(
x
)
3
−
3
S
2
(
x
)
S
(
x
)
+
2
S
3
(
x
)
)
\begin{aligned} &\frac{1}{3!}\Big(f(U)-3f(A)+2f(A\cap B)\Big)\\ =&\frac{1}{6}\Big(S(x)^3-3S_2(x)S(x)+2S_3(x)\Big) \end{aligned}
=3!1(f(U)−3f(A)+2f(A∩B))61(S(x)3−3S2(x)S(x)+2S3(x))
#include <bits/stdc++.h>
using namespace std;
const int N = 4e5 + 50;
const double Pi = acos (-1.0);
const int mov = 2e4, lim = 1 << 17;
struct node {
double x, y;
node operator + (node r) { return {x + r.x, y + r.y}; }
node operator - (node r) { return {x - r.x, y - r.y}; }
node operator * (node r) { return {x * r.x - y * r.y, x * r.y + y * r.x}; }
}s[N], s2[N], s3[N], f1[N], f2[N];
int n, a[N], bit = 17, rev[N];
void fft (node* f, int flag) {
for (int i = 0; i < lim; i ++)
if (rev[i] > i) swap (f[i], f[rev[i]]);
for (int k = 1; k < lim; k <<= 1) {
node wn = {cos(Pi / k), flag * sin(Pi / k)};
for (int i = 0; i < lim; i += k << 1) {
node w = {1, 0};
for (int j = 0; j < k; j ++, w = w * wn) {
node x = f[i + j], y = w * f[i + j + k];
f[i + j] = x + y;
f[i + j + k] = x - y;
}
}
}
if (flag == 1) return ;
for (int i = 0; i < lim; i ++) f[i].x /= lim;
}
int main () {
cin >> n;
for (int i = 1; i <= n; i ++)
cin >> a[i], a[i] += mov;
for (int i = 1; i <= n; i ++) {
s[a[i]] = {1, 0};
s2[a[i] * 2] = {1, 0};
s3[a[i] * 3] = {1, 0};
}
for (int i = 0; i < lim; i ++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
fft (s, 1), fft (s2, 1);
for (int i = 0; i < lim; i ++) {
f1[i] = s[i] * s[i] * s[i];
f2[i] = s2[i] * s[i];
}
fft (f1, -1), fft (f2, -1);
for (int i = 0; i < lim; i ++) {
int cnt = (f1[i].x - 3 * f2[i].x + 2 * s3[i].x) / 6 + 0.5;
if (cnt > 0) cout << i - mov*3 << " : " << cnt << endl;
}
return 0;
}
例题3:https://ac.nowcoder.com/acm/contest/26012/C
T
T
T 串在匹配
S
S
S 串具有容错范围
k
k
k,也就是说
T
[
j
]
T[j]
T[j] 匹配
S
[
i
]
S[i]
S[i],当
S
S
S 串周围
[
i
−
k
,
i
+
k
]
[i-k,i+k]
[i−k,i+k] 有
T
[
j
]
T[j]
T[j],则说明
T
T
T 串与
S
[
i
…
i
+
∣
T
∣
−
1
]
S[i\dots i+|T|-1]
S[i…i+∣T∣−1] 串匹配。
思路: 依次考虑每一个字符
′
A
′
,
′
T
′
,
′
G
′
,
′
C
′
'A','T','G','C'
′A′,′T′,′G′,′C′,下面以
′
A
′
'A'
′A′ 为例
先将
S
S
S 串进行扩展,以
S
S
S 串中每一个字符
′
A
′
'A'
′A′ 为中心,半径为
k
k
k 的范围全部改成
′
A
′
'A'
′A′,构成
S
′
S'
S′ 串。
那么,考虑
S
S
S 串中的每一个起点
i
i
i ,记
f
(
i
,
′
A
′
)
f(i,'A')
f(i,′A′) 为
∑
j
=
1
∣
T
∣
[
T
[
j
]
=
′
A
′
]
⋅
[
S
′
[
i
+
j
−
1
]
=
′
A
′
]
\sum_{j=1}^{|T|}[T[j]='A']\cdot[S'[i+j-1]='A']
j=1∑∣T∣[T[j]=′A′]⋅[S′[i+j−1]=′A′]
定义
T
(
x
)
T(x)
T(x) 和
S
(
x
)
S(x)
S(x) ,这里需要将
S
′
S'
S′ 串进行翻转
T
(
x
)
=
∑
j
=
1
∣
T
∣
[
T
[
j
]
=
′
A
′
]
x
j
S
(
x
)
=
∑
i
=
1
∣
S
∣
[
S
′
[
n
−
i
]
=
′
A
′
]
x
i
T(x)=\sum_{j=1}^{|T|}[T[j]='A']x^j\\ S(x)=\sum_{i=1}^{|S|}[S'[n-i]='A']x^i
T(x)=j=1∑∣T∣[T[j]=′A′]xjS(x)=i=1∑∣S∣[S′[n−i]=′A′]xi
那么
[
x
u
]
(
T
(
x
)
×
S
(
x
)
)
[x^u]\Big(T(x)\times S(x)\Big)
[xu](T(x)×S(x)) 为
∑
j
=
1
∣
T
∣
[
T
[
j
]
=
′
A
′
]
⋅
[
S
′
[
n
−
(
u
−
j
)
]
=
′
A
′
]
\sum_{j=1}^{|T|}[T[j]='A']\cdot [S'[n-(u-j)]='A']
j=1∑∣T∣[T[j]=′A′]⋅[S′[n−(u−j)]=′A′]
由于
i
+
j
−
1
=
n
−
(
u
−
j
)
⇒
u
=
n
−
i
+
1
i+j-1=n-(u-j)\Rightarrow u=n-i+1
i+j−1=n−(u−j)⇒u=n−i+1
f
(
i
,
′
A
′
)
=
[
x
n
−
i
+
1
]
(
T
(
x
)
×
S
(
x
)
)
f(i,'A')=[x^{n-i+1}]\Big(T(x)\times S(x)\Big)
f(i,′A′)=[xn−i+1](T(x)×S(x))
最后判断满足下面条件的
i
i
i 的数量
f
(
i
,
′
A
′
)
+
f
(
i
,
′
T
′
)
+
f
(
i
,
′
G
′
)
+
f
(
i
,
′
C
′
)
=
∣
T
∣
f(i,'A') + f(i,'T') + f(i,'G') + f(i,'C')=|T|
f(i,′A′)+f(i,′T′)+f(i,′G′)+f(i,′C′)=∣T∣
#include <bits/stdc++.h>
using namespace std;
const int N = 6e5 + 50;
const double Pi = acos (-1.0);
const char ch[] = {'A', 'G', 'T', 'C'};
char s[N], t[N];
int n, m, k, a[N];
int lim, bit, rev[N];
struct node {
double x, y;
node (double _x=0, double _y=0): x(_x), y(_y) {}
node operator + (node r) { return node (x + r.x, y + r.y); }
node operator - (node r) { return node (x - r.x, y - r.y); }
node operator * (node r) { return node (x*r.x-y*r.y, x*r.y + y*r.x); }
}f[N], g[N];
void fft (node* ff, int inv) {
for (int i = 0; i < lim; i ++)
if (rev[i] > i) swap (ff[i], ff[rev[i]]);
for (int k = 1; k < lim; k <<= 1) {
node wn = node (cos (Pi / k), inv * sin (Pi / k));
for (int i = 0; i < lim; i += k << 1) {
node w = node (1, 0);
for (int j = 0; j < k; j ++, w = w * wn) {
node x = ff[i + j], y = w * ff[i + j + k];
ff[i + j] = x + y;
ff[i + j + k] = x - y;
}
}
}
if (inv == 1) return ;
for (int i = 0; i < lim; i ++) ff[i].x /= lim;
}
void extend (node* ff) {
static int st[N];
for (int i = 0; i < n; i ++) st[i] = ff[i].x + 0.5;
int cnt = 0, r = 0, l = 0;
for (int i = 0; i < n; i ++) {
while (r - i <= k && r < n) cnt += st[r ++];
if (cnt) ff[i] = node (1, 0);
cnt -= st[i];
}
cnt = 0, l = n - 1;
for (int i = n - 1; i >= 0; i --) {
while (i - l <= k && l >= 0) cnt += st[l --];
if (cnt) ff[i] = node (1, 0);
cnt -= st[i];
}
}
void print (node* ff) {
for (int i = 0; i < lim; i ++)
cout << int(ff[i].x + 0.5) << endl;
}
int main () {
cin >> n >> m >> k >> s >> t;
for (lim = 1; lim < n + m; lim <<= 1) bit ++;
for (int i = 0; i < lim; i ++)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
for (int i = 0; i < 4; i ++) {
for (int j = 0; j < n; j ++) f[j] = node (s[j] == ch[i]);
for (int j = n; j < lim; j ++) f[j] = node ();
for (int j = 0; j < m; j ++) g[m - j - 1] = node (t[j] == ch[i]);
for (int j = m; j < lim; j ++) g[j] = node ();
extend(f);
fft (f, 1), fft (g, 1);
for (int j = 0; j < lim; j ++) f[j] = f[j] * g[j];
fft (f, -1);
for (int j = m - 1; j < lim; j ++) a[j] += f[j].x + 0.5;
}
int ans = 0;
for (int i = m - 1; i < lim; i ++) ans += a[i] == m;
cout << ans << endl;
return 0;
}
拉格朗日插值
拉格朗日插值定理: n n n 个点值 ( x i , y i ) (x_i,y_i) (xi,yi),满足 x i ≠ x j , ( i ≠ j ) x_i\not=x_j,(i\not=j) xi=xj,(i=j),它们唯一确定一个 n − 1 n-1 n−1 次多项式 f ( x ) f(x) f(x),如下
f ( x ) = ∑ i = 1 n y i ∏ j ≠ i x − x j x i − x j f(x)=\sum_{i=1}^ny_i\prod_{j\not=i}\frac{x-x_j}{x_i-x_j} f(x)=i=1∑nyij=i∏xi−xjx−xj
构造一个
f
(
x
)
f(x)
f(x),
n
−
1
n-1
n−1 次多项式,满足
f
(
x
i
)
=
y
i
f(x_i)=y_i
f(xi)=yi。
f
(
x
)
=
f
1
(
x
)
+
f
2
(
x
)
+
⋯
+
f
i
(
x
)
+
⋯
+
f
n
(
x
)
f(x)=f_1(x)+f_2(x)+\cdots+f_i(x)+\cdots+f_n(x)
f(x)=f1(x)+f2(x)+⋯+fi(x)+⋯+fn(x)
对于每一个
f
i
(
x
)
f_i(x)
fi(x),要满足
f
i
(
x
k
)
=
{
y
k
if
i
=
k
0
if
i
≠
k
f_i(x_k)=\begin{cases} y_k &\text{if} &i=k\\ 0 &\text{if} &i\not= k\\ \end{cases}
fi(xk)={yk0ififi=ki=k
那么
f
i
(
x
)
=
A
(
x
−
x
1
)
⋯
(
x
−
x
i
−
1
)
(
x
−
x
i
+
1
)
⋯
(
x
−
x
n
)
=
A
∏
j
≠
i
(
x
−
x
j
)
\begin{aligned} f_i(x)&=A(x-x_1)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)\\ &=A\prod_{j\not= i}(x-x_j) \end{aligned}
fi(x)=A(x−x1)⋯(x−xi−1)(x−xi+1)⋯(x−xn)=Aj=i∏(x−xj)
由于
f
i
(
x
i
)
=
y
i
f_i(x_i)=y_i
fi(xi)=yi,有
f
i
(
x
i
)
=
A
∏
j
≠
i
(
x
i
−
x
j
)
=
y
i
⇒
A
=
y
i
∏
j
≠
i
(
x
i
−
x
j
)
\begin{aligned} &f_i(x_i) =A\prod_{j\not= i}(x_i-x_j)=y_i\\ \Rightarrow &A=\frac{y_i}{\prod_{j\not= i}(x_i-x_j)} \end{aligned}
⇒fi(xi)=Aj=i∏(xi−xj)=yiA=∏j=i(xi−xj)yi
因此
f
(
x
)
=
f
1
(
x
)
+
f
2
(
x
)
+
⋯
+
f
n
(
x
)
=
∑
i
=
1
n
y
i
∏
j
≠
i
(
x
i
−
x
j
)
∏
j
≠
i
(
x
−
x
j
)
=
∑
i
=
1
n
y
i
∏
j
≠
i
x
−
x
j
x
i
−
x
j
\begin{aligned} f(x)&=f_1(x)+f_2(x)+\cdots+f_n(x)\\ &=\sum_{i=1}^n\frac{y_i}{\prod_{j\not= i}(x_i-x_j)}\prod_{j\not= i}(x-x_j)\\ &=\sum_{i=1}^ny_i\prod_{j\not=i}\frac{x-x_j}{x_i-x_j} \end{aligned}
f(x)=f1(x)+f2(x)+⋯+fn(x)=i=1∑n∏j=i(xi−xj)yij=i∏(x−xj)=i=1∑nyij=i∏xi−xjx−xj
它和
C
R
T
CRT
CRT 的关系?
- 中国剩余定理(
C
R
T
CRT
CRT):
{ x ≡ a 1 ( m o d m 1 ) x ≡ a 2 ( m o d m 2 ) ⋮ x ≡ a i ( m o d m i ) ⋮ x ≡ a n ( m o d m n ) \begin{cases} x &\equiv &a_1&(\mod m_1)\\ x &\equiv &a_2&(\mod m_2)\\ & \vdots\\ x &\equiv &a_i&(\mod m_i)\\ & \vdots\\ x &\equiv &a_n&(\mod m_n) \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧xxxx≡≡⋮≡⋮≡a1a2aian(modm1)(modm2)(modmi)(modmn)
其中 ( m i , m j ) = 1 , ( i ≠ j ) (m_i,m_j)=1,(i\not=j) (mi,mj)=1,(i=j),如何构造 x x x
x = x 1 + x 2 + ⋯ + x i + ⋯ + x n x=x_1+x_2+\cdots+x_i+\cdots+x_n x=x1+x2+⋯+xi+⋯+xn
对于每一个 x i x_i xi,要满足
x i m o d m k = { a i if i = k 0 if i ≠ k x_i\mod m_k=\begin{cases} a_i &\text{if} &i=k\\ 0 &\text{if} &i\not= k\\ \end{cases} ximodmk={ai0ififi=ki=k
那么
x i = A m 1 ⋯ m i − 1 m i + 1 ⋯ m n = A ∏ j ≠ i m j \begin{aligned} x_i&=Am_1\cdots m_{i-1}m_{i+1}\cdots m_n\\ &=A\prod_{j\not= i}m_j \end{aligned} xi=Am1⋯mi−1mi+1⋯mn=Aj=i∏mj
由于 x i m o d m i = a i x_i\mod m_i=a_i ximodmi=ai,有
( x i = A ∏ j ≠ i x j ) ≡ a i ( m o d m i ) ⇒ A ≡ a i ( ∏ j ≠ i x j ) − 1 ( m o d m i ) ⇒ A ≡ a i t i ( m o d m i ) \begin{aligned} &(x_i =A\prod_{j\not= i}x_j)\equiv a_i&(\mod m_i)\\ \Rightarrow &A\equiv a_i(\prod_{j\not= i}x_j)^{-1}&(\mod m_i)\\ \Rightarrow &A\equiv a_it_i&(\mod m_i)\\ \end{aligned} ⇒⇒(xi=Aj=i∏xj)≡aiA≡ai(j=i∏xj)−1A≡aiti(modmi)(modmi)(modmi)
因此
x ≡ x 1 + x 2 + ⋯ + x n ≡ ∑ i = 1 n a i t i ∏ j ≠ i m j ( m o d ∏ i = 1 n m i ) \begin{aligned} x&\equiv x_1+x_2+\cdots+x_n\\ &\equiv \sum_{i=1}^n a_it_i\prod_{j\not= i}m_j&(\mod \prod_{i=1}^nm_i) \end{aligned} x≡x1+x2+⋯+xn≡i=1∑naitij=i∏mj(modi=1∏nmi)
暴力实现:先计算 g ( x ) = ∏ i = 1 n ( x − x i ) g(x)=\prod_{i=1}^n(x-x_i) g(x)=∏i=1n(x−xi),再求 n n n 次 g ( x ) / ( x − x i ) g(x)/(x-x_i) g(x)/(x−xi),按照权值将它们相加。 O ( n 2 ) O(n^2) O(n2)
f ( x ) = ∑ i = 1 n y i ∏ j ≠ i ( x i − x j ) g ( x ) ( x − x i ) f(x)=\sum_{i=1}^n\frac{y_i}{\prod_{j\not=i}(x_i-x_j)}\frac{g(x)}{(x-x_i)} f(x)=i=1∑n∏j=i(xi−xj)yi(x−xi)g(x)
特殊情况1:只需要求一个 f ( x 0 ) f(x_0) f(x0)。 O ( n 2 ) O(n^2) O(n2)
f ( x 0 ) = ∑ i = 1 n y i ∏ j ≠ i x 0 − x j x i − x j f(x_0)=\sum_{i=1}^ny_i\prod_{j\not=i}\frac{x_0-x_j}{x_i-x_j} f(x0)=i=1∑nyij=i∏xi−xjx0−xj
特殊情况2:只需要求一个 f ( x 0 ) f(x_0) f(x0),且 x i = i x_i=i xi=i。 O ( n ) O(n) O(n)
f ( x 0 ) = ∑ i = 1 n y i ∏ j ≠ i ( x 0 − j ) ( − 1 ) n − i ( n − i ) ! ( i − 1 ) ! f(x_0)=\sum_{i=1}^ny_i\frac{\prod_{j\not=i}(x_0-j)}{(-1)^{n-i}(n-i)!(i-1)!} f(x0)=i=1∑nyi(−1)n−i(n−i)!(i−1)!∏j=i(x0−j)
预处理 ∏ j ≠ i ( x 0 − j ) \prod_{j\not=i}(x_0-j) ∏j=i(x0−j) 的前缀积和后缀积,拼在一起。
例题4:https://ac.nowcoder.com/acm/contest/26012/D
通过观察,
f
(
n
)
=
∑
i
=
1
n
i
k
f(n)=\sum_{i=1}^ni^k
f(n)=∑i=1nik 是
k
+
1
k+1
k+1 次多项式。例如:
k
=
0
f
(
n
)
=
n
k
=
1
f
(
n
)
=
n
(
n
+
1
)
2
k
=
2
f
(
n
)
=
n
(
n
+
1
)
(
2
n
+
1
)
6
⋮
\begin{aligned} &k=0 &f(n)&=n \\ &k=1 &f(n)&=\frac{n(n+1)}{2} \\ &k=2 &f(n)&=\frac{n(n+1)(2n+1)}{6}\\ & & \vdots \end{aligned}
k=0k=1k=2f(n)f(n)f(n)⋮=n=2n(n+1)=6n(n+1)(2n+1)
结合第二类斯特林数的公式,有
n
m
=
∑
k
=
0
m
{
m
k
}
n
k
‾
n^m=\sum_{k=0}^m\begin{Bmatrix} m\\ k \end{Bmatrix}n^{\underline{k}}
nm=k=0∑m{mk}nk
其中
n
m
n^m
nm 的含义为把
m
m
m 个互不相同的球放入
n
n
n 个不同的盒子的方案数。换一种理解方式为枚举有
k
k
k 个盒子一定要放,将
m
m
m 个球都放进这
k
k
k 个盒子。
∑
k
=
0
m
(
n
k
)
k
!
{
m
k
}
=
∑
k
=
0
m
{
m
k
}
n
k
‾
\sum_{k=0}^m\binom{n}{k}k!\begin{Bmatrix} m\\ k \end{Bmatrix}=\sum_{k=0}^m\begin{Bmatrix} m\\ k \end{Bmatrix}n^{\underline{k}}
k=0∑m(kn)k!{mk}=k=0∑m{mk}nk
那么
f
(
n
)
=
∑
i
=
1
n
i
k
=
∑
i
=
1
n
∑
j
=
1
k
{
k
j
}
i
j
‾
=
∑
j
=
1
k
{
k
j
}
∑
i
=
1
n
i
j
‾
f(n)=\sum_{i=1}^ni^k=\sum_{i=1}^n\sum_{j=1}^k\begin{Bmatrix} k\\ j \end{Bmatrix}i^{\underline{j}}=\sum_{j=1}^k\begin{Bmatrix} k\\ j \end{Bmatrix}\sum_{i=1}^ni^{\underline{j}}
f(n)=i=1∑nik=i=1∑nj=1∑k{kj}ij=j=1∑k{kj}i=1∑nij
因为
∑
i
=
1
n
i
j
‾
=
j
(
j
−
1
)
⋯
1
+
(
j
+
1
)
j
⋯
2
+
⋯
+
n
(
n
−
1
)
⋯
(
n
−
j
+
1
)
=
(
j
+
1
)
j
⋯
1
−
j
(
j
−
1
)
⋯
0
j
+
1
+
(
j
+
2
)
(
j
+
1
)
⋯
2
−
(
j
+
1
)
j
⋯
1
j
+
1
+
⋯
+
(
n
+
1
)
n
⋯
(
n
−
j
+
1
)
−
n
(
n
−
1
)
⋯
(
n
−
j
)
j
+
1
=
(
n
+
1
)
n
⋯
(
n
−
j
+
1
)
j
+
1
=
(
n
+
1
)
j
+
1
‾
j
+
1
\begin{aligned} \sum_{i=1}^ni^{\underline{j}}&=j(j-1)\cdots 1+(j+1)j\cdots 2+\cdots+n(n-1)\cdots(n-j+1)\\ &=\frac{(j+1)j\cdots 1-j(j-1)\cdots 0}{j+1}+\frac{(j+2)(j+1)\cdots 2-(j+1)j\cdots 1}{j+1}\\ & \ \ \ \ +\cdots+\frac{(n+1)n\cdots(n-j+1)-n(n-1)\cdots(n-j)}{j+1}\\ &=\frac{(n+1)n\cdots(n-j+1)}{j+1}\\ &=\frac{(n+1)^{\underline{j+1}}}{j+1} \end{aligned}
i=1∑nij=j(j−1)⋯1+(j+1)j⋯2+⋯+n(n−1)⋯(n−j+1)=j+1(j+1)j⋯1−j(j−1)⋯0+j+1(j+2)(j+1)⋯2−(j+1)j⋯1 +⋯+j+1(n+1)n⋯(n−j+1)−n(n−1)⋯(n−j)=j+1(n+1)n⋯(n−j+1)=j+1(n+1)j+1
由此可以得知
f
(
n
)
f(n)
f(n) 是
k
+
1
k+1
k+1 次多项式。那么我们需要
k
+
2
k+2
k+2 个点,利用快速幂加前缀和
O
(
k
l
o
g
k
)
O(klogk)
O(klogk) 得到
(
0
,
f
(
0
)
)
,
(
1
,
f
(
1
)
)
,
⋯
(
k
+
1
,
f
(
k
+
1
)
)
\Big(0,f(0)\Big),\Big(1,f(1)\Big),\cdots \Big(k+1,f(k+1)\Big)
(0,f(0)),(1,f(1)),⋯(k+1,f(k+1))。再利用拉格朗日插值
O
(
k
)
O(k)
O(k) 得到
f
(
n
0
)
f(n_0)
f(n0)。
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N = 1e3 + 50;
const int MOD = 9999991;
int n, m, f[N], g[N], x[N];
int pre[N], suf[N], fac, fnv[N];
int fpow (int a, int p) {
int res = 1;
for (; p; p >>= 1, a = (ll)a * a % MOD)
if (p & 1) res = (ll)res * a % MOD;
return res;
}
int lagrange (int n, int* x, int* y, int k) {
if (k <= n) return y[k];
int ans = 0;
pre[0] = (k - x[0]) % MOD, suf[n + 1] = 1;
for (int i = 1; i <= n; i ++) pre[i] = pre[i - 1] * (ll)(k - x[i]) % MOD;
for (int i = n; i >= 0; i --) suf[i] = suf[i + 1] * (ll)(k - x[i]) % MOD;
for (int i = 0; i <= n; i ++) {
int up = (ll)(i == 0 ? 1 : pre[i - 1]) * suf[i + 1] % MOD;
int down = (ll)fnv[i] * (n - i & 1 ? -1 : 1) * fnv[n - i] % MOD;
ans = (ans + (ll)y[i] * up % MOD * down % MOD) % MOD;
}
return (ans % MOD + MOD) % MOD;
}
int main () {
int T; scanf ("%d", &T);
fac = 1;
for (int i = 1; i < N; i ++) fac = (ll)fac * i % MOD;
fnv[N - 1] = fpow (fac, MOD - 2);
for (int i = N - 1; i > 0; i --)
fnv[i - 1] = fnv[i] * (ll)i % MOD;
while (T --) {
scanf ("%d%d", &n, &m);
for (int i = 0; i <= n; i ++) {
scanf ("%d", &f[i]);
x[i] = i;
}
f[n + 1] = lagrange(n, x, f, n + 1);
g[0] = f[0], x[n + 1] = n + 1;
for (int i = 1; i <= n + 1; i ++)
g[i] = (g[i - 1] + f[i]) % MOD;
while (m --) {
int lef, rig; scanf ("%d%d", &lef, &rig);
int ans = lagrange (n + 1, x, g, lef - 1);
ans = lagrange (n + 1, x, g, rig) - ans;
printf("%d\n", (ans % MOD + MOD) % MOD);
}
}
return 0;
}
例题5:https://ac.nowcoder.com/acm/contest/26012/E
求
∑
i
=
L
R
f
(
i
)
\sum_{i=L}^Rf(i)
∑i=LRf(i),那么可以联想到
S
(
R
)
−
S
(
L
−
1
)
S(R)-S(L-1)
S(R)−S(L−1),其中
S
(
x
)
=
f
(
1
)
+
f
(
2
)
+
⋯
+
f
(
x
)
=
∑
i
=
1
x
f
(
i
)
S(x)=f(1)+f(2)+\cdots+f(x)=\sum_{i=1}^xf(i)
S(x)=f(1)+f(2)+⋯+f(x)=i=1∑xf(i)
是
n
+
1
n+1
n+1 次多项式,类似于积分,那么我们只需要
n
+
2
n+2
n+2 个点值就可以得到
S
(
x
)
S(x)
S(x)。并且
O
(
n
)
O(n)
O(n) 计算
S
(
R
)
S(R)
S(R) 和
S
(
L
−
1
)
S(L-1)
S(L−1)。
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N = 1e3 + 50;
const int MOD = 9999991;
int n, m, f[N], g[N], x[N];
int pre[N], suf[N], fac, fnv[N];
int fpow (int a, int p) {
int res = 1;
for (; p; p >>= 1, a = (ll)a * a % MOD)
if (p & 1) res = (ll)res * a % MOD;
return res;
}
int lagrange (int n, int* x, int* y, int k) {
if (k <= n) return y[k];
int ans = 0;
pre[0] = (k - x[0]) % MOD, suf[n + 1] = 1;
for (int i = 1; i <= n; i ++) pre[i] = pre[i - 1] * (ll)(k - x[i]) % MOD;
for (int i = n; i >= 0; i --) suf[i] = suf[i + 1] * (ll)(k - x[i]) % MOD;
for (int i = 0; i <= n; i ++) {
int up = (ll)(i == 0 ? 1 : pre[i - 1]) * suf[i + 1] % MOD;
int down = (ll)fnv[i] * (n - i & 1 ? -1 : 1) * fnv[n - i] % MOD;
ans = (ans + (ll)y[i] * up % MOD * down % MOD) % MOD;
}
return (ans % MOD + MOD) % MOD;
}
int main () {
int T; scanf ("%d", &T);
fac = 1;
for (int i = 1; i < N; i ++) fac = (ll)fac * i % MOD;
fnv[N - 1] = fpow (fac, MOD - 2);
for (int i = N - 1; i > 0; i --)
fnv[i - 1] = fnv[i] * (ll)i % MOD;
while (T --) {
scanf ("%d%d", &n, &m);
for (int i = 0; i <= n; i ++) {
scanf ("%d", &f[i]);
x[i] = i;
}
f[n + 1] = lagrange(n, x, f, n + 1);
g[0] = f[0], x[n + 1] = n + 1;
for (int i = 1; i <= n + 1; i ++)
g[i] = (g[i - 1] + f[i]) % MOD;
while (m --) {
int lef, rig; scanf ("%d%d", &lef, &rig);
int ans = lagrange (n + 1, x, g, lef - 1);
ans = lagrange (n + 1, x, g, rig) - ans;
printf("%d\n", (ans % MOD + MOD) % MOD);
}
}
return 0;
}
例题6:https://ac.nowcoder.com/acm/contest/26012/F
定义
f
(
i
,
x
)
f(i,x)
f(i,x) 为
b
i
=
x
b_i=x
bi=x 的方案数,
[
i
,
n
]
[i,n]
[i,n] 已经确立,那么
f
(
i
,
x
)
=
∑
j
=
a
i
+
1
x
f
(
i
+
1
,
j
)
f(i,x)=\sum_{j=a_{i+1}}^xf(i+1,j)
f(i,x)=j=ai+1∑xf(i+1,j)
边界情况
a
n
≤
x
≤
R
,
f
(
n
,
x
)
=
1
a_n\le x\le R,\ f(n,x)=1
an≤x≤R, f(n,x)=1,其中
a
i
:
=
m
a
x
{
a
i
,
a
i
+
1
,
⋯
,
a
n
}
a_i:=max\{a_i,a_{i+1},\cdots,a_n\}
ai:=max{ai,ai+1,⋯,an}。根据经验,可以知道
f
(
i
,
x
)
f(i,x)
f(i,x) 是
n
−
i
n-i
n−i 次多项式,每次求出
f
(
i
+
1
,
x
)
f(i+1,x)
f(i+1,x) ,再通过拉格朗日插值求出
f
(
i
,
x
)
f(i,x)
f(i,x)。
S
(
i
+
1
,
x
)
=
∑
j
=
1
x
f
(
i
+
1
,
j
)
⇒
f
(
i
,
x
)
=
S
(
i
+
1
,
x
)
−
S
(
i
+
1
,
a
i
+
1
−
1
)
S(i+1,x)=\sum_{j=1}^xf(i+1,j)\\ \Rightarrow f(i,x)=S(i+1,x)-S(i+1,a_{i+1}-1)
S(i+1,x)=j=1∑xf(i+1,j)⇒f(i,x)=S(i+1,x)−S(i+1,ai+1−1)
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N = 5e3 + 50;
const int MOD = 1e9 + 7;
int f[N], pre[N], suf[N];
int n, r, a[N], fnv[N], x[N];
int fpow (int a, int p) {
int res = 1;
for (; p; p >>= 1, a = (ll)a * a % MOD)
if (p & 1) res = (ll)res * a % MOD;
return res;
}
int lagrange (int n, int* x, int* y, int k) {
int ans = 0;
pre[0] = (k - x[0]), suf[n + 1] = 1;
for (int i = 1; i <= n; i ++) pre[i] = (ll)pre[i - 1] * (k - x[i]) % MOD;
for (int i = n; i >= 0; i --) suf[i] = (ll)suf[i + 1] * (k - x[i]) % MOD;
for (int i = 0; i <= n; i ++) {
int up = (ll)(i == 0 ? 1 : pre[i - 1]) * suf[i + 1] % MOD;
int down = (ll)fnv[i] * fnv[n - i] * (n - i & 1 ? -1 : 1) % MOD;
ans = (ans + (ll)y[i] * down % MOD * up % MOD) % MOD;
}
return (ans % MOD + MOD) % MOD;
}
int main () {
scanf ("%d%d", &n, &r);
for (int i = 1; i <= n; i ++)
scanf ("%d", &a[i]);
for (int i = n - 1; i >= 1; i --)
a[i] = max (a[i], a[i + 1]);
for (int i = 0; i <= n; i ++) f[i] = 1;
for (int i = 0; i < N; i ++) x[i] = i;
int fac = 1;
for (int i = 1; i < N; i ++) fac = (ll)fac * i % MOD;
fnv[N - 1] = fpow (fac, MOD - 2);
for (int i = N - 1; i >= 0; i --) fnv[i - 1] = (ll)fnv[i] * i % MOD;
for (int i = n - 1; i >= 0; i --) {
for (int j = 1; j <= n; j ++) {
f[j] = (f[j] + f[j - 1]) % MOD;
}
int temp = lagrange (n - i, x, f, a[i + 1] - 1);
for (int j = 0; j <= n; j ++) {
f[j] = (f[j] - temp) % MOD;
}
}
printf("%d\n", lagrange(n, x, f, r));
return 0;
}