파울하버 공식

1. 개요

파울하버 공식 Faulhaber's Formula
대상	거듭제곱 합 $S_p(n) = 1^p + 2^p + \cdots + n^p$
형태	$n$ 에 대한 $p+1$ 차 다항식, 계수는 베르누이 수
최고차 두 항	$\dfrac{n^{p+1}}{p+1} + \dfrac{n^p}{2}$
홀수 $p$	$u = n(n+1)/2$ 의 다항식 (파울하버 1631, 증명 야코비 1834)
통일 공식	야코프 베르누이 『추측술』 (1713)
대칭	$S_p(-1-n) = (-1)^{p+1}S_p(n)$

$1 + 2 + \cdots + n = n(n+1)/2$ 를 아는 사람은 많다. 그다음 칸, 그 다음 칸을 400년 전에 17제곱까지 손으로 채운 사람이 있었다.

파울하버 공식(Faulhaber’s formula)은 거듭제곱 합

$$S_p(n) = \sum_{k=1}^{n} k^p$$

을 $n$ 의 다항식으로 쓰는 공식이다. 오늘날 표준형은 베르누이 수 $B_j$ 를 써서

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

로 쓴다. 여기서 $B_j^{+}$ 는 $B_1 = +\tfrac12$ 관례의 베르누이 수다. 처음 몇 개는

$$S_1 = \tfrac{n^2}{2} + \tfrac{n}{2},\quad S_2 = \tfrac{n^3}{3} + \tfrac{n^2}{2} + \tfrac{n}{6},\quad S_3 = \tfrac{n^4}{4} + \tfrac{n^3}{2} + \tfrac{n^2}{4}.$$

이름은 17세기 울름의 계산가 요한 파울하버에게서 왔다. 그는 1631년 『Academia Algebrae』에서 17제곱 합까지의 공식을 내놓았지만 베르누이 수라는 일반 구조는 보지 못했고, 그 통일은 야코프 베르누이의 『추측술』(1713, 사후 출간) 몫이었다. 그래서 위 공식을 「베르누이 공식」이라 부르는 문헌도 많다. 한편 파울하버가 정말로 알아챈 것은 다른 구조였고, 그게 이 문서 둘째 절의 주제다.

베르누이 수 문서가 위 공식의 정확 산술 검산(10제곱 합, 32자리)을 이미 다뤘으니, 이 문서는 $B_1$ 부호 관례 문제, 파울하버의 $u$-다항식, 공식이 가진 구조(대칭·인수·선형 계수), 다른 유도 경로, 그리고 이 다항식을 부동소수점으로 평가하면 무슨 일이 생기는지를 다룬다.

2. $B_1$ 의 부호 — 끝이 $n$ 이냐 $n-1$ 이냐

베르누이 수는 $B_1$ 만 관례가 갈린다. 생성함수 $t/(e^t - 1)$ 에서 나오는 쪽은 $B_1 = -\tfrac12$ 이고, $t/(1 - e^{-t})$ 에서 나오는 쪽은 $B_1 = +\tfrac12$ 다. 나머지 $B_j$ 는 같다. 이 차이가 거듭제곱 합에서는 합의 끝점으로 나타난다.

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

즉 $B_1 = -\tfrac12$ 관례로 위 식을 그대로 쓰면 $n$ 까지가 아니라 $n-1$ 까지의 합이 나온다. 두 결과의 차이가 정확히 $n^p$, 곧 $j = 1$ 항의 부호 차이 $\binom{p+1}{1}\cdot 1 \cdot n^p/(p+1)$ 다. 베르누이 수 문서는 $B_1 = -\tfrac12$ 관례에 $(-1)^j$ 를 곱해 같은 결과를 냈다($(-1)^j B_j^{-} = B_j^{+}$, 홀수 $j \ge 3$ 에서는 $B_j = 0$ 이라 부호가 의미 없다). 라이브러리마다 관례가 다르니(SymPy 는 버전에 따라 바뀌기까지 했다) 공식을 코드로 옮길 때 $n = 1$ 하나만 넣어 봐도 바로 잡힌다.¹

3. 파울하버가 본 것 — 홀수 제곱은 $u = n(n+1)/2$ 의 다항식

$u = n(n+1)/2 = S_1$ 으로 두면 놀랍게도

$$S_3 = u^2,\qquad S_5 = \frac{4u^3 - u^2}{3},\qquad S_7 = \frac{6u^4 - 4u^3 + u^2}{3}$$

이다. $n = 1, \dots, 19$ 에서 정확 산술로 모두 확인된다. 일반적으로 홀수 $p = 2m-1$ 에 대해 $S_p$ 는 $u$ 의 $m$ 차 다항식이고, $p \ge 3$ 이면 $u^2$ 로 나누어떨어진다. 짝수 제곱은 $S_2 = \tfrac{2n+1}{3}u$ 로 시작해 $S_{2m}/S_2$ 가 $u$ 의 다항식이 된다. 예를 들어

