伽马函数从定义到工程实践：收敛域、解析延拓与数值计算要点

“伽马函数”这四个字在不少人的数学记忆里停留得很浅：教材用一页纸给出一个反常积分，顺手推出 Γ(n+1)=n!，然后翻篇。等到真正动手做数值计算、推导概率分布的归一化常数、或者要把一个反常积分写成闭式解时，才会发现那一页纸远远不够用——为什么收敛域偏偏是 Re(z)>0，为什么负数处会炸，为什么代码库里明明有 gamma 函数却几乎所有人都改用 lgamma，这些问题在课本上常常被“容易验证”“显然成立”几个字带过去。

这篇想做的事情很具体：把伽马函数从定义到性质的整条链条摊开，把教材省略掉的环节补上——收敛域是怎么一段一段算出来的、解析延拓每一步在做什么、几个常用公式各自的有效边界在哪里、以及在把这些公式落进代码时会在哪些地方被咬一口。内容面向两类人：一类是正在学数学分析或复变函数、想把这块知识补扎实的同学；另一类是做统计、物理、工程计算，需要在程序里正确调用这个函数的从业者。定义和推导会写清楚，但不会为了严谨把可读性牺牲掉；工程细节会给到能直接照着敲的参数和写法，也会说明为什么那样写而不是另一样。

1. 从阶乘插值说起：伽马函数的真实来历

1.1 0.5 的阶乘该等于多少

先抛开积分式子不谈，问一个更朴素的问题。1!=1，2!=2，3!=6，4!=24，把这几个点画在坐标纸上，能不能用一条光滑曲线连起来，然后问这条曲线在 x=0.5 处取多少？

问题一问出来，麻烦就来了：它没有唯一答案。给定平面上有限个点，穿过它们的光滑函数有无穷多条。你可以让曲线在相邻两点之间随便鼓包或凹陷，只要保持光滑就行，而这些改动完全不影响它在整数点上的取值。所以“把阶乘连续化”这件事，光靠“穿过这些点”是远远不够的，必须额外加约束。

这就是很多人第一次接触伽马函数时最容易漏掉的逻辑：伽马函数不是“阶乘的唯一连续推广”，它是“在若干附加条件下唯一的推广”。附加条件是什么、为什么非得是这几条，比那个积分定义本身重要得多。搞清楚这一点，后面所有性质才有落脚点；搞不清楚，就只能把 Γ(1/2)=√π 当口诀背。

常见的那组条件是三条：f(1)=1；f(x+1)=x·f(x)，也就是递推结构保持不变；以及 ln f 在正半轴上是凸函数。前两条保证它和阶乘在整数点咬合、递归关系照搬，第三条负责砍掉那些“能取值但乱抖”的候选。

1.2 为什么是“对数凸”而不是“光滑”或“单调”

这里值得停一下，因为“为什么选凸性”是真正区分理解深度的地方。

光滑这个条件太弱。多项式插值、样条插值、无穷多种解析函数都能穿过那四个点并且光滑，加了等于没加。单调这个条件又太强且方向不对——阶乘本身在正半轴上确实是递增的，但伽马函数在 (0, 1.4616) 这一段其实是递减的，最小值约 0.8856，之后才转为递增。如果硬套“单调递增”，连阶乘本身都不能被正确推广。

对数凸则是刚好卡在中间的那个约束。它的意思是 ln f 满足凸性，等价于对任意 x₁、x₂ 和 λ∈[0,1]，有

$$\ln f\big((1-\lambda)x_1+\lambda x_2\big)\le (1-\lambda)\ln f(x_1)+\lambda \ln f(x_2)$$

把这个不等式指数化，就变成 f((1-λ)x₁+λx₂) ≤ f(x₁)^{1-λ} f(x₂)^λ。

阶乘满足这个吗？取中间点验证一下：f(3)f(1)=6·1=6，f(2)²=4，6≥4，成立。再取 f(4)f(2)=24·2=48，f(3)²=36，48≥36，也成立。这不是巧合，用柯西–施瓦茨或者 Hölder 不等式可以直接对积分定义证出来。

关键在于，对数凸这个条件足够严格，能把候选函数锁死到唯一一个；同时又足够宽松，能被积分定义满足。宽严之间的这个平衡点，才是伽马函数定义的核心设计。严格来说，这个唯一性结论叫 Bohr–Mollerup 定理，后面第 3 章会专门讲。

