我有自己的非常快的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)的方法呢？谢谢。
一个简单的立方近似，针对 x ∈ {-1, -½, 0, ½, 1} 的 Lagrange 插值多项式为：
double acos(x) {
return (-0.69813170079773212 * x * x - 0.87266462599716477) * x + 1.5707963267948966;
}
它的最大误差约为0.18弧度。
有多余的内存吗？使用查找表（如有必要，包含插值）会是最快的方式。
nVidia提供一些非常好的资源，展示了如何近似计算一些本来很昂贵的数学函数，例如：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
这非常接近！根据您所需的精度程度，这可能是一个很好的选择。
3.14159265358979 而不是 math.pi 的理由是什么？ - ideasman42我有自己的算法。它非常精确并且速度相当快。它是基于我构建的四次收敛理论工作的。这个方程式很有趣，你可以在这里看到它以及它对我的对数近似收敛速度的影响：https://www.desmos.com/calculator/yb04qt8jx4
这是我的反余弦代码：
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
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
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度）。这是一个图形（反余弦为黑色，三次多项式逼近为红色，上述函数为蓝色）进行比较：
上述系数已被选择以使整个定义域内的最大绝对误差最小化。如果您愿意允许端点处的误差更大，则可以使区间(-0.98, 0.98)上的误差变得更小。分子为5次，分母为2次的多项式与上述函数速度相当，但略微不准确。通过使用更高次数的多项式，可以提高精度，但会牺牲性能。
关于性能的说明：计算这两个多项式仍然非常便宜，并且您可以使用融合乘加指令。除法并不那么糟糕，因为您可以使用硬件倒数近似和乘法。与acos近似的误差相比，倒数近似的误差可以忽略不计。在2.6 GHz Skylake i7上，使用AVX，这种逼近法可以每6个周期大约执行8个反余弦操作。（这是吞吐量，延迟时间比6个周期长。）
另一种方法是使用复数。从de Moivre's formula可以得出:
ⅈx = cos(π/2*x) + ⅈ*sin(π/2*x)
令θ = π/2*x。那么x = 2θ/π，因此
如何在没有sin和cos的情况下计算ⅈ的幂次？从预先计算的2的幂表开始：
要计算任意值的ⅈx，请将指数近似为二进制分数，然后从表中相应地乘以对应值。
例如，要找到72° = 0.8π/2的sin和cos：
ⅈ0.8
≈ ⅈ205/256
= ⅈ0b11001101
= ⅈ1/2 * ⅈ1/4 * ⅈ1/32 * ⅈ1/64 * ⅈ1/256
= 0.3078496400415349 + 0.9514350209690084*ⅈ
要找到asin和acos，您可以使用此带有二分法的表：
例如，要查找asin(0.6)（3-4-5三角形中最小的角度）：一种快速且准确到约0.5度的反余弦函数实现可以基于观察，对于 x ∈ [0,1]，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
这是一个拥有许多选项的优秀网站： 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));
}
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
}}
很遗憾，我没有足够的声望来进行评论。 这是Nvidia函数的一个小修改，可以处理应该小于等于1的数字，同时尽可能保持性能。
这可能很重要，因为舍入误差可能导致应该为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);
}
math.h库中的sinf函数仅花费大约2.5倍于你的近似函数的时间。考虑到你的函数是内联的而库调用不是，这并没有太大的区别。我猜测，如果你添加范围缩减，使它可以像标准函数一样使用，你将会有完全相同的速度。 - Damon