Message ID  20200317144844.10751fabian.schriever@gtdgmbh.de 

State  Accepted 
Commit  9e8da7bd2138aaefcb746be3bcce2787c75a5849 
Headers  show 
Series 

Related  show 
On Mar 17 15:48, Fabian Schriever wrote: > 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. >  > newlib/libm/math/k_tan.c  29 +++++++++++++++++++++ > 1 file changed, 21 insertions(+), 8 deletions() Pushed. Thanks, Corinna  Corinna Vinschen Cygwin Maintainer Red Hat
diff git a/newlib/libm/math/k_tan.c b/newlib/libm/math/k_tan.c index 9f5b30760..4be82d599 100644  a/newlib/libm/math/k_tan.c +++ b/newlib/libm/math/k_tan.c @@ 84,14 +84,27 @@ T[] = { __int32_t ix,hx; GET_HIGH_WORD(hx,x); ix = hx&0x7fffffff; /* high word of x */  if(ix<0x3e300000) /* x < 2**28 */  {if((int)x==0) { /* generate inexact */  __uint32_t low;  GET_LOW_WORD(low,x);  if(((ixlow)(iy+1))==0) return one/fabs(x);  else return (iy==1)? x: one/x;  }  } + if(ix<0x3e300000) { /* x < 2**28 */ + if((int)x==0) { /* generate inexact */ + __uint32_t low; + GET_LOW_WORD(low,x); + if(((ixlow)(iy+1))==0) return one/fabs(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(hx<0) {x = x; y = y;} z = pio4x;