逼近反三角函数

Question

我需要在仅有正弦函数、余弦函数和基本固定小数点算术（没有浮点数）的环境中实现asin、acos和atan。

我已经有了一个相当好的平方根函数。

我能否使用这些工具来实现相对高效的反三角函数？

我不需要太高的精度（浮点数精度非常有限），基本逼近即可。

我已经最终决定采用表格查找方法，但我想知道是否有更简洁的选择（不需要编写几百行代码来实现基本的数学计算）。

编辑：

为了澄清事情：我需要在每秒35帧的情况下运行该函数数百次。

- Matěj Zábský

---

可能是三角函数如何工作？的重复问题。 - Jason S

提出的重复问题更多地涉及三角函数的工作原理（就像它的标题一样）。这是关于反三角函数的问题。 - Teepeemm

Answer 1

在固定点环境（S15.16）中，我成功地使用了CORDIC算法（请参阅维基百科的一般描述）来计算atan2（y，x），然后使用已知的函数恒等式涉及平方根导出asin（）和acos（）。

asin(x) = atan2 (x, sqrt ((1.0 + x) * (1.0 - x)))
acos(x) = atan2 (sqrt ((1.0 + x) * (1.0 - x)), x)

事实证明，在双倍精度（double）上找到有用的CORDIC迭代描述来计算 atan2() 比我想象的更难。下面的网站似乎包含了足够详细的描述，并讨论了两种备选方法，即多项式逼近和查找表：

http://ch.mathworks.com/examples/matlab-fixed-point-designer/615-calculate-fixed-point-arctangent

- njuffa

---

从维基百科来看，CORDIC甚至不使用三角函数（很棒！）；我想你所做的是一次搜索；考虑到sin()，cos()，牛顿-拉弗森或类似方法可能更好？（需要较少的迭代，尽管迭代的成本会有所不同。） - petrelharp

我建议研究CORDIC的原因是它只需要固定点算术。CORDIC最常见的用途可能是实现sin / cos，这也是我第一次了解它的方式（在1987年）。但是，还有很多其他函数也可以使用CORDIC计算，例如atan2。由于我没有任何关于使用CORDIC计算atan2的代码，因此我尝试找到一个具有足够详细信息的网站，以便有人可以基于它进行实现。我上面发布的链接是我在几分钟内通过搜索引擎找到的最好的页面。 - njuffa

Answer 2

以下代码应很容易调整为定点数格式。它使用分式逼近来计算标准化到 [0 1) 区间的反正切值（可以乘以 Pi/2 来获得真实的反正切值）。然后，您可以使用众所周知的公式从反正切值中获取反正弦/反余弦值。

normalized_atan(x) ~ (b x + x^2) / (1 + 2 b x + x^2)
 
where b = 0.596227

最大误差为0.1620度。

#include <stdint.h>
#include <math.h>

// Approximates atan(x) normalized to the [-1,1] range
// with a maximum error of 0.1620 degrees.

float norm_atan( float x )
{
    static const uint32_t sign_mask = 0x80000000;
    static const float b = 0.596227f;

    // Extract the sign bit
    uint32_t ux_s  = sign_mask & (uint32_t &)x;

    // Calculate the arctangent in the first quadrant
    float bx_a = ::fabs( b * x );
    float num = bx_a + x * x;
    float atan_1q = num / ( 1.f + bx_a + num );

    // Restore the sign bit
    uint32_t atan_2q = ux_s | (uint32_t &)atan_1q;
    return (float &)atan_2q;
}

// Approximates atan2(y, x) normalized to the [0,4) range
// with a maximum error of 0.1620 degrees

