
Date: Wed, 11 Apr 2018 02:52:11 +0200 From: Szabolcs Nagy <nsz@...t70.net> To: musl@...ts.openwall.com Subject: Re: catan errors * Rich Felker <dalias@...c.org> [20180410 19:23:28 0400]: > On Wed, Apr 11, 2018 at 01:08:50AM +0200, Szabolcs Nagy wrote: > > * Rich Felker <dalias@...c.org> [20180410 15:50:07 0400]: > > > t = y  1.0; > > > a = x2 + (t * t); > > >  if (a == 0.0) > > >  goto ovrf; > > > > why does a==0 imply x==0? > > OK, at least by my reasoning it's only algebraically true not > necessarily as doubles. > > > if x < sqrt(2)*0x1p538, x2 underflows to 0 in nearest rounding mode. > > > > to handle this correctly extra work would need to be done, so i think > > either way is fine (leaving the goto there or not are both wrong, but > > we dont guarantee correct complex functions yet) > > I'm fine with just moving the check back here. But the special case > (maybe not nearspecial cases with internal over/underflow though) > works fine just replacing the *I with CMPLX. See attached. > i think this makes sense for now. > Rich > diff git a/src/complex/catan.c b/src/complex/catan.c > index 39ce6cf..7dc2afe 100644 >  a/src/complex/catan.c > +++ b/src/complex/catan.c > @@ 91,29 +91,17 @@ double complex catan(double complex z) > x = creal(z); > y = cimag(z); > >  if (x == 0.0 && y > 1.0) >  goto ovrf; >  > x2 = x * x; > a = 1.0  x2  (y * y); >  if (a == 0.0) >  goto ovrf; > > t = 0.5 * atan2(2.0 * x, a); > w = _redupi(t); > > t = y  1.0; > a = x2 + (t * t); >  if (a == 0.0) >  goto ovrf; > > t = y + 1.0; > a = (x2 + t * t)/a; >  w = w + (0.25 * log(a)) * I; >  return w; >  > ovrf: >  // FIXME >  w = MAXNUM + MAXNUM * I; > + w = CMPLX(w, 0.25 * log(a)); > return w; > }
Powered by blists  more mailing lists