企业官网建设流程全解析

伽马函数的定义和性质，这两块内容我这些年被问过很多次，但问法差别挺大：做概率统计的人关心它和伽马分布、卡方分布的关系；写数值代码的人关心它为什么总在负整数上抛异常、为什么gamma(200)直接返回inf；学复分析的人则盯着它的定义域和解析延拓。我自己的体会是，这三拨人卡住的地方其实是同一处——没把定义看清楚，也没搞明白那几条性质各自成立的条件和适用场合。所以这篇按我平时带人的顺序来：先把定义从阶乘里拆出来，再把性质一条条过一遍，最后落到能跑起来的数值实现和我自己踩过的坑上。不管你是刚学完微积分，还是已经在项目里调过几十次lgamma，应该都能从里面捞到点东西。

1. 先把伽马函数当成"阶乘的连续延拓"来理解

1.1 阶乘的三个痛点：为什么非得造一个新函数

阶乘是个特别朴素的东西，n! = 1×2×…×n，小学奥数就在用。但真拿它去算东西，三个问题会立刻冒出来。第一，它只在非负整数上有定义，你问我 0.5 的阶乘是多少，这个问题在阶乘的语言体系里根本没有答案；第二，它长得太快，20! 就已经是 2.4×10^18，超过 64 位整数上限，100! 用双精度浮点直接溢出，做组合计数的时候很容易一算就废；第三，它只给出离散的点，做渐近分析、做积分估计、做连续优化的时候，你需要一条光滑曲线从这些点上穿过去，光有离散值没法求导、没法积分。

这三个痛点对应的正是三个真实需求。做概率的人要把二项分布推广到实数次试验，中间必然出现非整数的"阶乘"；做算法分析的人要估计 n! 的量级，需要一条能求对数的光滑曲线；做积分的人遇到形如 ∫₀^∞ t^{n}e^{-t}dt 的式子，希望有一个统一的符号把它写下来。伽马函数就是为这三件事同时准备的工具，在正整数上它精确等于 (n-1)!，在其他位置它给出一条光滑、可导、可积、可解析延拓的曲线。

这里要强调一个容易混淆的点：伽马函数满足的是 Γ(n) = (n-1)!，不是 Γ(n) = n!。这个 "-1" 的偏移是历史遗留，也是新手最容易在第一次写代码时搞错的地方。你写gamma(5)得到的是 24 而不是 120，很多 bug 就是从这儿来的。

1.2 一个积分怎么把离散阶乘铺成连续曲线

构造思路其实很自然，从一个熟得不能再熟的积分出发：

∫₀^∞ e^{-t} dt = 1

现在往里面塞一个 t：

∫₀^∞ t e^{-t} dt = 1

再来一次分部积分，∫ t² e^{-t} dt = 2，∫ t³ e^{-t} dt = 6，∫ t⁴ e^{-t} dt = 24。规律非常清楚：

∫₀^∞ t^{n} e^{-t} dt = n!

这个式子对每个非负整数 n 都成立。既然它成立，那么把离散的 n 换成连续的变量，整条曲线就出来了。标准写法是用 z-1 代替 n，也就是

Γ(z) = ∫₀^∞ t^{z-1} e^{-t} dt

为什么偏要用 z-1 而不是 z？纯粹是为了让递推关系写得干净。用 z-1 的话，分部积分一次得到 Γ(z+1) = zΓ(z)，形式和阶乘的 (n+1)! = (n+1)·n! 完全对齐。要是用 z，公式就变成 Γ(z+1) = (z+1)Γ(z)，到处多一个加一，写推导的时候会很难受。这种"为了公式美观而做的下标平移"在数学里很常见，理解动机之后就不用死记了。

还有个问题值得想一想：为什么核函数一定是 e^{-t}？换别的行不行？答案是，要让"分部积分一次就掉一个因子"这个性质成立，核函数必须满足 f'(t) = -f(t)，也就是 f(t) = Ce^{-t} 这一族。所以 e^{-t} 不是随便挑的，是被递推关系唯一确定下来的。这也是为什么伽马函数看起来"来路不明"却又处处出现——它不是人造的技巧，而是把"分部积分降幂"这个动作连续化之后的必然产物。

