베르누이 수

1. 개요

베르누이 수 Bernoulli Numbers
생성함수	$t/(e^t - 1) = \sum_{n\ge0} B_n t^n/n!$
처음 몇 개	$1,\ -\tfrac12,\ \tfrac16,\ 0,\ -\tfrac1{30},\ 0,\ \tfrac1{42},\ 0,\ -\tfrac1{30}$
홀수 번째	$n \ge 3$ 홀수면 $B_n = 0$
제타와의 관계	$\zeta(2n) = (-1)^{n+1}B_{2n}(2\pi)^{2n}/(2(2n)!)$
증가 속도	$\lvert B_{2n}\rvert \sim 4\sqrt{\pi n}\,(n/\pi e)^{2n}$
배정밀도 넘침	$B_{260}$ 부터 오버플로
쓰이는 곳	오일러-매클로린, 스털링 급수, tan·coth 전개, 거듭제곱 합

수치해석의 오차 전개를 열어 보면 꼭 이 숫자들이 계수로 앉아 있다.

베르누이 수(Bernoulli numbers) $B_n$ 은 생성함수

$$\frac{t}{e^t - 1} = \sum_{n=0}^{\infty} B_n \frac{t^n}{n!}$$

의 계수로 정의되는 유리수열이다. 야코프 베르누이가 유고집 『추측술』(1713)에서 거듭제곱 합 $1^p + 2^p + \cdots + n^p$ 의 일반 공식을 쓰다가 도입했고, 같은 수열을 세키 다카카즈가 비슷한 시기(1712년 유고)에 독립적으로 얻었다.

수치해석에서 베르누이 수를 「연구」할 일은 거의 없다. 대신 계수로 만날 일은 매우 많다. 사다리꼴 공식의 오차 전개, 오일러-매클로린 공식의 보정항, 감마 함수의 스털링 급수, 디감마 함수의 점근 전개, $\tan x$·$\coth x$ 의 테일러 계수가 전부 베르누이 수다. 이 문서는 그 공통 분모를 정리하고, 배정밀도에서 어떻게 계산하면 안전한가를 실측으로 보인다.

2. 정의와 값

생성함수에서 $t/(e^t-1) + t/2$ 가 우함수라는 사실을 쓰면 $B_1 = -1/2$ 을 빼고는 홀수 번째가 전부 0 이다. 짝수 번째는 부호가 번갈아 바뀐다.

$n$	$B_n$	소수
0	$1$	1
1	$-1/2$	−0.5
2	$1/6$	0.16667
4	$-1/30$	−0.03333
6	$1/42$	0.02381
8	$-1/30$	−0.03333
10	$5/66$	0.07576
12	$-691/2730$	−0.25311
14	$7/6$	1.16667
20	$-174611/330$	−529.12
30	$8615841276005/14322$	$6.016\times10^{8}$

처음 몇 개는 작아서 「얌전한 수열」처럼 보이지만 $B_{14}$ 에서 1 을 넘고 그 뒤로는 계승보다 빠르게 커진다. $B_{12}$ 의 분자 691 은 정수론에서 따로 유명한 수다¹.

부호 관례 주의. 위 생성함수는 $B_1 = -1/2$ 를 준다. $t e^t/(e^t-1)$ 을 쓰는 관례에서는 $B_1 = +1/2$ 이고 나머지는 같다. 거듭제곱 합 공식을 $\sum_{k=1}^{n}$ 으로 쓰느냐 $\sum_{k=0}^{n-1}$ 로 쓰느냐에 따라 편한 쪽이 갈린다. DLMF 와 대부분의 수치 라이브러리는 $-1/2$ 쪽이다. 다른 교재의 공식을 옮겨 올 때 $B_1$ 이 들어간 항의 부호부터 확인하는 것이 국룰이다.

점화식. 생성함수 양변에 $e^t - 1$ 을 곱하고 계수를 맞추면

$$\sum_{k=0}^{n} \binom{n+1}{k} B_k = 0 \qquad (n \ge 1)$$

이 나온다. $B_0 = 1$ 에서 시작해 차례로 풀면 된다.

폰 슈타우트-클라우센 정리. $B_{2n}$ 의 분모는 「$p-1$ 이 $2n$ 을 나누는 소수 $p$」 의 곱이고, 정확히는 $B_{2n} + \sum_{(p-1)\mid 2n} 1/p$ 가 정수다. $B_{12}$ 에서 $p = 2, 3, 5, 7, 13$ 이므로 분모가 $2\cdot3\cdot5\cdot7\cdot13 = 2730$ 이고, $-691/2730 + 1/2 + 1/3 + 1/5 + 1/7 + 1/13 = 1$ 로 정수가 된다. 정확 산술 구현의 검산용으로 쓸 만하다.

