Floating constants folded as the target computes them

build.py

@@ -51,8 +51,7 @@
     name="all",
     deps=(
         [".+compiler", "examples+all"]
-        + ([] if IS_WINDOWS else ["lang/cem/cpp.ansi+test"])
-        + (
+        + ([] if IS_WINDOWS else ["modules/src/flt_arith+test", "lang/cem/cpp.ansi+test"])
             []
             if IS_WINDOWS
             else [f"plat/{p}/tests" for p in TEST_PLATS if p in PLATS]

lang/cem/cemcom.ansi/arith.c

@@ -28,6 +28,7 @@
 #include "field.h"
 #include "mes.h"
 #include "cstoper.h"
+#include "fltcstoper.h"
 #include "ch3bin.h"
 #include "ch3.h"
 #include "error.h"
@@ -440,6 +441,7 @@ void int2float(struct expr** expp, struct type* tp)
 		exp->ex_type = tp;
 		exp->ex_class = Float;
 		flt_arith2flt((arith)exp->VL_VALUE, &(exp->FL_ARITH), uns);
+		fltcstround(exp);
 	}
 	else
 	{
@@ -491,11 +493,14 @@ void float2float(struct expr** expp, struct type* tp)
 	/*	The expression *expp, which is of some floating type, is
 	    converted to the floating type tp.
 	    There is no need for an explicit conversion operator
-	    if the expression is a constant.
+	    if the expression is a constant: it is rounded to tp.
 	*/
 
 	if (is_fp_cst(*expp))
+	{
 		(*expp)->ex_type = tp;
+		fltcstround(*expp);
+	}
 	else
 	{
 		fp_used = 1;

lang/cem/cemcom.ansi/cemcom.ansi.1

@@ -78,6 +78,11 @@ and the \fIn\fR parameter for the alignment of that type.
 Absence of \fIm\fR or \fIn\fR causes the default value to be retained.
 To specify that the bitfields should be right adjusted instead of the
 default left adjustment, specify \fBr\fR as \fIc\fR parameter.
+Floating constants are computed as the target computes them: rounded to
+IEEE single (4 bytes) and double (8 bytes) precision. For a target
+whose floating point is not IEEE, \fBm\fR as \fIc\fR parameter gives
+the number of mantissa bits of an 8-byte floating type in \fIm\fR
+(a 4-byte one has 24); there are no subnormal numbers then.
 .IP \fB\-w\fR
 suppress warning messages.
 .IP \fB\-s\fR

lang/cem/cemcom.ansi/expr.c

@@ -26,6 +26,7 @@
 #include "sizes.h"
 #include "level.h"
 #include "cstoper.h"
+#include "fltcstoper.h"
 #include "error.h"
 
 extern char options[];
@@ -260,6 +261,7 @@ void float2expr(struct expr* expr)
 	assert(flt_status != FLT_NOFLT);
 	if (flt_status == FLT_OVFL)
 		expr_warning(expr, "internal floating point overflow");
+	fltcstround(expr);
 }
 
 struct expr* intexpr(arith ivalue, int fund)

lang/cem/cemcom.ansi/fltcstoper.c

@@ -24,6 +24,33 @@
 
 extern int ResultKnown;
 
+int flt_dbl_bits; /* 0: IEEE */
+
+void fltcstround(struct expr* expr)
+{
+	/*	The constant expr is rounded to the precision of its type, as
+	    the target would have it: 4-byte and 8-byte types are IEEE
+	    single and double (with their subnormal numbers), unless -Vm
+	    gave the mantissa bits of a non-IEEE 8-byte type (no subnormal
+	    numbers then). Other sizes are left alone.
+	*/
+	switch ((int)expr->ex_type->tp_size)
+	{
+		case 4:
+			if (flt_dbl_bits)
+				flt_rnd(&expr->FL_ARITH, 24, -16384);
+			else
+				flt_rnd(&expr->FL_ARITH, 24, -126);
+			break;
+		case 8:
+			if (flt_dbl_bits)
+				flt_rnd(&expr->FL_ARITH, flt_dbl_bits, -16384);
+			else
+				flt_rnd(&expr->FL_ARITH, 53, -1022);
+			break;
+	}
+}
+
 void fltcstbin(struct expr** expp, int oper, struct expr* expr)
 {
 	/*	The operation oper is performed on the constant
@@ -113,6 +140,7 @@ void fltcstbin(struct expr** expp, int oper, struct expr* expr)
 	else
 	{
 		(*expp)->FL_ARITH = o1;
+		fltcstround(*expp);
 	}
 	(*expp)->ex_flags |= expr->ex_flags;
 	(*expp)->ex_flags &= ~EX_PARENS;

lang/cem/cemcom.ansi/fltcstoper.h

@@ -11,5 +11,10 @@
 struct expr;
 
 void fltcstbin(struct expr** expp, int oper, struct expr* expr);
+void fltcstround(struct expr* expr);
+
+/* Mantissa bits of an 8-byte floating type when the target's floating point
+   is not IEEE (-Vm), else 0. */
+extern int flt_dbl_bits;
 
 #endif /* FLTCSTOPER_H_ */

lang/cem/cemcom.ansi/options.c

@@ -18,6 +18,7 @@
 #include "arith.h"
 #include "sizes.h"
 #include "align.h"
+#include "fltcstoper.h"
 #include "error.h"
 
 char options[128]; /* one for every char	*/
@@ -195,6 +196,9 @@ void do_option(char* text)
 						if (algn != 0)
 							pointer_align = algn;
 						break;
+					case 'm': /* mantissa bits of a non-IEEE double */
+						flt_dbl_bits = sz;
+						break;
 					case 'r': /* adjust bitfields right	*/
 #ifndef NOBITFIELD
 						options['r'] = 1;

modules/src/flt_arith/build.py

@@ -1,4 +1,5 @@
-from build.c import clibrary
+from build.c import clibrary, cprogram
+from build.utils import test
 
 
 clibrary(
@@ -16,10 +17,22 @@
         "./flt_nrm.c",
         "./b64_sft.c",
         "./flt_umin.c",
+        "./flt_rnd.c",
         "./flt_chk.c",
         "./split.c",
         "./flt_misc.h",
     ],
     hdrs={"flt_arith.h": "./flt_arith.h", "flt_misc.h": "./flt_misc.h"},
     deps=["h", "modules/h"],
 )
+
+# flt_add, flt_sub, flt_mul and flt_div rounded with flt_rnd against the
+# host's IEEE double and float arithmetic.
+cprogram(
+    name="ieeetest",
+    srcs=["./ieeetest.c"],
+    deps=[".+flt_arith"],
+    ldflags=["-lm"],
+)
+
+test(name="test", command=".+ieeetest")

modules/src/flt_arith/flt_add.c

@@ -5,12 +5,33 @@
 
 /* $Id$ */
 
+#include <stdint.h>
 #include "flt_misc.h"
 
+/* Shifts the mantissa m right by n bits; returns whether any of the bits
+   shifted out was set. */
+static int sft_lost(struct flt_mantissa *m, int n)
+{
+	uint32_t mask;
+	int lost;
+
+	if (n > 63)
+		lost = (m->flt_h_32 | m->flt_l_32) != 0;
+	else if (n >= 32) {
+		mask = n == 32 ? 0 : (1UL << (n - 32)) - 1;
+		lost = m->flt_l_32 != 0 || (m->flt_h_32 & mask) != 0;
+	}
+	else
+		lost = (m->flt_l_32 & ((1UL << n) - 1)) != 0;
+	flt_b64_sft(m, n);
+	return lost;
+}
+
 void flt_add(flt_arith *e1, flt_arith *e2,flt_arith *res)
 {
 	flt_arith ce1, ce2;
 	int diff;
+	int lost = 0;	/* bits of the smaller operand shifted out */
 
 	flt_status = 0;
 	if ((e2->m1 | e2->m2) == 0L) {
@@ -31,11 +52,11 @@ void flt_add(flt_arith *e1, flt_arith *e2,flt_arith *res)
 	if (diff < 0) {
 		diff = -diff;
 		e2->flt_exp += diff;
-		flt_b64_sft(&(e2->flt_mantissa), diff);
+		lost = sft_lost(&(e2->flt_mantissa), diff);
 	}
 	else if (diff > 0) {
 		e1->flt_exp += diff;
-		flt_b64_sft(&(e1->flt_mantissa), diff);
+		lost = sft_lost(&(e1->flt_mantissa), diff);
 	}
 	if (e1->flt_sign != e2->flt_sign) {
 		/* e2 + e1 = e2 - (-e1) */
@@ -56,15 +77,24 @@ void flt_add(flt_arith *e1, flt_arith *e2,flt_arith *res)
                 	e2->m2 -= e1->m2;
 			*res = *e2;
         	}
+		if (lost) {
+			/* the exact difference is a little less */
+			if (res->m2-- == 0)
+				res->m1--;
+		}
 	}
 	else {
 		*res = *e2;
 		if (flt_b64_add(&res->flt_mantissa,&e1->flt_mantissa)) {/* addition carry */
-			flt_b64_sft(&res->flt_mantissa, 1);
+			lost |= sft_lost(&res->flt_mantissa, 1);
 			res->m1 |= 0x80000000L;	/* set max bit	*/
 			res->flt_exp++;		/* increase the exponent */
 		}
 	}
+	/* Bits that did not fit are kept as a sticky lowest bit (rounding to
+	   odd), so that flt_rnd() can round the result once more correctly. */
+	if (lost)
+		res->m2 |= 1;
 	flt_nrm(res);
 	flt_chk(res);
 }

modules/src/flt_arith/flt_arith.3

@@ -38,6 +38,8 @@ extern int	flt_status;
 .PP
 .B void flt_umin(flt_arith *e)
 .PP
+.B void flt_rnd(flt_arith *e, int bits, int emin)
+.PP
 .B void flt_modf(flt_arith *e, flt_arith *ipart, flt_arith *fpart)
 .PP
 .B int flt_cmp(flt_arith *e1, flt_arith *e2)
@@ -99,6 +101,27 @@ negates the number indicated by
 and stores the result indirectly through
 .IR e .
 .PP
+.B flt_rnd
+rounds the number indicated by
+.I e
+to
+.I bits
+bits of mantissa, to nearest and ties to even; below 2 to the power
+.I emin
+fewer bits are kept, as in a subnormal number.
+.BR flt_add ,
+.BR flt_sub ,
+.B flt_mul
+and
+.B flt_div
+keep the bits that do not fit in their 64-bit mantissa as a sticky lowest
+bit: with operands rounded by
+.B flt_rnd
+(at most 62 bits), the result rounded by
+.B flt_rnd
+to the same precision is the correctly rounded result, as the target
+computes it (53 and \-1022 for an IEEE double, 24 and \-126 for a float).
+.PP
 .B flt_modf
 splits the number indicated by
 .I e

modules/src/flt_arith/flt_arith.h

@@ -99,6 +99,16 @@ arith flt_flt2arith(flt_arith *e, int uns);
 /** Negates the number indicated by `e` and stores the result
  *  in `e`. */
 void flt_umin(flt_arith *e);
+/** Rounds the number indicated by `e` to `bits` bits of mantissa, to
+ *  nearest and ties to even. Below 2^`emin` fewer bits are kept, as in a
+ *  subnormal number. flt_add(), flt_sub(), flt_mul() and flt_div() keep
+ *  the bits that do not fit in the 64-bit mantissa as a sticky lowest
+ *  bit: with operands rounded by flt_rnd() (at most 62 bits), the result
+ *  rounded by flt_rnd() to the same precision is the correctly rounded
+ *  result, as the target computes it. E.g. bits 53, emin -1022 for an
+ *  IEEE double, 24 and -126 for a float.
+ */
+void flt_rnd(flt_arith *e, int bits, int emin);
 /** Shifts the mantissa m |n| bits left or right, depending on the sign
  *  of n. If n is negative, it is a left-shift; If n is positive,
  *  it is a right shift.

modules/src/flt_arith/flt_div.c

@@ -112,6 +112,14 @@ void flt_div(flt_arith *e1,flt_arith *e2,flt_arith *res)
 	}
 	res->m1 = result[0];
 	res->m2 = result[1];
+	/* a remainder as a sticky lowest bit (rounding to odd), so that
+	   flt_rnd() can round the quotient once more correctly */
+	for (j = 0; j < 8; j++) {
+		if (u[j]) {
+			res->m2 |= 1;
+			break;
+		}
+	}
 
 	flt_nrm(res);
 	flt_chk(res);

modules/src/flt_arith/flt_mul.c

@@ -62,14 +62,9 @@ void flt_mul(flt_arith *e1, flt_arith *e2, flt_arith *res)
 	 */
 	res->m1 = ((uint32_t)result[0] << 16) + result[1];
 	res->m2 = ((uint32_t)result[2] << 16) + result[3];
-	if (result[4] & 0x8000) {
-		if (++res->m2 == 0) {
-			res->m2 = 0;
-			if (++res->m1 == 0) {
-				res->m1 = 0x80000000;
-				res->flt_exp++;
-			}
-		}
-	}
+	/* the lower half as a sticky lowest bit (rounding to odd), so that
+	   flt_rnd() can round the product once more correctly */
+	if (result[4] | result[5] | result[6] | result[7])
+		res->m2 |= 1;
 	flt_chk(res);
 }

modules/src/flt_arith/flt_rnd.c

@@ -0,0 +1,51 @@
+/* Rounding to the precision of a target type. */
+
+#include <stdint.h>
+#include "flt_misc.h"
+
+void
+flt_rnd(flt_arith *e, int bits, int emin)
+{
+	/*	Round e to bits bits, to nearest, ties to even; below 2^emin
+		fewer bits are kept, as in a subnormal number.
+		The value is m * 2^(flt_exp - 63), m the 64-bit mantissa.
+	*/
+	uint64_t m = ((uint64_t)e->m1 << 32) | e->m2;
+	uint64_t keep, rest, half;
+	int drop;
+
+	if (m == 0)
+		return;
+	if (e->flt_exp < emin)
+		bits -= emin - e->flt_exp;
+	drop = 64 - bits;	/* the number of bits to round off */
+	if (drop <= 0)
+		return;
+	if (drop > 64) {
+		/* less than half the smallest unit */
+		keep = 0;
+	}
+	else if (drop == 64) {
+		/* m is at least half a unit; exactly half rounds to 0 */
+		keep = m > ((uint64_t)1 << 63);
+	}
+	else {
+		keep = m >> drop;
+		rest = m & (((uint64_t)1 << drop) - 1);
+		half = (uint64_t)1 << (drop - 1);
+		if (rest > half || (rest == half && (keep & 1)))
+			keep++;
+	}
+	if (keep == 0) {
+		e->m1 = e->m2 = 0;
+		e->flt_exp = 0;
+		e->flt_sign = 0;
+		return;
+	}
+	/* keep units of 2^(flt_exp - 63 + drop): normalize again */
+	e->m1 = (uint32_t)(keep >> 32);
+	e->m2 = (uint32_t)keep;
+	e->flt_exp += drop;
+	flt_nrm(e);
+	flt_chk(e);
+}