先看一下差分序列和斯特林數。https://riteme.github.io/blog/2016-11-29/delta-and-stirling.html
數學上，伯努利數 \(B_n\)的第一次發現與下述數列和的公式有關：$$\sum_{k=1} ^ {n} k ^ m = 1 ^ m + 2 ^ m + 3 ^ m + \dots + n ^ m$$其中\(m\)為固定的任意正整數。
這個數列的和的公式必定是變量為\(m\)，次數為\(n+1\)的多項式，稱為伯努利多項式。伯努利多項式的系數與伯努利數有密切關系如下。 $$\sum_{k=1} ^ {n} k ^ m = \frac {1} {m+1} \sum_{k=0} ^ m \dbinom{m+1}{k} B_k n ^{m+1-k}$$。
舉個例子： 取\(m=1\)，我們有\(1+2+\dots+n=\frac{1}{2}(B_0n ^ 2 + 2 B_1 ^ {+} n ^ 1 )=\frac{1}{2}(n ^ 2 + n)\)。
伯努利數 滿足條件\(B_0=1\)，且$$\sum_{k=0} ^ n \dbinom{n+1}{k} B_k = 0$$於是有遞歸形式定義伯努利數：$$B_n=-\frac{1}{n-1} \sum_{k=0} ^ {n-1} \dbinom{n+1}{k} B_k$$
可以計算前5項伯努利數。
\(B_0 = 1 , B_0 = -\frac{1}{2}，B_2 = \frac{1}{6}，B_3 = 0 ，B_4 = - \frac{1}{30}，B_5 = 0\)。
代碼：
伯努利數取模
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 10000;
const int mod = 1e9+7;
ll B[maxn]; //伯努利數
ll C[maxn][maxn]; //組合數
ll inv[maxn]; //逆元（計算伯努利數）
void init()
{
//預處理組合數
for (int i = 0; i < maxn; i++)
{
C[i][0] = C[i][i] = 1;
for (int k = 1; k < i; k++)
{
C[i][k] = (C[i - 1][k] % mod + C[i - 1][k - 1] % mod) % mod;
}
}
//預處理逆元
inv[1] = 1;
for (int i = 2; i < maxn; i++)
{
inv[i] = (mod - mod / i) * inv[mod % i] % mod;
}
//預處理伯努利數
B[0] = 1;
for (int i = 1; i < maxn; i++)
{
ll ans = 0;
if (i == maxn - 1)
break;
for (int k = 0; k < i; k++)
{
ans += C[i + 1][k] * B[k];
ans %= mod;
}
ans = (ans * (-inv[i + 1]) % mod + mod) % mod;
B[i] = ans;
}
}
int main(int argc, char const *argv[]) {
init();
std::cout << B[0] << '\n';
return 0;
}
計算伯努利數
將每個伯努利數的分子和分母分別存在\(B[k][0]\) 和 \(B[k][1]\)，負號在分母上。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 1000;
const int mod = 1e9+7;
ll B[maxn][maxn]; //伯努利數
ll C[maxn][maxn]; //組合數
ll inv[maxn]; //逆元（計算伯努利數）
ll lcm(ll a, ll b)
{
return a / __gcd(a, b) * b;
}
void bernoulli()
{
// 預處理組合數
for (int i = 0; i < maxn; i++)
{
C[i][0] = C[i][i] = 1;
for (int k = 1; k < i; k++)
{
C[i][k] = (C[i - 1][k] % mod + C[i - 1][k - 1] % mod) % mod;
}
}
// 預處理伯努利數 B[i][0] 是分子，B[i][1] 是分母
ll s[2], b[2], l, g;
B[0][0] = B[0][1] = 1;
for (int m = 1; m < maxn; m++)
{
s[0] = 1;
s[1] = 1;
for (int i = 1; i < m; i++)
{
b[0] = C[m + 1][i] * B[i][0];
b[1] = B[i][1];
l = lcm(s[1], b[1]);
s[0] = l / s[1] * s[0] + l / b[1] * b[0];
s[1] = l;
}
s[0] = -s[0];
if (s[0])
{
g = __gcd(s[0], s[1] * C[m + 1][m]);
B[m][0] = s[0] / g;
B[m][1] = s[1] * C[m + 1][m] / g;
}
else
B[m][0] = 0, B[m][1] = 1;
}
}
int main(int argc, char const *argv[]) {
bernoulli();
std::cout << B[1][0] << "/" <<B[1][1] << '\n';
std::cout << B[2][0] << "/" <<B[2][1] << '\n';
std::cout << B[4][0] << "/" <<B[4][1] << '\n';
return 0;
}
計算自然數冪和
1. 給你 \(n\)，\(k\)，求\(\sum_{i = 0}^{n}i^k\)，其中\(n\) 非常大。
解法一：直接遞推
利用差分的思想，先作差：
\((i + 1)^{k + 1} - i^{k + 1} = \binom{k + 1}{1}i^k + \binom{k + 1}{2} i^{k - 1} + ... + \binom{k + 1}{k}i + 1\)
再計算前綴和:
\((n + 1)^{k + 1} - 1 = \binom{k + 1}{1}\sum\limits_{i = 0}^{n}i^k + \binom{k + 1}{2}\sum\limits_{i = 0}^{n}i^{k-1} + ... + \binom{k + 1}{k}\sum\limits_{i = 0}^{n}i + n\)
整理一下就是：
\(\sum\limits_{i = 0}^{n}i^k = \frac{1}{k + 1} \left[(n + 1)^{k + 1} - \left(\binom{k + 1}{2}\sum\limits_{i = 0}^{n}i^{k-1} + ... + \binom{k + 1}{k}\sum\limits_{i = 0}^{n}i + n + 1\right)\right]\).
這樣做的復雜度是\(O(k^2)\)。
解法二：利用斯特林數。
前置知識（參考具體數學）
1.下降冪：\(x^{\underline{n}} = \frac{x!}{(x-n)!}\)
2.第一類斯特林數：\(\left[\begin{matrix} n \\ k \end{matrix}\right] = (n - 1)\left[\begin{matrix} n - 1 \\ k \end{matrix}\right] + \left[\begin{matrix} n - 1 \\ k - 1 \end{matrix}\right]\)
3.兩個恆等式：$\begin{aligned}
x^{\underline{n}} = \sum_{k = 0}^{n} \left[\begin{matrix} n \ k \end{matrix}\right] (-1)^{n - k} x^k \
\
\sum_{k = 0}^{n} \binom{k}{m} = \binom{n + 1}{m + 1}
\end{aligned}$
可以發現：
$\begin{aligned}
\binom{n}{k} = \frac{n^{\underline{k}}}{k!} = \frac{\sum\limits_{i = 0}^{k} \left[\begin{matrix} k \ i \end{matrix}\right] (-1)^{k - i} n^i}{k!} \
k! \binom{n}{k} = \sum_{i = 0}^{k} \left[\begin{matrix} k \ i \end{matrix}\right] (-1)^{k - i} n^i \
\
n^k = k! \binom{n}{k} - \sum_{i = 0}^{\color{red}{k - 1}} \left[\begin{matrix} k \ i \end{matrix}\right] (-1)^{k - i} n^i
\end{aligned}$
那么:
\(\begin{aligned} \sum_{i = 0}^{n} i^k &= \sum_{i = 0}^{n} \left( k! \binom{i}{k} - \sum_{j = 0}^{k - 1} \left[\begin{matrix} k \\ j \end{matrix}\right] (-1)^{k - j} i^j \right) \\ &= k! \sum_{i = 0}^{n} \binom{i}{k} - \sum_{i = 0}^{n} \sum_{j = 0}^{k - 1} \left[\begin{matrix} k \\ j \end{matrix}\right] (-1)^{k - j} i^j \\ &= k! \binom{n + 1}{k + 1} - \sum_{j = 0}^{k - 1} \left[\begin{matrix} k \\ j \end{matrix}\right] (-1)^{k - j} \color{red}{\sum_{i = 0}^{n} i^j} \\ &= \frac{(n + 1)^{\underline{k + 1}}}{k + 1} - \sum_{j = 0}^{k - 1} \left[\begin{matrix} k \\ j \end{matrix}\right] (-1)^{k - j} \color{red}{\sum_{i = 0}^{n} i^j} \end{aligned}\)