快速 Arc Cos 算法？答案

【问题标题】：Fast Arc Cos algorithm?快速 Arc Cos 算法？
【发布时间】：2011-03-23 18:40:11
【问题描述】：

我有自己的，非常快的 cos 函数：

float sine(float x)
{
    const float B = 4/pi;
    const float C = -4/(pi*pi);

    float y = B * x + C * x * abs(x);

    //  const float Q = 0.775;
    const float P = 0.225;

    y = P * (y * abs(y) - y) + y;   // Q * y + P * y * abs(y)


    return y;
}

float cosine(float x)
{
    return sine(x + (pi / 2));
}

但是现在当我分析时，我看到 acos() 正在杀死处理器。我不需要非常精确的。什么是计算acos（x）的快速方法谢谢。

【问题讨论】：

您非常快速的函数在 [-pi,pi] 中的平均误差为 16%，并且在该区间之外完全无法使用。在我的系统上，来自math.h 的标准sinf 只需要大约2.5 倍的时间。考虑到您的函数是内联的而 lib 调用不是，这实际上并没有太大区别。我的猜测是，如果您添加了范围缩减，因此它可以以与标准功能相同的方式使用，您将拥有完全相同的速度。
否，最大误差为 0.001 (1/10th %)。您是否忘记应用更正？ (y = P * bla...) 看原文和讨论：devmaster.net/forums/topic/4648-fast-and-accurate-sinecosine 其次，sin 和 cos 以 +-pi 为界是非常常见的情况，尤其是在图形和模拟中，这两者都经常需要一个快速的近似 sin/cos。
这是一个非常有趣的问题，感谢您的提问！

标签： c++ c algorithm math performance

【解决方案1】：

有备用内存吗？查找表（如果需要，可以使用插值）是最快的。

【讨论】：

我如何将它实现为 C 函数？
@Jex：边界检查你的论点（它必须在 -1 和 1 之间）。然后乘以 2 的幂，比如 64，得到范围 (-64, 64)。添加 64 使其成为非负数 (0, 128)。使用整数部分来索引查找表，如果需要，使用小数部分在两个最接近的条目之间进行插值。如果您不想插值，请尝试添加 64.5 并取舍，这与四舍五入相同。
查找表需要索引，这将需要浮点到整数的转换，这可能会影响性能。
@phkahler：浮点到整数的转换在 x86 上非常便宜，几乎与 FP 添加一样便宜，如您所见 in Agner Fog's latency/throughput/uop tables。对索引进行范围检查以确保它不会在表外建立索引可能同样昂贵。 int idx = x * 4096.0 在 Intel Haswell 上会有大约 9 个周期的延迟。到目前为止，最昂贵的部分是来自一个体面大小的表的缓存未命中。如果没有一堆不依赖于 acos 结果的并行计算，那么大表可能会更慢（尤其是缓存竞争）。

【解决方案2】：

一个简单的三次近似，x ∈ {-1, -½, 0, ½, 1} 的拉格朗日多项式是：

double acos(x) {
   return (-0.69813170079773212 * x * x - 0.87266462599716477) * x + 1.5707963267948966;
}

最大误差约为 0.18 rad。

【讨论】：

最大误差为 10.31 度。相当大，但在某些解决方案中可能就足够了。适用于计算速度比精度更重要的地方。可能四次近似会产生更高的精度并且仍然比原生 acos 更快？
确定这个公式没有错误？刚刚用 Wolfram Alpha 试了一下，它看起来不正确：wolframalpha.com/input/?i=y%3D%282%2F9*pixx-5*pi%2F18%29*x%2Bpi%2F2

【解决方案3】：

您可以采取的另一种方法是使用复数。来自de Moivre's formula，

ⅈ^x = cos(π/2*x) + ⅈ*sin(π/2*x)

令 θ = π/2*x。那么x = 2θ/π，所以

sin(θ) = ℑ(ⅈ^^2θ/π)
cos(θ) = ℜ(ⅈ^^2θ/π)

