如何实现arctan？

问题描述 投票：14回答：3

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

algorithm floating-point trigonometry bit processor
3个回答
7
投票

三角函数确实具有难看的实现，这些实现很hacky，并且有些麻烦。我认为在这里很难找到可以解释实际使用算法的人。

这里是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不必完成所有这些工作。

19
投票

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

  1. 用于将整个输入域映射到窄间隔的参数减少
  2. 窄区间（一次近似区间）上的核心近似的计算
  3. 基于自变量约简的中间结果的扩展以产生最终结果

参数归约通常基于众所周知的三角身份，可以在各种标准参考（例如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 n

），i = 1 ... 2 n，该表允许使用任意窄的一次近似间隔以牺牲额外的表存储为代价。这是时空之间的经典编程权衡。

核心区间上的近似值通常是足够程度的极小多项式近似值。由于浮点除法的高昂成本，有理逼近通常在现代硬件上没有竞争力，而且由于两个多项式的计算加上除法带来的误差，有理逼近也遭受附加的数值误差。

用于最小极大多项式逼近的系数通常使用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 ul的错误结果，这在各种情况下都足够了。我使用Remez算法的简单实现创建了一个自定义的minimax逼近，该算法对所有中间步骤都使用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);
}
6
投票

摘要：很难。另外，有时会在SO周围徘徊的Eric Postpischil和Stephen Canon也非常擅长。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.