1.3 收敛域这条线是怎么划出来的

写完定义式，紧接着必须问一句：这个积分什么时候收敛？把积分拆成 (0,1) 和 (1,∞) 两段看。

在 t 趋近 0 这一端，e^{-t} 趋近 1，被积函数的行为等价于 t^{z-1}。众所周知 ∫₀¹ t^{a} dt 收敛当且仅当 a > -1，套进来就是 Re(z-1) > -1，即 Re(z) > 0。这一端是真正卡脖子的地方，因为 t^{z-1} 在原点附近会爆。

在 t 趋近无穷这一端，e^{-t} 的衰减是指数级的，比任何多项式 t^{Re(z)-1} 都快，所以只要 Re(z) 有限，这一端永远收敛。

合起来，积分定义的收敛域就是 Re(z) > 0 这个右半平面。但伽马函数的实际定义域远不止于此，因为可以靠递推式往外推：

Γ(z) = Γ(z+1) / z

只要 z 不是 0 或负整数，右边就有意义，于是左半平面也被覆盖了。这种操作叫解析延拓。从 -1 往 -2 推，从 -2 往 -3 推，一路推下去，得到的函数在整个复平面上解析，唯独在 z = 0, -1, -2, … 这些点上爆掉，构成一列一阶极点。这些极点的位置不是"算不出来"的意外，而是延拓过程里分母逐个变零的必然结果。

注意：很多教程写"伽马函数定义域是 z>0"，那是只针对积分表达式本身说的。做复分析或者查特殊函数手册时，标准结论是整个复平面去掉非正整数。两句话都对，区别在于你在讨论积分式还是延拓后的函数，写文档和代码注释时最好分开说清楚，不然接手的人会懵。

2. 定义式的几种等价写法，各自适合什么场合

2.1 积分定义：推导解析结论时的首选

积分形式 Γ(z) = ∫₀^∞ t^{z-1}e^{-t}dt 是几乎所有教科书的第一定义，它的好处是直观、可操作。你要证递推式，分部积分；要证与贝塔函数的关系，换个二重积分的积分顺序；要估计量级，用放缩。绝大多数基础结论从这条式子出发都能在两三行内推完。

但它有个明显短板：一旦离开 Re(z)>0，这条式子就不成立了，不能直接拿来算负半轴的值。所以我在做理论推导时用它，在做数值计算时完全不用它——数值上用它是自找麻烦，无穷区间积分要么截断要么变换，精度控制非常磨人。

2.2 无穷乘积定义：看清零点与极点的最佳视角

另一条路是把伽马函数的倒数写成一个无穷乘积：

1/Γ(z) = z·e^{γz}·∏_{n=1}^{∞} (1 + z/n)·e^{-z/n}

这里的 γ 是欧拉-马歇罗尼常数，约等于 0.5772。这条式子的价值在于它把结构暴露得非常清楚：1/Γ(z) 是一个整函数，没有极点，所以 Γ(z) 的极点只可能来自这个整函数的零点。你看乘积里那个 z 因子，在 z=0 处给一个零点；每个 (1+z/n) 因子在 z=-n 处给一个零点。零点位置一目了然，正好是非正整数。

反过来说，Γ(z) 本身没有零点，因为 1/Γ(z) 永远不为零（乘积里每一项都不为零，指数因子也不为零）。这个结论在解方程和做渐近估计时经常用到，是个容易被忽略但很好用的事实。

这条定义在理论分析里非常强，收敛性证明、反射公式推导、对数凸性证明都靠它。但它的收敛速度很慢，写代码直接照抄这个无穷乘积，算到 z=10 你可能要累几十万项才勉强对得上，实用价值几乎为零。

2.3 递推关系：数值计算真正的主力

Γ(z+1) = zΓ(z)

这条式子看起来平平无奇，但它才是所有实用算法的骨架。原因很直接：它允许你把任意参数搬到收敛快、精度高的区间去算。比如要算 Γ(0.001)，直接代定义式是个灾难，但用递推式写成

Γ(0.001) = Γ(1.001) / 0.001

分母那个 0.001 提出来，剩下的 Γ(1.001) 就在 1 附近，任何近似方法在那儿都表现极好。这个"参数搬家"的操作是所有伽马函数库的核心套路，没有例外。

