diff --git a/build.py b/build.py
index 0e254e465..57b14cb63 100644
--- a/build.py
+++ b/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]
diff --git a/lang/cem/cemcom.ansi/arith.c b/lang/cem/cemcom.ansi/arith.c
index 36b8012b7..d0c2ccfd7 100644
--- a/lang/cem/cemcom.ansi/arith.c
+++ b/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;
diff --git a/lang/cem/cemcom.ansi/cemcom.ansi.1 b/lang/cem/cemcom.ansi/cemcom.ansi.1
index 3ca5de847..f3a6ea3c0 100644
--- a/lang/cem/cemcom.ansi/cemcom.ansi.1
+++ b/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
diff --git a/lang/cem/cemcom.ansi/expr.c b/lang/cem/cemcom.ansi/expr.c
index 7ce5661d5..cdf1f6e64 100644
--- a/lang/cem/cemcom.ansi/expr.c
+++ b/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)
diff --git a/lang/cem/cemcom.ansi/fltcstoper.c b/lang/cem/cemcom.ansi/fltcstoper.c
index ef6f8bccf..d35ba6820 100644
--- a/lang/cem/cemcom.ansi/fltcstoper.c
+++ b/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;
diff --git a/lang/cem/cemcom.ansi/fltcstoper.h b/lang/cem/cemcom.ansi/fltcstoper.h
index b48b5dfb1..700139bf6 100644
--- a/lang/cem/cemcom.ansi/fltcstoper.h
+++ b/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_ */
diff --git a/lang/cem/cemcom.ansi/options.c b/lang/cem/cemcom.ansi/options.c
index 4561837cd..349884a12 100644
--- a/lang/cem/cemcom.ansi/options.c
+++ b/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;
diff --git a/modules/src/flt_arith/build.py b/modules/src/flt_arith/build.py
index 4e2bf29c8..c2191d5c1 100644
--- a/modules/src/flt_arith/build.py
+++ b/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,6 +17,7 @@
         "./flt_nrm.c",
         "./b64_sft.c",
         "./flt_umin.c",
+        "./flt_rnd.c",
         "./flt_chk.c",
         "./split.c",
         "./flt_misc.h",
@@ -23,3 +25,14 @@
     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")
diff --git a/modules/src/flt_arith/flt_add.c b/modules/src/flt_arith/flt_add.c
index 21be5c2ed..602950f6e 100644
--- a/modules/src/flt_arith/flt_add.c
+++ b/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);
 }
diff --git a/modules/src/flt_arith/flt_arith.3 b/modules/src/flt_arith/flt_arith.3
index 8d09c656e..01f0a2e1c 100644
--- a/modules/src/flt_arith/flt_arith.3
+++ b/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
diff --git a/modules/src/flt_arith/flt_arith.h b/modules/src/flt_arith/flt_arith.h
index 976caf599..ac6b50a3b 100644
--- a/modules/src/flt_arith/flt_arith.h
+++ b/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.
diff --git a/modules/src/flt_arith/flt_div.c b/modules/src/flt_arith/flt_div.c
index 7929b2e7d..92f0ec612 100644
--- a/modules/src/flt_arith/flt_div.c
+++ b/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);
diff --git a/modules/src/flt_arith/flt_mul.c b/modules/src/flt_arith/flt_mul.c
index daf0c973f..ecee7ae9e 100644
--- a/modules/src/flt_arith/flt_mul.c
+++ b/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);
 }
diff --git a/modules/src/flt_arith/flt_rnd.c b/modules/src/flt_arith/flt_rnd.c
new file mode 100644
index 000000000..169ceef91
--- /dev/null
+++ b/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);
+}
diff --git a/modules/src/flt_arith/flt_str2fl.c b/modules/src/flt_arith/flt_str2fl.c
index f3dd7b7d6..d13343fa7 100644
--- a/modules/src/flt_arith/flt_str2fl.c
+++ b/modules/src/flt_arith/flt_str2fl.c
@@ -204,7 +204,12 @@ static void add_exponent(flt_arith *e, int exp)
 	divsz = exp / SMALLSZ;
 	modsz = exp % SMALLSZ;
 	if (!status) status = flt_status;