如果没有 sin 和 cos，你如何计算 ⅈ 的幂？从预先计算好的 2 次幂表开始：

ⅈ⁴ = 1
ⅈ² = -1
ⅈ¹ = ⅈ
ⅈ^1/2 = 0.7071067811865476 + 0.7071067811865475*ⅈ
ⅈ^1/4 = 0.9238795325112867 + 0.3826834323650898*ⅈ
ⅈ^1/8 = 0.9807852804032304 + 0.19509032201612825*ⅈ
ⅈ^1/16 = 0.9951847266721969 + 0.0980171403295606*ⅈ
ⅈ^1/32 = 0.9987954562051724 + 0.049067674327418015*ⅈ
ⅈ^1/64 = 0.9996988186962042 + 0.024541228522912288*ⅈ
ⅈ^1/128 = 0.9999247018391445 + 0.012271538285719925*ⅈ
ⅈ^1/256 = 0.9999811752826011 + 0.006135884649154475*ⅈ

要计算 ⅈ^x 的任意值，请将指数近似为二进制分数，然后将表中的相应值相乘。

例如，求 72° 的 sin 和 cos = 0.8π/2：

ⅈ^0.8 &大约; ⅈ^205/256 = ⅈ^0b11001101 = ⅈ^1/2 * ⅈ^1/4 * ⅈ^1/32 * ⅈ^1/64 * ⅈ ^1/256
= 0.3078496400415349 + 0.9514350209690084*ⅈ

sin(72°) & 约; 0.9514350209690084（“精确”值为 0.9510565162951535）
cos(72°) & 约; 0.3078496400415349（“精确”值为 0.30901699437494745）。

要查找 asin 和 acos，您可以将此表与二分法一起使用：

例如，求 asin(0.6)（3-4-5 三角形中的最小角）：

ⅈ⁰ = 1 + 0*ⅈ。 sin 太小了，所以把 x 增加 1/2。
ⅈ^1/2 = 0.7071067811865476 + 0.7071067811865475*ⅈ。罪过大，所以将 x 减少 1/4。
ⅈ^1/4 = 0.9238795325112867 + 0.3826834323650898*ⅈ。 sin 太小了，把 x 增加 1/8。
ⅈ^3/8 = 0.8314696123025452 + 0.5555702330196022*ⅈ。 sin 还是太小了，所以把 x 增加 1/16。
ⅈ^7/16 = 0.773010453362737 + 0.6343932841636455*ⅈ。罪过大，所以将 x 减少 1/32。
ⅈ^13/32 = 0.8032075314806449 + 0.5956993044924334*ⅈ。

每次增加 x 时，乘以 ⅈ 的相应幂。每次减少 x，除以相应的 ⅈ 次方。

如果我们停在这里，我们得到 acos(0.6) ≈ 13/32*π/2 = 0.6381360077604268（“精确”值是 0.6435011087932844。）

当然，准确性取决于迭代次数。对于快速而粗略的近似，使用 10 次迭代。对于“高精确度”，使用 50-60 次迭代。

【讨论】：

【解决方案4】：

我有自己的。它非常准确并且有点快。它的工作原理是我围绕四次收敛建立的定理。这真的很有趣，您可以在此处查看方程式以及它可以使我的自然对数近似收敛多快：https://www.desmos.com/calculator/yb04qt8jx4

这是我的 arccos 代码：

function acos(x)
    local a=1.43+0.59*x a=(a+(2+2*x)/a)/2
    local b=1.65-1.41*x b=(b+(2-2*x)/b)/2
    local c=0.88-0.77*x c=(c+(2-a)/c)/2
    return (8*(c+(2-a)/c)-(b+(2-2*x)/b))/6
end

其中很多只是平方根近似值。它也非常有效，除非你太接近取 0 的平方根。它的平均误差（不包括 x=0.99 到 1）为 0.0003。然而，问题在于，在 0.99 时它开始变得糟糕，而在 x=1 时，准确度的差异变为 0.05。当然，这可以通过在平方根上进行更多迭代来解决（哈哈，不），或者，如果 x>0.99 则使用一组不同的平方根线性化，但这会使代码变得又长又丑.