1.3 积分定义是被递推关系“逼”出来的

很多人以为 Γ(z)=∫₀^∞ t^{z-1}e^{-t}dt 是天上掉下来的，其实从递推关系倒推，这个式子是必然的产物。

要找 f 满足 f(x+1)=x·f(x)。看到右边那个“乘以 x”的结构，第一反应应该是什么？分部积分。因为分部积分最擅长的就是“把一个因子的指数降一，同时吐出这个指数”。

拿 I(a)=∫₀^∞ t^a e^{-t}dt 试一下：

$$I(a)=\Big[-t^a e^{-t}\Big]_0^\infty + a\int_0^\infty t^{a-1}e^{-t},dt = a,I(a-1)$$

边界项在两端都消失了——t→0 时 t^a→0，t→∞ 时 e^{-t} 把多项式压死。于是 I(a)=a·I(a-1)，正是想要的递推结构。

再取 a=0：I(0)=∫₀^∞ e^{-t}dt=1，正好对上 f(1)=1 这个归一化。

把指标平移一下，令 f(z)=I(z-1)，就得到 f(z+1)=z·f(z) 且 f(1)=1。所以那个积分定义不是凭空构造，而是“分部积分能产生乘以变量的效果”这一事实的直接后果。e^{-t} 负责在无穷远处保证收敛，t^{z-1} 负责提供那个乘子，两者的分工非常明确。

这也能解释为什么 e^{-x²} 那类形式做不了伽马函数：对它分部积分不会产生干净的单因子 x，只会拖出一堆系数。

2. 收敛域比定义本身更值得盯

2.1 Re(z)>0 是一段一段算出来的

积分定义写出来只是第一步，接下来必须回答“它在哪儿收敛”。做法是把积分在 t=1 处劈成两段，分别看。

先看 [1, ∞) 这一段。取 z=σ+iy，被积函数的模是

$$|t^{z-1}e^{-t}|=t^{\sigma-1}e^{-t}$$

当 t 很大时，e^{-t} 是指数级衰减，任何多项式级的 t^{σ-1} 都被它压制。所以这一段对任意复数 z 都绝对收敛，不需要任何附加条件。

再看 (0, 1] 这一段。这里 e^{-t} ≤ 1，所以

$$|t^{z-1}e^{-t}|\le t^{\sigma-1}$$

而 ∫₀¹ t^{σ-1}dt 在 σ>0 时等于 1/σ，收敛；σ=0 时变成 ∫₀¹ t^{-1}dt，对数发散；σ<0 时是更快的幂次发散。所以这一段的收敛条件精确地是 σ=Re(z)>0。

两段合起来，收敛域就是开右半平面 Re(z)>0，而且这个边界是紧的——在 Re(z)=0 上积分已经不再收敛。虚轴上的点 z=iy（y≠0）虽然被积函数的模是 t^{-1}，但相位因子 e^{iy ln t} 在 t→0 时振荡不趋于定值，极限不存在，同样按发散处理。

这里有个容易被忽略的实践含义：积分的收敛性来自两端的联合压制，缺一不可。t=0 附近的幂次可积性给出 Re(z)>0 这条线，t=∞ 处的指数衰减则无条件成立。如果你把 e^{-t} 换成别的衰减更慢的核函数，无穷远那一端就不再免费，收敛域会整体缩小。

2.2 解析延拓：拿递推关系当脚手架往外搭

定义只在右半平面活着，但伽马函数的实际使用几乎总要跨过这条线。负数、接近零的数，在统计和物理里出现得非常频繁。这时候就轮到解析延拓。

延拓的入口还是那条递推关系，只不过这次是反着用。把 Γ(z+1)=zΓ(z) 两边除以 z，得到

$$\Gamma(z)=\frac{\Gamma(z+1)}{z}$$

这个式子在 Re(z)>0 时当然成立，但它右边的 Γ(z+1) 的自变量比左边大 1，所以只要 z+1 落在收敛域里，右边就有意义。也就是说，任何满足 Re(z)>-1 的 z，只要不是 0，都能用这个式子算出一个值。

再推一步：Γ(z)=Γ(z+2)/[z(z+1)]，能把定义域推到 Re(z)>-2，去掉 0 和 -1。一般地，取足够大的整数 n 使得 Re(z+n+1)>0，就有

