有没有快速且准确的atan/arctan近似函数/算法？输入：x =（0，1]（最小x要求）。输出：弧度制
double FastArcTan(double x)
{
return M_PI_4*x - x*(fabs(x) - 1)*(0.2447 + 0.0663*fabs(x));
}
我在网上找到了这个函数，但它给出的最大误差为1.6弧度，这太大了。
快速准确的atan/arctan近似算法
7
我在网上找到了这个函数，但它给出的最大误差为1.6弧度，这太大了。
- Yingyan Wang
4
1您输入了错误的 Google 网址... - Eugene Sh.
1不确定那个公式是从哪里来的...你看过使用一阶逼近吗（如http://math.stackexchange.com/questions/128514/solving-the-arctan-of-an-angle-radians-by-hand中所指出的）...像x - (1.0/3) * (x ** 3) + (1.0/5) * (x ** 5)这样的公式不会完美（并且误差随着你接近1而增加），但我非常确定它将小于1.6弧度。 - Foon
1如果您愿意使用慢速、准确的反正切算法预计算一些值，那么您可以使用分段多项式插值。分段多项式插值的好处是您可以获得几乎任意数量的精度；这仅取决于您愿意拥有多少预计算值。您可以通过代数方法确定需要存储的值的数量以获得所需的精度。 - Matthew Pope
OP的公式来自于《Efficient Approximations for the Arctangent Function》。该公式的最大绝对误差为0.0015弧度（0.086º）。 - micycle
4个回答
11
OP报告了最大误差为1.6弧度（92度），这与下面测试OP代码的结果不一致，最大误差输入x：0到1范围约为0.0015弧度。我怀疑是代码错误或超出了0...1的范围。难道OP想表达的是1.6毫弧度吗？也许一个更准确的a*x^5 + b*x^3 + c*x对于OP来说仍然足够快速。在我的机器上，平均精度提高了约4倍，最差情况下提高了2倍。它使用了一个最优的三项式，如@Foon和@Matthew Pope所建议的。
#include <math.h>
#include <stdio.h>
#ifndef M_PI_4
#define M_PI_4 (3.1415926535897932384626433832795/4.0)
#endif
double FastArcTan(double x) {
return M_PI_4*x - x*(fabs(x) - 1)*(0.2447 + 0.0663*fabs(x));
}
#define A 0.0776509570923569
#define B -0.287434475393028
#define C (M_PI_4 - A - B)
#define FMT "% 16.8f"
double Fast2ArcTan(double x) {
double xx = x * x;
return ((A*xx + B)*xx + C)*x;
}
int main() {
double mxe1 = 0, mxe2 = 0;
double err1 = 0, err2 = 0;
int n = 100;
for (int i=-n;i<=n; i++) {
double x = 1.0*i/n;
double y = atan(x);
double y_fast1 = FastArcTan(x);
double y_fast2 = Fast2ArcTan(x);
printf("%3d x:% .3f y:" FMT "y1:" FMT "y2:" FMT "\n", i, x, y, y_fast1, y_fast2);
if (fabs(y_fast1 - y) > mxe1 ) mxe1 = fabs(y_fast1 - y);
if (fabs(y_fast2 - y) > mxe2 ) mxe2 = fabs(y_fast2 - y);
err1 += (y_fast1 - y)*(y_fast1 - y);
err2 += (y_fast2 - y)*(y_fast2 - y);
}
printf("max error1: " FMT "sum sq1:" FMT "\n", mxe1, err1);
printf("max error2: " FMT "sum sq2:" FMT "\n", mxe2, err2);
}
输出
...
96 x: 0.960 y: 0.76499283y1: 0.76582280y2: 0.76438526
97 x: 0.970 y: 0.77017091y1: 0.77082844y2: 0.76967407
98 x: 0.980 y: 0.77529750y1: 0.77575981y2: 0.77493733
99 x: 0.990 y: 0.78037308y1: 0.78061652y2: 0.78017777
100 x: 1.000 y: 0.78539816y1: 0.78539816y2: 0.78539816
max error1: 0.00150847sum sq1: 0.00023062
max error2: 0.00084283sum sq2: 0.00004826
不清楚为什么 OP 的代码在“输入：x =（0，1］”时使用了 fabs()。
- chux - Reinstate Monica
7
1感谢您的帮助。我意识到我的一些输入值大于1，这就是为什么它给出了最大1.6弧度的错误。OP的代码使用fabs()，因为它适用于[-1,1]的输入范围，而(0,1]只是我自己输入数据的范围。 - Yingyan Wang
@YingyanWang 对于
x > 1 的常见范围缩减，可以使用类似 return M_PI_2 - FastArcTan(1/x) 的代码，对于 x < 1 也是如此。请参考 https://dev59.com/sWAg5IYBdhLWcg3w5eZk#23097989 - chux - Reinstate Monica你是如何计算系数A、B、C的？ - Anh Dũng Lê
已经过去5年了。嗯，如果我没记错的话，对于
A,B，我在Excel中使用最佳二次拟合法进行了数百个均匀间隔数据点的拟合。 C =（M_PI_4 - A - B）用于将曲线拟合到x = 1.0处。 - chux - Reinstate Monica如果atan(1)应该接近pi/4，误差约为7.04E-4，则返回0.07606631753154554x^5-0.2854342007925691x^3+0.9947660466584719*x。 - undefined
显示剩余2条评论
4
我根据Pade-Chebyshev逼近方法为我的需求开发了以下函数。代码是在Visual Studio（x86）下编写的，并且可以在整个源数据范围内工作。该函数的准确度仅略低于标准的atanf函数，但速度明显更快。
// Fast arctangent with single precission
_declspec(naked) float _vectorcall arctg(float x)
{
// When |x|<=1 following formula is used:
//
// a0 + a0*a1*x^2 + a0*a2*x^4 + a0*a3*x^6
// arctg x = x * ------------------------------------------------------- = x*P(x)/Q(x)
// a0 + a0*(a1+b0)*x^2 + a0*(a2+b1)*x^4 + a0*(a3+b2)*x^6
//
// The a0 constant is reduced, but it is needed to less the error
// and prevent overflow at large |x|.
// P(x) and Q(x) are 3th degree polinomials of x^2.
// When 1<|x|<62919776 (approx.) used formula is
//
// pi/2*|x|*x^6*Q(1/x)-x^6*P(1/x)
// arctg x = pi/2*sgn(x)-arctg(1/x) = --------------------------------
// x*x^6*Q(1/x)
//
// Here x^6*P(1/x) and x^6*Q(1/x) are 3th degree polinomials of x^2 too.
// When |x|>=62919776 used formula is arctg x = pi/2*sgn(x)
// To improve accuracy at |x|>1 the constant pi/2 is replaced by the sum (pi-3)/2 and 3/2.
//
static const float ct[14] = // Constants table
{
6.28740248E-17f, // a0*(a1+b0)=a0*c1
4.86816205E-17f, // a0*a1
2.24874633E-18f, // a0*(a3+b2)=a0*c3
4.02179944E-19f, // a0*a3
4.25772129E-17f, // a0
4.25772129E-17f, // a0
2.50182216E-17f, // a0*(a2+b1)=a0*c2
1.25756219E-17f, // a0*a2
0.0707963258f, // (pi-3)/2
1.5f, // 3/2
1.0f, // 1
-0.0707963258f, // -(pi-3)/2
-1.5f, // -3/2
3.95889818E15f // Threshold of x^2 when arctg(x)=pi/2*sgn(x)
};
_asm
{
vshufps xmm1,xmm0,xmm0,0 // xmm1 = x # x : x # x
mov edx,offset ct // edx contains the address of constants table
vmulps xmm2,xmm1,xmm1 // xmm2 = x^2 # x^2 : x^2 # x^2
vmovups xmm3,[edx+16] // xmm3 = a0*a2 # a0*c2 : a0 # a0
vmulps xmm4,xmm2,xmm2 // xmm4 = y^2 # y^2 : y^2 # y^2
vucomiss xmm2,[edx+40] // Compare y=x^2 to 1
ja arctg_big // Jump if |x|>1
vfmadd231ps xmm3,xmm2,[edx] // xmm3 ~ a3*y+a2 # c3*y+c2 : a1*y+1 # c1*y+1
vmovhlps xmm1,xmm1,xmm3 // xmm1 ~ a3*y+a2 # c3*y+c2
vfmadd231ps xmm3,xmm4,xmm1 // xmm3 ~ a3*y^3+a2*y^2+a1*y+1 # c3*y^3+c2*y^2+c1*y+1
vmovshdup xmm2,xmm3 // xmm2 = P; xmm3 = Q
vdivss xmm2,xmm2,xmm3 // xmm2 = P/Q
vmulss xmm0,xmm0,xmm2 // xmm0 = x*P/Q = arctg(x)
ret // Return
arctg_big: // When |x|>1 use formula pi/2*sgn(x)-arctg(1/x)
vfmadd213ps xmm3,xmm2,[edx] // xmm3 ~ a2*y+a3 # c2*y+c3 : y+a1 # y+c1
vmovmskpd eax,xmm1 // eax=3, if x<0, otherwise eax=0
vmovhlps xmm0,xmm0,xmm3 // xmm0 ~ a2*y+a3 # c2*y+c3
vfmadd213ps xmm3,xmm4,xmm0 // xmm3 ~ y^3+a1*y^2+a2*y+a3 # y^3+c1*y^2+c2*y+c3
vmovss xmm0,[edx+4*eax+32] // xmm0 = (pi-3)/2*sgn(x)
vucomiss xmm2,[edx+52] // Compare y=x^2 to threshold value
jnb arctg_end // The data is already in xmm0, if |x|>=62919776
vmovshdup xmm4,xmm3 // xmm4 = P; xmm3 = Q
vmulss xmm1,xmm1,xmm3 // xmm1 = x*Q
vfmsub132ss xmm0,xmm4,xmm1 // xmm0 = (pi-3)/2*|x|*Q-P
vdivss xmm0,xmm0,xmm1 // xmm0 = (pi-3)/2*sgn(x)-P/(x*Q)
arctg_end: // Add to result 3/2*sgn(x)
vaddss xmm0,xmm0,[edx+4*eax+36] // xmm0 = pi/2*sgn(x)-P/(x*Q)
ret // Return
}
}
- aenigma
22
1该函数在整个范围内提供了几乎完整的浮点精度。误差主要是由于中间计算结果的舍入误差所致。逼近的分数有理函数本身的相对误差不超过3*10^-9，这超过了浮点格式的精度。 - aenigma
1@aenigma：我让我的启发式优化器运行了十天。找到的最准确的解决方案在|x|=0.544157028f处具有最大误差4.02802 ulps。系数为：
6.28740778e-17f，4.86816768e-17f，2.24874860e-18f，4.02182452e-19f，4.25772063e-17f，4.25772063e-17f，2.50182596e-17f，1.25756368e-17f。 - njuffa
1@njuffa 谢谢。我也在寻找更准确的最大ulp解决方案。但是程序将长时间运行。如果我找到了什么，我会分享的。我最初提出的系数在RMS相对误差方面提供了接近最优解的解决方案。 - aenigma
1@kimstik 所使用的近似精度应该取决于应用的需求。但是，如果创建了足够统一的算法，则希望确保接近数据格式的完全准确性。对于具有尾数最低数字的绝对误差为一半的浮点类型数字，平均相对误差约为4E-8。因此，为了消除误差的系统成分，一个好的近似应该提供一个数量级更高的准确性。 - aenigma
2@njuffa 这是一个提供最大误差约为3个ULP的解决方案。系数为：5.10944634e-17f，3.95610345e-17f，1.82743967e-18f，3.26830821e-19f，3.46002967e-17f，3.46002967e-17f，2.03310142e-17f，1.02195567e-17f。但它给出了稍大的均方相对误差：2.32272978e-8对比2.22521892e-8。 - aenigma
显示剩余17条评论
2
我发现这个公式适用于 x [0,+inf]。
atan(x) = 8x/(3+sqrt(25+80/3*x^2))
对于 x>2，我添加了 +0.02 来修正。
atan(x) = 8x/(3+sqrt(25+80/3*x^2))
对于 x>2，我添加了 +0.02 来修正。
- Johnwhile
6
1平均相对逼近误差是多少？你一直在优化常数吗？如果选取了合适的逼近函数，那么系数的选择应该会有良好的结果。当x<0时，使用您的函数是否存在问题？ - aenigma
这个公式可以在这里找到：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6010512/ 我只在[0,+inf]范围内使用它，但它的准确性不是很高。 - Johnwhile
1谢谢提供链接！非常有趣。显然，这个近似并不是非常准确的，但就其准确性而言，它表现出良好的性能，也就是说，它有存在的合理性。如果对其中的系数进行优化，并且以高效的方式编写代码（例如使用汇编语言），这个近似可以稍微改进一下。 - aenigma
顺便说一下，这个近似函数是奇函数，所以对于任何x（包括x<0）都适用。 - aenigma
你所给出的函数是反正切的下界（公式1.4，1.5）。为了得到更好的近似值，你可以尝试取上界和下界系数的平均值。例如，可以从公式1.7得到以下近似公式：atan(x) = pi^2x/(4+sqrt((pi^2-4)sqrt(32)+(2pix)^2))。 - aenigma
根据您的公式，在公式1.5中，8的第一个系数可以被替换为4+π*sqrt(15)/3。 - aenigma
0
反正切函数在 z=x/(x+pi/3) 空间中可以很好地近似为多项式。以下系数在您的区间上产生最大误差为 4.07e-9。
z0 = 0.0;
z1 = 1.0471985201240248;
z2 = 1.047102380719068;
z3 = 0.6678082665364844;
z4 = -0.16312555173050677;
z5 = -0.33989938855441837;
z6 = -5.8773286456840355;
z7 = 17.33961311116544;
z8 = -49.09950369422959;
z9 = 82.10190353000648;
z10 = -60.11210171537528;
z11 = 14.183376881192657;
你也可以通过更高阶的多项式（19次）将误差降至1.23e-14。
以下是一些C#代码（引用Accord.net），可用于查看精度和系数数量之间的权衡。
using Accord.Math;
using Accord.Statistics.Models.Regression.Linear;
using System.Management;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
foreach(var degree in new[] { 3, 5, 7, 11, 15, 19, 23, 27 })
{
var N = 5000;
var ols = new OrdinaryLeastSquares() { UseIntercept = false };
var inp = new double[N][];
var outputs = new double[N];
for (int i = 0; i < N; ++i)
{
double x = (double)i / (double)N;
inp[i] = new double[degree];
inp[i][0] = x / (Math.PI /3 + x);
for (int j = 1; j < degree; ++j)
inp[i][j] = inp[i][0] * inp[i][j - 1];
outputs[i] = Math.Atan(x);
}
MultipleLinearRegression regression = ols.Learn(inp, outputs);
var r = new List<double>();
foreach (var w in regression.Weights)
Console.WriteLine(w);
var model = inp.Dot(regression.Weights);
var diff = 0.0;
for (int i = 0; i < N; ++i)
diff = Math.Max(Math.Abs(model[i] - outputs[i]), diff);
Console.WriteLine("++++++++++++++++++++++");
Console.WriteLine(diff);
Console.WriteLine("++++++++++++++++++++++");
}
Console.ReadLine();
- kolbe
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接