ベルヌーイ数 (Bernoulli number)

概要

ここでは、素数mod上で計算することを考える。

ベルヌーイ数を求める

ベルヌーイ数 \(B_n\) は、以下により再帰的に計算できる。

\(\displaystyle B_0 = 1\)

\(\displaystyle B_n = - \frac{1}{n+1} \sum_{k=0}^{n-1} {n+1 \choose k} B_k\)

この式を愚直に計算することで \(B_0\) から \(B_p\) までのベルヌーイ数は (逆数計算を除き) \(O(p^2)\) で計算できる。

\(1^p + 2^p + ... + n^p\) の値を求める

ベルヌーイ数を用いて、以下の式により \(\displaystyle \sum_{k=1}^n k^p\) を計算できる。

\(\displaystyle \sum_{k=1}^n k^p = \frac{1}{p+1} \sum_{j=0}^p (-1)^j {p+1 \choose j} B_j n^{p+1-j}\)

この式は (逆数計算を除き) \(O(p)\) で計算できる。

実装

# calculate bernoulli numbers B_0 ~ B_p on modulo mod
def bernoulli(p, mod):
    rev = [1]*(p+2)
    cur = 1
    for i in range(1, p+1):
        cur = cur * i % mod
        rev[i+1] = cur
    cur = pow(cur * (p+1), mod-2, mod)
    for i in range(p+1, -1, -1):
        rev[i] = (rev[i] * cur) % mod
        cur = cur * i % mod

    B = [0]*(p+1)
    B[0] = 1
    if p > 1:
        B[1] = (-rev[2]) % mod
    for i in range(2, p+1, 2):
        v = 0
        tmp = 1
        for j in range(i):
            v = (v + tmp * B[j]) % mod
            tmp = ((i+1-j) * rev[j+1]) * tmp % mod
        B[i] = (-rev[i+1] * v) % mod
    return B, rev

# calculate ∑_{k=1}^n k^p on modulo mod
def calc(n, p, mod):
    B, rev = bernoulli(p, mod)

    res = 0
    tmp = pow(n, p+1, mod)
    rev_n = pow(n, mod-2, mod)
    for i in range(0, p+1):
        res = (res + tmp * B[i]) % mod
        tmp = (-(p+1-i) * rev[i+1] % mod) * (rev_n * tmp % mod) % mod
    return (res * rev[p+1]) % mod


print(calc(10, 2, 10**9 + 7))
# => "385"
print(calc(100, 200, 10**9 + 7))
# => "710636539"

参考

戻る