如果您不太关心准确性，您可以只对每个平方根进行一次迭代，这仍应使您保持在 0.0162 范围内或就准确性而言：

function acos(x)
    local a=1.43+0.59*x a=(a+(2+2*x)/a)/2
    local b=1.65-1.41*x b=(b+(2-2*x)/b)/2
    local c=0.88-0.77*x c=(c+(2-a)/c)/2
    return 8/3*c-b/3
end

如果您没问题，您可以使用预先存在的平方根代码。它将摆脱在 x=1 处有点疯狂的方程：

function acos(x)
    local a = math.sqrt(2+2*x)
    local b = math.sqrt(2-2*x)
    local c = math.sqrt(2-a)
    return 8/3*d-b/3
end

不过，坦率地说，如果您真的时间紧迫，请记住您可以将 arccos 线性化为 3.14159-1.57079x，然后这样做：

function acos(x)
    return 1.57079-1.57079*x
end

无论如何，如果您想查看我的 arccos 近似方程列表，您可以转到 https://www.desmos.com/calculator/tcaty2sv8l 我知道我的近似值对于某些事情不是最好的，但是如果您正在做一些我的近似值会有用，请使用它们，但请尽量给我荣誉。

【讨论】：

这就是你对最后一个的意思吗？ 1.57079-1.57079*x.
对于任何使用 c# 的人来说，这可能是一个很好的第一行： if (x 1D || Double.IsNaN(x)) return Double.NaN;与.net框架acos功能保持一致：msdn.microsoft.com/en-us/library/…
你的第一个实现在 x = -1 时差很多，比如 0.5rad。

【解决方案5】：

nVidia has some great resources 展示了如何逼近其他非常昂贵的数学函数，例如：acos asin atan2 等等等等……

当执行速度（在合理范围内）比精度更重要时，这些算法会产生良好的结果。这是他们的 acos 函数：

// Absolute error <= 6.7e-5
float acos(float x) {
  float negate = float(x < 0);
  x = abs(x);
  float 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 * sqrt(1.0-x);
  ret = ret - 2 * negate * ret;
  return negate * 3.14159265358979 + ret;
}

以下是计算 acos(0.5) 时的结果：

nVidia:   result: 1.0471513828611643
math.h:   result: 1.0471975511965976

这很接近！根据您所需的精确度，这对您来说可能是一个不错的选择。

【讨论】：