-	flt_mul(e, (neg ? r_10pow : s10pow) + modsz, &x);
+	/* s10pow holds 10 ** i exactly: dividing by it is exact up to the
+	   sticky bit, multiplying by r_10pow is not */
+	if (neg)
+		flt_div(e, s10pow + modsz, &x);
+	else
+		flt_mul(e, s10pow + modsz, &x);
 	if (!status) status = flt_status;
 	while (divsz >= (int)BIGSZ) {
 		flt_mul(&x, neg ? &r_big_10pow[BIGSZ-1] : &big_10pow[BIGSZ-1],&x);
@@ -221,6 +226,7 @@ void flt_str2flt(const char *s, flt_arith *e)
 	int	c;
 	int		dotseen = 0;
 	int		digitseen = 0;
+	int		dropped = 0;	/* a non-zero digit did not fit */
 	int		exp = 0;
 
 	while (isspace(*s)) s++;
@@ -252,7 +258,10 @@ void flt_str2flt(const char *s, flt_arith *e)
 			a1.flt_l_32 = c - '0';
 			flt_b64_add(&(e->flt_mantissa), &a1);
 		}
-		else exp++;
+		else {
+			exp++;
+			if (c != '0') dropped = 1;
+		}
 		if (dotseen) exp--;
 	}
 	if (! digitseen) {
@@ -281,6 +290,8 @@ void flt_str2flt(const char *s, flt_arith *e)
 	e->flt_exp = 63;
 	flt_nrm(e);
 	add_exponent(e, exp);
+	/* the digits that did not fit as a sticky bit, for flt_rnd() */
+	if (dropped) e->m2 |= 1;
 	flt_chk(e);
 }
 
diff --git a/modules/src/flt_arith/ieeetest.c b/modules/src/flt_arith/ieeetest.c
new file mode 100644
index 000000000..3faba36ad
--- /dev/null
+++ b/modules/src/flt_arith/ieeetest.c
@@ -0,0 +1,173 @@
+/*
+ * flt_add, flt_sub, flt_mul and flt_div followed by flt_rnd must give what
+ * IEEE arithmetic gives: checked against the host's double and float on
+ * random operands, subnormal results and exact ties included.
+ *
+ *	ieeetest [count]
+ */
+
+#include <math.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "flt_arith.h"
+
+static uint64_t state = 0x2545F4914F6CDD1DULL;
+
+static uint64_t rnd(void)
+{
+	/* xorshift64*: the same numbers on every host */
+	state ^= state >> 12;
+	state ^= state << 25;
+	state ^= state >> 27;
+	return state * 0x2545F4914F6CDD1DULL;
+}
+
+static void d2flt(double d, flt_arith* e)
+{
+	int exp;
+	double f = frexp(fabs(d), &exp); /* d = f * 2^exp, f in [0.5, 1) */
+	uint64_t m = (uint64_t)ldexp(f, 64);
+
+	e->flt_sign = d < 0;
+	e->flt_exp = exp - 1;
+	e->flt_mantissa.flt_h_32 = (uint32_t)(m >> 32);
+	e->flt_mantissa.flt_l_32 = (uint32_t)m;
+	if (d == 0)
+		e->flt_exp = e->flt_sign = 0;
+}
+
+static double flt2d(flt_arith* e)
+{
+	uint64_t m = ((uint64_t)e->flt_mantissa.flt_h_32 << 32) | e->flt_mantissa.flt_l_32;
+	double d = ldexp((double)m, e->flt_exp - 63);
+
+	return e->flt_sign ? -d : d;
+}
+
+/* A random double with an exponent in [-emax, emax]; bits mantissa bits. */
+static double randnum(int emax, int bits)
+{
+	uint64_t m = rnd() >> (64 - bits) | (uint64_t)1 << (bits - 1);
+	double d = ldexp((double)m, (int)(rnd() % (2 * emax + 1)) - emax - bits + 1);
+
+	return (rnd() & 1) ? -d : d;
+}
+
+static long checked, failures;
+
+static void check(int isfloat, char op, double a, double b)
+{
+	volatile double dwant;
+	volatile float fwant;
+	double want, got;
+	flt_arith e1, e2, r;
+
+	if (isfloat)
+	{
+		volatile float fa = (float)a, fb = (float)b;
+
+		a = fa;
+		b = fb;
+		switch (op)
+		{
+			case '+': fwant = fa + fb; break;
+			case '-': fwant = fa - fb; break;
+			case '*': fwant = fa * fb; break;
+			default: fwant = fa / fb; break;
+		}
+		want = fwant;
+	}
+	else
+	{
+		switch (op)
+		{
+			case '+': dwant = a + b; break;
+			case '-': dwant = a - b; break;
+			case '*': dwant = a * b; break;
+			default: dwant = a / b; break;
+		}
+		want = dwant;
+	}
+	/* flt_arith has a wider range and no infinity */
+	if (!isfinite(a) || !isfinite(b) || !isfinite(want) || (op == '/' && b == 0))
+		return;
+
+	d2flt(a, &e1);
+	d2flt(b, &e2);
+	switch (op)
+	{
+		case '+': flt_add(&e1, &e2, &r); break;
+		case '-': flt_sub(&e1, &e2, &r); break;
+		case '*': flt_mul(&e1, &e2, &r); break;
+		default: flt_div(&e1, &e2, &r); break;
+	}
+	if (isfloat)
+		flt_rnd(&r, 24, -126);
+	else
+		flt_rnd(&r, 53, -1022);
+	got = flt2d(&r);
+	checked++;
+	/* flt_arith has no negative zero */
+	if (want == 0 ? got != 0 : got != want)
+	{
+		if (failures++ < 10)
+			printf("%s %a %c %a: want %a, got %a\n", isfloat ? "float" : "double", a, op, b,
+			    want, got);
+	}
+}
+
+int main(int argc, char* argv[])
+{
+	long count = argc > 1 ? atol(argv[1]) : 200000;
+	static const char ops[] = "+-*/";
+	long i;
+	int k, isfloat;
+
+	for (i = 0; i < count; i++)
+	{
+		for (isfloat = 0; isfloat <= 1; isfloat++)
+		{
+			int emax = isfloat ? 127 : 1023;
+			int bits = isfloat ? 24 : 53;
+			double a, b;
+
+			for (k = 0; k < 4; k++)
+			{
+				switch (rnd() % 5)
+				{
+					case 0: /* any magnitudes */
+						a = randnum(emax, bits);
+						b = randnum(emax, bits);
+						break;
+					case 1: /* close magnitudes: cancellation, carries */
+						a = randnum(4, bits);
+						b = randnum(4, bits);
+						break;
+					case 2: /* short mantissas: exact results and ties */
+						a = randnum(emax / 8, 1 + (int)(rnd() % 12));
+						b = randnum(emax / 8, 1 + (int)(rnd() % 12));
+						break;
+					case 3: /* a tiny, b near 1: subnormal products */
+						a = ldexp(randnum(0, bits), -(emax - 1) + (int)(rnd() % 8));
+						b = ldexp(randnum(0, bits), -(int)(rnd() % (bits + 4)));
+						break;
+					default: /* both tiny: subnormal sums */
+						a = ldexp(randnum(0, bits), -(emax - 1) + (int)(rnd() % 4));
+						b = ldexp(randnum(0, bits), -(emax - 1) + (int)(rnd() % 4));
+						break;
+				}
+				check(isfloat, ops[k], a, b);
+			}
+		}
+	}
+	if (failures)
+	{
+		printf("@@FAIL flt_arith: %ld of %ld results not as in IEEE arithmetic\n",
+		    failures, checked);
+		return 1;
+	}
+	printf("flt_arith: %ld results as in IEEE arithmetic\n", checked);
+	return 0;
+}
diff --git a/plat/linuxppc/emu/instructions.dat b/plat/linuxppc/emu/instructions.dat
index 5bf4411c6..e1aec6494 100644
--- a/plat/linuxppc/emu/instructions.dat
+++ b/plat/linuxppc/emu/instructions.dat
@@ -197,8 +197,12 @@
 <63--><FRT>.....<FRB><136----->R setcr1(R, cpu.fpr[FRT] = d2b(-fabs(fpr(FRB))));
 <63--><FRT><FRA><FRB><21------>R setcr1(R, cpu.fpr[FRT] = d2b(fpr(FRA) + fpr(FRB)));
 <63--><FRT><FRA><FRB><20------>R setcr1(R, cpu.fpr[FRT] = d2b(fpr(FRA) - fpr(FRB)));
