libs/newlib-cygwin
libs
/
newlib-cygwin
mirror of git://sourceware.org/git/newlib-cygwin.git
4
0
Fork 0
Code Issues Projects Releases Wiki Activity

Fix for k_tan.c specific inputs 

This fix for k_tan.c is a copy from fdlibm version 5.3 (see also
http://www.netlib.org/fdlibm/readme), adjusted to use the macros
available in newlib (SET_LOW_WORD).

This fix reduces the ULP error of the value shown in the fdlibm readme
(tan(1.7765241907548024E+269)) to 0.45 (thereby reducing the error by
1).

This issue only happens for large numbers that get reduced by the range
reduction to a value smaller in magnitude than 2^-28, that is also
reduced an uneven number of times. This seems rather unlikely given that
one ULP is (much) larger than 2^-28 for the values that may cause an
issue.  Although given the sheer number of values a double can
represent, it is still possible that there are more affected values,
finding them however will be quite hard, if not impossible.

We also took a look at how another library (libm in FreeBSD) handles the
issue: In FreeBSD the complete if branch which checks for values smaller
than 2^-28 (or rather 2^-27, another change done by FreeBSD) is moved
out of the kernel function and into the external function. This means
that the value that gets checked for this condition is the unreduced
value. Therefore the input value which caused a problem in the
fdlibm/newlib kernel tan will run through the full polynomial, including
the careful calculation of -1/(x+r). So the difference is really whether
r or y is used. r = y + p with p being the result of the polynomial with
1/3*x^3 being the largest (and magnitude defining) value. With x being
<2^-27 we therefore know that p is smaller than y (y has to be at least
the size of the value of x last mantissa bit divided by 2, which is at
least x*2^-51 for doubles) by enough to warrant saying that r ~ y.  So
we can conclude that the general implementation of this special case is
the same, FreeBSD simply has a different philosophy on when to handle
especially small numbers.
This commit is contained in:
Fabian Schriever 2020-03-17 15:48:44 +01:00 committed by Corinna Vinschen
parent 7947581905
commit 9e8da7bd21
1 changed files with 21 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
newlib/libm/math/k_tan.c
View File

@ -84,14 +84,27 @@ T[] = {
__int32_t ix,hx; __int32_t ix,hx;
GET_HIGH_WORD(hx,x); GET_HIGH_WORD(hx,x);
ix = hx&0x7fffffff; /* high word of |x| */ ix = hx&0x7fffffff; /* high word of |x| */
if(ix<0x3e300000) /* x < 2**-28 */ if(ix<0x3e300000) { /* x < 2**-28 */
{if((int)x==0) { /* generate inexact */ if((int)x==0) { /* generate inexact */
__uint32_t low; __uint32_t low;
GET_LOW_WORD(low,x); GET_LOW_WORD(low,x);
if(((ix|low)|(iy+1))==0) return one/fabs(x); if(((ix|low)|(iy+1))==0) return one/fabs(x);
else return (iy==1)? x: -one/x; else {
} if(iy==1)
} return x;
else {
double a, t;
z = w = x + y;
SET_LOW_WORD(z,0);
v = y - (z - x);
t = a = -one / w;
SET_LOW_WORD(t,0);
s = one + t * z;
return t + a * (s + t * v);
}
}
}
}
if(ix>=0x3FE59428) { /* |x|>=0.6744 */ if(ix>=0x3FE59428) { /* |x|>=0.6744 */
if(hx<0) {x = -x; y = -y;} if(hx<0) {x = -x; y = -y;}
z = pio4-x; z = pio4-x;