$$\Gamma(z)=\frac{\Gamma(z+n+1)}{z(z+1)(z+2)\cdots(z+n)}$$

右边是一个有定义的区域的函数除以若干线性因子。这就是解析延拓的全部操作：它不是另找一个函数来顶替，而是用原函数在已知区域的值，去定义新区域的值。全纯性由右边的表达式自动继承，因为在被挖掉的点之外，分母处处非零且解析。

说白了，递推关系在这里的角色不是“一条性质”，而是一台延拓机器。这个视角转换很关键：教材往往把递推公式和解析延拓分成两节讲，实际上它们是同一件事的正反两面。

2.3 负整数处的一阶极点与留数公式

延拓到负半轴之后，马上会遇到那些被挖掉的点：z=0, -1, -2, …。

以 z=-n 为例。看上面那个表达式，分子 Γ(z+n+1) 在 z=-n 处取值 Γ(1)=1，非零；分母 z(z+1)…(z+n) 里只有 (z+n) 这一项为零，其余因子在 z=-n 处都是非零常数，乘积为 (-n)(-n+1)…(-1)，也就是 (-1)^n·n!。

单零点在分母、非零值在分子，所以这是一阶极点。留数直接算：

$$\operatorname*{Res}{z=-n}\Gamma(z)=\lim{z\to -n}(z+n)\Gamma(z)=\frac{1}{(-1)^n n!}=\frac{(-1)^n}{n!}$$

代入几个值感受一下：z=0 处留数是 1，z=-1 处是 -1，z=-2 处是 1/2，z=-3 处是 -1/6。绝对值随 n 增大迅速衰减，说明极点附近的“爆炸强度”越来越弱，这和图像上负半轴振幅逐渐压向 0 的表现是一致的。

需要强调一点：写成 Γ(-n)=∞ 只是一种记号上的方便，严格说是不成立的写法。Γ 在那些点上根本没有有限值，也没有无穷值，它在那里就没有定义。习惯说法“伽马函数在负整数处取无穷”容易让人误以为它是个扩展实数系里的函数，实际上它是在那些点上有极点，性质完全不同。

顺带提一个有用的对偶视角：倒数 1/Γ(z) 是整函数，在 0,-1,-2,… 处各有一阶零点。这个事实在 Weierstrass 乘积表示下看得最清楚：

$$\frac{1}{\Gamma(z)}=z,e^{\gamma z}\prod_{n=1}^{\infty}\left(1+\frac{z}{n}\right)e^{-z/n}$$

其中 γ≈0.5772 是欧拉常数。这个式子一摆出来就能直接读出三件事：零点只在非正整数处；Γ 本身处处非零；以及 1/Γ 的增长是整函数里比较温和的那一类。做数值估计时想判断“Γ 会不会等于零”，从这个式子看答案是不可能。

3. 几条核心性质，各有各的有效边界

3.1 递推公式：既是最常用的性质，也是延拓的脚手架

Γ(z+1)=zΓ(z) 是全篇使用频率最高的一条。它有两个方向的用法。

正着用叫归约：把任意参数搬到 Re(z)∈(0,1] 这个小区间里去求值。工程实现里几乎所有的伽马函数库都是这个套路——先用递推把参数压进基准区间，再用 Lanczos 近似或者 Stirling 展开算出基准值，最后把沿途累积的因子乘回去。这个流程不是巧合，是因为近似公式在小区间上精度最好控制。

反着用叫延拓，就是上一节讲的那套。两个方向的数学依据一样，但工程意义完全不同：一个是省算力，一个是扩定义域。

高频陷阱在于整数点的取值偏移一位。Γ(n+1)=n!，所以 Γ(n)=(n-1)!，不是 n!。写代码的时候如果参数是从“序数”来的，脑子里想的是 n!，实际传进去的应该是 n+1。这个错误隐蔽性极强，因为数值本身看起来很正常，只是整体错了位置。我自己在小规模验证时习惯拿 Γ(5)=24、Γ(6)=120 这两个点当锚，先用它们把偏移量钉死，再往外推。

还有一个边界要注意：在 z=0 处不能用递推式从两个方向推，因为两边都是无穷，等式退化成 ∞=0·∞ 这种没有信息量的形式。

3.2 反射公式：负半轴取值的钥匙，也是精度陷阱

$$\Gamma(z)\Gamma(1-z)=\frac{\pi}{\sin(\pi z)}$$

