ベルヌーイ数を線形時間で求める

# はじめに

$n$ 番目のベルヌーイ数 $B_n$ を $O(n)$ で求める方法です。




$\displaystyle B_0=1, B_1=-\frac12$ で、 $3$ 以上の奇数に対し $B_n=0$ です。これらを除外し、以下では正の偶数に対して $B_n$ を求めることを考えます。




# 記法

以下では $\displaystyle F(x)=\sum_n \frac{a_n}{n!}x^n$ に対し $\displaystyle \left[\frac{x^n}{n!}\right]F(x)=a_n$ とします。


# 本題

ベルヌーイ数の定義は $\displaystyle B_n=\left[\frac{x^n}{n!} \right] \frac x{e^x-1} $ です。この定義から式変形をしていきます。


```math
\begin{align}
&\phantom{=} B_n\\
&= \left[\frac{x^n}{n!} \right] \frac x{e^x-1} \\
&=\left[\frac{x^n}{n!} \right] \frac{\log\left(1+\left(e^x-1\right)\right)}{e^x-1}\\
&=\left[\frac{x^n}{n!} \right] \sum_{k=0}^\infty \frac{(1-e^x)^k}{k+1}\\
&=\left[\frac{x^n}{n!} \right] \sum_{k=0}^n \frac{(1-e^x)^k}{k+1}
\end{align}
```

最後の式変形は $e^x-1$ の定数項が $0$ であるため $k > n$ なら $\displaystyle\left[ \frac{x^n}{n!} \right] (1-e^x)^k=0$ となるからです。

ここで、 $\displaystyle F(z)=\sum_{k=0}^n\frac{(1-z)^k}{k+1}$ とすれば答えは $\displaystyle \left[\frac{x^n}{n!} \right] F(e^x)$ となります。さらに、有理数列 $c$ を $\displaystyle F(z)=\sum_{k=0}^n c_k z^k$ とします。 $c$ が求まれば、

```math
\begin{align}
&\phantom{=} B_n\\
&= \left[\frac{x^n}{n!} \right] F(e^x) \\
&=\left[\frac{x^n}{n!} \right] \sum_{k=0}^n c_k e^{xk}\\
&=\sum_{k=0}^n \left[\frac{x^n}{n!} \right]\sum_{i=0}^\infty c_kk^i\frac{x^i}{i!}\\
&=\sum_{k=0}^nc_kk^n
\end{align}
```

からベルヌーイ数を求められます。各 $k$ に対する $k^n$ の列挙は線形篩を用いて $O(n)$ で可能です。詳しくは [こちらの記事](https://37zigen.com/linear-sieve/) をご覧ください。

これ以降は各 $c$ の値を求めることを考えます。

まず $c_0$ ですが、こちらは $\displaystyle c_0=F(0)=\sum_{k=0}^n\frac1{k+1}$ から求めることができます。 $1,2,\ldots,n+1$ の逆元を事前に求めておくことでこの値は $O(n)$ で求められます。

次に、 $F$ を微分して微分方程式を求めます。

```math
\begin{align}
&\phantom{=} F'(z)\\
&=-\sum_{k=0}^n \frac k{k+1}(1-z)^{k-1}\\
&=-\frac1{1-z}\sum_{k=0}^n \left(1-\frac1{k+1} \right)(1-z)^k\\
&=\frac1{z-1}\left(\sum_{k=0}^n(1-z)^k-\sum_{k=0}^n\frac{(1-z)^k}{1+k}\right) \\
&=\frac1{z-1}\left(\frac{1-(1-z)^{n+1}}{z}-F(z) \right)
\end{align}
```

上の式変形から $\displaystyle (z-1)F'(z)+F(z)=\frac{1-(1-z)^{n+1}}{z}=\sum_{k=0}^n \binom{n+1}{k+1}(-1)^kz^k$ が得られます。

この式の $z^i$ の係数を比較することを考えます。

$\displaystyle (z-1)F'(z)=\sum_{k=0}^n kc_kz^k-\sum_{k=1}^{n-1} (k+1)c_{k+1}z^k$ の $z^i$ の係数は $ic_i-(i+1)c_{i+1}$ 、 $F(z)$ の $z^i$ の係数は $c_i$ 、 $\displaystyle \sum_{k=0}^n \binom{n+1}{k+1}(-1)^kz^k$ の $z_i$ の係数は $\displaystyle (-1)^i\binom{n+1}{i+1}$ です。したがって、 $\displaystyle ic_i-(i+1)c_{i+1}+c_i=(-1)^i\binom{n+1}{i+1}$　、つまり $$c_{i+1}=c_i-\frac{(-1)^i}{i+1}\binom{n+1}{i+1} $$ という関係式が得られます。 $c_0$ は上で求めたので、この漸化式に従って $c_1,c_2,\ldots,c_n$ を求めることができます。逆元と二項係数は $O(n)$ で列挙できるので、この漸化式から全ての $c$ の値を $O(n)$ で求めることができます。

あとは上で確かめた関係式 $\displaystyle B_n=\sum_{k=0}^nc_kk^n$ を使うことでベルヌーイ数を $O(n)$ で求めることができました。

# 問題例 (ネタバレ含む)

[Exponent of PI (Xmas Contest 2023)](https://atcoder.jp/contests/xmascon23/tasks/xmascon23_e)

この問題の答えは $1\le a,b,c \le 250$ に対し $\displaystyle \frac{\zeta(2ab)\zeta(2bc)\zeta(2ac)}{\zeta(2abc)^2}$ となります。正の偶数 $n$ に対し $\displaystyle \zeta(n)=\frac{(-1)^{n/2+1}2^{n-1} B_n \pi^n}{n!}$ という関係式が成り立つので、ベルヌーイ数を求めることができればこの問題を解くことができます。 $2abc\le 2\times 250^3=	3.125\times10^7$ なので、ベルヌーイ数を線形時間で求めれば TL に間に合います。

[実装例(C++,1463ms)](https://atcoder.jp/contests/xmascon23/submissions/48926931)

# おわりに

今回のような $e^x-1$ のかたまりを作り、 $e^x-1$ の定数項が $0$ であることを利用して高次の項を落とすという問題は様々な問題に応用することができます。 [こちらのブログ](https://www.luogu.com.cn/blog/lmpp/post-zai-tan-binomial-sum) では今回の方法で $\displaystyle \sum_{i=0}^n \binom ni i^k$ や $\displaystyle \sum_{i=0}^n F_i i^k$ を高速に求める方法が載っています。



以上がベルヌーイ数を線形で求める方法とその応用例です。読んでいただきありがとうございました。