Re: mysterious remquo bug

[Greg Troxel <gdt%lexort.com@localhost>]

FreeBSD has a change, which resolves the sign of q bug (1 vs -1):
However, it still has q being 0 when it should be -2^n.

The patch at end causes
  my remquo test programs
  proj's regression tests
to all pass.

The first three hunks are from FreeBSD:

  https://cgit.freebsd.org/src/commit/lib/msun/src/s_remquo.c?id=1cbd288942b08217e99bf889e0967895d53af00c

which resolve the sign bug (hunks 1+3 I believe) and an unrelated
subnormal bug (that I doubt I'm hitting, hunk 2).

The fourth hunk is one I wrote, to get the sign right, even though it
was always a valid value modulo 2^k.  (C99 requires that it have the
same sign as x/y.)

Building without the last hunk, I see remquo returning values that I can
say are wrong per the specification  But the proj tests still pass.

Therefore, unless someone can point out why it's wrong, I'm going to
apply the first three hunks, calling it merging a fix that's been in
FreeBSD for 12 years.  I'll run atf without/with before committing.

I also plan to commit the last hunk, simplified some, again assuming no
contrary rationale appears on the list.  I am waiting to hear back from
a numerical expert regarding exactly what the C99 spec requires.


Index: lib/libm/src/s_remquo.c
===================================================================
RCS file: /cvsroot/src/lib/libm/src/s_remquo.c,v
retrieving revision 1.1
diff -u -p -r1.1 s_remquo.c
--- lib/libm/src/s_remquo.c	6 Feb 2011 01:53:38 -0000	1.1
+++ lib/libm/src/s_remquo.c	16 Sep 2024 00:41:12 -0000
@@ -50,7 +50,7 @@ remquo(double x, double y, int *quo)
 		goto fixup;	/* |x|<|y| return x or x-y */
 	    }
 	    if(lx==ly) {
-		*quo = 1;
+		*quo = (sxy ? -1 : 1);
 		return Zero[(u_int32_t)sx>>31];	/* |x|=|y| return x*0*/
 	    }
 	}
@@ -113,6 +113,7 @@ remquo(double x, double y, int *quo)
 
     /* convert back to floating value and restore the sign */
 	if((hx|lx)==0) {			/* return sign(x)*0 */
+	    q &= 0x7fffffff;
 	    *quo = (sxy ? -q : q);
 	    return Zero[(u_int32_t)sx>>31];
 	}
@@ -128,9 +129,9 @@ remquo(double x, double y, int *quo)
 		lx = (lx>>n)|((u_int32_t)hx<<(32-n));
 		hx >>= n;
 	    } else if (n<=31) {
-		lx = (hx<<(32-n))|(lx>>n); hx = sx;
+		lx = (hx<<(32-n))|(lx>>n); hx = 0;
 	    } else {
-		lx = hx>>(n-32); hx = sx;
+		lx = hx>>(n-32); hx = 0;
 	    }
 	}
 fixup:
@@ -149,5 +150,15 @@ fixup:
 	SET_HIGH_WORD(x,hx^sx);
 	q &= 0x7fffffff;
 	*quo = (sxy ? -q : q);
+	/*
+	 * If q is 0 and we need to return negative, we have to just
+	 * choose the largest negative number.
+	 */
+	if (q == 0) {
+	    if (sxy)
+	      *quo = 0x80000000 /* -8 */;
+	    else
+	      *quo = 0 /* 8 */;
+	}
 	return x;
 }