$$S_4 = \frac{2n+1}{3}\cdot\frac{u(6u - 1)}{5}.$$

파울하버는 이 $u$-형태로 17제곱까지 계산했고, 계수를 구하는 체계적 방법도 갖고 있었던 것으로 보인다. 그러나 「모든 홀수 $p$ 에서 성립한다」는 증명은 1834년 야코비가 처음 출판했다. 커누스(1993)는 파울하버의 원저를 다시 읽고 그가 남긴 계수표의 구조를 현대 언어로 재구성했다.

왜 $n(n+1)$ 인가는 다음 절의 대칭에서 바로 나온다.

4. 구조 — 대칭, 인수, 선형 계수

반사 대칭. 다항식 $S_p$ 를 음의 정수로 연장하면

$$S_p(-1-n) = (-1)^{p+1}S_p(n)$$

이 성립한다($p = 2, 3, 5, 8$ 에서 $n = 0, \dots, 9$ 로 확인). $p$ 가 홀수면 $S_p$ 는 $n = -\tfrac12$ 를 축으로 대칭이고, $-\tfrac12$ 에 대칭인 다항식은 $(n + \tfrac12)^2$ 의 다항식, 곧 $n(n+1) = (n+\tfrac12)^2 - \tfrac14$ 의 다항식이다. 파울하버의 $u$-형태가 여기서 나온다. 짝수 $p$ 면 반대칭이라 $(2n+1)$ 인수를 하나 품는다. 위 $S_2$, $S_4$ 의 $(2n+1)$ 이 그것이다.

인수. $S_p(0) = 0$ 이고 대칭에 의해 $S_p(-1) = 0$ 이므로 $n(n+1)$ 로 나누어떨어진다. 홀수 $p \ge 3$ 은 $n^2(n+1)^2$ 까지 나누어떨어진다.

계수에 대한 사실. 정확 산술로 펼친 $S_9$ 는

$$S_9 = \tfrac{1}{10}n^{10} + \tfrac12 n^9 + \tfrac34 n^8 - \tfrac{7}{10}n^6 + \tfrac12 n^4 - \tfrac{3}{20}n^2$$

다. 최고차 $\tfrac{1}{p+1}n^{p+1}$ 은 적분 $\int_0^n x^p\,dx$ 이고, 다음 항 $\tfrac12 n^p$ 는 사다리꼴의 끝점 보정이다. 그리고 $n$ 의 1차 계수가 정확히 $B_p$ 다 — $S_2$ 의 $\tfrac16 = B_2$, $S_4$ 의 $-\tfrac{1}{30} = B_4$, 홀수 $p$ 에서는 0. 다시 말해 $S_p'(0) = B_p$ 로 베르누이 수를 정의할 수도 있다. 이 「적분 + 끝점 보정 + 베르누이 수 보정」 구조는 오일러-매클로린 공식을 $f(x) = x^p$ 에 적용한 것과 같고, 다항식이라 보정항이 유한 개에서 끊겨 정확해진다.

5. 유도 경로들

• 망원합(파스칼 1654). $(k+1)^{p+1} - k^{p+1}$ 을 이항 전개해 $k = 1, \dots, n$ 으로 더하면 좌변이 망원되어 $S_p$ 를 $S_0, \dots, S_{p-1}$ 로 쓰는 점화식이 나온다. 손계산에는 이게 제일 빠르다. 에드워즈(1986)는 이 점화식을 하삼각 행렬로 쓰고 그 역행렬에서 베르누이 수가 튀어나오는 것을 보였다.
• 스털링 수. $k^p = \sum_j S(p, j)\,j!\binom{k}{j}$ ($S(p,j)$ 는 제2종 스털링 수)로 바꾸면 하키 스틱 항등식으로

$$S_p(n) = \sum_{j=0}^{p} S(p, j)\,j!\binom{n+1}{j+1}$$

이다($p = 3, 6, 10$, $n \le 14$ 로 확인). 계수가 전부 양의 정수라 아래의 부동소수점 문제가 없다. 이산 미적분에서 $\binom{k}{j}$ 가 거듭제곱 역할을 하는 「하강 계승」 관점이다.

• 생성함수. $\sum_{p} S_p(n)\,t^p/p! = \sum_{k=1}^n e^{kt} = e^t(e^{nt} - 1)/(e^t - 1)$ 에 $t/(e^t-1)$ 의 베르누이 전개를 곱하면 표준형이 바로 나온다.
• 보간. $S_p$ 가 $p+1$ 차 다항식이라는 것만 알면 $n = 0, 1, \dots, p+1$ 의 값으로 결정된다. 그 값들로 라그랑주 보간을 하면 매우 큰 $n$ 에 대한 $S_p(n) \bmod q$ 를 $O(p)$ 연산에 구할 수 있어 정수론·경시대회 쪽에서 표준 기법이다.

