lgamma

pub fn lgamma(comptime T: type, x: T) T

File

Code

pub fn lgamma(comptime T: type, x: T) T {
    if (T != f32 and T != f64) {
        @compileError("gamma not implemented for " ++ @typeName(T));
    }
    // common integer case first
    if (x == @trunc(x)) {
        // lgamma(-inf)  = +inf
        // lgamma(n)     = +inf for negative integers
        // lgamma(+-0.0) = +inf
        if (x <= 0) {
            return std.math.inf(T);
        }
        // lgamma(1) = +0.0
        // lgamma(2) = +0.0
        if (x < integer_result_table.len) {
            const i = @as(u8, @intFromFloat(x));
            return @log(@as(T, @floatCast(integer_result_table[i])));
        }
        // lgamma(+inf) = +inf
        if (std.math.isPositiveInf(x)) {
            return x;
        }
    }

    const abs = @abs(x);
    // perfect precision here
    if (abs < 0x1p-54) {
        return -@log(abs);
    }
    // obvious approach when overflow is not a problem
    const upper_bound = if (T == f64) 128 else 26;
    if (abs < upper_bound) {
        return @log(@abs(gamma(T, x)));
    }

    const log_base = @log(abs + lanczos_minus_half) - 1;
    const exponent = abs - 0.5;
    const log_series = @log(series(T, abs));
    const initial = exponent * log_base + log_series - lanczos;

    // use reflection formula for negatives
    if (x < 0) {
        const reflected = std.math.pi / (abs * sinpi(T, abs));
        return @log(@abs(reflected)) - initial;
    }
    return initial;
}