反过来，参数很大的时候也是同样的思路：把大参数一路减到 1 附近，把每次减掉的因子记到对数里，最后加回去。这套移位加递推的组合，后面第 4 节会展开写。

2.4 与贝塔函数的桥：工程积分里最常撞见

B(a,b) = Γ(a)Γ(b) / Γ(a+b)

这条关系式的实用频率高得惊人。凡是形如 ∫₀¹ t^{a-1}(1-t)^{b-1}dt 的积分，都可以直接换成三个伽马函数，很多原本要硬算的定积分就这么没了。举个具体例子：B(2,3) = Γ(2)Γ(3)/Γ(5) = 1×2/24 = 1/12，手算验证一遍 ∫₀¹ t(1-t)²dt，展开成 ∫(t - 2t² + t³)dt = 1/2 - 2/3 + 1/4 = 1/12，对上了。

在概率论里这条关系更是随处可见。贝塔分布的归一化常数是 1/B(α,β)，狄利克雷分布、t 分布、F 分布里也都藏着它。我在一个做统计可视化的模拟项目里就干过一件事：把贝塔分布的概率密度函数从"调库直接算"改成"手动用 lgamma 组合"，性能没什么变化，但对数值稳定性的理解清楚了很多，后面调参时不再盲猜。

3. 核心性质逐条拆解：哪些能直接拿来算，哪些有前提

3.1 特殊值：整数、半整数、1/2 这三族要背下来

先把最常用的几个值列清楚。整数族最简单，Γ(n) = (n-1)!，没什么好说的。真正有意思的是半整数族，起点是

Γ(1/2) = √π ≈ 1.7724538509

这个值怎么来的？把 z = 1/2 代进定义式，做换元 t = u²，得到 ∫₀^∞ u^{-1}e^{-u²}·2u du = 2∫₀^∞ e^{-u²}du，后面这个就是标准的高斯积分，等于 √π/2，乘 2 得 √π。整条推导只有两步，但结果出人意料地"不讲道理"——一个纯指数和幂次的积分，结果里冒出个圆周率。这种"跨常数"的现象是伽马函数最迷人的地方之一。

从 Γ(1/2) 出发，靠递推式可以批量生产所有半整数：

Γ(3/2) = (1/2)Γ(1/2) = √π/2 Γ(5/2) = (3/2)(1/2)Γ(1/2) = 3√π/4 Γ(7/2) = 15√π/8

一般式是 Γ(n + 1/2) = [(2n)! / (4^n·n!)]·√π。这个式子在做统计里的卡方分布、正态矩计算时高频出现，建议直接记住形式，省得每次现场推。

再看一个具体数：Γ(5.5) = Γ(11/2)。按公式，n=5 时系数是 (10)!/(4⁵·5!) = 3628800/(1024×120) = 3628800/122880 = 29.53125，乘 √π 得约 52.3428。这个值在纯手算时稍微麻烦，但有了模式就只是乘法。

3.2 反射公式：负半轴和 (0,1) 区间都能救

Γ(z)Γ(1-z) = π / sin(πz)

这是伽马函数最有用的性质之一，因为它把"远离原点"的参数和"靠近原点"的参数挂上了钩，同时还能处理负数。举两个用法。

求 Γ(-1/2)：取 z = -1/2，右边是 π/sin(-π/2) = -π，左边是 Γ(-1/2)Γ(3/2)，而 Γ(3/2) = √π/2，所以 Γ(-1/2) = -π/(√π/2) = -2√π ≈ -3.5449。整个过程只用了反射加递推，完全没碰积分定义。

求 Γ(0.1)：这个值在 (0,1) 里，直接算收敛但慢。用反射写成 Γ(0.1) = π/(sin(0.1π)·Γ(0.9))，把参数搬到 0.9，虽然改善有限但确实更快收敛。更极端的例子是 Γ(0.001)，用反射变成 π/(sin(0.001π)·Γ(0.999))，好处主要来自分母那个 sin 能算出精确的小量，避免直接计算时有效数字损失。

