arctan是如何实施的？

Question

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

Answer 1

x86处理器中FPATAN指令的实现通常是专有的.为了计算arctan或其他(逆)三角函数,常见算法遵循三个步骤:

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 ⁿ),i = 1 ... 2 ^n的表,它允许使用任意窄的主近似间隔,但需要额外的表存储.这是空间和时间之间的经典编程权衡.

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

最小极大多项式近似的系数通常使用Remez算法(http://en.wikipedia.org/wiki/Remez_algorithm)计算.像Maple和Mathematica这样的工具具有内置工具来计算这种近似值.通过确保所有系数都是可精确表示的机器编号,可以提高多项式近似的精度.我所知道的唯一具有内置功能的工具是Sollya(http://sollya.gforge.inria.fr/),它提供了一个fpminimax()功能.

多项式的评估通常使用Horner的方案(http://en.wikipedia.org/wiki/Horner%27s_method),该方案有效且准确,或者是Estrin方案的混合(http://en.wikipedia.org/wiki/ Estrin%27s_scheme)和Horner's.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++最好的,我知道,它只能作为一个专有扩展)是指定浮点常数作为直接表达底层的位模式,有效地避免了复杂的转换十六进制文字.

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

Answer 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不必完成所有这些工作.

Answer 3

总结:这很难.此外,Eric Postpischil和Stephen Canon,他们有时会在SO身边徘徊,他们非常擅长.

许多特殊功能的通常方法如下:

• 处理NaNs,无穷大和签名零作为特例.
• 如果数字太大而结果转向M_PI,则返回M_PI.调用此阈值M.
• 如果存在任何类型的参数减少标识,请使用它将参数置于更好的范围内.(这可能很棘手:对于sin和cos,这意味着你选择2pi 的精确值的倍数,以便你落在正确的范围内.)
• 分解[0,M)成有限的间隔.在每个区间使用相当高阶的arctan的Chebyshev近似.(这是离线完成的,它通常是你在这些实现中看到的所有魔法数字的来源.而且,人们可以使用Remez的交换算法略微收紧Chebyshev近似值,但我不知道这有什么帮助很多.)
• 找出参数所在的区间(使用ifs和stuff或仅使用表索引的技巧),并在该区间上评估Chebyshev序列.

这里特别需要一些属性:

• 该arctan实施应该是单调的; 也就是说,如果x < y,那么arctan(x) <= arctan(y).
• 该arctan实施应始终为1个ULP正确答案的内返回一个答案.请注意,这是一个相对误差范围.

评估Chebyshev系列并不是完全直截了当,以便这两个属性成立.两个doubles用于表示单个值的不同部分的技巧在这里很常见.然后可能会有一些案例表明实施是单调的.此外,接近零,泰勒近似arctan而不是切比雪夫近似---你是在一个相对误差范围之后并使用霍纳规则评估该系列应该起作用.

如果你正在寻找一个atan阅读的实现,fdlibm似乎不如目前在glibc中的那个讨厌.的参数减少似乎可以基于TRIG身份tan(a+b) = (tan(a) + tan(b)) / (1 - tan(a) tan(b)),使用0.5,1或1.5用于tan(a)为适当.