3. 어디서 튀어나오나

거듭제곱 합 (파울하버 공식).

$$\sum_{k=1}^{n} k^p = \frac{1}{p+1}\sum_{j=0}^{p} \binom{p+1}{j} (-1)^j B_j\, n^{p+1-j}$$

($B_1 = -1/2$ 관례). $p = 10$, $n = 1000$ 에 대해 정확 유리수 산술로 직접 합과 비교하면 32자리 정수 91409924241424243424241924242500 이 정확히 일치한다. 베르누이가 『추측술』에서 「1000 까지의 10제곱 합을 반 쿼터 시간에 구했다」고 자랑한 바로 그 계산이다.²

사다리꼴 공식의 오차. 오일러-매클로린 공식은 사다리꼴 근사 $T_h$ 와 적분의 차이를

$$T_h - \int_a^b f\,dx = \sum_{k=1}^{m} \frac{B_{2k}}{(2k)!}h^{2k}\bigl[f^{(2k-1)}(b) - f^{(2k-1)}(a)\bigr] + R_m$$

으로 준다. 오차가 $h^2$ 의 거듭제곱으로만 전개된다는 이 사실이 롬베르그 적분과 리처드슨 외삽의 근거이고, 주기함수에서는 괄호 안이 전부 0 이 되어 사다리꼴이 지수 수렴한다. 나머지항 $R_m$ 에는 주기 베르누이 함수 $\tilde B_{2m}(x) = B_{2m}(x - \lfloor x\rfloor)$ 가 들어간다. 여기서 $B_n(x)$ 는 베르누이 다항식으로, $B_n(0) = B_n$ 이다.

스털링 급수와 디감마 전개. $\ln\Gamma(x)$ 의 점근 전개 계수는 $B_{2k}/(2k(2k-1)x^{2k-1})$ 이고, 디감마 함수는 $\psi(x) \sim \ln x - 1/(2x) - \sum B_{2k}/(2k x^{2k})$ 이다. 둘 다 점근 급수라 항을 무한정 늘릴 수 없는데, 그 이유가 바로 아래의 증가 속도다.

삼각·쌍곡 함수. $x\coth x = \sum_{n} B_{2n}(2x)^{2n}/(2n)!$ 이고, $\tan x$ 의 계수도 $B_{2n}$ 으로 쓰인다. 반대로 $\tan x$ 의 테일러 계수에 $(2n-1)!$ 을 곱한 정수 $T_n = 1, 2, 16, 272, 7936, \dots$ 을 탄젠트 수라 부르고, $B_{2n} = (-1)^{n-1}\,2n\,T_n / (4^n(4^n-1))$ 이 성립한다. 아래 계산법에서 이 관계를 쓴다.

4. 증가 속도와 제타 함수

리만 제타 함수의 짝수점 값은 오일러가 처음 구했다.

$$\zeta(2n) = \frac{(-1)^{n+1} B_{2n} (2\pi)^{2n}}{2\,(2n)!}$$

$n = 1$ 이면 $\zeta(2) = \pi^2/6$ (바젤 문제), $n = 2$ 면 $\pi^4/90$ 이다. $n$ 이 커지면 $\zeta(2n) \to 1$ 이므로 거꾸로

$$\lvert B_{2n}\rvert = \frac{2\,(2n)!}{(2\pi)^{2n}}\zeta(2n) \approx \frac{2\,(2n)!}{(2\pi)^{2n}} \sim 4\sqrt{\pi n}\Bigl(\frac{n}{\pi e}\Bigr)^{2n}$$

이다. 마지막 근사의 상대오차는 $2n = 10, 20, 40$ 에서 각각 0.93%, 0.42%, 0.21% 로 줄어든다. 요점은 $\lvert B_{2n}\rvert$ 가 대략 $(2n)!/(2\pi)^{2n}$ 으로 자란다는 것이다. 그래서 오일러-매클로린 보정항의 크기는 $h$ 와 무관하게 $(h/2\pi)^{2k}(2k)!\,\lvert f^{(2k-1)}\rvert$ 꼴이 되고, $k$ 를 키우면 결국 계승이 이겨서 발산한다. 「보정항을 몇 개까지 쓰는 게 최선인가」 가 점근 급수의 최적 절단 문제로 넘어가는 지점이 여기다.