提示：反射公式在数值实现里的真正作用是处理 z < 0.5 的分支。很多开源实现就是靠它在 0.5 处一刀切，左边用反射、右边用连乘或级数，两段拼起来覆盖全轴。写自己的库时照这个思路做，代码会短很多。

3.3 倍乘公式：半整数族的批量算法

Γ(z)Γ(z+1/2) = 2^{1-2z}·√π·Γ(2z)

这条公式初看有点绕，但它的用途很实际：当你需要一次算两个相隔 1/2 的值时，可以合并成一次计算。比如统计里算 Beta(1/2, 1/2) 的归一化常数，就会出现 Γ(1/2)Γ(1/2+1/2) 这种组合，直接用倍乘化简成一项。

验证一下 z = 1/2 的情况：左边 Γ(1/2)Γ(1) = √π，右边 2⁰·√π·Γ(1) = √π，成立。z = 1 时：左边 Γ(1)Γ(3/2) = √π/2，右边 2^{-1}·√π·Γ(2) = √π/2，也成立。这种"代进去验一遍"的习惯建议保留，公式多的时候特别管用。

更一般的形式是多重倍乘公式，把 z, z+1/n, …, z+(n-1)/n 乘在一起，右边化简成 n 的幂次乘 Γ(nz)。这条在数论和某些积分恒等式里会碰到，日常工程用得少，但知道它存在，遇到奇怪的乘积式时能想起来。

3.4 对数凸性：为什么这个延拓是"唯一正确"的

这一条是理论上最漂亮的部分。问题是这样：满足 f(x+1) = x·f(x) 且 f(1) = 1 的函数有无穷多个，随便在 (0,1) 上画一条光滑曲线连到 1 就行了。伽马函数凭什么说自己是"标准答案"？

答案是再加一条限制：对数凸性，也就是 ln Γ(x) 是凸函数。加上这条之后，解就唯一了。这条定理的分量在于，它给出了一个完全公理化的刻画——你不需要知道积分式长什么样，只要这三条要求（递推、初值、对数凸性），函数就被钉死了。

为什么对数凸性而不是普通凸性？因为要做插值的话，"二阶导非负"这个条件太弱，能塞进去的曲线还是很多；而对数凸性更强，配上递推式就能通过 Hölder 不等式把中间值夹死，上下界一碰，唯一性就出来了。这套证明在标准教材里通常占一页左右，值得自己动手推一遍，比背结论有用得多。

对数凸性还有个实用推论：ln Γ 是凸函数，意味着它在 (0,∞) 上有唯一的极小值点。这个点在哪？由 ψ(x) = d/dx ln Γ(x) = 0 决定，数值大约在 x ≈ 1.4616 处，最小值约 0.8856。这个数字在做最优化和不等式估计时偶尔会用到，值得记一下。

3.5 极点、留数与渐进展开：远离原点时怎么办

极点部分前面说过，z = 0, -1, -2, … 全是一阶极点。每个极点的留数是

Res(Γ, -n) = (-1)^n / n!

这个结论从递推式一眼就能看出来：Γ(z) 在 z=0 附近的行为是 1/z 减去一个常数加高阶项，在 z=-1 附近是 -1/(z+1) 加常数，符号交替，系数是阶乘的倒数。

至于大参数的行为，靠的是渐进展开。最基础的一阶形式是

ln Γ(z) ≈ (z - 1/2)ln z - z + (1/2)ln(2π) + 1/(12z)

再加修正项就可以更精确，下一项是 -1/(360z³)，再下一项是 1/(1260z⁵)，符号交替，量级迅速衰减。这套展开在 z 大于 10 的时候已经准得离谱，z 到 100 的时候基本可以只用前两三项。

拿 Γ(1000) 试一下。ln 1000 = 6.907755279，(z-1/2)ln z = 999.5 × 6.907755279 ≈ 6904.3014，减 1000 得 5904.3014，加 0.5·ln(2π) = 0.9189385 得 5905.2203，加 1/12000 = 0.0000833 得 5905.2204。换算成常用对数：5905.2204 / 2.302585 ≈ 2564.61，也就是 Γ(1000) 约等于 4.02×10^2564。这个结果和阶乘的真实量级对得上，而整个计算只用了四次乘除和一次对数——注意，如果我直接算 999! 再取对数，中间过程早就溢出了。

