贝尔数(Bell Number)
贝尔数表示的是基数为 n n n的集合的划分成两两不想交的子集的划分方案数,记作 B n B_n Bn
贝尔数的前几项为: 1 , 1 , 2 , 5 , 15 , 52 , 203.......... 1,1,2,5,15,52,203.......... 1,1,2,5,15,52,203..........
1、组合意义推导
我们考虑第 n + 1 n+1 n+1个元素 b n + 1 b_{n+1} bn+1,那么可以对其分类讨论:
(1)、如果 b n + 1 b_{n+1} bn+1单独分为一个子集,则不用从 n n n个元素中选择元素出来和其组成一个子集,故此时的方案数为: C n 0 B n C_n^0B_n Cn0Bn
(2)、如果 b n + 1 b_{n+1} bn+1和一个元素组成一个子集,则需要从 n n n个元素中选出一个元素与之组成一个子集,故此时的方案数为: C n 1 B n − 1 C_n^1B_{n-1} Cn1Bn−1
(3)、如果 b n + 1 b_{n+1} bn+1和两个元素组成一个子集,则需要从 n n n个元素中选出两个元素与之组成一个子集,故此时的方案数为: C n 2 B n − 2 C_n^2B_{n-2} Cn2Bn−2
. . . . . . . . . . .......... ..........
以此类推,那么我们可以得到递推式:
B
n
+
1
=
∑
k
=
0
n
C
n
k
B
n
−
k
=
∑
k
=
0
n
C
n
k
B
k
\begin{aligned} B_{n+1}&=\sum_{k=0}^{n}C_n^{k}B_{n-k}\\ &=\sum_{k=0}^{n}C_n^{k}B_k \end{aligned}
Bn+1=k=0∑nCnkBn−k=k=0∑nCnkBk
另一种意义:
我们还可以考虑另外一种组合意义:
(1)、将 n n n个元素分成 1 1 1个子集,此时方案数: { n 1 } \begin{Bmatrix}n\\\\1\end{Bmatrix} ⎩ ⎨ ⎧n1⎭ ⎬ ⎫
(2)、将 n n n个元素分成 2 2 2个子集,此时方案数: { n 2 } \begin{Bmatrix}n\\\\2\end{Bmatrix} ⎩ ⎨ ⎧n2⎭ ⎬ ⎫
(3)、将 n n n个元素分成 3 3 3个子集,此时方案数: { n 3 } \begin{Bmatrix}n\\\\3\end{Bmatrix} ⎩ ⎨ ⎧n3⎭ ⎬ ⎫
. . . . . . . . ........ ........
以此类推,我们可以得到新的递推式:
B
n
=
∑
k
=
0
n
{
n
k
}
B_n=\sum_{k=0}^n\begin{Bmatrix}n\\\\k\end{Bmatrix}
Bn=k=0∑n⎩
⎨
⎧nk⎭
⎬
⎫
这种推导方式也推导出了,贝尔数与第二类斯特林数的关系
2、生成函数推导
设 b ( x ) 为序列【 B 1 , B 2 , B 3 . . . . B n 】的生成函数 则 b ( x ) = ∑ i = 0 ∞ B i x i i ! ∴ b ( x ) = ∑ i = 0 ∞ ∑ k = 0 i − 1 C i − 1 k B k x i i ! 改变 i 与 k 的求和顺序 ∴ b ( x ) = ∑ k = 0 ∞ B k ∑ i = k + 1 ∞ C i − 1 k x i i ! = ∑ k = 0 ∞ B k ∑ i = k + 1 ∞ ( i − 1 ) ! ( i − k − 1 ) ! k ! x i i ! = ∑ k = 0 ∞ B k ∑ i = k + 1 ∞ 1 i ( i − k − 1 ) ! k ! x i 对 b ( x ) 求关于 x 的导数 ∴ b ‘ ( x ) = ∑ k = 0 ∞ B k ∑ i = k + 1 ∞ 1 ( i − k − 1 ) ! k ! x i − 1 = ∑ k = 0 ∞ B k k ! ∑ i = k + 1 ∞ 1 ( i − k − 1 ) ! x i − 1 = ∑ k = 0 ∞ B k k ! ∑ i = 0 ∞ 1 i ! x i − 1 x k + 1 = ∑ k = 0 ∞ B k x k k ! ∑ i = 0 ∞ 1 i ! x i = b ( x ) e x ∴ b ‘ ( x ) b ( x ) = e x 又 ∵ l n ( f ( x ) ) = f ‘ ( x ) f ( x ) ∴ l n ( b ( x ) ) = e x ∴ b ( x ) = e e x + c , c 为常数 代入计算得 C = − 1 ∴ b ( x ) = e e x − 1 设b(x)为序列【B_1,B_2,B_3....B_n】的生成函数 \\ 则b(x)=\sum_{i=0}^{\infty}B_i\frac{x^i}{i!}\\ \therefore b(x)=\sum_{i=0}^{\infty}\sum_{k=0}^{i-1}C_{i-1}^{k}B_k\frac{x^i}{i!}\\改变i与k的求和顺序\\ \begin{aligned} \therefore b(x)&=\sum_{k=0}^{\infty}B_k\sum_{i=k+1}^{\infty}C_{i-1}^{k}\frac{x^i}{i!}\\ &=\sum_{k=0}^{\infty}B_k\sum_{i=k+1}^{\infty}\frac{(i-1)!}{(i-k-1)!k!}\frac{x^i}{i!}\\ &=\sum_{k=0}^{\infty}B_k\sum_{i=k+1}^{\infty}\frac{1}{i(i-k-1)!k!}x^i\\ \end{aligned} \\ 对b(x)求关于x的导数\\ \begin{aligned} \therefore b^{`}(x)&=\sum_{k=0}^{\infty}B_k\sum_{i=k+1}^{\infty}\frac{1}{(i-k-1)!k!}x^{i-1}\\ &=\sum_{k=0}^{\infty}\frac{B_k}{k!}\sum_{i=k+1}^{\infty}\frac{1}{(i-k-1)!}x^{i-1}\\ &=\sum_{k=0}^{\infty}\frac{B_k}{k!}\sum_{i=0}^{\infty}\frac{1}{i!}x^{i-1}x^{k+1}\\ &=\sum_{k=0}^{\infty}\frac{B_kx^k}{k!}\sum_{i=0}^{\infty}\frac{1}{i!}x^i\\ &=b(x)e^x\\ \end{aligned}\\ \therefore \frac{b^{`}(x)}{b(x)}=e^x \\ 又\because ln(f(x))=\frac{f^{`}(x)}{f(x)}\\ \therefore ln(b(x))=e^x \\ \therefore b(x)=e^{e^x+c},c为常数 \\代入计算得C=-1\\ \therefore b(x)=e^{e^x-1} 设b(x)为序列【B1,B2,B3....Bn】的生成函数则b(x)=i=0∑∞Bii!xi∴b(x)=i=0∑∞k=0∑i−1Ci−1kBki!xi改变i与k的求和顺序∴b(x)=k=0∑∞Bki=k+1∑∞Ci−1ki!xi=k=0∑∞Bki=k+1∑∞(i−k−1)!k!(i−1)!i!xi=k=0∑∞Bki=k+1∑∞i(i−k−1)!k!1xi对b(x)求关于x的导数∴b‘(x)=k=0∑∞Bki=k+1∑∞(i−k−1)!k!1xi−1=k=0∑∞k!Bki=k+1∑∞(i−k−1)!1xi−1=k=0∑∞k!Bki=0∑∞i!1xi−1xk+1=k=0∑∞k!Bkxki=0∑∞i!1xi=b(x)ex∴b(x)b‘(x)=ex又∵ln(f(x))=f(x)f‘(x)∴ln(b(x))=ex∴b(x)=eex+c,c为常数代入计算得C=−1∴b(x)=eex−1
3、性质
引理:
B n + m = ∑ i = 0 n C n i B i ∑ j = 0 m j n − i { m j } B_{n+m}=\sum_{i=0}^{n}C_n^iB_i\sum_{j=0}^mj^{n-i}\begin{Bmatrix}m\\\\j\end{Bmatrix} Bn+m=i=0∑nCniBij=0∑mjn−i⎩ ⎨ ⎧mj⎭ ⎬ ⎫
组合意义证明:
把 n + m n+m n+m个元素分成 n n n个元素和 m m m个元素;首先枚举 m m m个元素划分成 j j j个集合,此时有 { m j } \begin{Bmatrix}m\\\\j\end{Bmatrix} ⎩ ⎨ ⎧mj⎭ ⎬ ⎫;再枚举 n n n个元素中选 i i i个出来任意划分 成子集,此时有 C n i B i C_n^iB_i CniBi;还剩下 n − i n-i n−i个元素,那么这 n − i n-i n−i个元素必然在 j j j个集合中,此时有 j n − i j^{n-i} jn−i,然后根据乘法原理将其相乘,便是最后的结果
3.1、Dobinski公式:
B n = 1 e ∑ k = 0 ∞ k n k ! B_n=\frac{1}{e}\sum_{k=0}^{\infty}\frac{k^n}{k!} Bn=e1k=0∑∞k!kn
期望为 n n n的泊松分数的 n n n次矩
3.2、Touchard同余
若对 ∀ p ∈ p r i m e s , 有 B p + n ≡ B n + B n + 1 ( m o d p ) 若对\forall p\in primes,有B_{p+n}\equiv B_n+B_{n+1}(mod\space p) 若对∀p∈primes,有Bp+n≡Bn+Bn+1(mod p)
证明:
B
n
+
p
=
∑
i
=
0
n
C
n
i
B
i
∑
j
=
0
p
j
n
−
i
{
p
j
}
=
∑
j
=
0
p
{
p
j
}
∑
i
=
0
n
C
n
i
B
i
j
n
−
i
又
∵
对于
∀
j
∈
(
1
,
p
)
且
p
∈
p
r
i
m
e
s
,
有
{
p
j
}
≡
0
(
m
o
d
p
)
∴
B
n
+
p
m
o
d
p
=
∑
j
=
0
p
{
p
j
}
∑
i
=
0
n
C
n
i
B
i
j
n
−
i
m
o
d
p
=
(
∑
j
=
0
p
{
p
j
}
m
o
d
p
)
(
∑
i
=
0
n
C
n
i
B
i
j
n
−
i
m
o
d
p
)
m
o
d
p
=
(
{
p
1
}
+
0
+
{
p
p
}
)
(
∑
i
=
0
n
C
n
i
B
i
j
n
−
i
m
o
d
p
)
m
o
d
p
=
(
{
p
1
}
(
∑
i
=
0
n
C
n
i
B
i
1
n
−
i
m
o
d
p
)
m
o
d
p
+
{
p
p
}
(
∑
i
=
0
n
C
n
i
B
i
p
n
−
i
m
o
d
p
)
m
o
d
p
)
m
o
d
p
=
(
∑
i
=
0
n
C
n
i
B
i
1
n
−
i
+
∑
i
=
0
n
C
n
i
B
i
p
n
−
i
)
m
o
d
p
=
(
∑
i
=
0
n
C
n
i
B
i
+
∑
i
=
0
n
C
n
i
B
i
p
n
−
i
)
m
o
d
p
=
(
B
n
+
1
+
∑
i
=
0
n
C
n
i
B
i
p
n
−
i
)
m
o
d
p
又
∵
对
∀
i
>
0
,
p
i
≡
0
(
m
o
d
p
)
∴
∑
i
=
0
n
C
n
i
B
i
p
n
−
i
m
o
d
p
=
C
n
n
B
n
p
n
−
n
m
o
d
p
=
B
n
m
o
d
p
∴
B
n
+
p
≡
B
n
+
1
+
B
n
(
m
o
d
p
)
故:结论得证
\begin{aligned} B_{n+p}&=\sum_{i=0}^{n}C_n^iB_i\sum_{j=0}^pj^{n-i}\begin{Bmatrix}p\\\\j\end{Bmatrix}\\ &=\sum_{j=0}^p\begin{Bmatrix}p\\\\j\end{Bmatrix}\sum_{i=0}^{n}C_n^iB_ij^{n-i}\\ 又\because 对于\forall &j\in(1,p)且p\in primes,有 \begin{Bmatrix}p\\\\j\end{Bmatrix}\equiv0(mod\space p) \\ \therefore B_{n+p}\space mod\space p&=\sum_{j=0}^p\begin{Bmatrix}p\\\\j\end{Bmatrix}\sum_{i=0}^{n}C_n^iB_ij^{n-i}\space mod\space p\\ &=(\sum_{j=0}^p\begin{Bmatrix}p\\\\j\end{Bmatrix}\space mod\space p)(\sum_{i=0}^{n}C_n^iB_ij^{n-i}\space mod\space p)\space mod\space p\\ &=(\begin{Bmatrix}p\\\\1\end{Bmatrix}+0+\begin{Bmatrix}p\\\\p\end{Bmatrix})(\sum_{i=0}^{n}C_n^iB_ij^{n-i}\space mod\space p)\space mod\space p\\ &=(\begin{Bmatrix}p\\\\1\end{Bmatrix}(\sum_{i=0}^{n}C_n^iB_i1^{n-i}\space mod\space p)\space mod\space p+\begin{Bmatrix}p\\\\p\end{Bmatrix}(\sum_{i=0}^{n}C_n^iB_ip^{n-i}\space mod\space p)\space mod\space p)\space mod\space p\\ &=(\sum_{i=0}^{n}C_n^iB_i1^{n-i}+\sum_{i=0}^{n}C_n^iB_ip^{n-i})\space mod\space p\\ &=(\sum_{i=0}^{n}C_n^iB_i+\sum_{i=0}^{n}C_n^iB_ip^{n-i})\space mod\space p\\ &=(B_{n+1}+\sum_{i=0}^{n}C_n^iB_ip^{n-i})\space mod\space p\\ 又\because 对\forall i>0,&p^i\equiv 0(mod \space p )\\ \therefore \sum_{i=0}^{n}C_n^i&B_ip^{n-i}\space mod\space p=C_n^nB_np^{n-n}\space mod\space p=B_n\space mod\space p\\ \therefore B_{n+p}&\equiv B_{n+1}+B_n(\space mod\space p)\\故:结论得证 \end{aligned}
Bn+p又∵对于∀∴Bn+p mod p又∵对∀i>0,∴i=0∑nCni∴Bn+p故:结论得证=i=0∑nCniBij=0∑pjn−i⎩
⎨
⎧pj⎭
⎬
⎫=j=0∑p⎩
⎨
⎧pj⎭
⎬
⎫i=0∑nCniBijn−ij∈(1,p)且p∈primes,有⎩
⎨
⎧pj⎭
⎬
⎫≡0(mod p)=j=0∑p⎩
⎨
⎧pj⎭
⎬
⎫i=0∑nCniBijn−i mod p=(j=0∑p⎩
⎨
⎧pj⎭
⎬
⎫ mod p)(i=0∑nCniBijn−i mod p) mod p=(⎩
⎨
⎧p1⎭
⎬
⎫+0+⎩
⎨
⎧pp⎭
⎬
⎫)(i=0∑nCniBijn−i mod p) mod p=(⎩
⎨
⎧p1⎭
⎬
⎫(i=0∑nCniBi1n−i mod p) mod p+⎩
⎨
⎧pp⎭
⎬
⎫(i=0∑nCniBipn−i mod p) mod p) mod p=(i=0∑nCniBi1n−i+i=0∑nCniBipn−i) mod p=(i=0∑nCniBi+i=0∑nCniBipn−i) mod p=(Bn+1+i=0∑nCniBipn−i) mod ppi≡0(mod p)Bipn−i mod p=CnnBnpn−n mod p=Bn mod p≡Bn+1+Bn( mod p)
TIPS:
证明:
对于
∀
j
∈
(
1
,
p
)
且
p
∈
p
r
i
m
e
s
,
有
{
p
j
}
≡
0
(
m
o
d
p
)
对于\forall j\in(1,p)且p\in primes,有 \begin{Bmatrix}p\\\\j\end{Bmatrix}\equiv0(mod\space p)
对于∀j∈(1,p)且p∈primes,有⎩
⎨
⎧pj⎭
⎬
⎫≡0(mod p)
由第二类斯特林数的通项公式得:
{
p
j
}
=
∑
k
=
0
m
(
−
1
)
k
k
!
(
j
−
k
)
p
(
j
−
k
)
!
∴
{
p
j
}
m
o
d
p
=
(
∑
k
=
0
j
(
−
1
)
k
k
!
(
j
−
k
)
p
(
j
−
k
)
!
)
m
o
d
p
=
(
∑
k
=
0
j
(
(
−
1
)
k
k
!
m
o
d
p
)
(
(
j
−
k
)
p
(
j
−
k
)
!
m
o
d
p
)
m
o
d
p
)
m
o
d
p
由费马小定理可得
:
对于
∀
p
∈
p
r
i
m
e
s
且
∀
a
满足
g
c
d
(
a
,
p
)
=
1
,都有
a
p
−
1
≡
1
(
m
o
d
p
)
∴
(
j
−
k
)
p
(
j
−
k
)
!
m
o
d
p
=
(
j
−
k
)
(
j
−
k
)
!
m
o
d
p
=
1
(
j
−
k
−
1
)
!
m
o
d
p
∴
{
p
j
}
m
o
d
p
=
(
∑
k
=
0
j
(
−
1
)
k
k
!
1
(
j
−
k
−
1
)
!
)
m
o
d
p
这个时候不妨考虑生成函数
A
(
x
)
=
∑
i
=
0
∞
(
∑
k
=
0
i
(
−
1
)
k
k
!
1
(
j
−
k
−
1
)
!
)
x
i
通过观察发现
A
(
x
)
其实是卷积的形式
不妨设
B
(
x
)
=
∑
i
=
0
∞
(
−
1
)
i
i
!
x
i
,
C
(
x
)
=
∑
i
=
0
∞
1
(
i
−
1
)
!
)
x
i
则
C
(
x
)
=
∑
i
=
0
∞
1
(
i
−
1
)
!
)
x
i
=
x
∑
i
=
0
∞
1
(
i
−
1
)
!
x
i
−
1
=
x
∑
i
=
0
∞
1
i
!
x
i
由
e
x
的泰勒展开可得
e
x
=
∑
i
=
0
∞
1
i
!
x
i
,
e
−
x
=
(
−
1
)
i
i
!
x
i
那么可得
A
(
x
)
=
B
(
x
)
C
(
x
)
=
e
−
x
x
e
x
=
x
∴
A
(
x
)
=
x
∴
对于
i
>
1
,
∑
k
=
0
i
(
−
1
)
k
k
!
1
(
j
−
k
−
1
)
!
=
0
∴
对于
∀
j
∈
(
1
,
p
)
且
p
∈
p
r
i
m
e
s
,
有
{
p
j
}
≡
0
(
m
o
d
p
)
结论得证
由第二类斯特林数的通项公式得:\begin{Bmatrix}p\\\\j\end{Bmatrix}=\sum_{k=0}^m\frac{(-1)^k}{k!}\frac{(j-k)^p}{(j-k)!}\\ \begin{aligned} \therefore \begin{Bmatrix}p\\\\j\end{Bmatrix}\space mod\space p&=(\sum_{k=0}^j\frac{(-1)^k}{k!}\frac{(j-k)^p}{(j-k)!})\space mod\space p\\ &=(\sum_{k=0}^j(\frac{(-1)^k}{k!}\space mod\space p)(\frac{(j-k)^p}{(j-k)!}\space mod\space p)\space mod\space p)\space mod\space p\\ \end{aligned} \\ 由费马小定理可得:对于\forall p\in primes且\forall a 满足gcd(a,p)=1,都有a^{p-1}\equiv 1(\space mod\space p)\\ \therefore \frac{(j-k)^p}{(j-k)!}\space mod\space p=\frac{(j-k)}{(j-k)!}\space mod\space p=\frac{1}{(j-k-1)!}\space mod\space p\\ \therefore \begin{Bmatrix}p\\\\j\end{Bmatrix}\space mod\space p=(\sum_{k=0}^j\frac{(-1)^k}{k!}\frac{1}{(j-k-1)!})\space mod\space p\\ 这个时候不妨考虑生成函数A(x)=\sum_{i=0}^{\infty}(\sum_{k=0}^i\frac{(-1)^k}{k!}\frac{1}{(j-k-1)!})x^i\\ 通过观察发现A(x)其实是卷积的形式\\ 不妨设B(x)=\sum_{i=0}^{\infty}\frac{(-1)^i}{i!}x^i,C(x)=\sum_{i=0}^{\infty}\frac{1}{(i-1)!})x^i \\ 则C(x)=\sum_{i=0}^{\infty}\frac{1}{(i-1)!})x^i=x\sum_{i=0}^{\infty}\frac{1}{(i-1)!}x^{i-1}=x\sum_{i=0}^{\infty}\frac{1}{i!}x^i\\ 由e^x的泰勒展开可得e^x=\sum_{i=0}^{\infty}\frac{1}{i!}x^i,e^{-x}=\frac{(-1)^i}{i!}x^i \\ 那么可得A(x)=B(x)C(x)=e^{-x}xe^x=x \\ \therefore A(x)=x\\ \therefore 对于i>1,\sum_{k=0}^i\frac{(-1)^k}{k!}\frac{1}{(j-k-1)!}=0 \\ \therefore对于\forall j\in(1,p)且p\in primes,有 \begin{Bmatrix}p\\\\j\end{Bmatrix}\equiv0(mod\space p) \\结论得证
由第二类斯特林数的通项公式得:⎩
⎨
⎧pj⎭
⎬
⎫=k=0∑mk!(−1)k(j−k)!(j−k)p∴⎩
⎨
⎧pj⎭
⎬
⎫ mod p=(k=0∑jk!(−1)k(j−k)!(j−k)p) mod p=(k=0∑j(k!(−1)k mod p)((j−k)!(j−k)p mod p) mod p) mod p由费马小定理可得:对于∀p∈primes且∀a满足gcd(a,p)=1,都有ap−1≡1( mod p)∴(j−k)!(j−k)p mod p=(j−k)!(j−k) mod p=(j−k−1)!1 mod p∴⎩
⎨
⎧pj⎭
⎬
⎫ mod p=(k=0∑jk!(−1)k(j−k−1)!1) mod p这个时候不妨考虑生成函数A(x)=i=0∑∞(k=0∑ik!(−1)k(j−k−1)!1)xi通过观察发现A(x)其实是卷积的形式不妨设B(x)=i=0∑∞i!(−1)ixi,C(x)=i=0∑∞(i−1)!1)xi则C(x)=i=0∑∞(i−1)!1)xi=xi=0∑∞(i−1)!1xi−1=xi=0∑∞i!1xi由ex的泰勒展开可得ex=i=0∑∞i!1xi,e−x=i!(−1)ixi那么可得A(x)=B(x)C(x)=e−xxex=x∴A(x)=x∴对于i>1,k=0∑ik!(−1)k(j−k−1)!1=0∴对于∀j∈(1,p)且p∈primes,有⎩
⎨
⎧pj⎭
⎬
⎫≡0(mod p)结论得证
4、贝尔三角形
构造三角矩阵
满足一下递推式:
a
n
,
m
=
{
1
,
n
=
1
,
m
=
1
a
n
−
1
,
n
−
1
,
n
>
1
,
m
=
1
a
n
,
m
−
1
+
a
n
−
1
,
m
−
1
,
n
>
1
,
m
>
1
a_{n,m}= \begin{cases} 1,n=1,m=1\\ a_{n-1,n-1},n>1,m=1\\ a_{n,m-1}+a_{n-1,m-1},n>1,m>1 \end{cases}
an,m=⎩
⎨
⎧1,n=1,m=1an−1,n−1,n>1,m=1an,m−1+an−1,m−1,n>1,m>1
1 1 2 2 3 5 5 7 10 15 15 20 27 37 52 52 67 87 114 151 203 203 255 322 409 523 674 877 \begin{aligned} & 1 \\ & 1\quad\qquad 2 \\ & 2\quad\qquad 3\quad\qquad 5 \\ & 5\quad\qquad 7\quad\qquad 10\qquad 15 \\ & 15\qquad 20\qquad 27\qquad 37\qquad 52 \\ & 52\qquad 67\qquad 87\qquad 114\qquad 151\qquad 203\\ & 203\qquad 255\qquad 322\qquad 409\qquad 523\qquad 674\qquad 877 \\ \end{aligned} 1122355710151520273752526787114151203203255322409523674877
每行的首项为贝尔数
5、快速求贝尔数在模意义下的第 n n n项
题目链接:https://vjudge.csgrandeur.cn/problem/HDU-4767
题意:求 B n m o d 95041567 B_n \space mod \space 95041567 Bn mod 95041567
思路:
首先,我们发现
95041567
95041567
95041567并不是一个素数;但是我们发现
95041567
95041567
95041567是多个不相同的素数的乘积,因此我们可以考虑CRT合并,所以我们可以列出以下方程:
95041567
=
Π
i
p
i
{
x
≡
a
1
(
m
o
d
p
1
)
x
≡
a
2
(
m
o
d
p
2
)
.
.
.
.
x
≡
a
n
(
m
o
d
p
k
)
95041567=\Pi_ip_i\\ \begin{cases} x\equiv a_1(\space mod \space p_1)\\ x\equiv a_2(\space mod \space p_2)\\ ....\\ x\equiv a_n(\space mod \space p_k)\\ \end{cases}
95041567=Πipi⎩
⎨
⎧x≡a1( mod p1)x≡a2( mod p2)....x≡an( mod pk)
因此我们将问题转化为了求出一组$a_i=B_n \space mod \space p_i\$
而对于每一组素数
p
i
p_i
pi,我们可以利用
T
o
u
c
h
a
r
d
同余
Touchard同余
Touchard同余
将
B
n
+
p
≡
B
n
+
1
+
B
n
(
m
o
d
p
)
改写成模意义下的矩阵相乘的形式
:
[
B
n
+
p
−
1
B
n
+
p
B
n
+
p
+
1
⋮
B
n
+
p
+
p
−
1
]
=
[
0
0
0
…
0
1
1
1
0
…
0
0
0
1
1
…
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
…
1
1
]
×
[
B
n
B
n
+
1
B
n
+
2
⋮
B
n
+
p
−
1
]
将 B_{n+p}\equiv B_{n+1}+B_n(\space mod\space p)改写成模意义下的矩阵相乘的形式: \\ \begin{bmatrix} B_{n+p-1}\\ B_{n+p}\\ B_{n+p+1}\\ \vdots \\ B_{n+p+p-1} \end{bmatrix}= \begin{bmatrix} &0 &0 &0 &\dots&0&1\\ &1 &1 &0 &\dots&0&0\\ &0 &1 &1 &\dots&0&0\\ &\vdots &\vdots &\vdots &\ddots &\vdots&\vdots\\ &0 &0 &0 &\dots&1&1\ \end{bmatrix}\times \begin{bmatrix} B_{n}\\ B_{n+1}\\ B_{n+2}\\ \vdots \\ B_{n+p-1} \end{bmatrix}
将Bn+p≡Bn+1+Bn( mod p)改写成模意义下的矩阵相乘的形式:
Bn+p−1Bn+pBn+p+1⋮Bn+p+p−1
=
010⋮0011⋮0001⋮0………⋱…000⋮1100⋮1
×
BnBn+1Bn+2⋮Bn+p−1
那么可以利用矩阵的递推推出第
n
n
n项与前几项的关系:
[
B
n
B
n
+
1
B
n
+
2
⋮
B
n
+
p
−
1
]
=
[
0
0
0
…
0
1
1
1
0
…
0
0
0
1
1
…
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
…
1
1
]
n
p
−
1
×
[
B
0
B
1
B
2
⋮
B
p
−
1
]
\begin{bmatrix} B_{n}\\ B_{n+1}\\ B_{n+2}\\ \vdots \\ B_{n+p-1} \end{bmatrix}= \begin{bmatrix} &0 &0 &0 &\dots&0&1\\ &1 &1 &0 &\dots&0&0\\ &0 &1 &1 &\dots&0&0\\ &\vdots &\vdots &\vdots &\ddots &\vdots&\vdots\\ &0 &0 &0 &\dots&1&1\ \end{bmatrix}^{\frac{n}{p-1}}\times \begin{bmatrix} B_{0}\\ B_{1}\\ B_{2}\\ \vdots \\ B_{p-1} \end{bmatrix}
BnBn+1Bn+2⋮Bn+p−1
=
010⋮0011⋮0001⋮0………⋱…000⋮1100⋮1
p−1n×
B0B1B2⋮Bp−1
我们发现中间那个矩阵可以利用矩阵快速幂求出,那么我们只需要预处理出前
p
i
p_i
pi个
B
i
B_i
Bi,然后利用矩阵乘法就可以得到
B
n
m
o
d
p
i
B_n\space mod \space p_i
Bn mod pi
现在我们已经知道一组 a i a_i ai,那么利用CRT合并得到的结果,便是最后的答案!
#include <iostream>
#include <algorithm>
#include <cmath>
#include <string.h>
#include <vector>
using namespace std;
#define int long long
#define endl "\n"
const int maxn = 1e2 + 10;
const int p = 95041567;
int bell[maxn][maxn], ans[maxn];
//矩阵快速幂
int cnt;
struct jz
{
int m[maxn][maxn];
jz()
{
memset(m, 0, sizeof(m));
}
} ot;
jz mul(jz a, jz b, int mod)
{
jz c;
for (int i = 1; i <= cnt; i++)
{
for (int j = 1; j <= cnt; j++)
{
c.m[i][j] = 0;
for (int k = 1; k <= cnt; k++)
c.m[i][j] = (c.m[i][j] + a.m[i][k] * b.m[k][j] % mod) % mod;
}
}
return c;
}
jz ksm(jz a, int b, int mod)
{
jz ans;
for (int i = 1; i <= cnt; i++)
for (int j = 1; j <= cnt; j++)
ans.m[i][j] = 1 * (i == j);
while (b)
{
if (b % 2)
ans = mul(ans, a, mod);
a = mul(a, a, mod);
b >>= 1;
}
return ans;
}
//模数分解
vector<int> mod;
vector<int> a;
void div(int n)
{
for (int i = 2; i * i <= n; i++)
{
if (n % i == 0)
{
mod.push_back(i);
while (n % i == 0)
{
n /= i;
}
}
}
if (n > 1)
{
mod.push_back(n);
}
}
//CRT合并
void exgcd(int a, int b, int &x, int &y)
{
if (b == 0)
{
x = 1;
y = 0;
return;
}
exgcd(b, a % b, x, y);
int temp = x;
x = y;
y = temp - a / b * y;
}
int CRT(vector<int> a, vector<int> mod)
{
int n = mod.size();
int M = 1, ans = 0;
for (int i = 0; i < n; i++)
{
M *= mod[i];
}
for (int i = 0; i < n; i++)
{
int m = M / mod[i];
int b, y;
exgcd(m, mod[i], b, y);
ans = (ans + a[i] * m * b % M) % M;
}
return (ans % M + M) % M;
}
//贝尔三角形递推前几项贝尔数
void init()
{
bell[1][1] = 1;
for (int i = 2; i < maxn; i++)
{
for (int j = 1; j < maxn; j++)
{
if (j == 1)
{
bell[i][j] = bell[i - 1][i - 1];
}
else
{
bell[i][j] = (bell[i][j - 1] + bell[i - 1][j - 1]) % p;
}
}
}
}
signed main()
{
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
div(95041567);
init();
int t;
cin >> t;
while (t--)
{
int n;
cin >> n;
a.clear();
for (auto mo : mod)
{
cnt = mo;
jz temp;
temp.m[1][mo] = 1;
for (int i = 2; i <= mo; i++)
{
for (int j = 1; j <= mo; j++)
{
if (i == j || i == j + 1)
{
temp.m[i][j] = 1;
}
}
}
ot = ksm(temp, n / (mo - 1), mo);
memset(ans, 0, sizeof(ans));
for (int i = 1; i <= mo; i++)
{
for (int j = 1; j <= mo; j++)
{
ans[i] = (ans[i] + ot.m[i][j] * bell[j][1] % mo) % mo;
}
}
a.push_back(ans[n % (mo - 1) + 1]);
}
int ans = CRT(a, mod);
cout << ans << endl;
}
return 0;
}