algorithm - arctan 是如何实现的？

Question

该库的许多实现深入到所有弧函数的 FPATAN 指令。 FPATAN 是如何实现的？假设我们有 1 位符号，M 位尾数和 N 位指数，得到这个数反正切的算法是什么？应该有这样的算法，因为 FPU 做到了。

Answer 1

FPATAN 指令在 x86 处理器中的实现通常是专有的。要计算反正切或其他(逆)三角函数，常用算法遵循三个步骤:

用于将完整输入域映射到窄区间的参数缩减

在窄区间(主近似区间)上计算核心近似值

基于参数缩减的中间结果扩展以产生最终结果

参数减少通常基于众所周知的三角恒等式，可以在各种标准引用资料中查找，例如 MathWorld ( http://mathworld.wolfram.com/InverseTangent.html )。对于 arctan 的计算，常用的恒等式是

反正切 (-x) = -反正切 (x)

arctan (1/x) = 0.5 * pi - arctan(x) [x > 0]

arctan (x) = arctan(c) + arctan((x - c)/(1 + x*c))

请注意，最后一个恒等式有助于构建值 arctan(i/2n), i = 1...2n 的表，这允许使用任意窄的主近似区间，但会增加表存储空间。这是空间和时间之间的经典编程权衡。

核心区间上的近似通常是具有足够度数的极小极大多项式近似。由于浮点除法的高成本，有理近似在现代硬件上通常没有竞争力，并且由于两个多项式的计算加上除法造成的误差，还存在额外的数值误差。

极小极大多项式近似的系数通常使用 Remez 算法 ( http://en.wikipedia.org/wiki/Remez_algorithm) 计算。 Maple 和 Mathematica 等工具具有计算此类近似值的内置工具。通过确保所有系数都是可准确表示的机器数，可以提高多项式近似的准确性。我所知道的唯一具有内置工具的工具是 Sollya ( http://sollya.gforge.inria.fr/ )，它提供了一个 fpminimax()功能。

多项式的计算通常使用高效且准确的霍纳方案 ( http://en.wikipedia.org/wiki/Horner%27s_method)，或埃斯特林方案 ( http://en.wikipedia.org/wiki/Estrin%27s_scheme) 和霍纳方案的混合。 Estrin 的方案允许人们充分利用超标量处理器提供的指令级并行性，对整体指令数的影响很小，并且经常(但并非总是)对准确性产生良性影响。

由于舍入步骤的数量减少并提供了一些防止减法抵消的保护，因此使用 FMA(融合乘法加法)可以提高任一评估方案的准确性和性能。 FMA 可以在许多处理器上找到，包括 GPU 和最新的 x86 CPU。在标准 C 和标准 C++ 中，FMA 操作公开为 fma()标准库函数，但是它需要在不提供硬件支持的平台上进行模拟，这使得它在这些平台上运行缓慢。

从编程的角度来看，在将近似和参数减少所需的浮点常量从文本表示转换为机器表示时，人们希望避免转换错误的风险。 ASCII 到浮点的转换例程因包含棘手的错误(例如 http://www.exploringbinary.com/php-hangs-on-numeric-value-2-2250738585072011e-308/ )而臭名昭著。标准 C 提供的一种机制( 不是 我所知道的 C++ 最好的，它仅作为专有扩展可用)是将浮点常量指定为直接表达底层位模式的十六进制文字，有效避免复杂的转换。

下面是用于计算 double arctan() 的 C 代码，它演示了上面提到的许多设计原则和技术。这种快速构建的代码缺乏其他答案中指出的实现的复杂性，但应提供小于 2 ulp 错误的结果，这在各种情况下可能就足够了。我使用 Remez 算法的简单实现创建了一个自定义的极小极大近似，该算法对所有中间步骤使用 1024 位浮点算法。我希望使用 Sollya 或类似的工具会产生数值上优越的近似值。

double my_atan (double x)
{
    double a, z, p, r, s, q, o;
    /* argument reduction: 
       arctan (-x) = -arctan(x); 
       arctan (1/x) = 1/2 * pi - arctan (x), when x > 0
    */
    z = fabs (x);
    a = (z > 1.0) ? 1.0 / z : z;
    /* evaluate minimax polynomial approximation */
    s = a * a; // a**2
    q = s * s; // a**4
    o = q * q; // a**8
    /* use Estrin's scheme for low-order terms */
    p = fma (fma (fma (-0x1.53e1d2a25ff34p-16, s, 0x1.d3b63dbb65af4p-13), q,
                  fma (-0x1.312788dde0801p-10, s, 0x1.f9690c82492dbp-9)), o,
             fma (fma (-0x1.2cf5aabc7cef3p-7, s, 0x1.162b0b2a3bfcep-6), q, 
                  fma (-0x1.a7256feb6fc5cp-6, s, 0x1.171560ce4a483p-5)));
    /* use Horner's scheme for high-order terms */
    p = fma (fma (fma (fma (fma (fma (fma (fma (fma (fma (fma (fma (p, s,
        -0x1.4f44d841450e1p-5), s,
         0x1.7ee3d3f36bb94p-5), s, 
        -0x1.ad32ae04a9fd1p-5), s,  
         0x1.e17813d66954fp-5), s, 
        -0x1.11089ca9a5bcdp-4), s,  
         0x1.3b12b2db51738p-4), s,
        -0x1.745d022f8dc5cp-4), s,
         0x1.c71c709dfe927p-4), s,
        -0x1.2492491fa1744p-3), s,
         0x1.99999999840d2p-3), s,
        -0x1.555555555544cp-2) * s, a, a);
    /* back substitution based on argument reduction */
    r = (z > 1.0) ? (0x1.921fb54442d18p+0 - p) : p;
    return copysign (r, x);
}