Module Name: src
Committed By: isaki
Date: Sat May 11 12:52:43 UTC 2013
Modified Files:
src/sys/arch/m68k/fpe: fpu_rem.c
Log Message:
Revise the algorithm after Step3.
almost written by Y.Sugahara and minor modify by me.
This works for all input of FMOD/FREM and of course solves PR/47810.
To generate a diff of this commit:
cvs rdiff -u -r1.15 -r1.16 src/sys/arch/m68k/fpe/fpu_rem.c
Please note that diffs are not public domain; they are subject to the
copyright notices on the relevant files.
Modified files:
Index: src/sys/arch/m68k/fpe/fpu_rem.c
diff -u src/sys/arch/m68k/fpe/fpu_rem.c:1.15 src/sys/arch/m68k/fpe/fpu_rem.c:1.16
--- src/sys/arch/m68k/fpe/fpu_rem.c:1.15 Sun May 5 13:25:20 2013
+++ src/sys/arch/m68k/fpe/fpu_rem.c Sat May 11 12:52:42 2013
@@ -1,4 +1,4 @@
-/* $NetBSD: fpu_rem.c,v 1.15 2013/05/05 13:25:20 isaki Exp $ */
+/* $NetBSD: fpu_rem.c,v 1.16 2013/05/11 12:52:42 isaki Exp $ */
/*
* Copyright (c) 1995 Ken Nakata
@@ -32,7 +32,7 @@
*/
#include <sys/cdefs.h>
-__KERNEL_RCSID(0, "$NetBSD: fpu_rem.c,v 1.15 2013/05/05 13:25:20 isaki Exp $");
+__KERNEL_RCSID(0, "$NetBSD: fpu_rem.c,v 1.16 2013/05/11 12:52:42 isaki Exp $");
#include <sys/types.h>
#include <sys/signal.h>
@@ -56,51 +56,68 @@ __KERNEL_RCSID(0, "$NetBSD: fpu_rem.c,v
* endif
*
* Step 3. Perform MOD(X,Y)
- * 3.1 If R = Y, go to Step 9.
+ * 3.1 If R = Y, then { Q := Q + 1, R := 0, go to Step 8. }
* 3.2 If R > Y, then { R := R - Y, Q := Q + 1}
* 3.3 If j = 0, go to Step 4.
* 3.4 k := k + 1, j := j - 1, Q := 2Q, R := 2R. Go to
* Step 3.1.
*
- * Step 4. At this point, R = X - QY = MOD(X,Y). Set
- * Last_Subtract := false (used in Step 7 below). If
- * MOD is requested, go to Step 6.
- *
- * Step 5. R = MOD(X,Y), but REM(X,Y) is requested.
- * 5.1 If R < Y/2, then R = MOD(X,Y) = REM(X,Y). Go to
- * Step 6.
- * 5.2 If R > Y/2, then { set Last_Subtract := true,
- * Q := Q + 1, Y := signY*Y }. Go to Step 6.
- * 5.3 This is the tricky case of R = Y/2. If Q is odd,
- * then { Q := Q + 1, signX := -signX }.
- *
- * Step 6. R := signX*R.
- *
- * Step 7. If Last_Subtract = true, R := R - Y.
- *
- * Step 8. Return signQ, last 7 bits of Q, and R as required.
- *
- * Step 9. At this point, R = 2^(-j)*X - Q Y = Y. Thus,
- * X = 2^(j)*(Q+1)Y. set Q := 2^(j)*(Q+1),
- * R := 0. Return signQ, last 7 bits of Q, and R.
+ * Step 4. R := signX*R.
+ *
+ * Step 5. If MOD is requested, go to Step .
+ *
+ * Step 6. Now, R = MOD(X,Y), convert to REM(X,Y) is requested.
+ * Do banker's rounding.
+ * If abs(R) > Y/2
+ * || (abs(R) == Y/2 && Q % 2 == 1) then
+ * { Q := Q + 1, R := R - signX * Y }.
+ *
+ * Step 7. Return signQ, last 7 bits of Q, and R as required.
*/
static struct fpn * __fpu_modrem(struct fpemu *fe, int is_mod);
+static int abscmp3(const struct fpn *a, const struct fpn *b);
+
+/* Absolute FORTRAN Compare */
+static int
+abscmp3(const struct fpn *a, const struct fpn *b)
+{
+ int i;
+
+ if (a->fp_exp < b->fp_exp) {
+ return -1;
+ } else if (a->fp_exp > b->fp_exp) {
+ return 1;
+ } else {
+ for (i = 0; i < 3; i++) {
+ if (a->fp_mant[i] < b->fp_mant[i])
+ return -1;
+ else if (a->fp_mant[i] > b->fp_mant[i])
+ return 1;
+ }
+ }
+ return 0;
+}
static struct fpn *
__fpu_modrem(struct fpemu *fe, int is_mod)
{
static struct fpn X, Y;
struct fpn *x, *y, *r;
- struct fpn r_bkup;
uint32_t signX, signY, signQ;
int j, k, l, q;
- int Last_Subtract;
+ int cmp;
+
+ if (ISNAN(&fe->fe_f1) || ISNAN(&fe->fe_f2))
+ return fpu_newnan(fe);
+ if (ISINF(&fe->fe_f1) || ISZERO(&fe->fe_f2))
+ return fpu_newnan(fe);
CPYFPN(&X, &fe->fe_f1);
CPYFPN(&Y, &fe->fe_f2);
x = &X;
y = &Y;
+ q = 0;
r = &fe->fe_f2;
/*
@@ -111,12 +128,17 @@ __fpu_modrem(struct fpemu *fe, int is_mo
signQ = (signX ^ signY);
x->fp_sign = y->fp_sign = 0;
+ /* Special treatment that just return input value but Q is necessary */
+ if (ISZERO(x) || ISINF(y)) {
+ r = &fe->fe_f1;
+ goto Step7;
+ }
+
/*
* Step 2
*/
l = x->fp_exp - y->fp_exp;
k = 0;
- q = 0;
CPYFPN(r, x);
if (l >= 0) {
r->fp_exp -= l;
@@ -125,21 +147,20 @@ __fpu_modrem(struct fpemu *fe, int is_mo
/*
* Step 3
*/
- while (y->fp_exp != r->fp_exp ||
- y->fp_mant[0] != r->fp_mant[0] ||
- y->fp_mant[1] != r->fp_mant[1] ||
- y->fp_mant[2] != r->fp_mant[2]) {
+ for (;;) {
+ cmp = abscmp3(r, y);
+
+ /* Step 3.1 */
+ if (cmp == 0)
+ break;
/* Step 3.2 */
- CPYFPN(&r_bkup, r);
- CPYFPN(&fe->fe_f1, r);
- CPYFPN(&fe->fe_f2, y);
- fe->fe_f2.fp_sign = 1;
- r = fpu_add(fe);
- if (r->fp_sign == 0) {
+ if (cmp > 0) {
+ CPYFPN(&fe->fe_f1, r);
+ CPYFPN(&fe->fe_f2, y);
+ fe->fe_f2.fp_sign = 1;
+ r = fpu_add(fe);
q++;
- } else {
- CPYFPN(r, &r_bkup);
}
/* Step 3.3 */
@@ -152,75 +173,48 @@ __fpu_modrem(struct fpemu *fe, int is_mo
q += q;
r->fp_exp++;
}
- /* Step 9 */
- goto Step9;
+ /* R == Y */
+ q++;
+ r->fp_class = FPC_ZERO;
+ goto Step7;
}
Step4:
- Last_Subtract = 0;
- if (is_mod)
- goto Step6;
+ r->fp_sign = signX;
/*
* Step 5
*/
- /* Step 5.1 */
- if (r->fp_exp + 1 < y->fp_exp ||
- (r->fp_exp + 1 == y->fp_exp &&
- (r->fp_mant[0] < y->fp_mant[0] ||
- r->fp_mant[1] < y->fp_mant[1] ||
- r->fp_mant[2] < y->fp_mant[2]))) {
- /* if r < y/2 */
- goto Step6;
- }
- /* Step 5.2 */
- if (r->fp_exp + 1 != y->fp_exp ||
- r->fp_mant[0] != y->fp_mant[0] ||
- r->fp_mant[1] != y->fp_mant[1] ||
- r->fp_mant[2] != y->fp_mant[2]) {
- /* if (!(r < y/2) && !(r == y/2)) */
- Last_Subtract = 1;
- q++;
- y->fp_sign = signY;
- } else {
- /* Step 5.3 */
- /* r == y/2 */
- if (q % 2) {
- q++;
- signX = !signX;
- }
- }
-
- Step6:
- r->fp_sign = signX;
+ if (is_mod)
+ goto Step7;
/*
- * Step 7
+ * Step 6
*/
- if (Last_Subtract) {
+ /* y = y / 2 */
+ y->fp_exp--;
+ /* abscmp3 ignore sign */
+ cmp = abscmp3(r, y);
+ /* revert y */
+ y->fp_exp++;
+
+ if (cmp > 0 || (cmp == 0 && q % 2)) {
+ q++;
CPYFPN(&fe->fe_f1, r);
CPYFPN(&fe->fe_f2, y);
- fe->fe_f2.fp_sign = !y->fp_sign;
+ fe->fe_f2.fp_sign = !signX;
r = fpu_add(fe);
}
+
/*
- * Step 8
+ * Step 7
*/
+ Step7:
q &= 0x7f;
q |= (signQ << 7);
fe->fe_fpframe->fpf_fpsr =
fe->fe_fpsr =
(fe->fe_fpsr & ~FPSR_QTT) | (q << 16);
return r;
-
- Step9:
- fe->fe_f1.fp_class = FPC_ZERO;
- q++;
- q &= 0x7f;
- q |= (signQ << 7);
- fe->fe_fpframe->fpf_fpsr =
- fe->fe_fpsr =
- (fe->fe_fpsr & ~FPSR_QTT) | (q << 16);
- return &fe->fe_f1;
}
struct fpn *