这条公式的作用是把 z 和 1-z 配成一对，已知一半求另一半。比如算出 Γ(0.3) 之后，Γ(0.7) 可以直接用 π/[sin(0.3π)·Γ(0.3)] 拿到，不需要重新做一次数值积分。

但它有个非常实际的弱点：sin(πz) 在 z 接近整数时趋近于零，倒数关系会把输入的微小误差放大成灾难。举个例子，如果 z=1+10⁻⁸，那么 sin(πz)≈-π·10⁻⁸，反射公式右边大约是 -10⁸/π。此时如果 z 本身的浮点表示有相对误差 10⁻¹⁶，经过这个除法之后结果的相对误差会被放大到和 1 差不多的量级，也就是结果的有效数字全部丢光。

实际处理办法是：参数距离整数小于某个阈值（比如 10⁻⁸）时绕开反射公式，改用递推公式从邻近的安全点推过去。递推的误差传播是乘法型的，不会出现这种病态的放大。

另一个常见误解是“反射公式可以用来求 Γ(n) 的值”。不行。z 取整数时 sin(πz)=0，右边无定义，公式直接失效，这对应的正是负整数处的极点。

3.3 倍增公式：处理半整数点的批量求值

$$\Gamma(z)\Gamma!\left(z+\tfrac12\right)=2^{1-2z}\sqrt{\pi},\Gamma(2z)$$

这个也叫 Legendre 倍角公式，形式上像是反射公式的表亲，但用途不一样。反射公式处理的是 z 与 1-z 的对称，倍增公式处理的是步长为 1/2 的递推。

什么场景会用到？凡是出现“参数相差半个整数”的地方。统计学里的卡方分布、物理学里的相空间积分、组合数学里含半整数的 Beta 函数，都很容易撞上 Γ(n/2) 这类形式。有了倍增公式，一个步长为 2 的参数序列可以整体降到步长为 1，计算量直接减半。

验算一下 z=1 的情形：左边是 Γ(1)Γ(3/2)=1·(√π/2)=√π/2；右边是 2^{-1}·√π·Γ(2)=√π/2。对得上。做任何公式记忆之前，先代两个小整数验一遍，是最省事的防错手段。

3.4 对数凸性与 Bohr–Mollerup：唯一性的来源

回到开头那个问题：阶乘的连续推广凭什么只有一个？答案就在这里。

定理的陈述是：若 f 在 (0,∞) 上取正值、满足 f(1)=1、满足 f(x+1)=x·f(x)、并且 ln f 是凸函数，则 f 恒等于 Γ。

证明思路分两段。先证明 Γ 满足这三个条件——前两条前面已经推过，第三条的凸性来自积分定义上的 Hölder 不等式：把 Γ((1-λ)x+λy) 写成 ∫ t^{(1-λ)x+λy-1}e^{-t}dt，然后把 t^{...} 拆成两个幂的乘积，套 Hölder，就能得到那个不等式。再证明满足条件的函数唯一——用递推把参数夹在两个整数之间，配合凸性给出的上下界，两边同时收紧，函数值被迫等于 Γ。

这个定理的实际价值不在于会背结论，而在于它告诉你唯一性是有条件的。如果把“对数凸”这一条拿掉，能穿过同样整数点、满足同样递推的函数有无穷多个，其中最出名的一个反例是 Hadamard 构造的那类带振荡修正项的版本。所以下次看到“伽马函数是阶乘的唯一推广”这种说法，心里应该自动补上后半句：“在正性、归一化、递推和对数凸这四条同时成立的前提下”。

3.5 Digamma：取对数再求导，结构反而更清爽

