2024.11.26

· · 个人记录

今天上课突然发现有人作业写了 std::tgamma 这个函数,去查了一下,好像可以 O(1) 计算 n-1 的阶乘,C++11 引入的,ICPC 可用

在 GCC 中,std::tgamma 的实现通常依赖于底层的数学库(libm)。在 GNU libm 库中,tgamma 函数的实现使用了一些经典的方法,如 Lanczos 逼近法反射公式。以下是具体的实现细节概述:

关于 tgamma

1. 基本结构

tgamma 函数的实现可以分为几个部分:

  1. 特殊值处理

    • 检测输入是否为负整数或零(伽马函数在这些点未定义)。
    • 对极端值(如 NaNInf)返回适当的结果。
  2. 正数的伽马值计算

    • 对于大于一定范围的正数,直接使用 斯特林近似
    • 对于中等范围的正数,使用 Lanczos 逼近法
  3. 负数和小数的处理

    • 使用 反射公式: [ \Gamma(1-x)\Gamma(x) = \frac{\pi}{\sin(\pi x)} ] 将负数或小数转化为正数计算。

2. 特殊值处理

在伽马函数中,以下特殊情况需要单独处理:

这部分通常用简单的分支语句处理。

3. 中等范围正数:Lanczos 逼近法

在 GNU 的数学库中,tgamma 的核心计算部分大多基于 Lanczos 逼近法。这是一个快速且高精度的伽马函数近似算法。

Lanczos 逼近公式 Lanczos 方法通过一个展开式逼近伽马函数: [ \Gamma(x) \approx \sqrt{2\pi} \cdot \left(x + g - \frac{1}{2}\right)^{x-\frac{1}{2}} \cdot e^{-(x+g-\frac{1}{2})} \cdot A(x) ] 其中:

这种方法在保证精度的同时计算复杂度较低。

4. 大范围正数:斯特林近似

对于较大的 ( x ),伽马函数值增长迅速,因此采用 斯特林近似: [ \Gamma(x) \sim \sqrt{2\pi} x^{x-\frac{1}{2}} e^{-x} ] 该公式适用于 ( x > 10 ) 或更大时,计算简单且精度足够。

5. 负数和小数:反射公式

对于负数输入 ( x ) 或小数输入 ( x < 1 ),使用反射公式: [ \Gamma(x) = \frac{\pi}{\sin(\pi x) \cdot \Gamma(1-x)} ] 这将计算转化为对 ( 1-x ) 的正数处理,避免直接处理负数。

6. GCC 的具体实现(简化伪代码)

以下是 tgammalibm 中的一个典型实现结构:

double tgamma(double x) {
    if (x <= 0 && x == (int)x) {
        // 负整数或零,伽马函数未定义
        errno = EDOM;
        return HUGE_VAL;
    }

    if (x < 0) {
        // 负数:使用反射公式
        double sin_pi_x = sin(M_PI * x);
        if (sin_pi_x == 0) {
            errno = ERANGE;
            return HUGE_VAL;
        }
        return M_PI / (sin_pi_x * tgamma(1 - x));
    }

    if (x > THRESHOLD) {
        // 大数:斯特林近似
        return stirling_approximation(x);
    }

    // 中等范围:Lanczos 逼近
    return lanczos_approximation(x);
}

7. 性能优化

总结

GCC 中的 std::tgamma 主要依赖 Lanczos 逼近法反射公式,并结合 斯特林近似 来高效处理伽马函数的不同输入范围。它的实现经过高度优化,既保证了高精度又能适应各种数值计算场景。如果需要具体实现代码,可以查阅 GNU libm 的源码,通常位于 glibc 仓库的 math 模块中。

关于 midpoint

GPT 结束,然后刚想去逛逛菜园,发现他们在讨论 std::midpoint,更好的计算区间中点,向第一个参数靠拢,可以避免溢出以及精度问题。

更逆天的是还能 midpoint(r,l)

关于 bits/////////stdc++.h

帮舍友调代码的时候发现了一个骚操作,原来编译器早就想到了可能的错误。

GPT FS 的结果

后续测试了 Linux 下的编译,无论多少个 / 都能过编,反向的话只有两条的情况 \\\\ 可以。

讨论区 可能有用的东西

关于重载[]运算符

芝士++ 大量空格可以用来卡 getchar() 快读

快读快写为什么会负优化?