4. 数值实现：从零写一个靠得住的 lgamma

4.1 先算对数：为什么不直接算 Γ

前面提过 Γ(172) 左右就跨过双精度浮点的上限（约 1.8×10^308），而实际工程里参数动辄几百上千，直接算必然溢出。所以几乎所有严肃的实现都不提供"Γ 本身"，而是提供 ln Γ，业界习惯叫lgamma。

这里有个非常实际的取舍：既然 Γ 会溢出，那我能不能先算 ln Γ，再 exp 回去？能，但要注意 exp(lgamma(x)) 依然会在同样的地方溢出——因为问题出在结果本身太大，不是算法有问题。所以正确的做法是，能在对数域里做完的计算，就一直待在对数域里，到最后一步再决定要不要取指数。

举个我在一个统计计算 Demo 里碰到的真实场景：算贝塔分布的对数密度，公式是

ln p(x) = (α-1)ln x + (β-1)ln(1-x) - ln B(α,β)

其中 ln B(α,β) = lgamma(α) + lgamma(β) - lgamma(α+β)。当 α、β 都取到 500 这个量级时，三个 lgamma 各自都是几千的量级，但相减之后结果很小。这个"大数相减得小数"的结构是典型的精度陷阱，后面第 5 节会专门讲怎么处理。

4.2 一套能用的连乘近似，附完整代码

工程上最流行的实数实现是一类连乘近似，思路是把 Γ(z+1) 拆成 √(2π)·t^{z+1/2}·e^{-t}·A(z) 的形式，其中 A(z) 是一个有理式的部分分式和。下面这套系数在公开资料里流传很广，9 个系数覆盖双精度需求足够：

import math def my_lgamma(z): """仅对实数 z 实现 log|Gamma(z)|，覆盖正负全轴。""" # 经典 9 系数连乘近似 g = 7 c = [0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7] if z < 0.5: # 反射公式，处理左半轴和小参数 # Gamma(z)Gamma(1-z) = pi / sin(pi z) return (math.log(math.pi) - math.log(abs(math.sin(math.pi * z))) - my_lgamma(1.0 - z)) z -= 1.0 x = c[0] for i in range(1, g + 2): x += c[i] / (z + i) t = z + g + 0.5 return (0.5 * math.log(2 * math.pi) + (z + 0.5) * math.log(t) - t + math.log(x))

这段代码有几个细节值得单独说。第一，函数内部把 z 减了 1，因为这套系数拟合的是 Γ(z+1) 而不是 Γ(z)，这个 1 的偏移和前面定义的偏移是两回事，别搞混。第二，循环上界是 g+2，因为 g=7 的配置需要 9 个系数，容易写成 g+1 导致最后一个系数没被用上。第三，反射分支用的是 log(π) - log|sin(πz)| 而不是 log(π/sin)，这样在 sin 接近零时不会先产生一个巨大的中间值。

跑一下验证：my_lgamma(5)应该等于 ln(24) ≈ 3.1780538，my_lgamma(0.5)应该等于 ln(√π) ≈ 0.5723649，my_lgamma(-0.5)应该等于 ln(2√π) ≈ 1.2655121（因为 Γ(-0.5) 是负数，这里取的是绝对值对数）。三个都对得上，说明实现没写歪。

4.3 大参数用渐进，小参数用移位，两头都要顾

连乘近似在 z 处于 [0.5, 10] 这个区间表现最好，相对误差能压到 1e-15 量级。参数再大一点，它也不会崩，但会开始吃力：t = z+g+0.5 里那个固定加数在大 z 时贡献越来越小，等于白白浪费了精度。

稳妥的工程做法是分段：

渐进展开那段可以直接写成：

def stirling_lgamma(z): """z >= 10 时使用，返回 log Gamma(z)""" inv = 1.0 / z inv2 = inv * inv # 1/(12z) - 1/(360z^3) + 1/(1260z^5) - 1/(1680z^7) series = inv / 12.0 - inv * inv2 / 360.0 \ + inv * inv2 * inv2 / 1260.0 \ - inv * inv2 * inv2 * inv2 / 1680.0 return ((z - 0.5) * math.log(z) - z + 0.5 * math.log(2 * math.pi) + series)

