先看一下差分序列和斯特林数。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] << '
';
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] << '
';
std::cout << B[2][0] << "/" <<B[2][1] << '
';
std::cout << B[4][0] << "/" <<B[4][1] << '
';
return 0;
}
计算自然数幂和
1. 给你 (n)，(k)，求(sum_{i = 0}^{n}i^k)，其中(n) 非常大。
解法一：直接递推
利用差分的思想，先作差：
((i + 1)^{k + 1} - i^{k + 1} = inom{k + 1}{1}i^k + inom{k + 1}{2} i^{k - 1} + ... + inom{k + 1}{k}i + 1)
再计算前缀和:
((n + 1)^{k + 1} - 1 = inom{k + 1}{1}sumlimits_{i = 0}^{n}i^k + inom{k + 1}{2}sumlimits_{i = 0}^{n}i^{k-1} + ... + inom{k + 1}{k}sumlimits_{i = 0}^{n}i + n)
整理一下就是：
(sumlimits_{i = 0}^{n}i^k = frac{1}{k + 1} left[(n + 1)^{k + 1} - left(inom{k + 1}{2}sumlimits_{i = 0}^{n}i^{k-1} + ... + inom{k + 1}{k}sumlimits_{i = 0}^{n}i + n + 1 ight) ight]).
这样做的复杂度是(O(k^2))。
解法二：利用斯特林数。
前置知识（参考具体数学）
1.下降幂：(x^{underline{n}} = frac{x!}{(x-n)!})
2.第一类斯特林数：(left[egin{matrix} n \ k end{matrix} ight] = (n - 1)left[egin{matrix} n - 1 \ k end{matrix} ight] + left[egin{matrix} n - 1 \ k - 1 end{matrix} ight])
3.两个恒等式：$egin{aligned}
x^{underline{n}} = sum_{k = 0}^{n} left[egin{matrix} n k end{matrix}
ight] (-1)^{n - k} x^k
sum_{k = 0}^{n} inom{k}{m} = inom{n + 1}{m + 1}
end{aligned}$
可以发现：
$egin{aligned}
inom{n}{k} = frac{n^{underline{k}}}{k!} = frac{sumlimits_{i = 0}^{k} left[egin{matrix} k i end{matrix}
ight] (-1)^{k - i} n^i}{k!}
k! inom{n}{k} = sum_{i = 0}^{k} left[egin{matrix} k i end{matrix}
ight] (-1)^{k - i} n^i
n^k = k! inom{n}{k} - sum_{i = 0}^{color{red}{k - 1}} left[egin{matrix} k i end{matrix}
ight] (-1)^{k - i} n^i
end{aligned}$
那么:
(egin{aligned} sum_{i = 0}^{n} i^k &= sum_{i = 0}^{n} left( k! inom{i}{k} - sum_{j = 0}^{k - 1} left[egin{matrix} k \ j end{matrix} ight] (-1)^{k - j} i^j ight) \ &= k! sum_{i = 0}^{n} inom{i}{k} - sum_{i = 0}^{n} sum_{j = 0}^{k - 1} left[egin{matrix} k \ j end{matrix} ight] (-1)^{k - j} i^j \ &= k! inom{n + 1}{k + 1} - sum_{j = 0}^{k - 1} left[egin{matrix} k \ j end{matrix} ight] (-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[egin{matrix} k \ j end{matrix} ight] (-1)^{k - j} color{red}{sum_{i = 0}^{n} i^j} end{aligned})