自然数幂和与伯努利数

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

红色的项是之前已经求出的,，所以我们只需要(O(k^2)) 预处理出第一类斯特林数即可在(O(k^2))递推计算。

解法三：利用拉格朗日插值。。。(不写了...)

相关题目：

http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=1864

http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1228

http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1258

http://codeforces.com/problemset/problem/622/F