배정밀도에서는 $\lvert B_{258}\rvert \approx 1.34\times10^{306}$ 이 마지막으로 표현되고 $B_{260}$ 부터 오버플로한다. 그 너머가 필요하면 $\ln\lvert B_{2n}\rvert$ 를 위 식(로그 감마 + $\ln\zeta$)으로 직접 계산한다.

5. 계산법 — 무엇이 안정한가

여러 알고리즘이 있고, 정확 유리수 산술에서는 모두 맞는 답을 준다. 부동소수점에서는 얘기가 다르다. 각 방법을 배정밀도로 돌려 정확값과 비교한 실측이다.

방법	$B_{10}$ 상대오차	$B_{20}$	$B_{40}$	$B_{200}$
점화식 $\sum\binom{n+1}{k}B_k = 0$	$9\times10^{-15}$	$4\times10^{-16}$	$1\times10^{-15}$	$2\times10^{-15}$
아키야마-다니가와 알고리즘	$6\times10^{-9}$	12	$10^{19}$	—

표준 점화식은 의외로 멀쩡하다. $n = 200$ 까지 상대오차가 $10^{-14}$ 대에 머문다. 「점화식은 불안정하다」는 말이 돌지만 적어도 배정밀도로 오버플로 직전까지는 실측상 그렇지 않다. 다만 $B_n$ 하나에 $O(n)$ 번, 전체 $O(n^2)$ 번의 연산이 든다.

아키야마-다니가와 알고리즘은 $a_{0,m} = 1/(m+1)$ 에서 시작해 $a_{j,m} = (j+1)(a_{j,m} - a_{j,m+1})$ 꼴의 차분을 반복하는, 짧고 우아한 표 계산이다. 정수(유리수) 산술용으로는 깔끔하지만 부동소수점에서는 차분이 반올림 오차를 단계마다 증폭시켜, $B_{20}$ 에서 이미 부호도 크기도 틀린 값($5986.6$, 정답 $-529.1$)이 나온다. 블로그·코딩 테스트 풀이에 자주 보이는 구현이라 그대로 float 로 옮겼다가 당하는 일이 많다³.

브렌트-하비(2011)의 탄젠트 수 알고리즘은 정수만으로 $T_1, \dots, T_N$ 을 $O(N^2)$ 번의 덧셈·곱셈으로 만든 뒤 위 관계식으로 $B_{2n}$ 을 얻는다. 모든 중간값이 양의 정수라 소거가 없고, 임의 정밀도에서 $B_{2n}$ 을 대량으로 뽑을 때 쓴다. 매우 큰 $n$ 하나만 필요하면 $\zeta(2n)$ 을 오일러 곱으로 고정밀 계산하고 분모는 폰 슈타우트-클라우센으로 정확히 정하는 「제타 + 분모」 방식이 가장 빠르다(분자의 자릿수만큼만 정밀도를 확보하면 반올림으로 정확한 정수가 나온다).

실무 결론은 단순하다. 수치 코드에서 쓰는 베르누이 수는 $B_{30}$ 정도까지가 전부이고, 그 값은 상수 테이블로 박는다. 스털링 급수나 디감마 전개도 $B_{14}$ 근처에서 끊는 것이 보통이다. 계산해서 쓸 일이 생기면 점화식이나 탄젠트 수, 아니면 정확 산술이다. 아키야마-다니가와를 float 로 돌리지만 않으면 된다.

6. 관련 문서

• 오일러-매클로린 공식 · 점근 급수 · 롬베르그 적분 · 수치적분
• 감마 함수 · 디감마 함수 · 리만 제타 함수
• 부동소수점 연산 · 리처드슨 외삽 · 테일러 급수
• 파울하버 공식 · 탄젠트 수 · 쿠머 합동식

7. Footnotes

1. 라마누잔 타우 함수의 합동식 $\tau(n) \equiv \sigma_{11}(n) \pmod{691}$ 에 나오는 그 691 이다. 수치해석 하다가 갑자기 모듈러 형식으로 빠지면 현생이 위험하니 여기까지만. ↩

2. 원문은 「intra semi-quadrantem horae」, 7분 30초 안에 했다는 뜻이다. 300년 전 수학자도 계산 시간을 벤치마크로 자랑했다. 지금 우리는 같은 걸 sum(k**10 for k in range(1001)) 로 1밀리초에 한다. ↩

3. 「정수 수열을 정수 알고리즘으로 구한다」는 원래 맥락을 떼어 내고 double 로 옮긴 순간 알고리즘은 같아도 다른 물건이 된다. 안정성은 알고리즘이 아니라 알고리즘 + 산술의 성질이다. ↩