为什么不干脆全用渐进？因为 z 小于 10 的时候展开式发散得很快，z=1 时误差能到百分之一以上，完全不可用。这就是为什么必须有分段，而不能"一套公式打天下"。

还有一个更稳的组合策略：先把大参数一路用递推减到 [0.5, 10] 区间，用连乘算完，再把减掉的那部分对数加回来。这样做的好处是全程只用一套核心算法，误差来源单一，好排查。代价是循环次数随 z 线性增长，z=10000 就要循环一万次，性能上不划算。所以实际选型要看场景：单次调用、精度优先就分段；批量调用、参数分布集中就预先选好分支。

4.4 复数参数与精度校验

复数版本在物理和信号处理里很常见，实现上就是把实数分支换成级数加反射，代码量会翻几倍。如果只是偶尔用一下，直接调成熟库省事得多，scipy.special.gamma、mpmath.gamma都是现成的。我的建议是：除非你要在嵌入式或者受限环境里跑，否则别自己造复数轮子。

校验方法我一般用三招。第一招是硬编码几个已知值做断言，Γ(1)=1、Γ(0.5)=√π、Γ(5)=24、Γ(-1.5)=4√π/3≈2.3633。第二招是用递推关系做自洽检查，随机取 z，验证 Γ(z+1) 和 zΓ(z) 的差异在容差范围内，这招能抓出绝大多数分支写错的问题。第三招是和成熟库对比，取 1000 个对数均匀分布的点，算最大相对误差，正常应该在 1e-14 以下。

注意：做对比测试时，取点方式很关键。如果只在 [0.5, 10] 里均匀取点，你会觉得自己的实现完美无缺；一定要覆盖负半轴、极小值附近（1.46 附近）、以及远离原点的大参数区，才能真正暴露问题。我吃过这个亏，测了半天全绿，结果上线后用户传进来一个 1e-3 的参数直接崩了。

5. 常见问题与排查技巧实录

5.1 溢出与下溢：不只是"太大了"这么简单

最常见的报错就是返回inf。判断标准前面算过，Γ(x) 大约在 171.6 处跨过双精度上限。但很多人不知道另一半：小参数向下溢出同样会发生，而且更隐蔽。Γ(0.001) 大约是 999 左右，还好；但如果算的是贝塔函数的倒数或者某些归一化常数的倒数，很容易出现 1e-400 这种数，直接塌成 0，后面所有乘法全变成 0，最后结果看起来"什么也没算出来"，也没有任何报错。

处理办法统一：全程用对数。需要乘法的地方改成加法，需要除法的地方改成减法，需要比较大小的地方直接比对数（对数单调，比较结果不变）。只有在最终输出给用户的时候才取一次 exp，并且在外层做好范围判断，超出就返回一个明确的哨兵值而不是硬算。

5.2 极点附近的报错：负整数要和"接近负整数"区别对待

传进来 -3 这个值，正确的行为是报错或者返回 inf，因为在 -3 处函数真的有极点。但传进来 -2.9999999，那就应该正常算出一个极大的值（正负取决于从哪边逼近）。这两种情况的代码路径完全不同，混在一起处理就会出问题。

实现上的建议是：在入口处判断z == round(z) and z <= 0，命中就走极点分支；否则一律走正常计算。用==判断浮点数相等在别的地方是坏习惯，但在这里恰好是合适的，因为用户如果真想传一个精确的负整数，传进来的就是这样。更保险的做法是加一个极小的容差，比如abs(z - round(z)) < 1e-12，命中后提示"参数过于接近极点，结果不可靠"。

还有一个在概率计算里特别常见的坑：算 Dirichlet 分布或者贝塔分布时，参数里有 α=1、β=1 这种值，然后某个中间式子里出现了 Γ(α+β-n) 且 α+β-n 恰好是 0 或负整数。数学上这说明该配置在这个点上确实退化，代码里应该提前拦截，而不是让它在深处炸掉，留一堆看不懂的调用栈。

5.3 精度灾难：大数相减的典型症状与解法

