先看一下差分序列和斯特林数。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}\)