6. 수치 계산 — 다항식으로 평가하면 위험하다

$S_p(n)$ 을 계산하는 방법은 (1) 직접 합, (2) 펼친 다항식을 호너 방법으로 평가, (3) 정수 산술 셋이다. 배정밀도에서 (2) 를 쓰면 언제 깨지나 실측해 보자.

$p$	$n$	다항식의 최대 계수	다항식 평가 상대오차
10	1000	1	$2.0\times10^{-16}$
30	1000	$7.8\times10^{9}$	$2.2\times10^{-16}$
40	3	$2.5\times10^{17}$	$4.8\times10^{-13}$
40	2	$2.5\times10^{17}$	$1.1\times10^{-8}$
40	1	$2.5\times10^{17}$	5.0 (답 1 대신 6)
60	2	$2.8\times10^{35}$	2025

$n$ 이 크면 최고차 항이 지배해 문제가 없다. 그러나 $p$ 가 크고 $n$ 이 작으면 계수들(베르누이 수가 들어간 거대한 유리수)이 부호를 바꿔 가며 서로 상쇄해야 작은 답이 나오는데, 배정밀도는 그 상쇄를 감당하지 못한다. $S_{40}(1) = 1$ 을 6 으로, $S_{60}(2) = 1 + 2^{60}$ 을 2000 배 넘게 틀린다. 베르누이 수가 $|B_{2k}| \sim 4\sqrt{\pi k}\,(k/\pi e)^{2k}$ 로 계승보다 빨리 커진다는 것(베르누이 수의 증가 속도 절)이 그대로 계수 크기로 나타난 것이다.

실무 결론은 간단하다. 부동소수점이 필요하면 직접 합(같은 표의 모든 경우에서 상대오차 $4\times10^{-16}$ 이하)이나 스털링 수 형태를 쓰고, 정확한 값이 필요하면 정수 산술을 쓴다. 펼친 파울하버 다항식은 기호 계산과 증명용이다.²

7. 시뮬레이션 코드에서 만나는 곳

• 연산 횟수 세기. 가우스 소거법의 소거 단계는 $k$ 번째 피벗에서 나눗셈 $(n-k)$ 번, 곱셈·덧셈 $2(n-k)^2$ 번을 하므로 총 연산이 $\sum_{k}\bigl[(n-k) + 2(n-k)^2\bigr] = 2S_2(n-1) + S_1(n-1) = \tfrac23 n^3 - \tfrac12 n^2 - \tfrac16 n$ 이다. $n = 100$ 에서 정확히 661 650 이고, 흔히 쓰는 근사 $\tfrac23 n^3 = 666\,667$ 과 0.8% 차이다. 「LU 는 $\tfrac23 n^3$」의 출처가 $S_2$ 의 최고차항이다.³
• 등간격 격자 위의 모멘트. 격자점 $x_k = kh$ 에서 $\sum x_k^p$ 를 닫힌꼴로 쓰면 수치적분 공식(뉴턴-코츠)의 다항식 정확도를 손으로 검산하거나, 균일 격자의 관성 모멘트·분산을 루프 없이 계산할 수 있다. 사다리꼴 규칙이 $x^p$ 를 얼마나 틀리는지가 정확히 위 구조절의 「끝점 보정 + 베르누이 보정」 항이다.
• 차분 연산자 시험. $S_p$ 는 후진 차분 $\nabla S_p(n) = S_p(n) - S_p(n-1) = n^p$ 를 만족하는 다항식이므로, 이산 적분(누적합) 루틴의 정확도 시험용 정답으로 쓰기 좋다.

8. 관련 문서

• 베르누이 수 · 오일러-매클로린 공식 · 리만 제타 함수 · 점근 급수
• 수치적분 · 롬베르그 적분 · 가우스 소거법 · 부동소수점 연산
• 스털링 수 · 탄젠트 수 · 보간과 근사 · 호너 방법

9. Footnotes

1. $B_1$ 의 부호 전쟁은 아직도 끝나지 않았다. 「$+\tfrac12$ 가 자연스럽다」 파와 「생성함수가 $t/(e^t-1)$ 인데 무슨 소리냐」 파가 위키 편집 분쟁까지 벌인다. 결론: 둘 다 맞고, 코드에서는 $n = 1$ 을 넣어 봐라. ↩

2. 「닫힌꼴이 있으니 루프보다 빠르고 정확하겠지」라는 직관이 정면으로 배신당하는 사례. 닫힌꼴은 정확한 산술에서 정확한 것이지, 부동소수점에서 정확하다는 보장은 아무도 안 했다. ↩

3. 연산 횟수를 $n^3$ 의 최고차항만 남기고 부르는 건 국룰이지만, $n = 10$ 짜리 작은 블록을 수백만 개 푸는 코드에서는 저차항 때문에 $\tfrac23 n^3$ 근사가 8% 넘게 어긋난다. 작은 행렬의 세계에서는 저차항도 현생이다. ↩