## 1. Faulhaber's Formula 概述

Faulhaber's formula 用于计算 **前 n 个正整数的 p 次幂之和**：
$$S_p(n) = \sum_{k=1}^n k^p = 1^p + 2^p + 3^p + \cdots + n^p$$

公式的一般形式为：
$$S_p(n) = \frac{1}{p+1} \sum_{j=0}^p \binom{p+1}{j} B_j \cdot n^{p+1-j}$$

其中 $B_j$ 是伯努利数（Bernoulli numbers）。

### 特例

**p = 0**: $S_0(n) = n$

**p = 1**: $S_1(n) = \frac{n(n+1)}{2}$（三角形数）

**p = 2**: $S_2(n) = \frac{n(n+1)(2n+1)}{6}$（平方金字塔数）

**p = 3**: $S_3(n) = \left[\frac{n(n+1)}{2}\right]^2$（有趣的事实：等于 $S_1(n)^2$）

## 2. 历史背景

- **Johann Faulhaber**（1580-1635）：德国数学家，计算了 $p=1$ 到 $p=17$ 的公式，但未发现通式
- **Jacob Bernoulli**（1654-1705）：独立研究并发现了通用公式，引入了伯努利数
- 公式因此以 Faulhaber 命名，但现代形式主要归功于 Bernoulli

## 3. 伯努利数（Bernoulli Numbers）详解

### 3.1 定义

伯努利数由生成函数定义：
$$\frac{x}{e^x - 1} = \sum_{k=0}^\infty B_k \frac{x^k}{k!}$$

### 3.2 递推关系

更实用的计算方法是**递推公式**：
$$\sum_{k=0}^{m} \binom{m+1}{k} B_k = 0 \quad \text{对于} \ m \ge 1$$

由此可得：
$$B_m = -\frac{1}{m+1} \sum_{k=0}^{m-1} \binom{m+1}{k} B_k$$

### 3.3 前几个伯努利数

| n | $B_n$ | 说明 |
|---|-------|------|
| 0 | 1 | $B_0 = 1$ |
| 1 | -1/2 | |
| 2 | 1/6 | |
| 3 | 0 | （所有奇数大于1的都是0）|
| 4 | -1/30 | |
| 5 | 0 | |
| 6 | 1/42 | |
| 7 | 0 | |
| 8 | -1/30 | |
| 9 | 0 | |
| 10 | 5/66 | |
| 12 | -691/2730 | 这是第一个分子大于1的|
| 14 | 7/6 | |

**重要性质**：
- 对于所有奇数 $n > 1$，有 $B_n = 0$
- $B_1 = -1/2$ 是唯一的非零奇数项
- 符号交替变化：$B_2 > 0, B_4 < 0, B_6 > 0, \ldots$

### 3.4 计算示例

**计算 $B_2$：**

使用递推式（m=2）：
$$\binom{3}{0}B_0 + \binom{3}{1}B_1 + \binom{3}{2}B_2 = 0$$
$$1 \cdot 1 + 3 \cdot \left(-\frac{1}{2}\right) + 3 \cdot B_2 = 0$$
$$1 - \frac{3}{2} + 3B_2 = 0$$
$$3B_2 = \frac{1}{2}$$
$$B_2 = \frac{1}{6}$$

**计算 $B_4$：**

使用递推式（m=4）：
$$\sum_{k=0}^4 \binom{5}{k} B_k = 0$$
$$1 \cdot 1 + 5 \cdot \left(-\frac{1}{2}\right) + 10 \cdot \frac{1}{6} + 10 \cdot 0 + 5 \cdot B_4 = 0$$
$$1 - \frac{5}{2} + \frac{10}{6} + 5B_4 = 0$$
$$-\frac{3}{2} + \frac{5}{3} + 5B_4 = 0$$
$$-\frac{9}{6} + \frac{10}{6} + 5B_4 = 0$$
$$\frac{1}{6} + 5B_4 = 0$$
$$B_4 = -\frac{1}{30}$$

