Skip site navigation (1)Skip section navigation (2)
Date:      Sat, 28 Jul 2012 21:29:42 -0500
From:      Stephen Montgomery-Smith <stephen@missouri.edu>
To:        Bruce Evans <brde@optusnet.com.au>
Cc:        freebsd-bugs@freebsd.org, FreeBSD-gnats-submit@freebsd.org, Stephen Montgomery-Smith <stephen@freebsd.org>
Subject:   Re: bin/170206: complex arcsinh, log, etc.
Message-ID:  <5014A016.5090205@missouri.edu>
In-Reply-To: <50145DBF.4090609@missouri.edu>
References:  <201207270247.q6R2lkeR021134@wilberforce.math.missouri.edu> <20120727233939.A7820@besplex.bde.org> <5012A96E.9090400@missouri.edu> <20120728142915.K909@besplex.bde.org> <50137C24.1060004@missouri.edu> <20120728171345.T1911@besplex.bde.org> <50140956.1030603@missouri.edu> <50141018.3040203@missouri.edu> <50145DBF.4090609@missouri.edu>

next in thread | previous in thread | raw e-mail | index | archive | help
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--



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?5014A016.5090205>