与 nvidia 网站上“参考实现”中的评论相比，绝对错误不是
有什么理由使用3.14159265358979 而不是math.pi ？
@ideasman42 Python 与我的回答无关。我的目标只是指出 nVidia 的文档将近似方法作为一种资源。所以为了清楚起见，我编辑了我的答案。但是要回答您的问题：我的猜测是这些数字被选择为可以很好地协同工作。因此，在大多数情况下，使用 math.pi 似乎并没有太大的不同，直到您遇到错误阈值会恶化的边缘情况。
问的原因是有一个非常小的差异导致了这个问题：当将它移植到任何其他语言时 - 是否应该使用 pi 的常量？或者这是一个故意稍微修改的 pi，调整为更好地使用近似值？ （当然可以不用太麻烦测试）
float(x

【解决方案6】：

您可以使用多项式as suggested by dan04 来逼近反余弦，但多项式在 -1 和 1 附近是一个非常糟糕的逼近，其中反余弦的导数趋于无穷大。当您增加多项式的次数时，您会快速达到收益递减，并且仍然很难在端点周围获得良好的近似值。在这种情况下，有理函数（两个多项式的商）可以提供更好的近似值。

acos(x) ≈ π/2 + (ax + bx³) / (1 + cx² + dx⁴)

在哪里

a = -0.939115566365855
b =  0.9217841528914573
c = -1.2845906244690837
d =  0.295624144969963174

在区间 (-1, 1) 上的最大绝对误差为 0.017 弧度（0.96 度）。这是a plot（黑色为反余弦，红色为三次多项式近似，蓝色为上述函数）进行比较：

已选择上述系数以最小化整个域的最大绝对误差。如果您愿意在端点处允许更大的误差，则区间（-0.98, 0.98）上的误差可以变得更小。 5 次分子和 2 次分母与上述函数一样快，但精度稍差。以性能为代价，您可以通过使用更高次多项式来提高准确性。

关于性能的说明：计算两个多项式仍然非常便宜，您可以使用融合乘加指令。除法还不错，因为您可以使用硬件倒数逼近和乘法。与 acos 近似中的误差相比，倒数近似中的误差可以忽略不计。在 2.6 GHz Skylake i7 上，这个近似值可以使用 AVX 每 6 个周期执行大约 8 个反余弦。（即吞吐量，延迟大于6个周期。）

【讨论】：

这些系数有来源吗？
@Gokul 它们是由这个脚本计算的：github.com/ruuda/convector/blob/…

【解决方案7】：

这是一个很棒的网站，有很多选择： https://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/arcsin/onlyelem.html

我个人使用以下代码进行 Chebyshev-Pade 商近似：

double arccos(double x) {
const double pi = 3.141592653;
    return pi / 2 - (.5689111419 - .2644381021*x - .4212611542*(2*x - 1)*(2*x - 1)
         + .1475622352*(2*x - 1)*(2*x - 1)*(2*x - 1))
         / (2.006022274 - 2.343685222*x + .3316406750*(2*x - 1)*(2*x - 1) +
             .02607135626*(2*x - 1)*(2*x - 1)*(2*x - 1));
}

【讨论】：

这在 x = -1 时相差甚远，例如 0.5rad。不可用。

【解决方案8】：

快速反余弦实现，精确到大约 0.5 度，可以基于 observation，对于 [0,1] 中的 x，acos(x) ≈ √(2*(1-x))。一个额外的比例因子提高了接近零的精度。最佳因子可以通过简单的二分搜索找到。负参数按照 acos (-x) = π - acos (x) 处理。

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>

// Approximate acos(a) with relative error < 5.15e-3
// This uses an idea from Robert Harley's posting in comp.arch.arithmetic on 1996/07/12
// https://groups.google.com/forum/#!original/comp.arch.arithmetic/wqCPkCCXqWs/T9qCkHtGE2YJ
float fast_acos (float a)
{
    const float PI = 3.14159265f;
    const float C  = 0.10501094f;
    float r, s, t, u;
    t = (a < 0) ? (-a) : a;  // handle negative arguments
    u = 1.0f - t;
    s = sqrtf (u + u);
    r = C * u * s + s;  // or fmaf (C * u, s, s) if FMA support in hardware
    if (a < 0) r = PI - r;  // handle negative arguments
    return r;
}

float uint_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof(r));
    return r;
}

int main (void)
{
    double maxrelerr = 0.0;
    uint32_t a = 0;
    do {
        float x = uint_as_float (a);
        float r = fast_acos (x);
        double xx = (double)x;
        double res = (double)r;
        double ref = acos (xx);
        double relerr = (res - ref) / ref;
        if (fabs (relerr) > maxrelerr) {
            maxrelerr = fabs (relerr);
            printf ("xx=% 15.8e  res=% 15.8e  ref=% 15.8e  rel.err=% 15.8e\n",
                    xx, res, ref, relerr);
        }
        a++;
    } while (a);
    printf ("maximum relative error = %15.8e\n", maxrelerr);
    return EXIT_SUCCESS;
}

上述测试脚手架的输出应该类似于：

