On 07/28/2012 04:46 PM, Stephen Montgomery-Smith wrote:
Here are some diffs to catrig.c so that it completely passes Peter Jeremy's program www.rulingia.com/~peter/ctest.c. That is, it seems to get all the infs and nans correct.
And I think I messed up these diffs as well. Can we try this instead?
--- catrig.c-old1 2012-07-28 15:00:46.000000000 -0500 +++ catrig.c 2012-07-28 21:10:20.000000000 -0500 @@ -89,7 +89,7 @@ S = hypot(x,y-1); A = 0.5*(R + S); - if (A < 10) { + if (A < 10 && isfinite(A)) { fp = f(x,1+y,&fpuf); fm = f(x,1-y,&fmuf); if (fpuf == 1 && fmuf == 1) { @@ -108,9 +108,24 @@ } else { *rx = log1p(fp + fm + sqrt((fp+fm)*(A+1))); } - } else + } else if (isinf(y)) + *rx = y; + else *rx = log(A + sqrt(A*A-1)); + if (!isfinite(y)) { + *B_good = 0; + if (isfinite(x)) *A2my2 = 0; + else if (isnan(x)) *A2my2 = x; + else *A2my2 = y; + return; + } + if (isnan(x) && y == 0) { + *B_good = 0; + *A2my2 = 1; + return; + } + *B = y/A; /* = 0.5*(R - S) */ *B_good = 1; @@ -147,7 +162,7 @@ x = fabs(x); y = fabs(y); - if (cabs(z) > 1e20) { + if (cabs(z) > 1e20 && isfinite(x) && isfinite(y)) { if (huge+x>one) { /* set inexact flag. */ if (sx == 0) return clog(2*z); if (sx == 1) return -clog(-2*z); @@ -164,8 +179,8 @@ else ry = atan2(y,A2my2); - if (sx == 1) rx = -rx; - if (sy == 1) ry = -ry; + if (sx == 1) rx = copysign(rx, -1); + if (sy == 1) ry = copysign(ry, -1); return cpack(rx,ry); } @@ -206,7 +221,7 @@ x = fabs(x); y = fabs(y); - if (cabs(z) > 1e20) { + if (cabs(z) > 1e20 && isfinite(x) && isfinite(y)) { if (huge+x>one) { /* set inexact flag. */ w = clog(2*z); if (signbit(cimag(w)) == 0) @@ -232,7 +247,7 @@ rx = atan2(A2my2,-x); } - if (sy==0) ry = -ry; + if (sy==0) ry = copysign(ry, -1); return cpack(rx,ry); } @@ -271,11 +286,23 @@ if (huge+x>one) /* set inexact flag. */ return z; - if (cabs(z) > 1e20) - if (huge+x>one) { /* set inexact flag. */ - if (signbit(x) == 0) - return cpack(0,M_PI_2); - return cpack(0,-M_PI_2); + if (isinf(x) && isnan(y)) + return cpack(copysign(0, x), y); + + if (isnan(x) && isinf(y)) + return cpack(copysign(0, x), copysign(M_PI_2,y)); + + if (x == 0 && isnan(y)) + return cpack(x, y); + + if (isnan(x) || isnan(y)) + return clog(z); + + if (cabs(z) > 1e20) { + if (isinf(x) || isinf(y)) + return cpack(copysign(0,x),copysign(M_PI_2,y)); + if (huge+x>one) /* set inexact flag. */ + return cpack(copysign(0,x),copysign(M_PI_2,y)); } if (fabs(y) < 1e-100) { @@ -290,13 +317,17 @@ if (hp < 0.5 || hm < 0.5) rx = 0.25*(log(hp/hm)); + else if (x == 0) + rx = x; else if (x > 0) rx = 0.25*log1p(4*x/hm); else rx = -0.25*log1p(-4*x/hp); if (x==1 || x==-1) { - if (signbit(y) == 0) + if (y==0) + ry = y; + else if (signbit(y) == 0) ry = atan2(2, -y)/2; else ry = atan2(-2, y)/2;
_______________________________________________ freebsd-bugs@freebsd.org mailing list http://lists.freebsd.org/mailman/listinfo/freebsd-bugs To unsubscribe, send any mail to "freebsd-bugs-unsubscr...@freebsd.org"