如何解决仅使用单精度浮点在[0，pi]上近似余弦
我目前正在研究余弦的近似值。由于最终的目标设备是使用32位浮点ALU / LU进行的自行开发，并且具有针对C的专用编译器，因此我无法使用c库数学函数（cosf等）。我的目标是编写各种在准确性和指令/周期数方面不同的方法。
我已经尝试了很多不同的近似算法，从fdlibm，taylor扩展，pade近似，使用maple的remez算法等等。
但是，只要我仅使用浮点精度实现它们，就会大大降低精度。并确保：我知道以双精度，更高的精度根本没有问题...
现在，我有一些近似值，精确到pi / 2（发生最大误差的范围）附近几千ulp，我感到我受到单精度转换的限制。
要解决主题参数减少：输入为弧度。我认为参数减法会由于除法/乘法而导致更多的精度损失。由于我的总输入范围仅为0..pi，因此我决定将参数减为0..pi / 2。
因此，我的问题是：有人知道高精度的余弦函数的单个精度近似值吗？是否有用于优化单精度近似值的算法？您是否知道内置cosf函数是否在内部以单精度或双精度计算值？ 〜
float ua_cos_v2(float x)
{
float output;
float myPi = 3.1415927410125732421875f;
if (x < 0) x = -x;
int quad = (int32_t)(x*0.63661977236f);//quad = x/(pi/2) = x*2/pi
if (x<1.58f && x> 1.57f) //exclude approximation around pi/2
{
output = -(x - 1.57079637050628662109375f) - 2.0e-12f*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f) + 0.16666667163372039794921875f*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f) + 2.0e-13f*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)+ 0.000198412701138295233249664306640625f*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f);
output -= 4.37E-08f;
}
else {
float param_x;
int param_quad = -1;
switch (quad)
{
case 0:
param_x = x;
break;
case 1:
param_x = myPi - x;
param_quad = 1;
break;
case 2:
param_x = x - myPi;
break;
case 3:
param_x = 2 * myPi - x;
break;
}
float c1 = 1.0f,c2 = -0.5f,c3 = 0.0416666679084300994873046875f,c4 = -0.001388888922519981861114501953125f,c5 = 0.00002480158218531869351863861083984375f,c6 = -2.75569362884198199026286602020263671875E-7f,c7 = 2.08583283978214240050874650478363037109375E-9f,c8 = -1.10807162057025010426514199934899806976318359375E-11f;
float _x2 = param_x * param_x;
output = c1 + _x2*(c2 + _x2*(c3 + _x2*(c4 + _x2*(c5 + _x2*(c6 + _x2*(c7
+ _x2* c8))))));
if (param_quad == 1 || param_quad == 0)
output = -output;
}
return output;
}
〜
如果我忘记了任何信息，请随时询问！
预先感谢
解决方法
当然，仅使用本机精度运算，就可以在[0，π]上以任何期望的误差范围> = 0.5 ulp计算余弦。但是，目标越接近正确舍入的函数，就需要越多的前期设计工作和运行时的计算工作。
超越函数的实现通常包括参数减少，核心近似，最终修正以抵消参数减少。在参数减少涉及减法的情况下，需要通过显式或隐式使用更高的精度来避免灾难性的取消。隐式技术可以设计为仅依赖于本机精度计算，例如，在使用IEEE-754 1.57079637e+0f - 4.37113883e-8f（单精度）时，可以将π常量拆分为未评估的和，例如binary32。
当硬件提供融合乘加（FMA）操作时，通过本机精度计算实现高精度要容易得多。 OP没有指定他们的目标平台是否提供此操作，因此我将首先展示一种非常简单的方法，该方法仅依靠乘法和加法即可提供中等精度（最大误差float映射到IEEE-754 binary32格式。
以下内容基于Colin Wallace的博客文章，标题为“使用Chebyshev多项式将sin（x）近似为5 ULP”，在撰写本文时不在线。我最初是检索它here的，而Google目前保留了缓存的副本here。他们提议通过使用sin（x）/（x *（x *（x²-π²）））的x²中的多项式来近似[-π，π]上的正弦，然后将其乘以x *（x²-π²）。更精确地计算a²-b²的标准技巧是将其重写为（a-b）*（a + b）。将π表示为两个浮点数pi_high和pi_low的未求和，可以避免减法期间的灾难性抵消，这将计算x²-π²转换为((x - pi_hi) - pi_lo) * ((x + pi_hi) + pi_lo)。
理想情况下，多项式核逼近应使用极小极大值逼近，它使 min 最大化 max 最大误差。我在这里做了。可以使用Maple或Mathematics等各种标准工具，也可以根据Remez算法创建自己的代码。
对于[0，PI]上的余弦计算，我们可以利用cos（t）= sin（π/ 2-t）的事实。将x =（π/ 2-t）代入x *（x-π/ 2）*（x +π/ 2）得到（π/ 2-t）*（3π/ 2-t）*（-π/ 2 -t）。可以像以前那样将常量分为高低部分（或者使用另一种常见的习惯用法，分为头部和尾部）。
/* Approximate cosine on [0,PI] with maximum error of 4.704174 ulp */
float cosine (float x)
{
const float half_pi_hi = 1.57079637e+0f; // 0x1.921fb6p+0
const float half_pi_lo = -4.37113883e-8f; // -0x1.777a5cp-25
const float three_half_pi_hi = 4.71238899e+0f; // 0x1.2d97c8p+2
const float three_half_pi_lo = -1.19248806e-8f; // -0x1.99bc5cp-27
float p,s,hpmx,thpmx,nhpmx;
/* cos(x) = sin (pi/2 - x) = sin (hpmx) */
hpmx = (half_pi_hi - x) + half_pi_lo; // pi/2-x
thpmx = (three_half_pi_hi - x) + three_half_pi_lo; // 3*pi/2 - x
nhpmx = (-half_pi_hi - x) - half_pi_lo; // -pi/2 - x
/* P(hpmx*hpmx) ~= sin (hpmx) / (hpmx * (hpmx * hpmx - pi * pi)) */
s = hpmx * hpmx;
p = 1.32729383e-10f;
p = p * s - 2.33177868e-8f;
p = p * s + 2.52223435e-6f;
p = p * s - 1.73503853e-4f;
p = p * s + 6.62087463e-3f;
p = p * s - 1.01321176e-1f;
return hpmx * nhpmx * thpmx * p;
}
下面，我将展示一种经典方法，该方法首先在记录象限时将自变量简化为[-π/ 4，π/ 4]。然后，该象限告诉我们是否需要在此主近似区间上计算正弦或余弦的多项式近似，以及是否需要翻转最终结果的符号。此代码假定目标平台支持IEEE-754指定的FMA操作，并且已通过标准C函数fmaf()进行了单精度映射。
该代码非常简单，除了具有用于舍入象限的舍入模式的浮点到整数转换外，该舍入模式用于计算象限，该转换通过“幻数加法”方法执行并与乘法相结合2 /π（等于除以π/ 2）。最大误差小于1.5 ulps。
/* compute cosine on [0,PI] with maximum error of 1.429027 ulp */
float my_cosf (float a)
{
const float half_pi_hi = 1.57079637e+0f; // 0x1.921fb6p+0
const float half_pi_lo = -4.37113883e-8f; // -0x1.777a5cp-25
float c,j,r,sa,t;
int i;
/* subtract closest multiple of pi/2 giving reduced argument and quadrant */
j = fmaf (a,6.36619747e-1f,12582912.f) - 12582912.f; // 2/pi,1.5 * 2**23
a = fmaf (j,-half_pi_hi,a);
a = fmaf (j,-half_pi_lo,a);
/* phase shift of pi/2 (one quadrant) for cosine */
i = (int)j;
i = i + 1;
sa = a * a;
/* Approximate cosine on [-PI/4,+PI/4] with maximum error of 0.87444 ulp */
c = 2.44677067e-5f; // 0x1.9a8000p-16
c = fmaf (c,-1.38877297e-3f); // -0x1.6c0efap-10
c = fmaf (c,4.16666567e-2f); // 0x1.555550p-5
c = fmaf (c,-5.00000000e-1f); // -0x1.000000p-1
c = fmaf (c,1.00000000e+0f); // 1.00000000p+0
/* Approximate sine on [-PI/4,+PI/4] with maximum error of 0.64196 ulp */
s = 2.86567956e-6f; // 0x1.80a000p-19
s = fmaf (s,-1.98559923e-4f); // -0x1.a0690cp-13
s = fmaf (s,8.33338592e-3f); // 0x1.111182p-7
s = fmaf (s,-1.66666672e-1f); // -0x1.555556p-3
t = a * sa;
s = fmaf (s,t,a);
/* select sine approximation or cosine approximation based on quadrant */
r = (i & 1) ? c : s;
/* adjust sign based on quadrant */
r = (i & 2) ? (0.0f - r) : r;
return r;
}
事实证明，在这种特殊情况下，使用FMA在准确性方面仅提供了很小的好处。如果我用fmaf(a,b,c)替换对((a)*(b)+(c))的调用，则最大错误最小增加到1.451367 ulps，即保持在1.5 ulps以下。
我看到@njuffa有很好的方法，但想提出另一种方法：
- 角度可能最初是度数，而不是弧度，并且可以利用它。
- 不取决于
float是IEEE。 - fma可能是weak，所以不要使用它。
使用整数数学执行范围缩小，然后通过自调整泰勒级数找到答案。
#include <assert.h>
static float my_sinf_helper(float xx,float term,unsigned n) {
if (term + 1.0f == 1.0f) {
return term;
}
return term - my_sinf_helper(xx,xx * term / ((n + 1) * (n + 2)),n + 2);
}
static float my_cosf_helper(float xx,unsigned n) {
if (term + 1.0f == 1.0f) {
return term;
}
return term - xx * my_cosf_helper(xx,term / ((n + 1) * (n + 2)),n + 2);
}
// valid for [-pi/4 + pi/4]
static float my_sinf_primary(float x) {
return x * my_sinf_helper(x * x,1.0,1);
}
// valid for [-pi/4 + pi/4]
static float my_cosf_primary(float x) {
return my_cosf_helper(x * x,0);
}
#define MY_PIf 3.1415926535897932384626433832795f
#define D2Rf(d) ((d)*(MY_PIf/180))
float my_cosdf(float x) {
if (x < 0) {x = -x;}
unsigned long long ux = (unsigned long long) x;
x -= (float) ux;
unsigned ux_primary = ux % 360u;
int uxq = ux_primary%90;
if (uxq >= 45) uxq -= 90;
x += uxq;
switch (ux_primary/45) {
case 7: //
case 0: return my_cosf_primary(D2Rf(x));
case 1: //
case 2: return -my_sinf_primary(D2Rf(x));
case 3: //
case 4: return -my_cosf_primary(D2Rf(x));
case 5: //
case 6: return my_sinf_primary(D2Rf(x));
}
assert(0);
return 0;
}
测试代码
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define DBL_FMT "%+24.17e"
typedef struct {
double x,y0,y1,adiff;
unsigned n;
} test;
test worst = {0};
int my_cosd_test(float x) {
test t;
t.x = x;
t.y0 = cos(x*acos(-1)/180);
t.y1 = my_cosdf(x);
t.adiff = fabs(t.y1 - t.y0);
if (t.adiff > worst.adiff) {
t.n = worst.n + 1;
printf("n:%3u x:" DBL_FMT " y0:" DBL_FMT " y1:" DBL_FMT " d:" DBL_FMT "\n",//
t.n,t.x,t.y0,t.y1,t.adiff);
fflush(stdout);
worst = t;
if (t.n > 100)
exit(-1);
}
return t.adiff != 0.0;
}
float rand_float_finite(void) {
union {
float f;
unsigned char uc[sizeof(float)];
} u;
do {
for (size_t i = 0; i < sizeof u.uc / sizeof u.uc[0]; i++) {
u.uc[i] = (unsigned char) rand();
}
} while (!isfinite(u.f) || fabs(u.f) > 5000);
return u.f;
}
int my_cosd_tests(unsigned n) {
my_cosd_test(0.0);
for (unsigned i = 0; i < n; i++) {
my_cosd_test(rand_float_finite());
}
return 0;
}
int main(void) {
my_cosd_tests(1000000);
}
最差的投放错误：+ 8.2e-08。最大递归深度注释：6。
n: 14 x:+3.64442993164062500e+03 y0:+7.14107074054115110e-01 y1:+7.14107155799865723e-01 d:+8.17457506130381262e-08
我将在以后再进行评论。我确实看到了更广泛的测试，达到了大约9e-08最坏情况的错误以及x > about 1e10的一些TBD问题。
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。