$$\psi(z)=\frac{d}{dz}\ln\Gamma(z)=\frac{\Gamma'(z)}{\Gamma(z)}$$

这个函数叫双伽马函数，是伽马函数的对数导数。它出现的频率比想象中高——最大似然估计里对伽马分布、贝塔分布求导的时候，链式法则一定会把它拽出来。

几条常用关系值得记：

为什么工程上更愿意用 ψ 而不是直接对 Γ 求导？因为 Γ 的增长太快，直接做数值微分会在分母上炸掉；而 ln Γ 的增长是线性的（Stirling 公式给出 lnΓ(z)≈z ln z - z），数值上非常温顺。求导之后得到的 ψ 增长只有 ln z 量级，误差控制好得多。

再往下一阶是 ψ'(z)，叫三伽马函数，在统计里做 Fisher 信息矩阵时会碰到；更高阶统称多伽马函数。它们都保留了“取对数再求导”带来的数值稳定性。

4. 手算 Γ(1/2) 的完整路径

4.1 极坐标代换那一步为什么必须这么换

Γ(1/2)=√π 这个结论几乎人人都知道，但完整的推导链条里有一个关键转折，值得单独说清楚。

直接算：Γ(1/2)=∫₀^∞ t^{-1/2}e^{-t}dt。第一步换元，令 t=u²，则 dt=2u du，t^{-1/2}=1/u，代入得

$$\Gamma!\left(\tfrac12\right)=\int_0^\infty \frac{1}{u}e^{-u^2}\cdot 2u,du=2\int_0^\infty e^{-u^2}du$$

到这一步，问题转化成了高斯积分。麻烦在于 ∫₀^∞ e^{-u²}du 没有初等原函数，硬算是不可能的。这时候必须换策略：把积分平方，变成二重积分。

$$\left(\int_0^\infty e^{-u^2}du\right)^2=\int_0^\infty!!\int_0^\infty e^{-(x^2+y^2)}dx,dy$$

平方之后出现的 x²+y² 是关键信号——它在提示极坐标。换成 r、θ 之后

$$\int_0^{\pi/2}!!\int_0^\infty e^{-r^2}r,dr,d\theta=\frac{\pi}{2}\cdot\frac{1}{2}=\frac{\pi}{4}$$

所以高斯积分本身等于 √π/2，代回去得到 Γ(1/2)=√π。

整套推导里最关键的一步是令 t=u²。这个替换不是为了好看，而是因为只有把 e^{-t} 变成 e^{-u²}，平方之后才凑得出 x²+y² 这个极坐标的识别标志。如果换成别的替换，比如 t=v³，平方后得到的是 e^{-(v³+w³)}，形不成半径平方，走不下去。选替换的时候盯着“平方后能不能合并成一个半径”，方向就不会错。

4.2 从 Γ(1/2) 推向全部半整数点

拿到基准值之后，剩下的半整数点全靠递推往上推：

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

通式可以写成

$$\Gamma!\left(n+\tfrac12\right)=\frac{(2n)!}{4^n,n!}\sqrt{\pi},\qquad n=0,1,2,\dots$$

拿 n=4 验一遍：右边是 8!·√π/(4⁴·4!)=40320√π/6144=6.5625√π≈11.6317。用递推从 Γ(1/2) 一路乘上去：1.7725×0.5×1.5×2.5×3.5=11.6317，对得上。手算验证一遍很有必要，因为通式里的 (2n)!/(4ⁿn!) 这种结构很容易在记忆时把因子 4 写成 2。

顺带注意，这张表不是单调的：Γ(1/2) 到 Γ(3/2) 是下降的，之后才转为上升。这对应前面说的最小值点位于 x≈1.4616 处。

4.3 负半轴的取值与符号交替规律

负半轴没有独立的公式，全靠递推反推。反推的写法是 Γ(z)=Γ(z+1)/z。

$$\Gamma!\left(-\tfrac12\right)=\frac{\Gamma(1/2)}{-1/2}=-2\sqrt{\pi}\approx -3.5449$$

$$\Gamma!\left(-\tfrac32\right)=\frac{\Gamma(-1/2)}{-3/2}=\frac{4\sqrt{\pi}}{3}\approx 2.3633$$

$$\Gamma!\left(-\tfrac52\right)=\frac{\Gamma(-3/2)}{-5/2}=-\frac{8\sqrt{\pi}}{15}\approx -0.9453$$

看符号：负、正、负，交替出现。规律是 Γ 在区间 (-n-1, -n) 上的符号为 (-1)^{n+1}。原因很好理解——每往左推一格就多除一个负数，符号翻转一次。

绝对值则在递减：3.5449 → 2.3633 → 0.9453，越往左越接近 0。这和图像上负半轴的振荡逐渐衰减是一致的，也解释了为什么远负半轴的 Γ 值在数值上基本可以忽略。

这里有个符号陷阱值得单独提醒：用反射公式算负半轴时，sin(πz) 的符号由 z 的小数部分决定，很容易搞反。比如算 Γ(-0.5)，反射公式给的是 π/[sin(-π/2)·Γ(1.5)]=π/[-1·0.8862]=-3.5449，负