Implementation:ClickHouse ClickHouse Musl Lgammal
| Knowledge Sources | |
|---|---|
| Domains | Mathematics, Portability |
| Last Updated | 2026-02-08 00:00 GMT |
Overview
Long double precision implementation of the logarithm of the gamma function with sign tracking.
Description
This implementation computes `lgammal_r(x, sg)` which returns the natural logarithm of the absolute value of the gamma function and sets the sign via a pointer. It uses rational polynomial approximations in different regions: argument reduction for small values, polynomial approximations around the minimum, rational approximations for the primary interval, asymptotic series for large values, and reflection formula for negative arguments. The code is derived from OpenBSD and Sun Microsystems implementations optimized for 64-bit and 80-bit long double formats.
Usage
Use this function when computing gamma-related statistical or probabilistic calculations requiring extended precision, particularly for values where standard double precision is insufficient or when working with very large or very small arguments.
Code Reference
Source Location
- Repository: ClickHouse
- File: base/glibc-compatibility/musl/lgammal.c
- Lines: 1-330
Signature
long double lgammal_r(long double x, int *sg);
Import
#include <math.h>
// Or directly:
// #include "base/glibc-compatibility/musl/lgammal.c"
I/O Contract
| Input Range | Output | Sign Parameter | Special Cases |
|---|---|---|---|
| x > 0 | log(Gamma(x)) | *sg = 1 | Regular computation |
| x < 0, non-integer | log(abs(Gamma(x))) | *sg = ±1 based on sin(πx) | Uses reflection formula |
| x = 0 | +∞ | *sg = ±1 | Pole at zero |
| x = negative integer | +∞ | *sg = ±1 | Pole at negative integers |
| x = 1 or 2 | 0 | *sg = 1 | Gamma(1) = Gamma(2) = 1 |
| x = NaN | NaN | *sg unchanged | Propagates NaN |
| x = ±∞ | +∞ | *sg = 1 | Gamma increases without bound |
Usage Examples
// Basic usage - compute log(Gamma(x))
int sign;
long double x = 5.5L;
long double result = lgammal_r(x, &sign);
// result ≈ 3.957, sign = 1
long double gamma_x = expl(result) * sign;
// Handle negative values
x = -3.7L;
result = lgammal_r(x, &sign);
// Computes log(|Gamma(-3.7)|) and sets sign appropriately
// Statistical computation - Beta function
long double log_beta(long double a, long double b) {
int sign_a, sign_b, sign_sum;
long double lgamma_a = lgammal_r(a, &sign_a);
long double lgamma_b = lgammal_r(b, &sign_b);
long double lgamma_ab = lgammal_r(a + b, &sign_sum);
return lgamma_a + lgamma_b - lgamma_ab;
}
// Binomial coefficient in log space
long double log_binom(long double n, long double k) {
int s1, s2, s3;
long double result = lgammal_r(n + 1.0L, &s1)
- lgammal_r(k + 1.0L, &s2)
- lgammal_r(n - k + 1.0L, &s3);
return result;
}
// Check for overflow before exponentiating
x = 100.0L;
result = lgammal_r(x, &sign);
if (result < LDBL_MAX_EXP * logl(2.0L)) {
long double gamma_val = sign * expl(result);
printf("Gamma(%Lf) = %Le\n", x, gamma_val);
} else {
printf("Gamma(%Lf) would overflow\n", x);
}