用伽马函数表示伯恩斯坦基函数

如果你接触过 Bezier 曲线，几乎一定见过伯恩斯坦基函数。它看起来像一个“公式细节”，但其实背后有一个很实用的问题：

同样的数学式子，怎么写得更不容易算崩？

先说结论

伯恩斯坦基函数的标准形式是

\[ B_r^d(t)=\binom{d}{r}t^r(1-t)^{d-r},\quad r=0,1,\ldots,d,\; t\in[0,1]. \]

关键在组合数 \(\binom{d}{r}\)。它常见写法是

\[ \binom{d}{r}=\frac{d!}{r!(d-r)!}. \]

问题来了：当阶数 \(d\) 很大时，阶乘会非常大，直接计算容易溢出或者损失精度。

更稳妥的做法是：用伽马函数，特别是用 gammaln（对数伽马函数）来算。

直觉理解：它在“分配权重”

可以把 \(B_r^d(t)\) 理解成一组权重函数：

1. 每个 \(r\) 对应一个权重。
2. 这些权重在 \(t\in[0,1]\) 内非负。
3. 所有权重加起来等于 1。

这也是为什么它非常适合曲线插值和形状控制。

数学推导：阶乘如何变成伽马函数

伽马函数有一个经典关系：

\[ \Gamma(n+1)=n!,\quad n\in\mathbb{N}. \]

把它代入组合数：

\[ \binom{d}{r}=\frac{\Gamma(d+1)}{\Gamma(r+1)\Gamma(d-r+1)}. \]

再代回伯恩斯坦基函数：

\[ B_r^d(t)=\frac{\Gamma(d+1)}{\Gamma(r+1)\Gamma(d-r+1)}t^r(1-t)^{d-r}. \]

到这一步还是“直接算伽马函数”。在数值计算里，我们通常再往前走一步，用对数形式：

\[ \mathrm{gammaln}(x)=\log\Gamma(x). \]

于是

\[ \binom{d}{r}=\exp\left(\mathrm{gammaln}(d+1)-\mathrm{gammaln}(r+1)-\mathrm{gammaln}(d-r+1)\right). \]

最终得到更稳定的伯恩斯坦基函数写法：

\[ B_r^d(t)=\exp\left(\mathrm{gammaln}(d+1)-\mathrm{gammaln}(r+1)-\mathrm{gammaln}(d-r+1)\right)t^r(1-t)^{d-r}. \]

代码实现（MATLAB）

本文使用的实现如下：

% Standard Bernstein Basis Polynomial
function v = bern(r, d_, t)
    v = exp(gammaln(d_+1)-gammaln(r+1)-gammaln(d_-r+1)) .* t.^r .* (1-t).^(d_-r);
end

为什么 gammaln 更推荐

1. 避免“大数爆炸” 直接算 \(d!\) 很快就会变成超大数，浮点数容易顶不住；gammaln 先在对数域处理，稳定很多。

2. 降低精度损失 组合数可能极大或极小。对数域里先做加减，再 exp，比“先乘出大数再相除”更安全。

3. 更适合高阶曲线曲面计算 在高阶 Bezier/Bernstein 场景中，更不容易出现 Inf、NaN。

4. 数学上更通用 伽马函数是阶乘的连续推广，在一些扩展模型中更方便。