Von Staudt–Clausenの定理でベルヌーイ数の分母を求める

今回はベルヌーイ数に関する話題です。まずベルヌーイ数になじみが無い場合はこちらのyoutubeが超おすすめです。

[べき乗和の公式と「6」② - 明日話したくなる「数」のお話 #19 【ベルヌーイ数】](https://www.youtube.com/watch?v=wvtm1QQymHc)

### ベルヌーイ数を求める

まず上のビデオでも紹介されているベルヌーイ数の再帰的な定義は以下になります。

```math 
\binom{k+1}{k}B_k = k+1 -\sum_{i=0}^{k-1}\binom{k+1}{i}B_i
```
これをそのままプログラムにします$k=10$まで求めます。
```python 
import sympy as sp
from functools import lru_cache
@lru_cache(maxsize=None)
def bern(k):
  if k == 0: return sp.Rational(1,1)
  ret = k + 1 - sum(sp.binomial(k+1,i)*bern(i) for i in range(k))
  return ret / sp.binomial(k+1,k)

for k in range(10+1):
  print(f"{k}: {bern(k)}")
# 0: 1
# 1: 1/2
# 2: 1/6
# 3: 0
# 4: -1/30
# 5: 0
# 6 : 1/42
# 7: 0
# 8: -1/30
# 9: 0
# 10: 5/66
```
ベルヌーイ数はべき乗和の公式に使える等、話題には事欠かないようですが、今回はその分母に注目します。

まず次の例題を考えます。

### 【例題１】k=2024のときBkの分母を求めよ

上のプログラムで$bern(2024)$を求めようとすると全然終わりません。そこでVon Staudt–Clausenの定理の登場となります。[Von Staudt–Clausen theorem (Wikipedia) 英語](https://en.wikipedia.org/wiki/Von_Staudt%E2%80%93Clausen_theorem)によると、Von Staudt–Clausenの定理は

```math 
B_{2n} + \sum_{(p-1)|2n}\frac{1}{p} \in \mathbb{Z}
```

ですが何のことか分かりませんね（笑）。私もこのビデオを見てやっと分かりました。

[ベルヌーイ数とゼータ関数「6」③ - 明日話したくなる「数」のお話 #20](https://www.youtube.com/watch?v=tl_O47Xc6MQ)

結論だけを言うと **2nの約数に１を足して素数になるものの積が$B_{2n}$の分母になる** ということです。

```math
\begin{align}
たとえば&2n=10のとき \\
約数は&[1,2,5,10] \\
１を足すと&[2,3,6,11] \\
そのうち素数は&[2,3,11] \\
その積は&66 \\
となりB_{10}の&分母が求まります。
\end{align}
```

本当にそうなるか$2n \le 40$で見てみると、確かに同じ値になっていることが確認されました。

|2n |$B_{2n}$の分子 |$B_{2n}$の分母|素数の(2nの約数+1)|その積|
|:---:|:---:|:---:|:---:|:---:|
|2|1|6|[2, 3]|6|
|4|-1|30|[2, 3, 5]|30|
|6|1|42|[2, 3, 7]|42|
|8|-1|30|[2, 3, 5]|30|
|10|5|66|[2, 3, 11]|66|
|12|-691|2730|[2, 3, 5, 7, 13]|2730|
|14|7|6|[2, 3]|6|
|16|-3617|510|[2, 3, 5, 17]|510|
|18|43867|798|[2, 3, 7, 19]|798|
|20|-174611|330|[2, 3, 5, 11]|330|
|22|854513|138|[2, 3, 23]|138|
|24|-236364091|2730|[2, 3, 5, 7, 13]|2730|
|26|8553103|6|[2, 3]|6|
|28|-23749461029|870|[2, 3, 5, 29]|870|
|30|8615841276005|14322|[2, 3, 7, 11, 31]|14322|
|32|-7709321041217|510|[2, 3, 5, 17]|510|
|34|2577687858367|6|[2, 3]|6|
|36|-26315271553053477373|1919190|[2, 3, 5, 7, 13, 19, 37]|1919190|
|38|2929993913841559|6|[2, 3]|6|
|40|-261082718496449122051|13530|[2, 3, 5, 11, 41]|13530|

このアルゴリズムを使って例題２を解きます。

```python
import sympy as sp
import numpy as np
k = 2024
print(f"B({k}).denominator = {np.prod([d+1 for d in sp.divisors(k) if sp.isprime(d+1)])}")
# B(2024).denominator = 2923791510
```

### 【例題３】k<10**8のとき(Bkの分母)=131070となるkの個数を求めよ

ベルヌーイ数の分母は簡単に求めることができるようになったので、すべてのkについてチェックすれば良いのですが、さすがに$10^8$は厳しいのでちょっと工夫します。

1. $13170=2 \times 3 \times 5 \times 17 \times257$なので$B_k$の分母が$131070$ということは、$k$の約数にそれぞれ$1$を足したものが素数であれば$[2,3,5,17,257]$のどれかになり、それ以外を含まない
1. したがって、$k$の約数には$[1,2,4,16,256]$が含まれるので、$k$はそれらの最小公倍数$256$の倍数である

この1と2を満たせば$B_k$の分母は必ず$13170$になるので、この２点チェックすれば十分となります。

```python 
import sympy as sp
import math
def isBk(k, bk):
  for d in sp.divisors(k):
    if sp.isprime(d+1) and d+1 not in sp.factorint(bk).keys():
      return False
  return True

bk = 131070
N = 10**8
print(f"{bk} = prod{list(sp.factorint(bk).keys())}")
k0 = math.lcm(*([p-1 for p in sp.factorint(bk).keys()]))
print(f"k = ({k0} x i)")
print(f"Answer: {sum(isBk(k, bk) for k in range(k0, N+1, k0))}")
# 131070 = prod[2, 3, 5, 17, 257]
# k = (256 x i)
# Answer: 28575
```

（開発環境：Google Colab）

この考え方はProject Euler、Problem 545: Faulhaber's Formulasを解くのに役に立ちます