From dcd015359ef2cc2267c716c2660b57b1a386d9ef Mon Sep 17 00:00:00 2001 From: Maarten Boot <11140490+maarten-boot@users.noreply.github.com> Date: Fri, 30 Oct 2026 20:34:47 +0100 Subject: [PATCH 1/4] linuxppc emu: fmul is A-form (FRA x FRC); add fadds/fsubs/fmuls/fdivs fmul was decoded as an X-form instruction with FRB, but it is A-form with FRA and FRC, as the assembler (mach/powerpc/as/mach3.c) encodes it: a multiply with FRC other than f0 stopped the emulator with "unimplemented instruction", and one that matched computed FRA * FRB. The single-precision fadds, fsubs, fmuls and fdivs (opcode 59) were missing. They compute in double and round to float, which is exact for these operations. No test multiplied or did float arithmetic at run time on linuxppc before tests/plat/floats/constfold. Co-Authored-By: Claude Opus 5.5 --- plat/linuxppc/emu/instructions.dat | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) 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-->.....<136----->R setcr1(R, cpu.fpr[FRT] = d2b(-fabs(fpr(FRB)))); <63--><21------>R setcr1(R, cpu.fpr[FRT] = d2b(fpr(FRA) + fpr(FRB))); <63--><20------>R setcr1(R, cpu.fpr[FRT] = d2b(fpr(FRA) - fpr(FRB))); -<63--><25------>R setcr1(R, cpu.fpr[FRT] = d2b(fpr(FRA) * fpr(FRB))); +<63-->.....<25->R setcr1(R, cpu.fpr[FRT] = d2b(fpr(FRA) * fpr(FRC))); <63--><18------>R setcr1(R, cpu.fpr[FRT] = d2b(fpr(FRA) / fpr(FRB))); +<59-->.....<21->R setcr1(R, cpu.fpr[FRT] = d2b((float)(fpr(FRA) + fpr(FRB)))); +<59-->.....<20->R setcr1(R, cpu.fpr[FRT] = d2b((float)(fpr(FRA) - fpr(FRB)))); +<59-->.....<25->R setcr1(R, cpu.fpr[FRT] = d2b((float)(fpr(FRA) * fpr(FRC)))); +<59-->.....<18->R setcr1(R, cpu.fpr[FRT] = d2b((float)(fpr(FRA) / fpr(FRB)))); # Floating point comparisons From 3312e5d9bd14d2d987a19fc5903440c65c6b1872 Mon Sep 17 00:00:00 2001 From: Maarten Boot <11140490+maarten-boot@users.noreply.github.com> Date: Fri, 30 Oct 2026 17:31:07 +0100 Subject: [PATCH 2/4] flt_arith: sticky bit in add, sub, mul and div; flt_rnd() The front ends fold float constants with flt_arith, on 64-bit mantissas, and the back end rounds the result to the target type only at the end. To make that rounding exact, the operations now keep what does not fit as a sticky lowest bit (rounding to odd): the bits shifted out of the smaller operand of flt_add/flt_sub, the lower half of flt_mul's 128-bit product (was rounded half-up at 64 bits), flt_div's remainder (was dropped). flt_rnd(e, bits, emin) rounds to bits bits, to nearest and ties to even, with fewer bits below 2^emin as in a subnormal number. With operands rounded by flt_rnd, the result rounded by flt_rnd to the same precision is what IEEE arithmetic gives. Test: modules/src/flt_arith/ieeetest.c, part of all: 1,580,135 random double and float operations (subnormal results, exact results and ties included) against the host's IEEE arithmetic. The old add, mul and div give 3, 24 and 38 wrong results. Co-Authored-By: Claude Opus 5.5 --- build.py | 1 + modules/src/flt_arith/build.py | 15 ++- modules/src/flt_arith/flt_add.c | 36 ++++++- modules/src/flt_arith/flt_arith.3 | 23 ++++ modules/src/flt_arith/flt_arith.h | 10 ++ modules/src/flt_arith/flt_div.c | 8 ++ modules/src/flt_arith/flt_mul.c | 13 +-- modules/src/flt_arith/flt_rnd.c | 51 +++++++++ modules/src/flt_arith/ieeetest.c | 173 ++++++++++++++++++++++++++++++ 9 files changed, 317 insertions(+), 13 deletions(-) create mode 100644 modules/src/flt_arith/flt_rnd.c create mode 100644 modules/src/flt_arith/ieeetest.c diff --git a/build.py b/build.py index 2c94443d9..2a2110593 100644 --- a/build.py +++ b/build.py @@ -51,6 +51,7 @@ name="all", deps=( [".+compiler", "examples+all"] + + ([] if IS_WINDOWS else ["modules/src/flt_arith+test"]) + ( [] if IS_WINDOWS 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 #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 +#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/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 +#include +#include +#include +#include +#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; +} From 56beb86fcb82b7e225af380f468c922421e5fb82 Mon Sep 17 00:00:00 2001 From: Maarten Boot <11140490+maarten-boot@users.noreply.github.com> Date: Fri, 30 Oct 2026 17:44:45 +0100 Subject: [PATCH 3/4] cemcom.ansi: round float constants to their type, as the target does Float constants were kept and folded with flt_arith's 64-bit mantissa and rounded to the target type only when written out, so a folded constant could differ from the same computation at run time: 11-39% of folded + - * / results were 1 ulp off (0.1 + 0.2 gave the double nearest 0.3, not 0.30000000000000004). A float constant also kept 64 bits: (double)0.1f was the double 0.1, and 0.1f == 0.1 was true. Now every float constant is rounded to its type with flt_rnd(): a literal, a folded result, an integer converted to a floating type, and a conversion between floating types (which only changed the type). A 4-byte type is IEEE single and an 8-byte type IEEE double, subnormal numbers included. For a target whose floating point is not IEEE, the new -Vm gives the mantissa bits of an 8-byte type; pdpv7 uses -Vm56. flt_str2flt scaled the digits by multiplying with reciprocal powers of ten, which is not exact; with flt_mul's sticky bit even "1.0" came out as more than 1, and an exact tie like 9007199254740993.0 rounded the wrong way. Negative exponents now divide by the exact powers 10^0 to 10^27, and digits that do not fit in 64 bits set the sticky bit instead of being dropped. Co-Authored-By: Claude Opus 5.5 --- lang/cem/cemcom.ansi/arith.c | 7 ++++++- lang/cem/cemcom.ansi/cemcom.ansi.1 | 5 +++++ lang/cem/cemcom.ansi/expr.c | 2 ++ lang/cem/cemcom.ansi/fltcstoper.c | 28 ++++++++++++++++++++++++++++ lang/cem/cemcom.ansi/fltcstoper.h | 5 +++++ lang/cem/cemcom.ansi/options.c | 4 ++++ modules/src/flt_arith/flt_str2fl.c | 15 +++++++++++++-- plat/pdpv7/descr | 2 ++ 8 files changed, 65 insertions(+), 3 deletions(-) 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/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/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 From 5039af8bb9d987fe04bb69dbc76f22c3bec9ff9b Mon Sep 17 00:00:00 2001 From: Maarten Boot <11140490+maarten-boot@users.noreply.github.com> Date: Fri, 30 Oct 2026 18:49:07 +0100 Subject: [PATCH 4/4] tests/plat: constfold, floating constants folded as the target computes Folded floating constants compared with the same expressions computed at run time from volatile operands, and with exact values: + - * / of double and float, float precision of float constants ((double)0.1f, 0.1f == 0.1), 0.1 + 0.2 == 0.30000000000000004, ties rounded to even (9007199254740993.0, 16777217.0f, (float)16777217L) and subnormal results. In the floats set. With the compiler before the two previous commits, 20 of its 26 checks fail. Co-Authored-By: Claude Opus 5.5 --- tests/plat/build.py | 1 + tests/plat/floats/constfold,e,.c | 59 ++++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+) create mode 100644 tests/plat/floats/constfold,e,.c diff --git a/tests/plat/build.py b/tests/plat/build.py index 30bdb962b..74000b8b1 100644 --- a/tests/plat/build.py +++ b/tests/plat/build.py @@ -51,6 +51,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(); +}