bin/170206: complex arcsinh, log, etc.

[Stephen Montgomery-Smith]

The following reply was made to PR bin/170206; it has been noted by GNATS.

From: Stephen Montgomery-Smith <stephen at missouri.edu>
To: Bruce Evans <brde at optusnet.com.au>
Cc: freebsd-bugs at freebsd.org, FreeBSD-gnats-submit at freebsd.org,
        Stephen Montgomery-Smith <stephen at freebsd.org>
Subject: Re: bin/170206: complex arcsinh, log, etc.
Date: Sat, 28 Jul 2012 21:29:42 -0500

 This is a multi-part message in MIME format.
 --------------050105020207030702040007
 Content-Type: text/plain; charset=ISO-8859-1; format=flowed
 Content-Transfer-Encoding: 7bit
 
 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?
 
 
 
 
 --------------050105020207030702040007
 Content-Type: text/x-diff;
  name="ca.diff"
 Content-Transfer-Encoding: 7bit
 Content-Disposition: attachment;
  filename="ca.diff"
 
 --- 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;
 
 --------------050105020207030702040007--