-<63--><FRT><FRA><FRB><25------>R setcr1(R, cpu.fpr[FRT] = d2b(fpr(FRA) * fpr(FRB)));
+<63--><FRT><FRA>.....<FRC><25->R setcr1(R, cpu.fpr[FRT] = d2b(fpr(FRA) * fpr(FRC)));
 <63--><FRT><FRA><FRB><18------>R setcr1(R, cpu.fpr[FRT] = d2b(fpr(FRA) / fpr(FRB)));
+<59--><FRT><FRA><FRB>.....<21->R setcr1(R, cpu.fpr[FRT] = d2b((float)(fpr(FRA) + fpr(FRB))));
+<59--><FRT><FRA><FRB>.....<20->R setcr1(R, cpu.fpr[FRT] = d2b((float)(fpr(FRA) - fpr(FRB))));
+<59--><FRT><FRA>.....<FRC><25->R setcr1(R, cpu.fpr[FRT] = d2b((float)(fpr(FRA) * fpr(FRC))));
+<59--><FRT><FRA><FRB>.....<18->R setcr1(R, cpu.fpr[FRT] = d2b((float)(fpr(FRA) / fpr(FRB))));
 
 # Floating point comparisons
 
diff --git a/plat/pdpv7/descr b/plat/pdpv7/descr
index 1fcc21161..f430c6b3c 100644
--- a/plat/pdpv7/descr
+++ b/plat/pdpv7/descr
@@ -21,6 +21,8 @@ var CPP_F=-D__unix
 var ALIGN=-a0:2 -a1:2 -a2:2 -a3:2
 
 var C_INCLUDES=-I{PLATFORMDIR}/include -I{EM}/share/ack/include/ansi