xx= 0.00000000e+000  res= 1.56272149e+000  ref= 1.57079633e+000  rel.err=-5.14060021e-003
xx= 2.98023259e-008  res= 1.56272137e+000  ref= 1.57079630e+000  rel.err=-5.14065723e-003
xx= 8.94069672e-008  res= 1.56272125e+000  ref= 1.57079624e+000  rel.err=-5.14069537e-003
xx=-2.98023259e-008  res= 1.57887137e+000  ref= 1.57079636e+000  rel.err= 5.14071269e-003
xx=-8.94069672e-008  res= 1.57887149e+000  ref= 1.57079642e+000  rel.err= 5.14075044e-003
maximum relative error = 5.14075044e-003

【讨论】：

【解决方案9】：

如果您使用的是 Microsoft VC++，这里有一个内联 __asm x87 FPU 代码版本，没有所有 CRT 填充、错误检查等，与您可以找到的最早的经典 ASM 代码不同，它使用 FMUL 而不是较慢的 FDIV .它与 Microsoft VC++ 2005 Express/Pro 一起编译/工作，我出于各种原因一直坚持使用它。

使用“__declspec(naked)/__fastcall”设置函数、正确提取参数、处理堆栈有点棘手，所以不适合胆小的人。如果它无法在您的版本上编译并出现错误，请不要打扰，除非您有经验。或者问我，我可以用稍微友好的 __asm{} 块重写它。如果它是循环中函数的关键部分，我会手动内联它，以便在需要时进一步提高性能。

extern float __fastcall fs_acos(float x);
extern double __fastcall fs_Acos(double x);

// ACOS(x)- Computes the arccosine of ST(0)
// Allowable range: -1<=x<=+1
// Derivative Formulas: acos(x) = atan(sqrt((1 - x * x)/(x * x))) OR
// acos(x) = atan2(sqrt(1 - x * x), x)
// e.g. acos(-1.0) = 3.1415927

__declspec(naked) float __fastcall fs_acos(float x) { __asm {
    FLD   DWORD PTR [ESP+4] ;// Load/Push parameter 'x' to FPU stack
    FLD1            ;// Load 1.0
    FADD  ST, ST(1) ;// Compute 1.0 + 'x'
    FLD1            ;// Load 1.0
    FSUB  ST, ST(2) ;// Compute 1.0 - 'x'
    FMULP ST(1), ST ;// Compute (1-x) * (1+x)
    FSQRT           ;// Compute sqrt(result)
    FXCH  ST(1)
    FPATAN          ;// Compute arctangent of result / 'x' (ST1/ST0)
    RET 4
}}

__declspec(naked) double __fastcall fs_Acos(double x) { __asm { //
    FLD   QWORD PTR [ESP+4] ;// Load/Push parameter 'x' to FPU stack
    FLD1            ;// Load 1.0
    FADD  ST, ST(1) ;// Compute (1.0 + 'x')
    FLD1            ;// Load 1.0
    FSUB  ST, ST(2) ;// Compute (1.0 - 'x')
    FMULP ST(1), ST ;// Compute (1-x) * (1+x)
    FSQRT           ;// Compute sqrt((1-x) * (1+x))
    FXCH  ST(1) 
    FPATAN          ;// Compute arctangent of result / 'x' (ST1/ST0)
    RET 8
}}

【讨论】：

我怀疑 FPU 会比 SSE 指令更快，而且它对于 x64 目标不可用，因为 MSVC 不允许此类目标的内联 asm 块

【解决方案10】：

很遗憾，我没有足够的声誉来发表评论。这是对 Nvidia 函数的一个小修改，它处理数字应该

这可能很重要，因为舍入误差会导致应该是 1.0 的数字（哦，稍微）大于 1.0。


double safer_acos(double x) {
  double negate = double(x < 0);
  x = abs(x);
  x -= double(x>1.0)*(x-1.0); // <- equivalent to min(1.0,x), but faster
  double 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 * sqrt(1.0-x);
  ret = ret - 2 * negate * ret;
  return negate * 3.14159265358979 + ret;

  // In a single line (no gain using gcc)
  //return negate * 3.14159265358979 + (((((-0.0187293*x)+ 0.0742610)*x - 0.2121144)*x + 1.5707288)* sqrt(1.0-x))*(1.0-2.0*negate);

}

【讨论】：