## 4. 使用 Faulhaber's Formula 计算具体例子

### 示例：计算 $S_4(n) = \sum_{k=1}^n k^4$

使用公式：
$$S_4(n) = \frac{1}{5} \sum_{j=0}^4 \binom{5}{j} B_j \cdot n^{5-j}$$

展开：
$$S_4(n) = \frac{1}{5}\left[\binom{5}{0}B_0 n^5 + \binom{5}{1}B_1 n^4 + \binom{5}{2}B_2 n^3 + \binom{5}{3}B_3 n^2 + \binom{5}{4}B_4 n\right]$$

代入伯努利数：
$$S_4(n) = \frac{1}{5}\left[1 \cdot 1 \cdot n^5 + 5 \cdot \left(-\frac{1}{2}\right) n^4 + 10 \cdot \frac{1}{6} n^3 + 10 \cdot 0 \cdot n^2 + 5 \cdot \left(-\frac{1}{30}\right) n\right]$$

化简：
$$S_4(n) = \frac{1}{5}\left[n^5 - \frac{5}{2}n^4 + \frac{10}{6}n^3 - \frac{1}{6}n\right]$$

$$= \frac{1}{5}\left[n^5 - \frac{5}{2}n^4 + \frac{5}{3}n^3 - \frac{1}{6}n\right]$$

$$= \frac{n^5}{5} - \frac{n^4}{2} + \frac{n^3}{3} - \frac{n}{30}$$

通分后：
$$S_4(n) = \frac{6n^5 - 15n^4 + 10n^3 - n}{30}$$

$$= \frac{n(n+1)(2n+1)(3n^2+3n-1)}{30}$$

## 5. 算法实现思路

### 计算伯努利数的 Python 代码框架：

```python
def bernoulli_numbers(n):
    """
    计算前n个伯努利数 B_0 到 B_n
    """
    B = [0] * (n + 1)
    B[0] = 1  # B_0 = 1
    
    for m in range(1, n + 1):
        sum_val = 0
        for k in range(m):
            # 二项式系数 binomial(m+1, k)
            binom = binomial_coefficient(m + 1, k)
            sum_val += binom * B[k]
        B[m] = -sum_val / (m + 1)
    
    return B
```

### 计算幂和：

```python
def faulhaber_sum(p, n):
    """
    使用Faulhaber公式计算 sum_{k=1}^n k^p
    """
    B = bernoulli_numbers(p)
    result = 0
    
    for j in range(p + 1):
        binom = binomial_coefficient(p + 1, j)
        term = binom * B[j] * (n ** (p + 1 - j))
        result += term
    
    return result / (p + 1)
```

## 6. 重要关系与推广

### 6.1 与黎曼ζ函数的关系

对于偶数伯努利数：
$$B_{2k} = (-1)^{k+1} \frac{2(2k)!}{(2\pi)^{2k}} \zeta(2k)$$

例如：
- $\zeta(2) = \sum_{n=1}^\infty \frac{1}{n^2} = \frac{\pi^2}{6}$
- $B_2 = \frac{1}{6}$

### 6.2 与多项式的关系

$S_p(n)$ 是关于 $n$ 的 $p+1$ 次多项式，且：
- 常数项为 0
- 系数为有理数
- 首项系数为 $\frac{1}{p+1}$

### 6.3 生成函数视角

幂和也有生成函数表示：
$$\sum_{p=0}^\infty S_p(n) \frac{t^p}{p!} = \frac{t}{e^t - 1} \cdot \frac{e^{(n+1)t} - 1}{t}$$

这正是伯努利数生成函数与几何级数的乘积。

## 7. 计算复杂度与优化

- 直接计算：$O(n \cdot p)$（暴力求和）
- Faulhaber公式：$O(p^2)$（计算伯努利数）+ $O(p)$（求值）
- 对于**固定p，多次查询不同n**：预处理伯努利数后每次查询 $O(p)$
- 对于**大p**（如 $p > 10^6$）：需要更高级算法和模运算