float norm_atan2( float y, float x )
{
    static const uint32_t sign_mask = 0x80000000;
    static const float b = 0.596227f;

    // Extract the sign bits
    uint32_t ux_s  = sign_mask & (uint32_t &)x;
    uint32_t uy_s  = sign_mask & (uint32_t &)y;

    // Determine the quadrant offset
    float q = (float)( ( ~ux_s & uy_s ) >> 29 | ux_s >> 30 ); 

    // Calculate the arctangent in the first quadrant
    float bxy_a = ::fabs( b * x * y );
    float num = bxy_a + y * y;
    float atan_1q =  num / ( x * x + bxy_a + num );

    // Translate it to the proper quadrant
    uint32_t uatan_2q = (ux_s ^ uy_s) | (uint32_t &)atan_1q;
    return q + (float &)uatan_2q;
} 

如果您需要更高的精度，可以使用第三阶有理函数：

normalized_atan(x) ~ ( c x + x^2 + x^3) / ( 1 + (c + 1) x + (c + 1) x^2 + x^3)
 
where c = (1 + sqrt(17)) / 8

这个技术的最大近似误差为0.00811度。

- rcor

Answer 3

您可能需要使用近似值：使用无限级数，直到解决方案足够接近为止。

例如：
arcsin(z) = Sigma((2n!)/((2^2n)*(n!)^2)*((z^(2n+1))/(2n+1)))其中n在[0,无穷大)。

- amit

Answer 4

http://en.wikipedia.org/wiki/Inverse_trigonometric_functions#Expression_as_definite_integrals

您可以使用您的平方根函数进行数值积分，通过无限级数逼近：

- nulvinge

Answer 5

您是否需要对 arcsin(x) 函数进行高精度计算？如果不需要，您可以在N个节点上计算 arcsin 并将值保存在内存中。我建议使用线性逼近法。如果 x = A*x_(N) + (1-A)*x_(N+1)，那么 x = A*arcsin(x_(N)) + (1-A)*arcsin(x_(N+1))，其中 arcsin(x_(N)) 是已知的。

- Eugeny89

---

是的，这就是我在原帖中提到的表查找。我不认为有理由在运行时计算它，我会将这些值直接嵌入程序中，因此实际的asin计算只是两个值之间的插值。 - Matěj Zábský

Answer 6

我在这个类似的问题中提交我的答案。

nVidia拥有一些我用于自己的工具，例如：acos、asin、atan2等等...

这些算法产生足够精确的结果。以下是一个直接使用它们代码的Python示例：

import math
def nVidia_acos(x):
    negate = float(x<0)
    x=abs(x)
    ret = -0.0187293
    ret = ret * x
    ret = ret + 0.0742610
    ret = ret * x
    ret = ret - 0.2121144
    ret = ret * x
    ret = ret + 1.5707288
    ret = ret * math.sqrt(1.0-x)
    ret = ret - 2 * negate * ret
    return negate * 3.14159265358979 + ret

以下是用于比较的结果：

nVidia_acos(0.5)  result: 1.0471513828611643
math.acos(0.5)    result: 1.0471975511965976

这很接近了！将结果乘以57.29577951即可得到以度为单位的结果，这也是他们“度数”公式的一部分。

- Fnord

Answer 7

也许可以采用类似牛顿-拉普森法的智能暴力方法。

因此，为了解决asin()问题，您可以在sin()上使用最陡下降法。

- aero

---

你可以从一个小的查找表开始选择起点，以加快计算速度。 - Karoly Horvath

Answer 8

使用多项式逼近。最小二乘拟合是最简单的（Microsoft Excel有它），而Chebyshev逼近更准确。

这个问题以前已经讨论过：三角函数是如何工作的？

- Jason S

Answer 9

只有连续函数可以用多项式逼近。而arcsin(x)在点x=1处不连续，arccos(x)也是如此。但是，在这种情况下，将范围缩小到区间1，sqrt(1/2)可以避免这种情况。我们有arcsin(x)=pi/2-arccos(x)，arccos(x)=pi/2-arcsin(x)。您可以使用matlab进行极小值逼近。仅在范围[0，sqrt(1/2)]内进行逼近（如果请求该arcsin的角度大于sqrt(1/2)，则找到cos(x)）。arctangent函数仅适用于x<1。arctan(x)=pi/2-arctan(1/x)。

- florin