From 3738a96e052603403e085e9a1024289ba3e09188 Mon Sep 17 00:00:00 2001
From: nsz <nsz@port70.net>
Date: Tue, 8 May 2012 00:22:56 +0200
Subject: [PATCH] math: fix remquo.c when x==-y and a subnormal remainder bug
 as well backported fix from freebsd:
 http://svnweb.FreeBSD.org/base?view=revision&revision=233973

---
 src/math/remquo.c  | 7 ++++---
 src/math/remquof.c | 3 ++-
 src/math/remquol.c | 3 ++-
 3 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/src/math/remquo.c b/src/math/remquo.c
index e92984ed..085466e6 100644
--- a/src/math/remquo.c
+++ b/src/math/remquo.c
@@ -44,7 +44,7 @@ double remquo(double x, double y, int *quo)
 			goto fixup;
 		}
 		if (lx == ly) {            /* |x| = |y| return x*0 */
-			*quo = 1;
+			*quo = sxy ? -1 : 1;
 			return Zero[(uint32_t)sx>>31];
 		}
 	}
@@ -127,6 +127,7 @@ double 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[(uint32_t)sx>>31];
 	}
@@ -144,10 +145,10 @@ double remquo(double x, double y, int *quo)
 			hx >>= n;
 		} else if (n <= 31) {
 			lx = (hx<<(32-n))|(lx>>n);
-			hx = sx;
+			hx = 0;
 		} else {
 			lx = hx>>(n-32);
-			hx = sx;
+			hx = 0;
 		}
 	}
 fixup:
diff --git a/src/math/remquof.c b/src/math/remquof.c
index 11569ce8..536a050a 100644
--- a/src/math/remquof.c
+++ b/src/math/remquof.c
@@ -41,7 +41,7 @@ float remquof(float x, float y, int *quo)
 		q = 0;
 		goto fixup;
 	} else if(hx==hy) {  /* |x| = |y| return x*0*/
-		*quo = 1;
+		*quo = sxy ? -1 : 1;
 		return Zero[(uint32_t)sx>>31];
 	}
 
@@ -92,6 +92,7 @@ float remquof(float x, float y, int *quo)
 
 	/* convert back to floating value and restore the sign */
 	if (hx == 0) {                             /* return sign(x)*0 */
+		q &= 0x7fffffff;
 		*quo = sxy ? -q : q;
 		return Zero[(uint32_t)sx>>31];
 	}
diff --git a/src/math/remquol.c b/src/math/remquol.c
index 721231b4..a2e11728 100644
--- a/src/math/remquol.c
+++ b/src/math/remquol.c
@@ -94,7 +94,7 @@ long double remquol(long double x, long double y, int *quo)
 			goto fixup;       /* |x|<|y| return x or x-y */
 		}
 		if (ux.bits.manh == uy.bits.manh && ux.bits.manl == uy.bits.manl) {
-			*quo = 1;
+			*quo = sxy ? -1 : 1;
 			return Zero[sx];  /* |x|=|y| return x*0*/
 		}
 	}
@@ -152,6 +152,7 @@ long double remquol(long double x, long 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[sx];
 	}
-- 
2.25.1