前面提到的 lgamma(500) + lgamma(500) - lgamma(1000) 就是经典案例。三个数各自都是 2000 多，相减得到的结果在几十这个量级，有效数字白白丢了两位。参数再大一点，比如 α=β=1e6，三个 lgamma 都在 1e7 量级，相减之后的有效位数只剩个位数，结果基本不可信。

标准解法是利用渐进展开做抵消。把 lgamma(α)、lgamma(β)、lgamma(α+β) 分别展开，那些 1e7 量级的主项在代数上会精确抵消，剩下的修正项量级很小，直接用双精度算就行。具体展开形式在讲大参数性质那一节里给过，把三份代进去做符号整理即可。这个过程手推一次会有点长，但推完之后你对"为什么库函数能在极端参数下给出准确结果"就彻底明白了。

如果不想手推，还有个偷懒但有效的办法：用高精度库（比如mpmath设 50 位精度）算一遍作为基准，和双精度结果对比，如果差异超过预期，就说明你的实现里存在抵消问题，需要专门处理。这个方法我用来定位过好几个隐藏很深的精度 bug。

5.4 一张速查表顶十次翻书

6. 这些性质最后都落到哪些场景里

6.1 概率分布族：一整个家族都靠它撑着

伽马分布、卡方分布、贝塔分布、t 分布、狄利克雷分布，这几个的密度函数里全都有伽马函数。以卡方分布为例，自由度 k 的密度是

f(x) = x^{k/2-1}·e^{-x/2} / (2^{k/2}·Γ(k/2))

自由度取奇数的时候，k/2 就是半整数，这时候 Γ(k/2) 的值就是前面那个半整数族的公式，直接套。取偶数时是整数阶乘。这就是为什么统计教材里总要单独讲一遍半整数阶乘——它不是为了好看，是真要用。

实际写代码的时候，密度函数一律在对数域里拼：(k/2-1)*log(x) - x/2 - (k/2)*log(2) - lgamma(k/2)。这套写法在 k 很大的时候依然稳定，而直接按公式算分子分母，k 过百就废了。

6.2 组合与算法分析：把离散阶乘变成连续函数

做算法复杂度分析时，n! 的估计是基本操作。传统的做法是取对数然后用渐进展开，而这套展开正是伽马函数大参数性质的直接应用。有了连续版本之后，很多原本只能对整数讨论的东西可以求导了，比如"n 取多大时某个组合数超过阈值"这类问题，把 n 当连续变量处理后直接解方程就行，不用逐个试。

还有一个不那么显眼但很实用的技巧：形如 ∫₀^∞ x^{a-1}e^{-bx}dx 的积分，换元 t = bx 之后直接等于 Γ(a)/b^a。这类积分在算期望、算矩、算归一化常数时天天出现，记住这个换算可以省掉大量分部积分。举个具体的：∫₀^∞ x²e^{-2x}dx = Γ(3)/2³ = 2/8 = 0.25，两步搞定。

6.3 几何与物理：高维体积里的那个分母

n 维单位球的体积公式是 V_n = π^{n/2} / Γ(n/2 + 1)。n=2 时，π/Γ(2) = π，正确；n=3 时，π^{1.5}/Γ(2.5) = 5.5683/1.3293 ≈ 4.1888，正好是 4π/3，正确。这个公式的美妙之处在于它能一路推广到非整数维，做分数维分析的时候直接可用。

物理和工程里的积分变换、某些特殊函数的级数展开、热传导方程的解，也都会冒出伽马函数。遇到的时候不用慌，绝大多数情况只需要两件事：知道它在整数和半整数上的取值，以及会用递推式把参数搬到好算的区间。

我个人在实际操作中的体会是，伽马函数真正难的不是那些公式，而是"什么时候该用哪条"。做题时容易陷入"背了一堆公式但不知道用哪个"的状态，破局的办法只有一个：把它当成工具箱而不是知识清单。你手上要解决的具体问题是什么形状——是负参数，是半整数，是大参数，还是要跟贝塔函数打交道——想清楚这一点，对应的那条性质自然就跳出来了。另外强烈建议自己动手写一遍 lgamma 的实现，哪怕只用一次，写完之后再回头看那些性质，感受会完全不一样：每一条公式都不再是纸上的符号，而是你在代码某个分支里真真切切写过的那几行。