AuroraMiddleware/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2024-11-24 14:00:30 +00:00
Code Issues Packages Projects Releases Wiki Activity

math: Fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny

Browse Source 
Checked on x86_64-linux-gnu and i686-linux-gnu.
This commit is contained in:
Paul Zimmermann 2020-08-07 16:14:53 -03:00 committed by Adhemerval Zanella
parent 1cfb471528
commit b7dd366dbe
2 changed files with 18 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
math/auto-libm-test-in
View File

@ -5748,6 +5748,8 @@ j0 0x1p16382
j0 0x1p16383
# the next value generates larger error bounds on x86_64 (binary32)
j0 0x2.602774p+0 xfail-rounding:ibm128-libgcc
# the next value exercises the flt-32 code path for x >= 2^127
j0 0x8.2f4ecp+124
j1 -1.0
j1 0.0

17
sysdeps/ieee754/flt-32/e_j0f.c
View File

@ -55,7 +55,22 @@ __ieee754_j0f(float x)
z = -__cosf(x+x);
if ((s*c)<zero) cc = z/ss;
else ss = z/cc;
}
} else {
/* We subtract (exactly) a value x0 such that
cos(x0)+sin(x0) is very near to 0, and use the identity
sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
sin(x) + cos(x) with extra accuracy. */
float x0 = 0xe.d4108p+124f;
float y = x - x0; /* exact */
/* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
z = __sinf (y);
float eps = 0x1.5f263ep-24f;
/* cos(x0) ~ -sin(x0) + eps */
z += eps * __cosf (x);
/* now z ~ (sin(x)-cos(x))*cos(x0) */
float cosx0 = -0xb.504f3p-4f;
cc = z / cosx0;
}
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)