+# PDP-11 floating point: double has 56 mantissa bits, no subnormal numbers
+var CC_ALIGN=-Vm56
 
 name be
 	from .m.g
diff --git a/tests/plat/build.py b/tests/plat/build.py
index 3d08d0377..8bd5c387a 100644
--- a/tests/plat/build.py
+++ b/tests/plat/build.py
@@ -61,6 +61,7 @@
         "tests/plat/core/xorv,e,.e",
     ],
     "floats": [
+        "tests/plat/floats/constfold,e,.c",
         "tests/plat/floats/doublecmp,e,.c",
         "tests/plat/floats/from_d_to_si,e,.c",
         "tests/plat/floats/from_d_to_ui,e,.c",
diff --git a/tests/plat/floats/constfold,e,.c b/tests/plat/floats/constfold,e,.c
new file mode 100644
index 000000000..7285bb230
--- /dev/null
+++ b/tests/plat/floats/constfold,e,.c
@@ -0,0 +1,59 @@
+#include "test.h"
+
+/* Floating constants folded by the compiler must have the value the target
+ * computes: each operand rounded to its type, each result rounded to its
+ * type, to nearest and ties to even (IEEE). */
+
+/* Operands in globals to defeat constant folding. (Not volatile: the
+ * powerpc ncg table has no pattern for the nop cemcom.ansi emits after
+ * a volatile access.) */
+double d243 = 2.43, d622 = 622.780, d483 = 483.667, d807 = 807.214;
+double d866 = 866.963, d603 = 603.873, d975 = 975.296, d16 = 16.764;
+double dp1 = 0.1, dp2 = 0.2, d3 = 3.0;
+float f243 = 2.43f, f622 = 622.780f, f24 = 24.41f, f58 = 58.848f;
+float f975 = 975.296f, f16 = 16.764f, fp1 = 0.1f;
+
+/* The folded values, stored so that comparing them is not folded. */
+double d;
+float f;
+
+/* Bypasses the CRT, so there's no stdio. */
+void _m_a_i_n(void)
+{
+    /* folded as computed at run time; these were 1 ulp off */
+    d = 0.1 + 0.2;           ASSERT(d == dp1 + dp2);
+    d = 2.43 + 622.780;      ASSERT(d == d243 + d622);
+    d = 975.296 + 16.764;    ASSERT(d == d975 + d16);
+    d = 866.963 - 603.873;   ASSERT(d == d866 - d603);
+    d = 2.43 * 622.780;      ASSERT(d == d243 * d622);
+    d = 2.43 / 622.780;      ASSERT(d == d243 / d622);
+    d = 483.667 / 807.214;   ASSERT(d == d483 / d807);
+    f = 24.41f + 58.848f;    ASSERT(f == f24 + f58);
+    f = 975.296f - 16.764f;  ASSERT(f == f975 - f16);
+    f = 2.43f * 622.780f;    ASSERT(f == f243 * f622);
+    f = 2.43f / 622.780f;    ASSERT(f == f243 / f622);
+
+    /* a float constant has float precision */
+    d = (double)0.1f;        ASSERT(d == (double)fp1);
+    d = 0.1f * 3.0;          ASSERT(d == fp1 * d3);
+    d = 0.1f == 0.1;         ASSERT(d == 0.0);
+    d = 0.5f == 0.5;         ASSERT(d == 1.0);
+
+    /* exact values */
+    d = 0.1 + 0.2;           ASSERT(d == 0.30000000000000004);
+                             ASSERT(d != 0.3);
+
+    /* ties round to even */
+    d = 9007199254740993.0;  ASSERT(d == 9007199254740992.0);
+    d = 9007199254740995.0;  ASSERT(d == 9007199254740996.0);
+    f = 16777217.0f;         ASSERT(f == 16777216.0f);
+    f = (float)16777217L;    ASSERT(f == 16777216.0f);
+    f = (float)16777219L;    ASSERT(f == 16777220.0f);
+
+    /* subnormal results: half the smallest one is a tie, to 0 */
+    d = 4.9406564584124654e-324 / 2.0;  ASSERT(d == 0.0);
+    d = 4.9406564584124654e-324 * 1.5;  ASSERT(d == 9.8813129168249309e-324);
+    d = 2.2250738585072014e-308 / 4.0;  ASSERT(d == 5.5626846462680035e-309);
+
+    finished();
+}