arctan 是如何实现的？答案

【问题标题】：How is arctan implemented?arctan 是如何实现的？
【发布时间】：2014-05-27 16:44:56
【问题描述】：

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

【问题讨论】：

标签： algorithm floating-point trigonometry bit processor

【解决方案1】：

在 x86 处理器中 FPATAN 指令的实现通常是专有的。要计算 arctan 或其他（反）三角函数，常用算法遵循三个步骤：

将整个输入域映射到窄区间的参数缩减
在窄区间（初级逼近区间）上计算核心逼近
根据参数缩减扩展中间结果以生成最终结果

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

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

请注意，最后一个恒等式有助于构建值表 arctan(i/2ⁿ), i = 1...2ⁿ，其中允许以额外的表存储为代价使用任意窄的主近似间隔。这是空间和时间之间的经典编程权衡。

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

极小多项式近似的系数通常使用 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)，或 Estrin 方案 (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 提供的一种机制（not C++ 我知道，它只能作为专有扩展使用）是将浮点常量指定为直接表示底层位模式的十六进制文字，有效避免复杂的转换。

下面是计算双精度 arctan() 的 C 代码，它演示了上面提到的许多设计原则和技术。这种快速构建的代码缺乏其他答案中指出的实现的复杂性，但应该提供少于 2 ulps 错误的结果，这在各种情况下可能就足够了。我使用 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);
}

【讨论】：

出于好奇，在任何情况下，使用弧度进行三角函数计算比使用整数个细分可以实现更好的精度吗？当然，使用以度、象限或整圆为单位测量的角度会更容易和更精确地减少模量。我知道为什么弧度在微积分中很有用，但是不能精确表示一个完整圆的角度单位的数量似乎很恶心。
一些平台提供sinpi() 和cospi() 函数，它们接受的参数是pi 的倍数，这使得参数减少变得容易。否则，sin、cos、tan 的精确参数归约是困难的，并且无论使用弧度还是度数，本质上都需要多精度中间计算。规范参考是：Mary H. Payne 和 Robert N. Hanek，三角函数的弧度缩减，ACM SIGNUM 通讯，卷。 18，没有。 1，1983 年 1 月，第 19 - 24 页
度参数约简的配套论文是：Mary H. Payne 和 Robert N. Hanek，三角函数的度数约简，ACM SIGNUM 通讯，卷。 18. 没有。 2，1983 年 4 月，第 18 - 19 页
为什么在度数情况下需要多精度归约？可以肯定的是，在倍数的情况下更容易，但是 fpmod(x, 360.0) 被指定为对所有 x 值绝对精确，不是吗？顺便说一句，我不确定在使用弧度时超精确的参数缩减有多大用处。如果尝试使用Math.Sin(x*2.0*Math.Pi) 计算 sin(2πx)，则以 2.0*Math.Pi 为模执行参数缩减比以 2π 为模执行的结果会更准确。
@chux 我同意按程度减少三角函数参数很容易。不幸的是，当一个人说错话时，没有办法纠正评论（除了在宽限期内）。不过，我建议使用remquo (angle,90.0) 而不是fmod()。

【解决方案2】：

三角函数确实有非常丑陋的实现，它们很笨拙并且做了大量的摆弄。我认为在这里很难找到能够解释实际使用的算法的人。

这是一个 atan2 实现：https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/e_atan2.c;h=a287ca6656b210c77367eec3c46d72f18476d61d;hb=HEAD

编辑：实际上我找到了这个：http://www.netlib.org/fdlibm/e_atan2.c，它更容易理解，但可能因此更慢（？）。

FPU 在某些电路中完成所有这些工作，因此 CPU 不必完成所有这些工作。

【讨论】：

非常感谢。在第一个链接上，它还包括 mpatan.h 和 mpatan.c，其中有 atan 的实现——正是我想要的。
并非所有 FPU 都在硬件中执行此操作。可能有些架构没有三角指令。 SSE 也不支持三角函数，所以 MSVC 2013 在矢量化代码时必须实现一个软件
x86 CPU 中的 FPATAN 指令通常通过微码实现，即存储在处理器内部 ROM 中的小程序。虽然此类程序可能使用可见 ISA 中不可用的特殊操作，但通常不涉及特殊电路。
second implementation of atan2 要短很多，因为它使用了atan。

【解决方案3】：

总结：这很难。此外，有时在 SO 周围闲逛的 Eric Postpischil 和 Stephen Canon 也非常擅长。

许多特殊功能的常用方法如下：

将 NaN、无穷大和有符号零作为特殊情况处理。
如果数字太大以至于结果舍入为M_PI，则返回M_PI。将此阈值称为M。
如果有任何类型的参数缩减标识，请使用它来将参数带入更好的范围。（这可能很棘手：对于 sin 和 cos，这意味着您选择 2pi 的精确值的倍数，以便落在正确的范围内.)
将[0,M) 分解为有限多个区间。在每个间隔上使用Chebyshev approximation 到相当高阶的反正切。（这是离线完成的，它通常是您在这些实现中看到的所有幻数的来源。此外，可以使用 Remez 的交换算法稍微收紧切比雪夫近似值，但我不知道这有什么帮助.)
找出参数在哪个区间（使用ifs 之类的东西，或者只是使用表索引的技巧），并在该区间评估切比雪夫级数。

这里有一些特别需要的属性：

arctan 实现应该是单调的；也就是说，如果x < y，那么arctan(x) <= arctan(y)。
arctan 实现应始终在正确答案的 1 ulp 内返回答案。请注意，这是一个相对误差范围。

评估切比雪夫级数以使这两个属性成立并不完全简单。使用两个doubles 来表示单个值的不同部分的技巧在这里很常见。然后可能有一些案例表明实现是单调的。此外，在接近零的情况下，arctan 的泰勒近似值而不是切比雪夫近似值——你在一个相对误差范围之后，使用霍纳规则评估系列应该可以工作。

如果您正在寻找要阅读的atan 实现，fdlibm 似乎没有 glibc 中的当前那么讨厌。参数缩减似乎基于三角标识tan(a+b) = (tan(a) + tan(b)) / (1 - tan(a) tan(b))，酌情使用0.5、1 或1.5 表示tan(a)。

【讨论】：

既然我们正在讨论这个主题，我也许应该在另一个问题中问这个问题，使用 Padé 逼近而不是多项式逼近的一个很好的理由是当要逼近的函数（例如反正切）趋于朝着 +/-inf 的有限极限。显然，大于 1 的多项式近似在那里永远不会有任何好处。现在我的问题是，因为无论如何我们都在进行参数缩减，并且近似值只用于，比如 [0 ... 0.5]，那么上述原因（我听说过的唯一一个）应该没那么重要，应该吗？
@PascalCuoq：我希望 k 度的切比雪夫近似和总度（分子度 + 分母度） k 的 Pade-Chebyshev 近似在近似表现良好的函数方面大致相同在一个紧凑的区间。在没有这样的参数减少方案的情况下，我猜你需要正确地获得学位的差异。（我只需要编写特殊函数的低质量实现，因此在某些情况下可能有更微妙的理由使用有理逼近而不是多项式逼近——我不知道。）
有理近似很少有竞争力。浮点除法比 FADD、FMUL 或 FMA 昂贵得多。此外，您必须处理两个多项式的误差加上除法的误差。在大多数情况下，您需要直接多项式或表加多项式。就多项式而言，您可能希望针对目标精度优化系数，例如Sollya 的 fpminimax() 函数提供的近似值。如果 FMA 可用，它将有助于保持较小的评估误差。 Estrin 的方案有助于提高超标量架构的性能。