diff --git a/doc/architectural/x86-long-double.md b/doc/architectural/x86-long-double.md new file mode 100644 index 00000000000..2901e915998 --- /dev/null +++ b/doc/architectural/x86-long-double.md @@ -0,0 +1,182 @@ +\page x86-long-double Modelling x86 80-bit extended `long double` + +\author Michael Tautschnig + +# Summary + +On x86 hardware, the `long double` C type is the x87 80-bit extended +precision floating-point format, even though it is stored in 12 bytes +(i386) or 16 bytes (x86_64) of memory. Most other modern targets use +either an IEEE 754 `binary64` (e.g. AArch64 macOS, Windows MSVC) or an +IEEE 754 `binary128` (e.g. AArch64 Linux, ppc64le) layout for `long +double`. CBMC must model the byte layout of `long double` faithfully +because programs do read `long double` bytes via union-based bit +twiddling -- most prominently the macOS SDK's `` inline helpers +for `signbit`, `__inline_isnan`, etc. + +This document describes how CBMC encodes the format and where the code +paths live. + +# Hardware reference + +The 80-bit extended format is documented in: + +- Intel® 64 and IA-32 Architectures Software Developer's Manual, vol. 1, + §4.2.2 ("Floating-Point Data Types"). +- The IA-32 ABI for Linux ("System V Application Binary Interface ‒ + i386 Architecture Processor Supplement"), §3.1.2. + +The relevant facts are: + +- The value occupies 80 bits. Storage padding extends it to 96 bits on + i386 (`alignof == 4`) and 128 bits on x86_64 (`alignof == 16`). +- Unlike IEEE binary formats, the leading "integer bit" of the + significand (the J-bit) is **explicit**: it is stored in the + encoding, rather than being implicit-1 as in IEEE. +- Exponent width is 15 bits with bias 16383. + +The on-disk byte layout, low to high (LSB first): + +``` +bits 0 .. 62 : explicit fraction (62 bits, no implicit leading 1) +bit 63 : explicit integer bit (J-bit) +bits 64 .. 78 : biased 15-bit exponent +bit 79 : sign bit +bits 80 .. 127 : storage padding (zero on Linux/macOS) +``` + +For example, on real macOS-15 x86_64 hardware (Xcode 16.4, Apple clang +17.0.0) the 16-byte little-endian dump of `1.0L` is: + +``` +00 00 00 00 00 00 00 80 ff 3f 00 00 00 00 00 00 +``` + +with bytes 0-7 holding the mantissa `0x8000000000000000` (only the +explicit integer bit set), bytes 8-9 holding the biased exponent +`0x3FFF = 16383`, and bytes 10-15 being storage padding. Negation flips +bit 79 of the value, giving `... ff bf ...`. + +# Modelling + +`ieee_float_spect` (in `src/util/ieee_float.h`) records the format: + +- `f` -- fraction width (63 for x86 extended). +- `e` -- exponent width (15 for x86 extended). +- `x86_extended` -- when true, the encoding has an explicit integer bit + at position `f`, so the value width is `f + e + 2` (sign + exponent + + J-bit + fraction) rather than `f + e + 1`. +- `storage_width_bits` -- when non-zero, the storage container has more + bits than the value; the high bits are padding zeros. + +Three factory methods are provided for the x86 case: + +- `ieee_float_spect::x86_80()` -- the bare 80-bit value, used internally + for solver-internal arithmetic. +- `ieee_float_spect::x86_96()` -- 80-bit value in 96-bit storage (i386 + `long double`). +- `ieee_float_spect::x86_128()` -- 80-bit value in 128-bit storage + (x86_64 `long double` on Linux, macOS, FreeBSD). + +The selection happens in `long_double_type()` (in +`src/util/c_types.cpp`) based on `config.ansi_c.long_double_width` and +`config.ansi_c.arch`. + +# Encoding paths + +Floating-point values flow through three layers, and each must agree on +the byte layout: + +1. **Constants** -- `ieee_float_valuet::pack` / `unpack` in + `src/util/ieee_float.cpp` produce the `mp_integer` bit pattern that + the type checker hands off as a constant initialiser. +2. **Expression-level solver** -- `float_bvt` in + `src/solvers/floatbv/float_bv.cpp` lowers floating-point expressions + to bit-vector expressions before flattening. +3. **SAT-level solver** -- `float_utilst` in + `src/solvers/floatbv/float_utils.cpp` produces literals/clauses for + the SAT back-end directly. + +For the x86 extended layout each layer: + +- Places the sign bit at `value_width() - 1` (i.e. bit 79), **not** at + the top of the storage container. +- Places the biased exponent at bits `[f + 1, f + 1 + e)`. +- Places the explicit integer bit at bit `f` (== 63). +- Reads the integer bit out of the encoding rather than synthesising a + hidden bit during `unpack`. +- Zero-pads the high storage bits during `pack`. + +The `boolbvt::convert_rest` handler for `ID_sign` (in +`src/solvers/flattening/boolbv.cpp`) similarly reads the sign from +`value_width() - 1` for `floatbv` operands. + +# SMT2 back-ends + +The SMT-LIB FloatingPoint theory only covers IEEE 754 binary formats, +which all use an *implicit* leading integer bit. The x87 extended +format stores that bit *explicitly*, and the value lives in a 96- or +128-bit container whose high bits are padding. Encoding an +x86-extended value as `(_ FloatingPoint 15 64)` would discard the +explicit J-bit and the storage padding, so it is not a faithful +encoding. + +`smt2_convt::use_FPA_for_type(typet)` (in +`src/solvers/smt2/smt2_conv.cpp`) returns `false` for x86-extended +`floatbv` types regardless of the global `use_FPA_theory` flag, and +`true` otherwise (delegating to `use_FPA_theory`). All dispatch sites +in `smt2_conv.cpp` consult this helper, which means: + +- The SMT-LIB sort emitted for an x86-extended type is always + `(_ BitVec storage_width)` (128 on x86_64, 96 on i386). +- All FP operations on x86-extended operands lower through the + existing `convert_floatbv` path, which emits the `float_bv.*` helper + functions defined elsewhere in the same module. +- Other floating-point types (binary32, binary64, binary128) continue + to use the SMT-LIB FloatingPoint theory on solvers that support it + (currently Z3, CVC5, Bitwuzla, and the cprover-smt2 solver). + +The new SMT2 back-end (`--incremental-smt2-solver`, in +`src/solvers/smt2_incremental/`) does not yet implement the +byte-extract / concatenation lowerings the byte-layout regression +tests rely on; those tests are therefore tagged `no-new-smt`. + +# Why this matters + +Before this change, CBMC modelled `long double` on x86_64 as IEEE +`binary128` (`quadruple_precision`). The two formats agree on the +*value* of representable numbers but disagree on every individual *bit +position*: in `binary128` the sign sits at bit 127, the exponent at +bits 112-126, and there is no explicit integer bit. + +The macOS SDK (and similar code in BSD and embedded systems) reaches +into `long double` via a union to extract the sign or classify the +value. With the binary128 model the `__sexp` field of the SDK's +classification union (bytes 8-9) is always zero, so +`signbit(-1.0L)` evaluates to `0` and dependent control flow never +sees a negative value -- a soundness bug in any analysis that relies +on these classifications. + +# Testing + +The fix is exercised from several angles: + +- `regression/cbmc/long-double-x86-bytes` and + `regression/cbmc/long-double-i386-bytes` check the byte layout + against hardware-captured reference bytes. +- `regression/cbmc/long-double-signbit-x86` replicates the macOS SDK + signbit pattern. +- `regression/cbmc/long-double-roundtrip-x86` exercises + `double <-> long double` conversion, including the symbolic case. +- `regression/cbmc/long-double-pseudo-encodings` checks that + `__CPROVER_isnanld` / `__CPROVER_isinfld` / `__CPROVER_isnormalld` + reject the x86 pseudo-encodings (pseudo-NaN, pseudo-Infinity, + pseudo-denormal-with-non-zero-exponent), in agreement with modern + Intel/AMD hardware. +- `regression/cbmc-library/{expl,exp2l,logl,log2l,log10l,__builtin_powil}` + cover the long-double-only Schraudolph-style fast-math models in + `src/ansi-c/library/math.c`, which were updated in lockstep so that + the bit manipulations target the x86-extended layout when active. +- `unit/util/ieee_float.cpp` covers the `ieee_float_spect` factories + and `ieee_float_valuet::pack` / `unpack` against the same hardware + references. diff --git a/regression/cbmc/long-double-arithmetic-x86/main.i b/regression/cbmc/long-double-arithmetic-x86/main.i new file mode 100644 index 00000000000..49f8d01e96e --- /dev/null +++ b/regression/cbmc/long-double-arithmetic-x86/main.i @@ -0,0 +1,67 @@ +// Pure long-double FP arithmetic. No byte-level access to the +// underlying storage, so no bv<->FP boundary; this exercises the FP +// semantics on long double across all four CBMC back-ends. +// +// On x86_64 with x86 80-bit extended `long double`, the arithmetic +// semantics match IEEE binary79 for valid encodings (and our SAT +// backend encodes the layout precisely), so the assertions below hold +// under SAT, --z3, --cprover-smt2 and --incremental-smt2-solver. +// +// This complements the byte-layout tests, which exercise the FP<->bv +// boundary and are restricted to back-ends that fully support it. +// +// The test is hermetic (no system headers) and provided as `main.i` +// so CBMC parses it without invoking the host's preprocessor. It can +// therefore run on any host -- including FreeBSD/OpenBSD CI runners +// that don't ship with gcc and arm64 Linux runners that don't have +// x86_64 cross-compilation tooling -- as long as the test.desc fixes +// `--arch x86_64 --os linux` so CBMC models x86 80-bit extended +// regardless of the host's native long-double layout. + +int main(void) +{ + // 1. Exact addition / multiplication / division on representable + // values. + long double a = 1.0L; + long double b = 2.0L; + long double c = a + b; + __CPROVER_assert(c == 3.0L, "1.0L + 2.0L == 3.0L"); + + long double d = c * 2.0L; + __CPROVER_assert(d == 6.0L, "3.0L * 2.0L == 6.0L"); + + long double e = d / 3.0L; + __CPROVER_assert(e == 2.0L, "6.0L / 3.0L == 2.0L"); + + // 2. Sign preservation under negation. + long double f = 3.5L; + __CPROVER_assert(-f < 0.0L, "negation gives < 0"); + __CPROVER_assert(f > 0.0L, "positive is > 0"); + __CPROVER_assert(-(-f) == f, "double negation"); + + // 3. Long double has 64 bits of significand precision (well above + // the 53 bits of binary64), so 1.0L + 2^-60 must not equal + // 1.0L. This would round to 1.0 if the model were treating long + // double as `double`. + long double tiny = 1.0L; + for(int i = 0; i < 60; ++i) + tiny /= 2.0L; + long double almost_one = 1.0L + tiny; + __CPROVER_assert(almost_one != 1.0L, "long double has > 60 bits precision"); + + // 4. Reciprocal of a wide-range value: 1/x for very large x is + // positive and finite. At 2^-100 the result is comfortably + // *normal* in extended precision (whose subnormal threshold is + // around 2^-16382), so this exercises the wide exponent range + // rather than subnormal handling. TODO: add a separate test + // for genuine subnormal arithmetic once the pack/unpack of + // canonical x86 denormals is exercised end-to-end. + long double large = 1.0L; + for(int i = 0; i < 100; ++i) + large *= 2.0L; + long double recip = 1.0L / large; + __CPROVER_assert(recip > 0.0L, "1/2^100 > 0"); + __CPROVER_assert(recip < 1.0L, "1/2^100 < 1"); + + return 0; +} diff --git a/regression/cbmc/long-double-arithmetic-x86/test.desc b/regression/cbmc/long-double-arithmetic-x86/test.desc new file mode 100644 index 00000000000..28b2d8219cb --- /dev/null +++ b/regression/cbmc/long-double-arithmetic-x86/test.desc @@ -0,0 +1,22 @@ +CORE +main.i +--arch x86_64 --os linux --floatbv +^EXIT=0$ +^SIGNAL=0$ +^VERIFICATION SUCCESSFUL$ +-- +^warning: ignoring +-- +Pure FP arithmetic on x86 80-bit extended `long double` -- no byte-level +access, no FP<->bv boundary. Runs under all four back-ends (SAT, +--z3, --cprover-smt2, --incremental-smt2-solver), demonstrating that +long-double FP arithmetic is precise across the board. Byte-layout +testing is in the long-double-{x86-bytes,i386-bytes,signbit-x86,...} +tests, which are restricted to back-ends that support the FP<->bv +boundary. + +The source is provided as a pre-processed `main.i` so CBMC parses it +without invoking the host preprocessor; combined with `--arch x86_64 +--os linux`, this lets the test run on hosts whose native long-double +layout differs from x86 80-bit extended (e.g. ARM Linux, macOS +arm64) and on hosts that don't ship gcc (FreeBSD, OpenBSD). diff --git a/regression/cbmc/long-double-i386-bytes/main.i b/regression/cbmc/long-double-i386-bytes/main.i new file mode 100644 index 00000000000..38478434470 --- /dev/null +++ b/regression/cbmc/long-double-i386-bytes/main.i @@ -0,0 +1,63 @@ +// Same byte-layout invariants as long-double-x86-bytes but for i386 +// (32-bit x86), where `long double` is the 80-bit x87 extended-precision +// format stored in 12 bytes. The 80-bit value layout is identical to +// x86_64; only the storage container size differs. +// +// The test is hermetic (no system headers) and provided as `main.i` +// so CBMC parses it without invoking the host's preprocessor. It can +// therefore run on any host (in particular hosts without i386 +// cross-compilation tooling installed) provided the test.desc fixes +// `--i386-linux` so CBMC models the 12-byte layout regardless. + +typedef unsigned char u8; +typedef unsigned short u16; +typedef unsigned long long u64; + +typedef union +{ + long double ld; + u8 bytes[12]; + struct + { + u64 mantissa; + u16 sign_exp; + u16 pad; + } parts; +} ld_layout_t; + +int main(void) +{ + // 1. Storage size is 12 bytes on i386. + __CPROVER_assert(sizeof(long double) == 12, "12-byte storage"); + + // 2. The 80-bit value of 1.0L is mantissa = 0x8000000000000000, + // sign_exp = 0x3FFF (sign 0, biased exponent 0x3FFF = 16383). + ld_layout_t one; + one.ld = 1.0L; + __CPROVER_assert( + one.parts.mantissa == 0x8000000000000000ULL, "mantissa for 1.0L"); + __CPROVER_assert(one.parts.sign_exp == 0x3FFF, "sign_exp for 1.0L"); + + // 3. Negation flips the sign bit at the top of sign_exp. + ld_layout_t neg_one; + neg_one.ld = -1.0L; + __CPROVER_assert( + neg_one.parts.mantissa == 0x8000000000000000ULL, "mantissa for -1.0L"); + __CPROVER_assert(neg_one.parts.sign_exp == 0xBFFF, "sign_exp for -1.0L"); + + // 4. The byte view matches the 12-byte little-endian sequence. + __CPROVER_assert(one.bytes[0] == 0x00, "byte 0"); + __CPROVER_assert(one.bytes[1] == 0x00, "byte 1"); + __CPROVER_assert(one.bytes[2] == 0x00, "byte 2"); + __CPROVER_assert(one.bytes[3] == 0x00, "byte 3"); + __CPROVER_assert(one.bytes[4] == 0x00, "byte 4"); + __CPROVER_assert(one.bytes[5] == 0x00, "byte 5"); + __CPROVER_assert(one.bytes[6] == 0x00, "byte 6"); + __CPROVER_assert(one.bytes[7] == 0x80, "byte 7"); + __CPROVER_assert(one.bytes[8] == 0xFF, "byte 8"); + __CPROVER_assert(one.bytes[9] == 0x3F, "byte 9"); + __CPROVER_assert(one.bytes[10] == 0x00, "byte 10"); + __CPROVER_assert(one.bytes[11] == 0x00, "byte 11"); + + return 0; +} diff --git a/regression/cbmc/long-double-i386-bytes/test.desc b/regression/cbmc/long-double-i386-bytes/test.desc new file mode 100644 index 00000000000..de16bf989e6 --- /dev/null +++ b/regression/cbmc/long-double-i386-bytes/test.desc @@ -0,0 +1,17 @@ +CORE +main.i +--i386-linux --floatbv +^EXIT=0$ +^SIGNAL=0$ +^VERIFICATION SUCCESSFUL$ +-- +^warning: ignoring +-- +Verifies the byte layout of x86 80-bit extended `long double` on i386 +(12-byte storage container). The 80-bit value layout is identical to +x86_64; only the storage container size differs. + +The source is provided as a pre-processed `main.i` so CBMC parses it +without invoking the host preprocessor; combined with `--i386-linux`, +this lets the test run on any host (in particular hosts without i386 +multilib system headers and ARM hosts whose gcc lacks `-m32`). diff --git a/regression/cbmc/long-double-math-x86/main.i b/regression/cbmc/long-double-math-x86/main.i new file mode 100644 index 00000000000..e73e70b7cf5 --- /dev/null +++ b/regression/cbmc/long-double-math-x86/main.i @@ -0,0 +1,31 @@ +// Host-independent coverage for the x86 80-bit extended `long double` +// fast-math library models (expl / logl), which encode their +// Schraudolph-style approximation directly in the x86-extended bit +// layout (see src/ansi-c/library/math.c and the +// __cprover_x86ld_{to,from}_schraudolph helpers). +// +// The corresponding regression/cbmc-library/{expl,logl,...} tests run +// under the *host* architecture, so they only exercise the x86-extended +// path when CI happens to run on x86_64. This test pins `--arch x86_64` +// (via test.desc) and is provided pre-processed (`main.i`, no system +// headers), so the x86-extended models are verified regardless of the +// runner's native long-double layout. +// +// The asserted ranges mirror the cbmc-library expl / logl tests; they +// are wide because the models are deliberate approximations. + +long double expl(long double); +long double logl(long double); + +int main(void) +{ + long double e = expl(1.0l); + __CPROVER_assert(e > 2.713l && e < 2.886l, "expl(1.0L) approximates e"); + + // ln(e) ~ 1, with e written as a long double literal so the test needs + // no . + long double one = logl(2.718281828459045235360287471352662498l); + __CPROVER_assert(one > 0.942l && one < 1.002l, "logl(e) approximates 1"); + + return 0; +} diff --git a/regression/cbmc/long-double-math-x86/test.desc b/regression/cbmc/long-double-math-x86/test.desc new file mode 100644 index 00000000000..5aa61a25d7c --- /dev/null +++ b/regression/cbmc/long-double-math-x86/test.desc @@ -0,0 +1,34 @@ +CORE no-new-smt +main.i +--arch x86_64 --floatbv +^EXIT=0$ +^SIGNAL=0$ +^VERIFICATION SUCCESSFUL$ +-- +^warning: ignoring +-- +Host-independent coverage for the x86 80-bit extended `long double` +fast-math library models `expl` and `logl`, which encode their +Schraudolph-style approximation directly in the x86-extended bit layout +(see src/ansi-c/library/math.c and the __cprover_x86ld_{to,from}_schraudolph +helpers). The regression/cbmc-library/{expl,logl,...} tests run under +the *host* architecture, so they only exercise the x86-extended path +when CI happens to run on x86_64; pinning `--arch x86_64` here selects +the x86 80-bit extended `long double` on every Unix-like target +(Linux, macOS, *BSD; including ARM hosts via the forced x86_64 arch). + +We deliberately do NOT pin `--os linux`: instantiating the `expl` / +`logl` library models requires preprocessing the library, and `--os +linux` forces the GCC preprocessor, which is absent on the FreeBSD, +OpenBSD and Windows CI runners. Leaving the OS at the host default lets +CBMC use the host's native preprocessor (clang on the BSDs, cl on +Windows, gcc on Linux), so the test runs everywhere; the byte-layout +tests (long-double-x86-bytes etc.) cover the strict layout guarantee +without calling library functions. + +The new SMT2 backend (--incremental-smt2-solver) is excluded via +no-new-smt: it bit-blasts the Schraudolph models into a problem large +enough to exhaust memory during propositional reduction. + +The source is provided as a pre-processed `main.i` so CBMC parses the +test harness without invoking the host preprocessor. diff --git a/regression/cbmc/long-double-pseudo-encodings/main.i b/regression/cbmc/long-double-pseudo-encodings/main.i new file mode 100644 index 00000000000..b2fad358358 --- /dev/null +++ b/regression/cbmc/long-double-pseudo-encodings/main.i @@ -0,0 +1,84 @@ +// Modern x86 hardware treats x86 80-bit extended encodings with the +// explicit integer bit cleared as invalid wherever the value would +// otherwise be classified as normal/infinity/NaN: +// +// - integer bit = 0 with non-zero biased exponent: "pseudo-denormal" +// or "unsupported" -- raises FE_INVALID, not a normal value. +// - integer bit = 0 with all-ones exponent + zero fraction: +// "pseudo-infinity" -- raises FE_INVALID, not infinity. +// - integer bit = 0 with all-ones exponent + non-zero fraction: +// "pseudo-NaN" -- raises FE_INVALID, not a NaN. +// +// CBMC's classification helpers (isnan, isinf, isnormal) follow the +// hardware convention: they require the integer bit to be set before +// classifying a value as NaN/Inf/normal. This test constructs each +// pseudo-encoding via union and confirms the helpers reject them. +// +// The test is hermetic (no system headers) and provided as `main.i` +// so CBMC parses it without invoking the host's preprocessor. This +// avoids host-specific issues (e.g. mingw constructs CBMC +// cannot parse, BSDs lacking gcc, arm64 hosts lacking x86 cross +// tooling) provided the test.desc fixes `--arch x86_64 --os linux`. + +typedef unsigned long long u64; +typedef unsigned short u16; + +typedef union +{ + long double ld; + struct + { + u64 mantissa; + u16 sign_exp; + u16 pad[3]; + } parts; +} ld_layout_t; + +int main(void) +{ + // 1. A pseudo-NaN: exponent all-ones, fraction non-zero, but the + // explicit integer bit (top bit of mantissa) is 0. + ld_layout_t pseudo_nan; + pseudo_nan.parts.mantissa = 0x0000000000000001ULL; // int_bit=0, frac!=0 + pseudo_nan.parts.sign_exp = 0x7FFF; // exp=all-ones, sign=0 + pseudo_nan.parts.pad[0] = 0; + pseudo_nan.parts.pad[1] = 0; + pseudo_nan.parts.pad[2] = 0; + __CPROVER_assert(!__CPROVER_isnanld(pseudo_nan.ld), "pseudo-NaN is not NaN"); + + // 2. A pseudo-Infinity: exponent all-ones, fraction zero, integer + // bit 0. + ld_layout_t pseudo_inf; + pseudo_inf.parts.mantissa = 0x0000000000000000ULL; // int_bit=0, frac=0 + pseudo_inf.parts.sign_exp = 0x7FFF; // exp=all-ones, sign=0 + pseudo_inf.parts.pad[0] = 0; + pseudo_inf.parts.pad[1] = 0; + pseudo_inf.parts.pad[2] = 0; + __CPROVER_assert( + !__CPROVER_isinfld(pseudo_inf.ld), "pseudo-Inf is not infinity"); + + // 3. A pseudo-denormal: non-zero exponent, integer bit 0. Hardware + // treats this as a denormal-like; CBMC's `isnormal` rejects it. + ld_layout_t pseudo_denormal; + pseudo_denormal.parts.mantissa = 0x0000000000000001ULL; // int_bit=0 + pseudo_denormal.parts.sign_exp = 0x3FFF; // exp=16383 + pseudo_denormal.parts.pad[0] = 0; + pseudo_denormal.parts.pad[1] = 0; + pseudo_denormal.parts.pad[2] = 0; + __CPROVER_assert( + !__CPROVER_isnormalld(pseudo_denormal.ld), + "pseudo-denormal is not normal"); + + // 4. By contrast, a real normal value (1.0L) is classified normal, + // a real NaN is classified NaN, and a real infinity as inf. + long double one = 1.0L; + __CPROVER_assert(__CPROVER_isnormalld(one), "1.0L is normal"); + + long double nan_value = 0.0L / 0.0L; + __CPROVER_assert(__CPROVER_isnanld(nan_value), "0/0 is NaN"); + + long double inf_value = 1.0L / 0.0L; + __CPROVER_assert(__CPROVER_isinfld(inf_value), "1/0 is infinite"); + + return 0; +} diff --git a/regression/cbmc/long-double-pseudo-encodings/test.desc b/regression/cbmc/long-double-pseudo-encodings/test.desc new file mode 100644 index 00000000000..28b6f78b68a --- /dev/null +++ b/regression/cbmc/long-double-pseudo-encodings/test.desc @@ -0,0 +1,21 @@ +CORE no-new-smt +main.i +--arch x86_64 --os linux --floatbv --no-simplify +^EXIT=0$ +^SIGNAL=0$ +^VERIFICATION SUCCESSFUL$ +-- +^warning: ignoring +-- +Verifies that CBMC's __CPROVER_isnanld / __CPROVER_isinfld / +__CPROVER_isnormalld classification helpers reject the x86 80-bit +extended pseudo-encodings (pseudo-NaN, pseudo-Infinity, +pseudo-denormal-with-non-zero-exponent), in agreement with modern +Intel/AMD hardware which raises FE_INVALID for these patterns. The +new SMT2 backend (--incremental-smt2-solver) does not yet implement +the byte-extract lowering pattern this test relies on; the test is +therefore excluded from that profile via no-new-smt. + +The source is provided as a pre-processed `main.i` so CBMC parses it +without invoking the host preprocessor; combined with `--arch x86_64 +--os linux`, this lets the test run on any host. diff --git a/regression/cbmc/long-double-roundtrip-x86/main.i b/regression/cbmc/long-double-roundtrip-x86/main.i new file mode 100644 index 00000000000..83579db763d --- /dev/null +++ b/regression/cbmc/long-double-roundtrip-x86/main.i @@ -0,0 +1,58 @@ +// Round-trip casts between `double` and x86 80-bit extended `long +// double`. The 80-bit format has a strictly wider value range and more +// mantissa bits than `double`, so any finite double-precision value +// promoted to long double must round-trip exactly back to its original +// double bit pattern. +// +// The test is hermetic (no system headers) and provided as `main.i` so +// CBMC parses it without invoking the host's preprocessor. The +// test.desc fixes `--arch x86_64 --os linux` so CBMC models x86 80-bit +// extended regardless of the host's native long-double layout. + +typedef unsigned long long u64; + +int main(void) +{ + // 1. A handful of finite double constants round-trip exactly through + // `long double`. We compare the bit patterns directly via union + // aliasing. + union double_bits + { + double d; + u64 i; + }; + + double samples[] = { + 0.0, -0.0, 1.0, -1.0, 0.5, -0.5, 2.0, -2.0, 1024.0, 1.0 / 3.0}; + const int n = sizeof(samples) / sizeof(samples[0]); + + for(int i = 0; i < n; ++i) + { + long double promoted = (long double)samples[i]; + double demoted = (double)promoted; + union double_bits before; + union double_bits after; + before.d = samples[i]; + after.d = demoted; + __CPROVER_assert(before.i == after.i, "double -> long double -> double"); + } + + // 2. A nondet finite double round-trips bit-for-bit (excluding NaN + // and +/-infinity, which require comparison via bit patterns + // rather than ==). + double d; + __CPROVER_assume(d == d); // not NaN + __CPROVER_assume(d != 1.0 / 0.0 && d != -1.0 / 0.0); // not infinity + + long double promoted = (long double)d; + double demoted = (double)promoted; + __CPROVER_assert(d == demoted, "symbolic double round-trips"); + + // 3. Sign is preserved through a double -> long double promotion. + if(d < 0.0) + __CPROVER_assert(promoted < 0.0L, "sign preserved (negative)"); + if(d > 0.0) + __CPROVER_assert(promoted > 0.0L, "sign preserved (positive)"); + + return 0; +} diff --git a/regression/cbmc/long-double-roundtrip-x86/test.desc b/regression/cbmc/long-double-roundtrip-x86/test.desc new file mode 100644 index 00000000000..7923ee88b3a --- /dev/null +++ b/regression/cbmc/long-double-roundtrip-x86/test.desc @@ -0,0 +1,29 @@ +CORE no-new-smt +main.i +--arch x86_64 --os linux --floatbv +^EXIT=0$ +^SIGNAL=0$ +^VERIFICATION SUCCESSFUL$ +-- +^warning: ignoring +-- +Verifies that casting a finite `double` to `long double` and back yields +the original `double` bit-for-bit, and that sign is preserved across +the promotion. Exercises both the constant-folding path (compile-time) +and the symbolic path (nondet input). + +This test reads the bit pattern of a `double` via union aliasing +(`union { double d; uint64_t i; }`) AND casts that `double` to a +`long double` (x86 80-bit extended) and back. Under the legacy +cprover-smt2 backend in FPA mode (`use_FPA_theory == true`) this needs +the float<->bit-vector boundary support added earlier in this branch: +flattening an FPA float to its IEEE bits (smt2_convt::flatten2bv), the +single-argument `((_ to_fp e f) bv)` reinterpretation in CPROVER's smt2 +solver, and the mixed FPA/bit-vector float typecast lowering. + +The new SMT2 backend (`--incremental-smt2-solver`) has its own +pre-existing limitations and is excluded via no-new-smt. + +The source is provided as a pre-processed `main.i` so CBMC parses it +without invoking the host preprocessor; combined with `--arch x86_64 +--os linux`, this lets the test run on any host. diff --git a/regression/cbmc/long-double-signbit-x86/main.i b/regression/cbmc/long-double-signbit-x86/main.i new file mode 100644 index 00000000000..ffd5714aac0 --- /dev/null +++ b/regression/cbmc/long-double-signbit-x86/main.i @@ -0,0 +1,57 @@ +// Exercises the union-based signbit pattern that the macOS SDK uses for +// `long double`. See on Xcode 16: +// +// __header_always_inline int __inline_signbitl(long double __x) { +// union { long double __ld; struct{ unsigned long long __m; +// unsigned short __sexp; } __p; } __u; +// __u.__ld = __x; +// return (int)(__u.__p.__sexp >> 15); +// } +// +// Before this fix CBMC encoded `long double` on x86_64 as IEEE binary128 +// and so `__sexp` (bytes 8-9 of the storage) was always read as zero, +// making `signbit(-1.0L) == 0`. With x86 80-bit extended modelling the +// sign bit is at bit 79 (i.e. the top of the 16-bit `__sexp` field). +// +// The test is hermetic (no system headers) and provided as `main.i` so +// CBMC parses it without invoking the host's preprocessor. The +// test.desc fixes `--arch x86_64 --os linux` so CBMC models x86 80-bit +// extended regardless of the host's native long-double layout. + +typedef unsigned long long u64; +typedef unsigned short u16; + +static int my_signbitl(long double x) +{ + union + { + long double ld; + struct + { + u64 m; + u16 sexp; + } p; + } u; + u.ld = x; + return (int)(u.p.sexp >> 15); +} + +int main(void) +{ + // 1. Negative non-zero values have signbit == 1. + __CPROVER_assert(my_signbitl(-1.0L) == 1, "signbit(-1.0L) == 1"); + __CPROVER_assert(my_signbitl(-2.0L) == 1, "signbit(-2.0L) == 1"); + __CPROVER_assert(my_signbitl(-0.5L) == 1, "signbit(-0.5L) == 1"); + + // 2. Positive values and positive zero have signbit == 0. + __CPROVER_assert(my_signbitl(1.0L) == 0, "signbit(1.0L) == 0"); + __CPROVER_assert(my_signbitl(2.0L) == 0, "signbit(2.0L) == 0"); + __CPROVER_assert(my_signbitl(0.0L) == 0, "signbit(0.0L) == 0"); + + // 3. Negation flips the sign. + long double v = 3.5L; + __CPROVER_assert(my_signbitl(v) == 0, "signbit(v) == 0 for v > 0"); + __CPROVER_assert(my_signbitl(-v) == 1, "signbit(-v) == 1"); + + return 0; +} diff --git a/regression/cbmc/long-double-signbit-x86/test.desc b/regression/cbmc/long-double-signbit-x86/test.desc new file mode 100644 index 00000000000..af632ce3f3e --- /dev/null +++ b/regression/cbmc/long-double-signbit-x86/test.desc @@ -0,0 +1,19 @@ +CORE no-new-smt +main.i +--arch x86_64 --os linux --floatbv --no-simplify +^EXIT=0$ +^SIGNAL=0$ +^VERIFICATION SUCCESSFUL$ +-- +^warning: ignoring +-- +Replicates the signbit() pattern from the macOS SDK's for +`long double`. Without proper x86 80-bit extended layout this returned +0 for negative values, defeating the SDK's branch. The new SMT2 +backend (--incremental-smt2-solver) does not yet correctly lower the +sign-bit extraction for x86-extended values; the test is therefore +excluded from that profile via no-new-smt. + +The source is provided as a pre-processed `main.i` so CBMC parses it +without invoking the host preprocessor; combined with `--arch x86_64 +--os linux`, this lets the test run on any host. diff --git a/regression/cbmc/long-double-x86-bytes/main.i b/regression/cbmc/long-double-x86-bytes/main.i new file mode 100644 index 00000000000..b4affe3bdd5 --- /dev/null +++ b/regression/cbmc/long-double-x86-bytes/main.i @@ -0,0 +1,106 @@ +// On x86_64 (Linux, macOS, FreeBSD) `long double` is the 80-bit x87 +// extended-precision format stored in 16 bytes (12 on i386). The stored +// representation is, low to high (i.e. byte 0 is the LSB): +// +// bits 0 .. 62 : explicit fraction (62 bits, no implicit leading 1) +// bit 63 : explicit integer bit (J-bit) -- 1 for normal/inf/NaN, +// 0 for zero/denormal +// bits 64 .. 78 : biased 15-bit exponent +// bit 79 : sign bit +// bits 80 .. 127 : storage padding (zero on Linux/macOS) +// +// The bias is 16383, so a finite normal value v is encoded as +// (-1)^sign * 2^(exp - 16383) * (1.f) +// where the leading 1 is stored explicitly in bit 63. +// +// The ground-truth byte sequences below were captured on real macOS-15 +// x86_64 hardware (Xcode 16.4, Apple clang 17.0.0): +// +// 1.0L -> 00 00 00 00 00 00 00 80 ff 3f 00 00 00 00 00 00 +// -1.0L -> 00 00 00 00 00 00 00 80 ff bf 00 00 00 00 00 00 +// 2.0L -> 00 00 00 00 00 00 00 80 00 40 00 00 00 00 00 00 +// 0.5L -> 00 00 00 00 00 00 00 80 fe 3f 00 00 00 00 00 00 +// 0.0L -> 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 +// +// The test is hermetic (no system headers) and provided as `main.i` so +// CBMC parses it without invoking the host's preprocessor. The +// test.desc fixes `--arch x86_64 --os linux` so CBMC models x86 80-bit +// extended in 16-byte storage regardless of the host's native long +// double layout. + +typedef unsigned char u8; +typedef unsigned short u16; +typedef unsigned long long u64; + +typedef union +{ + long double ld; + u8 bytes[16]; + struct + { + u64 mantissa; + u16 sign_exp; + u16 pad[3]; + } parts; +} ld_layout_t; + +int main(void) +{ + // 1. Storage size is 16 bytes. + __CPROVER_assert(sizeof(long double) == 16, "16-byte storage"); + + // 2. The 80-bit value of 1.0L is mantissa = 0x8000000000000000, + // sign_exp = 0x3FFF (sign 0, biased exponent 0x3FFF = 16383). + ld_layout_t one; + one.ld = 1.0L; + __CPROVER_assert( + one.parts.mantissa == 0x8000000000000000ULL, "mantissa for 1.0L"); + __CPROVER_assert(one.parts.sign_exp == 0x3FFF, "sign_exp for 1.0L"); + + // 3. Negation flips the sign bit at the top of sign_exp. + ld_layout_t neg_one; + neg_one.ld = -1.0L; + __CPROVER_assert( + neg_one.parts.mantissa == 0x8000000000000000ULL, "mantissa for -1.0L"); + __CPROVER_assert(neg_one.parts.sign_exp == 0xBFFF, "sign_exp for -1.0L"); + + // 4. Scaling by 2 increments the biased exponent by 1. + ld_layout_t two; + two.ld = 2.0L; + __CPROVER_assert( + two.parts.mantissa == 0x8000000000000000ULL, "mantissa for 2.0L"); + __CPROVER_assert(two.parts.sign_exp == 0x4000, "sign_exp for 2.0L"); + + // 5. Scaling by 1/2 decrements the biased exponent by 1. + ld_layout_t half; + half.ld = 0.5L; + __CPROVER_assert( + half.parts.mantissa == 0x8000000000000000ULL, "mantissa for 0.5L"); + __CPROVER_assert(half.parts.sign_exp == 0x3FFE, "sign_exp for 0.5L"); + + // 6. Zero has all bits clear (in particular no explicit integer bit). + ld_layout_t zero; + zero.ld = 0.0L; + __CPROVER_assert(zero.parts.mantissa == 0ULL, "mantissa for 0.0L"); + __CPROVER_assert(zero.parts.sign_exp == 0, "sign_exp for 0.0L"); + + // 7. The byte view matches the captured native sequence for 1.0L. + __CPROVER_assert(one.bytes[0] == 0x00, "byte 0"); + __CPROVER_assert(one.bytes[1] == 0x00, "byte 1"); + __CPROVER_assert(one.bytes[2] == 0x00, "byte 2"); + __CPROVER_assert(one.bytes[3] == 0x00, "byte 3"); + __CPROVER_assert(one.bytes[4] == 0x00, "byte 4"); + __CPROVER_assert(one.bytes[5] == 0x00, "byte 5"); + __CPROVER_assert(one.bytes[6] == 0x00, "byte 6"); + __CPROVER_assert(one.bytes[7] == 0x80, "byte 7"); + __CPROVER_assert(one.bytes[8] == 0xFF, "byte 8"); + __CPROVER_assert(one.bytes[9] == 0x3F, "byte 9"); + __CPROVER_assert(one.bytes[10] == 0x00, "byte 10"); + __CPROVER_assert(one.bytes[11] == 0x00, "byte 11"); + __CPROVER_assert(one.bytes[12] == 0x00, "byte 12"); + __CPROVER_assert(one.bytes[13] == 0x00, "byte 13"); + __CPROVER_assert(one.bytes[14] == 0x00, "byte 14"); + __CPROVER_assert(one.bytes[15] == 0x00, "byte 15"); + + return 0; +} diff --git a/regression/cbmc/long-double-x86-bytes/test.desc b/regression/cbmc/long-double-x86-bytes/test.desc new file mode 100644 index 00000000000..7c63fc6e0e2 --- /dev/null +++ b/regression/cbmc/long-double-x86-bytes/test.desc @@ -0,0 +1,21 @@ +CORE no-new-smt +main.i +--arch x86_64 --os linux --floatbv --no-simplify +^EXIT=0$ +^SIGNAL=0$ +^VERIFICATION SUCCESSFUL$ +-- +^warning: ignoring +-- +Verifies the byte layout of x86 80-bit extended `long double` matches the +hardware reference on x86_64. Native byte sequences were captured on +macOS-15 x86_64 (Apple clang 17.0.0). The new SMT2 backend +(--incremental-smt2-solver) does not yet implement the byte-extract +lowering pattern this test relies on; the test is therefore excluded +from that profile via no-new-smt. + +The source is provided as a pre-processed `main.i` so CBMC parses it +without invoking the host preprocessor; combined with `--arch x86_64 +--os linux`, this lets the test run on any host (FreeBSD/OpenBSD CI +runners that don't ship gcc, ARM hosts that don't have x86_64 cross +tooling, etc.). diff --git a/regression/cbmc/union-double-bits-fpa/main.c b/regression/cbmc/union-double-bits-fpa/main.c new file mode 100644 index 00000000000..57cee8614e2 --- /dev/null +++ b/regression/cbmc/union-double-bits-fpa/main.c @@ -0,0 +1,33 @@ +#include + +// Exercises smt2_convt::flatten2bv on an FPA-encoded `double`: reading the +// bytes of a double through a union forces the SMT2 back-end to turn the +// floating-point value into its IEEE-754 interchange bit pattern. The +// SMT-LIB FloatingPoint theory has no float-to-bit-vector operation, so +// before this was supported the cprover-smt2 / FPA path hit an invariant +// here. `--no-simplify` keeps the constant from being folded away before +// it reaches the back-end, so the flattening code is actually exercised. +// +// The check is meaningful primarily under the cprover-smt2 profile (FPA +// theory); under the SAT back-end the value is bit-vector-encoded anyway, +// in which case the assertions simply hold. + +int main(void) +{ + union + { + double d; + uint64_t bits; + } u; + + u.d = 1.0; + __CPROVER_assert(u.bits == 0x3FF0000000000000ull, "IEEE-754 bits of 1.0"); + + u.d = -2.0; + __CPROVER_assert(u.bits == 0xC000000000000000ull, "IEEE-754 bits of -2.0"); + + u.d = 0.0; + __CPROVER_assert(u.bits == 0x0000000000000000ull, "IEEE-754 bits of 0.0"); + + return 0; +} diff --git a/regression/cbmc/union-double-bits-fpa/test.desc b/regression/cbmc/union-double-bits-fpa/test.desc new file mode 100644 index 00000000000..852372b49b8 --- /dev/null +++ b/regression/cbmc/union-double-bits-fpa/test.desc @@ -0,0 +1,23 @@ +CORE +main.c +--floatbv --no-simplify +^EXIT=0$ +^SIGNAL=0$ +^VERIFICATION SUCCESSFUL$ +-- +^warning: ignoring +-- +Reads the IEEE-754 bit pattern of an FPA-encoded `double` through a +union, which forces `smt2_convt::flatten2bv` to turn the floating-point +value into a bit-vector. The SMT-LIB FloatingPoint theory has no +float-to-bit-vector operation, so the cprover-smt2 / FPA path previously +hit an invariant here; the bit pattern of a constant is now emitted as a +literal bit-vector instead. + +`--no-simplify` keeps the constant union read from being folded away +before it reaches the SMT2 back-end, so the flattening code is actually +exercised. The test is meaningful primarily under the cprover-smt2 +profile (FPA theory); under the SAT back-end the value is bit-vector +encoded anyway and the assertions hold trivially. `double` is IEEE +binary64 on every supported target, so no architecture pinning is +needed. diff --git a/regression/smt2_solver/fp/bit-to-fp2.desc b/regression/smt2_solver/fp/bit-to-fp2.desc new file mode 100644 index 00000000000..8ff69739e02 --- /dev/null +++ b/regression/smt2_solver/fp/bit-to-fp2.desc @@ -0,0 +1,7 @@ +CORE +bit-to-fp2.smt2 + +^EXIT=0$ +^SIGNAL=0$ +^unsat$ +-- diff --git a/regression/smt2_solver/fp/bit-to-fp2.smt2 b/regression/smt2_solver/fp/bit-to-fp2.smt2 new file mode 100644 index 00000000000..c0ead6c2827 --- /dev/null +++ b/regression/smt2_solver/fp/bit-to-fp2.smt2 @@ -0,0 +1,15 @@ +; +; Reinterpret the 64-bit IEEE-754 interchange pattern of 1.0 +; (0x3FF0000000000000 = 4607182418800017408) as a double via the +; single-argument bit-pattern overload of to_fp, and check it equals 1.0. +; +(define-fun B0 () Bool (= + ((_ to_fp 11 53) (_ bv4607182418800017408 64)) + (fp #b0 #b01111111111 #b0000000000000000000000000000000000000000000000000000))) + +(assert (not B0)) + +; expected to be UNSAT, i.e., they are equal +(check-sat) + +(exit) diff --git a/src/ansi-c/cprover_library.cpp b/src/ansi-c/cprover_library.cpp index 80a18ad2b95..171494d0721 100644 --- a/src/ansi-c/cprover_library.cpp +++ b/src/ansi-c/cprover_library.cpp @@ -8,6 +8,7 @@ Author: Daniel Kroening, kroening@kroening.com #include "cprover_library.h" +#include #include #include #include @@ -23,6 +24,15 @@ static std::string get_cprover_library_text( { std::ostringstream library_text; + // Decide whether the target uses x86 80-bit extended `long double`. + // We delegate to `c_types::long_double_is_x86_extended()` (which is + // defined in terms of the type-system decision in + // `c_types.cpp::long_double_type()`) so this and the type system stay + // in sync by construction; the math.c models use the macro to dispatch + // between layout-specific encodings of the Schraudolph-style fast-math + // approximations for `expl`, `logl`, `powl` and friends. + const bool ld_is_x86_extended = long_double_is_x86_extended(); + library_text << "#line 1 \"\"\n" "#define " CPROVER_PREFIX "malloc_failure_mode " << std::to_string(config.ansi_c.malloc_failure_mode) @@ -36,7 +46,10 @@ static std::string get_cprover_library_text( config.ansi_c.malloc_failure_mode_assert_then_assume) << "\n" "#define " CPROVER_PREFIX "malloc_may_fail " - << std::to_string(config.ansi_c.malloc_may_fail) << "\n"; + << std::to_string(config.ansi_c.malloc_may_fail) + << "\n" + "#define " CPROVER_PREFIX "LDBL_IS_X86_EXTENDED " + << (ld_is_x86_extended ? "1" : "0") << "\n"; library_text << "#line 1 \"\"\n" diff --git a/src/ansi-c/library/math.c b/src/ansi-c/library/math.c index 32b6afa9fa1..a13d153a3af 100644 --- a/src/ansi-c/library/math.c +++ b/src/ansi-c/library/math.c @@ -2308,6 +2308,87 @@ float expf(float x) return u.f; } +/* FUNCTION: __cprover_x86ld_to_schraudolph */ + +#ifndef __CPROVER_STDINT_H_INCLUDED +# include +# define __CPROVER_STDINT_H_INCLUDED +#endif + +#if __CPROVER_LDBL_IS_X86_EXTENDED +// Decode an x86 80-bit extended `long double` into the 32-bit Schraudolph +// integer representation used by the long-double-only fast-math models. +// +// Layout of the 80-bit value (low to high): +// bits 0..62 : 63-bit explicit fraction +// bit 63 : explicit integer bit (J-bit, always 1 for normal numbers) +// bits 64..78 : 15-bit biased exponent +// bit 79 : sign bit +// bits 80+ : storage padding (zero on Linux/macOS/FreeBSD) +// +// The Schraudolph 32-bit integer for value x = 2^q * (1 + top16_frac/2^16) is +// (q + 16383) << 16 | top16_frac, with the sign bit (bit 31) cleared. The +// long-double fast-math callers (logl, log2l, log10l, powl, __builtin_powil) +// only invoke this on x > 0 (the negative-input cases short-circuit via NaN +// or special-value handling above), so we mask the sign bit out +// unconditionally rather than threading it through the result. +int32_t __cprover_x86ld_to_schraudolph(long double x) +{ + union + { + long double l; + unsigned char b[sizeof(long double)]; + } u = {.l = x}; + uint64_t low_bits = ((uint64_t)u.b[0]) | ((uint64_t)u.b[1] << 8) | + ((uint64_t)u.b[2] << 16) | ((uint64_t)u.b[3] << 24) | + ((uint64_t)u.b[4] << 32) | ((uint64_t)u.b[5] << 40) | + ((uint64_t)u.b[6] << 48) | ((uint64_t)u.b[7] << 56); + uint16_t high_bits = (uint16_t)((uint16_t)u.b[8] | ((uint16_t)u.b[9] << 8)); + return ( + int32_t)(((uint32_t)(high_bits & 0x7FFFu) << 16) | (uint32_t)((low_bits >> 47) & 0xFFFFu)); +} +#endif + +/* FUNCTION: __cprover_x86ld_from_schraudolph */ + +#ifndef __CPROVER_STDINT_H_INCLUDED +# include +# define __CPROVER_STDINT_H_INCLUDED +#endif + +#if __CPROVER_LDBL_IS_X86_EXTENDED +// Inverse of `__cprover_x86ld_to_schraudolph`: encode the 32-bit Schraudolph +// integer `s` as an x86 80-bit extended `long double` value. See the layout +// comment in `__cprover_x86ld_to_schraudolph` for the bit positions; the +// J-bit is set to 1 (Schraudolph results are always positive non-zero) and +// any storage padding above the 80-bit value is zero. +// +// The byte-by-byte writes are independent of host endianness; x86 itself is +// always little-endian, which is what the union read at the use site +// expects. +long double __cprover_x86ld_from_schraudolph(int32_t s) +{ + uint16_t high_bits = (uint16_t)(((uint32_t)s >> 16) & 0xFFFFu); + uint64_t low_bits = ((uint64_t)((uint32_t)s & 0xFFFFu) << 47) | (1ULL << 63); + union + { + long double l; + unsigned char b[sizeof(long double)]; + } u = {.b = {0}}; + u.b[0] = (unsigned char)(low_bits >> 0); + u.b[1] = (unsigned char)(low_bits >> 8); + u.b[2] = (unsigned char)(low_bits >> 16); + u.b[3] = (unsigned char)(low_bits >> 24); + u.b[4] = (unsigned char)(low_bits >> 32); + u.b[5] = (unsigned char)(low_bits >> 40); + u.b[6] = (unsigned char)(low_bits >> 48); + u.b[7] = (unsigned char)(low_bits >> 56); + u.b[8] = (unsigned char)(high_bits >> 0); + u.b[9] = (unsigned char)(high_bits >> 8); + return u.l; +} +#endif + /* FUNCTION: expl */ #ifndef __CPROVER_MATH_H_INCLUDED @@ -2359,6 +2440,7 @@ long double expl(long double x) return HUGE_VALL; # pragma CPROVER check pop } + // Schraudolph-style fast approximation of exp(x) ≈ 2^(x/ln(2)). // 16 is 32 - 1 sign bit - 15 exponent bits int32_t bias = (1 << 16) * ((1 << 14) - 1); int32_t exp_a_x = (int32_t)(x / M_LN2 * (long double)(1 << 16)) + bias; @@ -2370,11 +2452,17 @@ long double expl(long double x) __CPROVER_assume(result >= lower); __CPROVER_assume(result <= upper); -# ifndef _MSC_VER - _Static_assert +# if __CPROVER_LDBL_IS_X86_EXTENDED + // See `__cprover_x86ld_from_schraudolph` for the byte-level encoding; + // this is just a one-liner once the helper is in scope. + long double __cprover_x86ld_from_schraudolph(int32_t); + return __cprover_x86ld_from_schraudolph(result); # else +# ifndef _MSC_VER + _Static_assert +# else static_assert -# endif +# endif (sizeof(long double) % sizeof(int32_t) == 0, "bit width of long double is a multiple of bit width of int32_t"); union @@ -2382,12 +2470,13 @@ long double expl(long double x) long double l; int32_t i[sizeof(long double) / sizeof(int32_t)]; } u = {.i = {0}}; -# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ +# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ u.i[sizeof(long double) / sizeof(int32_t) - 1] = result; -# else +# else u.i[0] = result; -# endif +# endif return u.l; +# endif #endif } @@ -2576,11 +2665,27 @@ long double logl(long double x) #if LDBL_MAX_EXP == DBL_MAX_EXP return log(x); #else -# ifndef _MSC_VER - _Static_assert + // Schraudolph-style fast logarithm: reconstruct the 32-bit integer that + // would have been produced by the corresponding `expl` encoding, then + // invert the Schraudolph linear formula. + int32_t bias = (1 << 16) * ((1 << 14) - 1); + int32_t exp_c = __VERIFIER_nondet_int32_t(); + __CPROVER_assume(exp_c >= -5641 && exp_c <= 1); +# if __CPROVER_LDBL_IS_X86_EXTENDED + // x86 80-bit extended layout: see __cprover_x86ld_to_schraudolph for + // the byte-level decoding. This function is only called with x > 0 + // (the early-out checks above handle x <= 0), so the sign bit is + // always zero and the helper masks it out unconditionally. + int32_t __cprover_x86ld_to_schraudolph(long double); + int32_t exp_a_x = __cprover_x86ld_to_schraudolph(x); + return ((long double)exp_a_x - (long double)(bias + exp_c)) * M_LN2 / + (long double)(1 << 16); # else +# ifndef _MSC_VER + _Static_assert +# else static_assert -# endif +# endif (sizeof(long double) % sizeof(int32_t) == 0, "bit width of long double is a multiple of bit width of int32_t"); union @@ -2588,16 +2693,14 @@ long double logl(long double x) long double l; int32_t i[sizeof(long double) / sizeof(int32_t)]; } u = {x}; - int32_t bias = (1 << 16) * ((1 << 14) - 1); - int32_t exp_c = __VERIFIER_nondet_int32_t(); - __CPROVER_assume(exp_c >= -5641 && exp_c <= 1); -# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ +# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ return ((long double)u.i[sizeof(long double) / sizeof(int32_t) - 1] - (long double)(bias + exp_c)) * M_LN2 / (long double)(1 << 16); -# else +# else return ((long double)u.i[0] - (long double)(bias + exp_c)) * M_LN2 / (long double)(1 << 16); +# endif # endif #endif } @@ -2785,11 +2888,26 @@ long double log2l(long double x) #if LDBL_MAX_EXP == DBL_MAX_EXP return log2(x); #else -# ifndef _MSC_VER - _Static_assert + // Schraudolph-style fast log2: same Schraudolph integer reconstruction as + // logl, then divide by 2^16 (no M_LN2 factor). + int32_t bias = (1 << 16) * ((1 << 14) - 1); + int32_t exp_c = __VERIFIER_nondet_int32_t(); + __CPROVER_assume(exp_c >= -5641 && exp_c <= 1); +# if __CPROVER_LDBL_IS_X86_EXTENDED + // x86 80-bit extended layout: see __cprover_x86ld_to_schraudolph for + // the byte-level decoding. This function is only called with x > 0 + // (the early-out checks above handle x <= 0), so the sign bit is + // always zero and the helper masks it out unconditionally. + int32_t __cprover_x86ld_to_schraudolph(long double); + int32_t exp_a_x = __cprover_x86ld_to_schraudolph(x); + return ((long double)exp_a_x - (long double)(bias + exp_c)) / + (long double)(1 << 16); # else +# ifndef _MSC_VER + _Static_assert +# else static_assert -# endif +# endif (sizeof(long double) % sizeof(int32_t) == 0, "bit width of long double is a multiple of bit width of int32_t"); union @@ -2797,16 +2915,14 @@ long double log2l(long double x) long double l; int32_t i[sizeof(long double) / sizeof(int32_t)]; } u = {x}; - int32_t bias = (1 << 16) * ((1 << 14) - 1); - int32_t exp_c = __VERIFIER_nondet_int32_t(); - __CPROVER_assume(exp_c >= -5641 && exp_c <= 1); -# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ +# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ return ((long double)u.i[sizeof(long double) / sizeof(int32_t) - 1] - (long double)(bias + exp_c)) / (long double)(1 << 16); -# else +# else return ((long double)u.i[0] - (long double)(bias + exp_c)) / (long double)(1 << 16); +# endif # endif #endif } @@ -2999,11 +3115,26 @@ long double log10l(long double x) #if LDBL_MAX_EXP == DBL_MAX_EXP return log10(x); #else -# ifndef _MSC_VER - _Static_assert + // Schraudolph-style fast log10: same Schraudolph integer reconstruction + // as logl, then multiply by ln(2)/ln(10) / 2^16. + int32_t bias = (1 << 16) * ((1 << 14) - 1); + int32_t exp_c = __VERIFIER_nondet_int32_t(); + __CPROVER_assume(exp_c >= -5641 && exp_c <= 1); +# if __CPROVER_LDBL_IS_X86_EXTENDED + // x86 80-bit extended layout: see __cprover_x86ld_to_schraudolph for + // the byte-level decoding. This function is only called with x > 0 + // (the early-out checks above handle x <= 0), so the sign bit is + // always zero and the helper masks it out unconditionally. + int32_t __cprover_x86ld_to_schraudolph(long double); + int32_t exp_a_x = __cprover_x86ld_to_schraudolph(x); + return ((long double)exp_a_x - (long double)(bias + exp_c)) * + (M_LN2 / M_LN10) / (long double)(1 << 16); # else +# ifndef _MSC_VER + _Static_assert +# else static_assert -# endif +# endif (sizeof(long double) % sizeof(int32_t) == 0, "bit width of long double is a multiple of bit width of int32_t"); union @@ -3011,16 +3142,14 @@ long double log10l(long double x) long double l; int32_t i[sizeof(long double) / sizeof(int32_t)]; } u = {x}; - int32_t bias = (1 << 16) * ((1 << 14) - 1); - int32_t exp_c = __VERIFIER_nondet_int32_t(); - __CPROVER_assume(exp_c >= -5641 && exp_c <= 1); -# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ +# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ return ((long double)u.i[sizeof(long double) / sizeof(int32_t) - 1] - (long double)(bias + exp_c)) * (M_LN2 / M_LN10) / (long double)(1 << 16); -# else +# else return ((long double)u.i[0] - (long double)(bias + exp_c)) * (M_LN2 / M_LN10) / (long double)(1 << 16); +# endif # endif #endif } @@ -3431,11 +3560,29 @@ long double powl(long double x, long double y) #if LDBL_MAX_EXP == DBL_MAX_EXP return pow(x, y); #else -# ifndef _MSC_VER - _Static_assert + // Schraudolph-style fast pow: pow(x, y) ≈ exp(y * log(x)) is implemented + // by extracting the Schraudolph 32-bit integer encoding of x, multiplying + // by y, and writing the result back as a long double. See `expl` and + // `logl` for the per-layout encoding details. Only the read of the + // encoding and the write-back of the result differ between the x86 + // 80-bit extended and the IEEE binary128 layout; the overflow handling + // and arithmetic in between are shared. + int32_t bias = (1 << 16) * ((1 << 14) - 1); + int32_t exp_c = __VERIFIER_nondet_int32_t(); + __CPROVER_assume(exp_c >= -5641 && exp_c <= 1); + +# if __CPROVER_LDBL_IS_X86_EXTENDED + // This branch is only reached after the early-out checks above have + // eliminated x <= 0, so the helper can mask out the sign bit + // unconditionally. + int32_t __cprover_x86ld_to_schraudolph(long double); + int32_t exponent = __cprover_x86ld_to_schraudolph(x); # else +# ifndef _MSC_VER + _Static_assert +# else static_assert -# endif +# endif (sizeof(long double) % sizeof(int32_t) == 0, "bit width of long double is a multiple of bit width of int32_t"); union U @@ -3443,14 +3590,13 @@ long double powl(long double x, long double y) long double l; int32_t i[sizeof(long double) / sizeof(int32_t)]; } u = {x}; -# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ +# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ int32_t exponent = u.i[sizeof(long double) / sizeof(int32_t) - 1]; -# else +# else int32_t exponent = u.i[0]; +# endif # endif - int32_t bias = (1 << 16) * ((1 << 14) - 1); - int32_t exp_c = __VERIFIER_nondet_int32_t(); - __CPROVER_assume(exp_c >= -5641 && exp_c <= 1); + # pragma CPROVER check push # pragma CPROVER check disable "signed-overflow" long double mult_result = y * (long double)(exponent - (bias + exp_c)); @@ -3464,13 +3610,19 @@ long double powl(long double x, long double y) # pragma CPROVER check pop } int32_t result = (int32_t)mult_result + (bias + exp_c); + +# if __CPROVER_LDBL_IS_X86_EXTENDED + long double __cprover_x86ld_from_schraudolph(int32_t); + return __cprover_x86ld_from_schraudolph(result); +# else union U result_u = {.i = {0}}; -# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ +# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ result_u.i[sizeof(long double) / sizeof(int32_t) - 1] = result; -# else +# else result_u.i[0] = result; -# endif +# endif return result_u.l; +# endif #endif } @@ -3876,11 +4028,25 @@ long double __builtin_powil(long double x, int y) #if LDBL_MAX_EXP == DBL_MAX_EXP return __builtin_powi(x, y); #else -# ifndef _MSC_VER - _Static_assert + // Schraudolph-style fast pow with integer exponent. Same encoding as + // `powl`; see `expl` and `logl` for layout details. Only the read of + // the encoding and the write-back of the result differ between the x86 + // 80-bit extended and the IEEE binary128 layout. + int32_t bias = (1 << 16) * ((1 << 14) - 1); + int32_t exp_c = __VERIFIER_nondet_int32_t(); + __CPROVER_assume(exp_c >= -5641 && exp_c <= 1); + +# if __CPROVER_LDBL_IS_X86_EXTENDED + // This branch is only reached after the early-out checks above have + // eliminated x <= 0. + int32_t __cprover_x86ld_to_schraudolph(long double); + int32_t exponent = __cprover_x86ld_to_schraudolph(x); # else +# ifndef _MSC_VER + _Static_assert +# else static_assert -# endif +# endif (sizeof(long double) % sizeof(int32_t) == 0, "bit width of long double is a multiple of bit width of int32_t"); union U @@ -3888,14 +4054,13 @@ long double __builtin_powil(long double x, int y) long double l; int32_t i[sizeof(long double) / sizeof(int32_t)]; } u = {x}; -# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ +# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ int32_t exponent = u.i[sizeof(long double) / sizeof(int32_t) - 1]; -# else +# else int32_t exponent = u.i[0]; +# endif # endif - int32_t bias = (1 << 16) * ((1 << 14) - 1); - int32_t exp_c = __VERIFIER_nondet_int32_t(); - __CPROVER_assume(exp_c >= -5641 && exp_c <= 1); + # pragma CPROVER check push # pragma CPROVER check disable "signed-overflow" long double mult_result = @@ -3910,12 +4075,18 @@ long double __builtin_powil(long double x, int y) # pragma CPROVER check pop } int32_t result = (int32_t)mult_result + (bias + exp_c); + +# if __CPROVER_LDBL_IS_X86_EXTENDED + long double __cprover_x86ld_from_schraudolph(int32_t); + return __cprover_x86ld_from_schraudolph(result); +# else union U result_u = {.i = {0}}; -# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ +# if !defined(__BYTE_ORDER__) || __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ result_u.i[sizeof(long double) / sizeof(int32_t) - 1] = result; -# else +# else result_u.i[0] = result; -# endif +# endif return result_u.l; +# endif #endif } diff --git a/src/ansi-c/library_check.sh b/src/ansi-c/library_check.sh index 6883984ae26..4fb9e4e21ec 100755 --- a/src/ansi-c/library_check.sh +++ b/src/ansi-c/library_check.sh @@ -30,6 +30,12 @@ perl -p -i -e 's/^__CPROVER_jsa_synthesise\n//' __functions perl -p -i -e 's/^java::java.io.InputStream.read:\(\)I\n//' __functions perl -p -i -e 's/^__CPROVER_contracts_library\n//' __functions +# Internal long-double helpers shared between the x86-extended fast-math +# models (expl, logl, log2l, log10l, powl, __builtin_powil); they have no +# user-visible name and are exercised transitively by those tests. +perl -p -i -e 's/^__cprover_x86ld_from_schraudolph\n//' __functions +perl -p -i -e 's/^__cprover_x86ld_to_schraudolph\n//' __functions + # Some functions are implicitly covered by running on different operating # systems: perl -p -i -e 's/^_creat\n//' __functions # creat, macOS diff --git a/src/solvers/flattening/boolbv.cpp b/src/solvers/flattening/boolbv.cpp index bdeb176dd8f..d648c9d2f85 100644 --- a/src/solvers/flattening/boolbv.cpp +++ b/src/solvers/flattening/boolbv.cpp @@ -411,8 +411,15 @@ literalt boolbvt::convert_rest(const exprt &expr) const bvt &bv = convert_bv(op); CHECK_RETURN(!bv.empty()); const irep_idt type_id = op.type().id(); - if(type_id == ID_signedbv || type_id == ID_fixedbv || type_id == ID_floatbv) + if(type_id == ID_signedbv || type_id == ID_fixedbv) return bv_utils.sign_bit(bv); + if(type_id == ID_floatbv) + { + // For x86 80-bit extended in padded storage the sign bit lives at + // spec.value_width()-1, not at the top of the storage container. + const ieee_float_spect spec(to_floatbv_type(op.type())); + return bv[spec.value_width() - 1]; + } if(type_id == ID_unsignedbv) return const_literal(false); } diff --git a/src/solvers/floatbv/float_bv.cpp b/src/solvers/floatbv/float_bv.cpp index ee657069fe6..6233ca8e572 100644 --- a/src/solvers/floatbv/float_bv.cpp +++ b/src/solvers/floatbv/float_bv.cpp @@ -201,8 +201,12 @@ ieee_float_spect float_bvt::get_spec(const exprt &expr) exprt float_bvt::abs(const exprt &op, const ieee_float_spect &spec) { - // we mask away the sign bit, which is the most significant bit - const mp_integer v = power(2, spec.width() - 1) - 1; + // we clear the sign bit; for x86 extended in padded storage the sign + // bit lives at spec.value_width() - 1 rather than at the top of the + // storage container. + const mp_integer all_ones = power(2, spec.width()) - 1; + const mp_integer sign_bit_mask = power(2, spec.value_width() - 1); + const mp_integer v = all_ones - sign_bit_mask; const constant_exprt mask(integer2bvrep(v, spec.width()), op.type()); @@ -211,8 +215,10 @@ exprt float_bvt::abs(const exprt &op, const ieee_float_spect &spec) exprt float_bvt::negation(const exprt &op, const ieee_float_spect &spec) { - // we flip the sign bit with an xor - const mp_integer v = power(2, spec.width() - 1); + // we flip the sign bit with an xor; for x86 extended in padded storage + // the sign bit lives at spec.value_width() - 1 rather than at the top + // of the storage container. + const mp_integer v = power(2, spec.value_width() - 1); constant_exprt mask(integer2bvrep(v, spec.width()), op.type()); @@ -234,22 +240,57 @@ exprt float_bvt::is_equal( exprt isnan1=isnan(src1, spec); const or_exprt nan(isnan0, isnan1); - const equal_exprt bitwise_equal(src0, src1); + // Bitwise equality must ignore any storage padding above the value + // bits. For x86 80-bit extended `long double` in padded 96-/128-bit + // storage the padding bits above value_width() are not part of the + // numeric value, and for symbolic inputs they are unconstrained -- + // without this mask a == a could be spuriously false because the + // SAT solver could pick different padding for the two operands. + const exprt bitwise_equal = value_bits_equal(src0, src1, spec); return and_exprt( or_exprt(bitwise_equal, both_zero), not_exprt(nan)); } -exprt float_bvt::is_zero(const exprt &src) +exprt float_bvt::value_bits(const exprt &src, const ieee_float_spect &spec) { - // we mask away the sign bit, which is the most significant bit - const floatbv_typet &type=to_floatbv_type(src.type()); - std::size_t width=type.get_width(); + if(spec.value_width() == spec.width()) + return src; + // Drop the storage-padding bits above value_width(). This is only + // reachable for x86 80-bit extended `long double` (96-/128-bit + // storage container, 80 value bits); for every other format the + // early-out above applies. + return extractbits_exprt(src, std::size_t{0}, bv_typet(spec.value_width())); +} - const mp_integer v = power(2, width - 1) - 1; +exprt float_bvt::value_bits_equal( + const exprt &src0, + const exprt &src1, + const ieee_float_spect &spec) +{ + return equal_exprt(value_bits(src0, spec), value_bits(src1, spec)); +} - constant_exprt mask(integer2bvrep(v, width), src.type()); +exprt float_bvt::is_zero(const exprt &src) +{ + // We test that all value bits except the sign bit are zero. For x86 + // 80-bit extended `long double` in padded 96-/128-bit storage, the + // padding bits above value_width() are not part of the numeric + // value, and for symbolic inputs they are unconstrained. Building + // the mask over value_width()-1 bits (i.e. all of fraction + + // explicit-integer + exponent, with the sign bit cleared and any + // storage padding cleared as well) ensures padding bits are masked + // out of the comparison. + const floatbv_typet &type = to_floatbv_type(src.type()); + const ieee_float_spect spec(type); + const std::size_t width = type.get_width(); + + // Mask of all-ones in bit positions [0, value_width()-1) and zero + // everywhere else (sign bit and any padding above). + const mp_integer mask_value = power(2, spec.value_width() - 1) - 1; + + constant_exprt mask(integer2bvrep(mask_value, width), src.type()); ieee_float_valuet z(type); z.make_zero(); @@ -305,9 +346,19 @@ void float_bvt::rounding_mode_bitst::get(const exprt &rm) exprt float_bvt::sign_bit(const exprt &op) { - const bitvector_typet &bv_type=to_bitvector_type(op.type()); - std::size_t width=bv_type.get_width(); - return extractbit_exprt(op, width-1); + // For floatbv operands, the sign bit lives at value_width()-1 (i.e. at + // the top of the encoded value, not the top of the storage container + // -- which matters for x86 80-bit extended in padded 96-/128-bit + // storage). For non-floatbv bit-vector operands (e.g. when this + // helper is reused from from_signed_integer), the sign is the top bit + // of the bit-vector. + if(op.type().id() == ID_floatbv) + { + const ieee_float_spect spec(to_floatbv_type(op.type())); + return extractbit_exprt(op, spec.value_width() - 1); + } + const bitvector_typet &bv_type = to_bitvector_type(op.type()); + return extractbit_exprt(op, bv_type.get_width() - 1); } exprt float_bvt::from_signed_integer( @@ -500,9 +551,16 @@ exprt float_bvt::isnormal( const exprt &src, const ieee_float_spect &spec) { - return and_exprt( - not_exprt(exponent_all_zeros(src, spec)), - not_exprt(exponent_all_ones(src, spec))); + exprt result = and_exprt( + not_exprt(exponent_all_zeros(src, spec)), + not_exprt(exponent_all_ones(src, spec))); + // For x86 80-bit extended, a normal value must also have the explicit + // integer bit set; integer-bit-cleared encodings with a non-zero + // exponent are pseudo-denormals (treated as invalid by modern x86 + // hardware). + if(spec.x86_extended) + result = and_exprt(result, extractbit_exprt(src, spec.f)); + return result; } /// Subtracts the exponents @@ -1040,9 +1098,14 @@ exprt float_bvt::isinf( const exprt &src, const ieee_float_spect &spec) { - return and_exprt( - exponent_all_ones(src, spec), - fraction_all_zeros(src, spec)); + exprt result = + and_exprt(exponent_all_ones(src, spec), fraction_all_zeros(src, spec)); + // For x86 80-bit extended, infinity requires the explicit integer + // bit to be 1; integer-bit-cleared encodings with all-ones exponent + // are pseudo-infinities and are invalid on modern x86 hardware. + if(spec.x86_extended) + result = and_exprt(result, extractbit_exprt(src, spec.f)); + return result; } exprt float_bvt::isfinite( @@ -1057,7 +1120,11 @@ exprt float_bvt::get_exponent( const exprt &src, const ieee_float_spect &spec) { - return extractbits_exprt(src, spec.f, unsignedbv_typet(spec.e)); + // For IEEE the exponent sits at bits [f, f+e); for x86 80-bit + // extended the explicit integer bit occupies bit f, so the exponent + // is shifted up by one to bits [f+1, f+1+e). + const std::size_t exp_lsb = spec.x86_extended ? spec.f + 1 : spec.f; + return extractbits_exprt(src, exp_lsb, unsignedbv_typet(spec.e)); } /// Gets the fraction without hidden bit in a floating-point bit-vector src @@ -1065,6 +1132,8 @@ exprt float_bvt::get_fraction( const exprt &src, const ieee_float_spect &spec) { + // The lower f bits of the encoding are the explicit fraction (no + // implicit/explicit integer bit) for both IEEE and x86 layouts. return extractbits_exprt(src, 0, unsignedbv_typet(spec.f)); } @@ -1072,8 +1141,14 @@ exprt float_bvt::isnan( const exprt &src, const ieee_float_spect &spec) { - return and_exprt(exponent_all_ones(src, spec), - not_exprt(fraction_all_zeros(src, spec))); + exprt result = and_exprt( + exponent_all_ones(src, spec), not_exprt(fraction_all_zeros(src, spec))); + // For x86 80-bit extended, a NaN requires the explicit integer bit + // to be 1; integer-bit-cleared encodings with all-ones exponent are + // pseudo-NaNs and are invalid on modern x86 hardware. + if(spec.x86_extended) + result = and_exprt(result, extractbit_exprt(src, spec.f)); + return result; } /// normalize fraction/exponent pair returns 'zero' if fraction is zero @@ -1547,11 +1622,22 @@ float_bvt::unbiased_floatt float_bvt::unpack( result.fraction=get_fraction(src, spec); - // add hidden bit - exprt hidden_bit=isnormal(src, spec); - result.fraction= - concatenation_exprt(hidden_bit, result.fraction, - unsignedbv_typet(spec.f+1)); + if(spec.x86_extended) + { + // The explicit integer bit lives at position spec.f in the encoding. + // Use it directly rather than synthesising a hidden bit. + exprt integer_bit = extractbit_exprt(src, spec.f); + result.fraction = concatenation_exprt( + integer_bit, result.fraction, unsignedbv_typet(spec.f + 1)); + } + else + { + // add hidden bit + exprt hidden_bit=isnormal(src, spec); + result.fraction= + concatenation_exprt(hidden_bit, result.fraction, + unsignedbv_typet(spec.f+1)); + } result.exponent=get_exponent(src, spec); @@ -1593,6 +1679,39 @@ exprt float_bvt::pack( const if_exprt exponent( infinity_or_NaN, from_integer(-1, src.exponent.type()), src.exponent); + if(spec.x86_extended) + { + // Layout (high to low): sign | exponent | integer_bit | fraction. + // The explicit integer (J) bit is 1 for normal, infinity, and NaN, + // and 0 for true zero and canonical denormals (where the biased + // exponent is zero). Setting it to 1 only when (NaN || Inf || + // exp != 0) avoids producing a pseudo-denormal pattern (J=1, + // exp=0, frac!=0) for results that round into the denormal range; + // ieee_float_valuet::unpack would otherwise route those through + // the normal-number branch and misinterpret the encoded value. + const equal_exprt exp_zero(exponent, from_integer(0, exponent.type())); + exprt integer_bit = or_exprt(infinity_or_NaN, not_exprt(exp_zero)); + + // Build the full spec.width()-bit vector in one concatenation. + const std::size_t value_bits = spec.f + spec.e + 2; + if(spec.width() > value_bits) + { + const std::size_t pad_bits = spec.width() - value_bits; + const constant_exprt pad( + integer2bvrep(0, pad_bits), bv_typet(pad_bits)); + return typecast_exprt( + concatenation_exprt( + {pad, sign_bit, exponent, integer_bit, fraction}, + bv_typet(spec.width())), + spec.to_type()); + } + return typecast_exprt( + concatenation_exprt( + {sign_bit, exponent, integer_bit, fraction}, + bv_typet(spec.width())), + spec.to_type()); + } + // stitch all three together return typecast_exprt( concatenation_exprt( diff --git a/src/solvers/floatbv/float_bv.h b/src/solvers/floatbv/float_bv.h index be1cb853809..8df0155ecf4 100644 --- a/src/solvers/floatbv/float_bv.h +++ b/src/solvers/floatbv/float_bv.h @@ -32,6 +32,18 @@ class float_bvt static exprt isnormal(const exprt &, const ieee_float_spect &); static exprt isfinite(const exprt &, const ieee_float_spect &); + /// Returns the lowest `spec.value_width()` bits of `src`, so any + /// storage padding above the value (used by x86 80-bit extended + /// `long double` in 96-/128-bit storage) is dropped before the + /// caller does a bit-wise comparison. For formats with no padding + /// this is a no-op and `src` is returned unchanged. + static exprt value_bits(const exprt &src, const ieee_float_spect &spec); + + /// Bitwise equality of two floatbv values that ignores any storage + /// padding above `spec.value_width()`. See `value_bits`. + static exprt + value_bits_equal(const exprt &, const exprt &, const ieee_float_spect &); + // add/sub exprt add_sub( bool subtract, diff --git a/src/solvers/floatbv/float_utils.cpp b/src/solvers/floatbv/float_utils.cpp index ba72bb09b21..e32ef7dc277 100644 --- a/src/solvers/floatbv/float_utils.cpp +++ b/src/solvers/floatbv/float_utils.cpp @@ -37,7 +37,10 @@ bvt float_utilst::from_signed_integer(const bvt &src) unbiased_floatt result; // we need to convert negative integers - result.sign=sign_bit(src); + // Note: src is an integer bitvector, not a float, so the sign is at + // the top bit of the integer (src.back()), not at the float sign + // position. + result.sign = src.back(); result.fraction=bv_utils.absolute_value(src); @@ -255,9 +258,15 @@ bvt float_utilst::conversion( literalt float_utilst::is_normal(const bvt &src) { - return prop.land( - !exponent_all_zeros(src), - !exponent_all_ones(src)); + literalt result = + prop.land(!exponent_all_zeros(src), !exponent_all_ones(src)); + // For x86 80-bit extended, a normal value must also have the explicit + // integer bit set; integer-bit-cleared encodings with a non-zero + // exponent are pseudo-denormals (treated as invalid by modern x86 + // hardware). + if(spec.x86_extended) + result = prop.land(result, src[spec.f]); + return result; } /// Subtracts the exponents @@ -709,17 +718,21 @@ bvt float_utilst::rem(const bvt &src1, const bvt &src2) bvt float_utilst::negate(const bvt &src) { PRECONDITION(!src.empty()); - bvt result=src; - literalt &sign_bit=result[result.size()-1]; - sign_bit=!sign_bit; + bvt result = src; + // For x86 80-bit extended in padded storage the sign bit lives at the + // top of the 80-bit value (spec.value_width()-1), not at the top of + // the 96-/128-bit storage container. + const std::size_t sign_pos = spec.value_width() - 1; + result[sign_pos] = !result[sign_pos]; return result; } bvt float_utilst::abs(const bvt &src) { PRECONDITION(!src.empty()); - bvt result=src; - result[result.size()-1]=const_literal(false); + bvt result = src; + const std::size_t sign_pos = spec.value_width() - 1; + result[sign_pos] = const_literal(false); return result; } @@ -747,7 +760,17 @@ literalt float_utilst::relation( if(rel==relt::LT || rel==relt::LE) { - literalt bitwise_equal=bv_utils.equal(src1, src2); + // Bitwise equality and unsigned comparison must ignore any storage + // padding above the value bits. For x86 80-bit extended `long + // double` in 96-/128-bit storage the padding bits above + // value_width() are not part of the numeric value, and for + // symbolic inputs they are unconstrained -- without this mask the + // SAT solver could pick different padding for the two operands + // and make a < b or a == b false for two values that are + // mathematically equal. + const bvt v1 = value_bits(src1); + const bvt v2 = value_bits(src2); + literalt bitwise_equal = bv_utils.equal(v1, v2); // signs different? trivial! Unless Zero. @@ -757,7 +780,7 @@ literalt float_utilst::relation( // as long as the signs match: compare like unsigned numbers // this works due to the BIAS - literalt less_than1=bv_utils.unsigned_less_than(src1, src2); + literalt less_than1 = bv_utils.unsigned_less_than(v1, v2); // if both are negative (and not the same), need to turn around! literalt less_than2= @@ -792,7 +815,8 @@ literalt float_utilst::relation( } else if(rel==relt::EQ) { - literalt bitwise_equal=bv_utils.equal(src1, src2); + // See the LT/LE case above for why padding is masked off. + literalt bitwise_equal = bv_utils.equal(value_bits(src1), value_bits(src2)); return prop.land( prop.lor(bitwise_equal, both_zero), @@ -807,32 +831,43 @@ literalt float_utilst::relation( literalt float_utilst::is_zero(const bvt &src) { PRECONDITION(!src.empty()); - bvt all_but_sign; - all_but_sign=src; - all_but_sign.resize(all_but_sign.size()-1); + // Check all value bits except the sign bit are zero. For x86 80-bit + // extended `long double` in padded storage we must first drop the + // padding bits above value_width(); padding bits are not part of the + // numeric value, and for symbolic inputs they are unconstrained, so + // including them would let the solver pick non-zero padding and + // spuriously make is_zero false. + bvt all_but_sign = src; + if(all_but_sign.size() > spec.value_width()) + all_but_sign.resize(spec.value_width()); + // The sign bit is now the top bit; drop it. + all_but_sign.pop_back(); return bv_utils.is_zero(all_but_sign); } literalt float_utilst::is_plus_inf(const bvt &src) { - bvt and_bv; - and_bv.push_back(!sign_bit(src)); - and_bv.push_back(exponent_all_ones(src)); - and_bv.push_back(fraction_all_zeros(src)); - return prop.land(and_bv); + return prop.land(!sign_bit(src), is_infinity(src)); } literalt float_utilst::is_infinity(const bvt &src) { - return prop.land( - exponent_all_ones(src), - fraction_all_zeros(src)); + literalt result = prop.land(exponent_all_ones(src), fraction_all_zeros(src)); + // For x86 80-bit extended, infinity requires the explicit integer + // bit to be 1; integer-bit-cleared encodings with all-ones exponent + // are pseudo-infinities and are invalid on modern x86 hardware. + if(spec.x86_extended) + result = prop.land(result, src[spec.f]); + return result; } /// Gets the unbiased exponent in a floating-point bit-vector bvt float_utilst::get_exponent(const bvt &src) { - return bv_utils.extract(src, spec.f, spec.f+spec.e-1); + // For x86 extended the explicit integer bit sits at spec.f, so the + // exponent starts one position higher. + const std::size_t exp_lsb = spec.x86_extended ? spec.f + 1 : spec.f; + return bv_utils.extract(src, exp_lsb, exp_lsb + spec.e - 1); } /// Gets the fraction without hidden bit in a floating-point bit-vector src @@ -843,43 +878,28 @@ bvt float_utilst::get_fraction(const bvt &src) literalt float_utilst::is_minus_inf(const bvt &src) { - bvt and_bv; - and_bv.push_back(sign_bit(src)); - and_bv.push_back(exponent_all_ones(src)); - and_bv.push_back(fraction_all_zeros(src)); - return prop.land(and_bv); + return prop.land(sign_bit(src), is_infinity(src)); } literalt float_utilst::is_NaN(const bvt &src) { - return prop.land(exponent_all_ones(src), - !fraction_all_zeros(src)); + literalt result = prop.land(exponent_all_ones(src), !fraction_all_zeros(src)); + // For x86 80-bit extended, a NaN requires the explicit integer bit + // to be 1; integer-bit-cleared encodings with all-ones exponent are + // pseudo-NaNs and are invalid on modern x86 hardware. + if(spec.x86_extended) + result = prop.land(result, src[spec.f]); + return result; } literalt float_utilst::exponent_all_ones(const bvt &src) { - bvt exponent=src; - - // removes the fractional part - exponent.erase(exponent.begin(), exponent.begin()+spec.f); - - // removes the sign - exponent.resize(spec.e); - - return bv_utils.is_all_ones(exponent); + return bv_utils.is_all_ones(get_exponent(src)); } literalt float_utilst::exponent_all_zeros(const bvt &src) { - bvt exponent=src; - - // removes the fractional part - exponent.erase(exponent.begin(), exponent.begin()+spec.f); - - // removes the sign - exponent.resize(spec.e); - - return bv_utils.is_zero(exponent); + return bv_utils.is_zero(get_exponent(src)); } literalt float_utilst::fraction_all_zeros(const bvt &src) @@ -1383,7 +1403,16 @@ float_utilst::unbiased_floatt float_utilst::unpack(const bvt &src) result.sign=sign_bit(src); result.fraction=get_fraction(src); - result.fraction.push_back(is_normal(src)); // add hidden bit + + if(spec.x86_extended) + { + // The explicit integer bit is at position spec.f in the encoding. + result.fraction.push_back(src[spec.f]); + } + else + { + result.fraction.push_back(is_normal(src)); // add hidden bit + } result.exponent=get_exponent(src); CHECK_RETURN(result.exponent.size() == spec.e); @@ -1413,8 +1442,8 @@ bvt float_utilst::pack(const biased_floatt &src) // do sign // we make this 'false' for NaN - result[result.size()-1]= - prop.lselect(src.NaN, const_literal(false), src.sign); + const std::size_t sign_pos = spec.value_width() - 1; + result[sign_pos] = prop.lselect(src.NaN, const_literal(false), src.sign); literalt infinity_or_NaN= prop.lor(src.NaN, src.infinity); @@ -1425,11 +1454,28 @@ bvt float_utilst::pack(const biased_floatt &src) result[0]=prop.lor(result[0], src.NaN); + // For x86 extended, set the explicit integer bit at position spec.f. + // It is 1 for normal, infinity, and NaN, and 0 for true zero and + // canonical denormals (where the biased exponent is zero). Setting + // it to 1 only when (NaN || Inf || exp != 0) avoids producing a + // pseudo-denormal pattern (J=1, exp=0, frac!=0) for results that + // round into the denormal range; ieee_float_valuet::unpack would + // otherwise route those through the normal-number branch and + // misinterpret the encoded value. + if(spec.x86_extended) + { + literalt exp_zero = bv_utils.is_zero(src.exponent); + result[spec.f] = prop.lor(infinity_or_NaN, !exp_zero); + } + // do exponent + const std::size_t exp_lsb = spec.x86_extended ? spec.f + 1 : spec.f; for(std::size_t i=0; i_ form keeps distinct + // floatbv types apart even when they share a storage width -- notably + // IEEE binary128 (f=112) and x86 80-bit extended long double (f=63), + // both 128 bits wide -- so their float_bv.* helper bodies do not get + // conflated. ieee_float_spect spec(to_floatbv_type(type)); - return "f"+std::to_string(spec.width())+"_"+std::to_string(spec.f); + return "f" + std::to_string(spec.width()) + "_" + std::to_string(spec.f); } else if(type.id() == ID_bv) { @@ -1194,7 +1219,17 @@ std::string smt2_convt::floatbv_suffix(const exprt &expr) const void smt2_convt::convert_floatbv(const exprt &expr) { - PRECONDITION(!use_FPA_theory); + // `convert_floatbv` lowers FP operations to their bit-vector encoding + // (via the `float_bv.` helper functions defined elsewhere). The + // helper's return sort matches `convert_type(expr.type())`, so an + // FPA-encoded float result is already produced as an FPA value by the + // helper body -- no reinterpretation is needed here. Operands that + // are FPA-encoded floats, however, must be flattened to their IEEE + // bit pattern (the bvfromfloat reinterpret) so the bit-vector helper + // body sees bit-vector inputs. This only matters for boundary- + // crossing floatbv_typecast (e.g. (long double)d, (double)long_double); + // for homogeneous ops the operands are already in the helper's + // encoding. if(expr.id()==ID_symbol) { @@ -1220,7 +1255,12 @@ void smt2_convt::convert_floatbv(const exprt &expr) for(const auto &op : expr.operands()) { out << ' '; - convert_expr(op); + // FPA-encoded float operands must be presented to the bit-vector + // helper as their IEEE bit pattern (the bvfromfloat reinterpret). + if(op.type().id() == ID_floatbv && use_FPA_for_type(op.type())) + convert_expr(typecast_exprt{op, bv_typet{boolbv_width(op.type())}}); + else + convert_expr(op); } out << ')'; @@ -1450,7 +1490,7 @@ void smt2_convt::convert_expr(const exprt &expr) else if(type.id() == ID_floatbv) { // this has no rounding mode - if(use_FPA_theory) + if(use_FPA_for_type(type)) { out << "(fp.neg "; convert_expr(unary_minus_expr.op()); @@ -1480,7 +1520,7 @@ void smt2_convt::convert_expr(const exprt &expr) if(op_type.id()==ID_floatbv) { - if(use_FPA_theory) + if(use_FPA_for_type(op_type)) { out << "(fp.isNegative "; convert_expr(sign_expr.op()); @@ -1663,7 +1703,7 @@ void smt2_convt::convert_expr(const exprt &expr) "operands of float equal and not equal expressions shall have same type"); // The FPA theory properly treats NaN and negative zero. - if(use_FPA_theory) + if(use_FPA_for_type(rel_expr.lhs().type())) { if(rel_expr.id() == ID_ieee_float_notequal) out << "(not "; @@ -2164,7 +2204,7 @@ void smt2_convt::convert_expr(const exprt &expr) } else if(type.id()==ID_floatbv) { - if(use_FPA_theory) + if(use_FPA_for_type(type)) { out << "(fp.abs "; convert_expr(abs_expr.op()); @@ -2186,7 +2226,7 @@ void smt2_convt::convert_expr(const exprt &expr) out << "false"; else if(op_type.id()==ID_floatbv) { - if(use_FPA_theory) + if(use_FPA_for_type(op_type)) { out << "(fp.isNaN "; convert_expr(isnan_expr.op()); @@ -2208,7 +2248,7 @@ void smt2_convt::convert_expr(const exprt &expr) out << "true"; else if(op_type.id()==ID_floatbv) { - if(use_FPA_theory) + if(use_FPA_for_type(op_type)) { out << "(and "; @@ -2238,7 +2278,7 @@ void smt2_convt::convert_expr(const exprt &expr) out << "false"; else if(op_type.id()==ID_floatbv) { - if(use_FPA_theory) + if(use_FPA_for_type(op_type)) { out << "(fp.isInfinite "; convert_expr(isinf_expr.op()); @@ -2260,7 +2300,7 @@ void smt2_convt::convert_expr(const exprt &expr) out << "true"; else if(op_type.id()==ID_floatbv) { - if(use_FPA_theory) + if(use_FPA_for_type(op_type)) { out << "(fp.isNormal "; convert_expr(isnormal_expr.op()); @@ -2839,7 +2879,7 @@ void smt2_convt::convert_typecast(const typecast_exprt &expr) } else if(src_type.id()==ID_floatbv) { - if(use_FPA_theory) + if(use_FPA_for_type(src_type)) { out << "(not (fp.isZero "; convert_expr(src); @@ -2958,7 +2998,7 @@ void smt2_convt::convert_typecast(const typecast_exprt &expr) { // this is _NOT_ a semantic conversion, but bit-wise - if(use_FPA_theory) + if(use_FPA_for_type(src_type)) { defined_expressionst::const_iterator it = defined_expressions.find(expr); @@ -3304,7 +3344,7 @@ void smt2_convt::convert_typecast(const typecast_exprt &expr) UNEXPECTEDCASE("Typecast bv -> float with wrong width"); } - if(use_FPA_theory) + if(use_FPA_for_type(dest_floatbv_type)) { out << "((_ to_fp " << dest_floatbv_type.get_e() << " " << dest_floatbv_type.get_f() + 1 << ") "; @@ -3394,8 +3434,9 @@ void smt2_convt::convert_floatbv_typecast(const floatbv_typecast_exprt &expr) const floatbv_typet &dst=to_floatbv_type(dest_type); - if(use_FPA_theory) + if(use_FPA_for_type(dst) && use_FPA_for_type(src_type)) { + // Both sides FPA-encoded: native FPA value conversion. out << "((_ to_fp " << dst.get_e() << " " << dst.get_f() + 1 << ") "; convert_rounding_mode_FPA(expr.op1()); @@ -3404,7 +3445,14 @@ void smt2_convt::convert_floatbv_typecast(const floatbv_typecast_exprt &expr) out << ")"; } else + { + // At least one side is bit-vector-encoded (x86 80-bit extended + // long double). Lower the value-preserving conversion through + // the bit-vector float_bv helper; convert_floatbv flattens any + // FPA-encoded operand to its IEEE bit pattern and unflattens an + // FPA-encoded result via (_ to_fp ...). convert_floatbv(expr); + } } else if(src_type.id()==ID_unsignedbv) { @@ -3425,7 +3473,7 @@ void smt2_convt::convert_floatbv_typecast(const floatbv_typecast_exprt &expr) const floatbv_typet &dst=to_floatbv_type(dest_type); - if(use_FPA_theory) + if(use_FPA_for_type(dst)) { out << "((_ to_fp_unsigned " << dst.get_e() << " " << dst.get_f() + 1 << ") "; @@ -3443,7 +3491,7 @@ void smt2_convt::convert_floatbv_typecast(const floatbv_typecast_exprt &expr) const floatbv_typet &dst=to_floatbv_type(dest_type); - if(use_FPA_theory) + if(use_FPA_for_type(dst)) { out << "((_ to_fp " << dst.get_e() << " " << dst.get_f() + 1 << ") "; @@ -3471,7 +3519,7 @@ void smt2_convt::convert_floatbv_typecast(const floatbv_typecast_exprt &expr) } else if(dest_type.id()==ID_signedbv) { - if(use_FPA_theory) + if(use_FPA_for_type(src_type)) { std::size_t dest_width=to_signedbv_type(dest_type).get_width(); out << "((_ fp.to_sbv " << dest_width << ") "; @@ -3485,7 +3533,7 @@ void smt2_convt::convert_floatbv_typecast(const floatbv_typecast_exprt &expr) } else if(dest_type.id()==ID_unsignedbv) { - if(use_FPA_theory) + if(use_FPA_for_type(src_type)) { std::size_t dest_width=to_unsignedbv_type(dest_type).get_width(); out << "((_ fp.to_ubv " << dest_width << ") "; @@ -3509,7 +3557,7 @@ void smt2_convt::convert_floatbv_round_to_integral( { PRECONDITION(expr.type().id() == ID_floatbv); - if(use_FPA_theory) + if(use_FPA_for_type(expr.type())) { out << "(fp.roundToIntegral "; convert_rounding_mode_FPA(expr.rounding_mode()); @@ -3687,7 +3735,7 @@ void smt2_convt::convert_constant(const constant_exprt &expr) const floatbv_typet &floatbv_type= to_floatbv_type(expr_type); - if(use_FPA_theory) + if(use_FPA_for_type(expr_type)) { /* CBMC stores floating point literals in the most computationally useful form; biased exponents and @@ -3943,7 +3991,7 @@ void smt2_convt::convert_relation(const binary_relation_exprt &expr) } else if(op_type.id()==ID_floatbv) { - if(use_FPA_theory) + if(use_FPA_for_type(op_type)) { out << "("; if(expr.id()==ID_le) @@ -4138,6 +4186,12 @@ void smt2_convt::convert_plus(const plus_exprt &expr) /// \par parameters: The expression representing the rounding mode. void smt2_convt::convert_rounding_mode_FPA(const exprt &expr) { + // Every caller has already gated on `use_FPA_for_type(...)` for the + // operand or result type before reaching this function, so when the + // global FPA flag is set the relevant type is also FPA-encoded. We + // assert against the global flag here because the caller-side gate + // determines that we are emitting native FPA code; the rounding mode + // itself is shared across all FP types and not specific to any one. PRECONDITION(use_FPA_theory); /* CProver uses the x86 numbering of the rounding-mode @@ -4208,7 +4262,8 @@ void smt2_convt::convert_floatbv_plus(const ieee_float_op_exprt &expr) (type.id() == ID_complex && to_complex_type(type).subtype().id() == ID_floatbv)); - if(use_FPA_theory) + if(use_FPA_for_type( + type.id() == ID_complex ? to_complex_type(type).subtype() : type)) { if(type.id()==ID_floatbv) { @@ -4334,7 +4389,7 @@ void smt2_convt::convert_floatbv_minus(const ieee_float_op_exprt &expr) expr.type().id() == ID_floatbv, "type of ieee floating point expression shall be floatbv"); - if(use_FPA_theory) + if(use_FPA_for_type(expr.type())) { out << "(fp.sub "; convert_rounding_mode_FPA(expr.rounding_mode()); @@ -4408,7 +4463,7 @@ void smt2_convt::convert_floatbv_div(const ieee_float_op_exprt &expr) expr.type().id() == ID_floatbv, "type of ieee floating point expression shall be floatbv"); - if(use_FPA_theory) + if(use_FPA_for_type(expr.type())) { out << "(fp.div "; convert_rounding_mode_FPA(expr.rounding_mode()); @@ -4503,7 +4558,7 @@ void smt2_convt::convert_floatbv_mult(const ieee_float_op_exprt &expr) expr.type().id() == ID_floatbv, "type of ieee floating point expression shall be floatbv"); - if(use_FPA_theory) + if(use_FPA_for_type(expr.type())) { out << "(fp.mul "; convert_rounding_mode_FPA(expr.rounding_mode()); @@ -4523,7 +4578,7 @@ void smt2_convt::convert_floatbv_rem(const binary_exprt &expr) expr.type().id() == ID_floatbv, "type of ieee floating point expression shall be floatbv"); - if(use_FPA_theory) + if(use_FPA_for_type(expr.type())) { // Note that these do not have a rounding mode out << "(fp.rem "; @@ -4546,7 +4601,7 @@ void smt2_convt::convert_floatbv_fma(const floatbv_fma_exprt &expr) expr.type().id() == ID_floatbv, "type of ieee floating point expression shall be floatbv"); - if(use_FPA_theory) + if(use_FPA_for_type(expr.type())) { out << "(fp.fma "; convert_rounding_mode_FPA(expr.rounding_mode()); @@ -5004,11 +5059,42 @@ void smt2_convt::flatten2bv(const exprt &expr) } else if(type.id()==ID_floatbv) { - INVARIANT( - !use_FPA_theory, - "floatbv expressions should be flattened when using FPA theory"); - - convert_expr(expr); + if(use_FPA_for_type(type)) + { + // The SMT-LIB FloatingPoint theory has no float-to-bit-vector + // operation, so an FPA-encoded value cannot be flattened with a + // single operator. Two shapes cover the cases that arise from + // type-punning (union / byte access): + // - a constant: the IEEE interchange bit pattern is the value's + // bit-vector representation, emitted as a literal; + // - a bit-pattern reinterpret `(float)bits` (a typecast from a + // same-width bit-vector): flattening recovers exactly those + // bits, i.e. the typecast operand. + if(expr.is_constant()) + { + const ieee_float_spect spec(to_floatbv_type(type)); + const mp_integer value = bvrep2integer( + to_constant_expr(expr).get_value(), spec.width(), false); + out << "(_ bv" << value << " " << spec.width() << ")"; + } + else if( + expr.id() == ID_typecast && + to_typecast_expr(expr).op().type().id() == ID_bv && + boolbv_width(to_typecast_expr(expr).op().type()) == boolbv_width(type)) + { + convert_expr(to_typecast_expr(expr).op()); + } + else + { + INVARIANT( + false, + "flatten2bv of a non-constant FPA float that is not a " + "bit-pattern reinterpret is unsupported: " + + id2string(expr.id())); + } + } + else + convert_expr(expr); } else convert_expr(expr); @@ -5643,8 +5729,20 @@ void smt2_convt::find_symbols(const exprt &expr) } } // clang-format off - else if(!use_FPA_theory && - expr.operands().size() >= 1 && + else if(expr.operands().size() >= 1 && + (!use_FPA_for_type(to_multi_ary_expr(expr).op0().type()) || + // A floatbv_typecast whose *result* is a bit-vector-encoded + // float (x86 80-bit extended) is lowered through the bit-vector + // float_bv helper regardless of the operand's encoding (see + // convert_floatbv_typecast). This covers conversions whose + // operand is not itself bit-vector-encoded -- an FPA-encoded + // float (e.g. double -> long double) or an integer (e.g. the + // Schraudolph exponent arithmetic, int -> long double). The + // reverse directions have a bit-vector-encoded operand and are + // already covered by the clause above. + (expr.id() == ID_floatbv_typecast && + expr.type().id() == ID_floatbv && + !use_FPA_for_type(expr.type()))) && (expr.id() == ID_floatbv_plus || expr.id() == ID_floatbv_minus || expr.id() == ID_floatbv_mult || @@ -5666,8 +5764,31 @@ void smt2_convt::find_symbols(const exprt &expr) expr.id() == ID_typecast || expr.id() == ID_abs) && to_multi_ary_expr(expr).op0().type().id() == ID_floatbv))) - // clang-format on { + // clang-format on + // Helper-function parameters and return value are bit-vectors. For + // an FPA-encoded float operand or result we use the bit-vector sort + // of the corresponding IEEE encoding and flatten/unflatten at the + // call site (see convert_floatbv); the float_bv body operates on the + // bit pattern regardless. + const auto bv_lowered_sort = [this](const typet &t) + { + if(t.id() == ID_floatbv && use_FPA_for_type(t)) + out << "(_ BitVec " << boolbv_width(t) << ")"; + else + convert_type(t); + }; + + // Pre-register the bit-vector reinterpretation of any FPA-encoded + // float operand so its bvfromfloat declaration is emitted before the + // call site references it. This is done per occurrence (independent + // of whether the shared helper body is defined below). + for(const auto &op : expr.operands()) + { + if(op.type().id() == ID_floatbv && use_FPA_for_type(op.type())) + find_symbols(typecast_exprt{op, bv_typet{boolbv_width(op.type())}}); + } + irep_idt function = convert_identifier("float_bv." + expr.id_string() + floatbv_suffix(expr)); @@ -5683,11 +5804,17 @@ void smt2_convt::find_symbols(const exprt &expr) if(i!=0) out << " "; out << "(op" << i << ' '; - convert_type(expr.operands()[i].type()); + bv_lowered_sort(expr.operands()[i].type()); out << ')'; } out << ") "; + // The return sort is the natural sort of the result type: a + // bit-vector for x86-extended (BV-encoded) results, the FPA sort + // for an FPA-encoded result. float_bv's lowering produces a + // value of the result type, which convert_expr renders in exactly + // that encoding (emitting (_ to_fp ...) over the packed bits for + // an FPA result), so no parameter-style flattening applies here. convert_type(expr.type()); // return type out << ' '; @@ -5697,6 +5824,32 @@ void smt2_convt::find_symbols(const exprt &expr) smt2_symbolt("op"+std::to_string(i), tmp1.operands()[i].type()); exprt tmp2=float_bv(tmp1); + + // For FPA-encoded float operands the parameter is declared as a + // bit-vector (see bv_lowered_sort) and the call site passes the + // flattened IEEE bit pattern. float_bv's lowering, however, emits + // bit operations on the original floatbv-typed `op_i` symbol; left + // as-is, convert_expr would try to flatten2bv that FPA value (which + // is forbidden). Retype those operand symbols to the matching + // bit-vector so the body operates on the bit pattern directly. + for(std::size_t i = 0; i < expr.operands().size(); i++) + { + const typet &op_type = expr.operands()[i].type(); + if(op_type.id() == ID_floatbv && use_FPA_for_type(op_type)) + { + const smt2_symbolt fpa_sym{"op" + std::to_string(i), op_type}; + const smt2_symbolt bv_sym{ + "op" + std::to_string(i), bv_typet{boolbv_width(op_type)}}; + for(auto it = tmp2.depth_begin(), itend = tmp2.depth_end(); + it != itend; + ++it) + { + if(*it == fpa_sym) + it.mutate() = bv_sym; + } + } + } + tmp2=letify(tmp2); CHECK_RETURN(!tmp2.is_nil()); @@ -5706,8 +5859,9 @@ void smt2_convt::find_symbols(const exprt &expr) } } else if( - use_FPA_theory && expr.id() == ID_typecast && + expr.id() == ID_typecast && to_typecast_expr(expr).op().type().id() == ID_floatbv && + use_FPA_for_type(to_typecast_expr(expr).op().type()) && expr.type().id() == ID_bv) { // This is _NOT_ a semantic conversion, but bit-wise. @@ -6078,7 +6232,7 @@ void smt2_convt::convert_type(const typet &type) { const floatbv_typet &floatbv_type=to_floatbv_type(type); - if(use_FPA_theory) + if(use_FPA_for_type(type)) out << "(_ FloatingPoint " << floatbv_type.get_e() << " " << floatbv_type.get_f() + 1 << ")"; diff --git a/src/solvers/smt2/smt2_conv.h b/src/solvers/smt2/smt2_conv.h index ca94493a6be..3730b8ff7d4 100644 --- a/src/solvers/smt2/smt2_conv.h +++ b/src/solvers/smt2/smt2_conv.h @@ -71,6 +71,15 @@ class smt2_convt : public stack_decision_proceduret bool use_lambda_for_array = false; bool emit_set_logic = true; + /// Whether to use the SMT-LIB FloatingPoint theory for the given type. + /// Returns `use_FPA_theory` for ordinary IEEE 754 binary FP types. + /// Returns `false` for x86 80-bit extended `long double`: that format has + /// an explicit integer (J) bit and is stored in a 96- or 128-bit container + /// with high-bit storage padding, so it does not correspond to any + /// SMT-LIB FloatingPoint sort. Such values are always emitted as bit + /// vectors and operations on them are lowered through `convert_floatbv`. + bool use_FPA_for_type(const typet &type) const; + exprt handle(const exprt &expr) override; void set_to(const exprt &expr, bool value) override; exprt get(const exprt &expr) const override; diff --git a/src/solvers/smt2/smt2_parser.cpp b/src/solvers/smt2/smt2_parser.cpp index 94ceb8b3306..bfef0c6977a 100644 --- a/src/solvers/smt2/smt2_parser.cpp +++ b/src/solvers/smt2/smt2_parser.cpp @@ -739,6 +739,38 @@ exprt smt2_parsert::function_application() auto rounding_mode = expression(); + // The SMT-LIB FloatingPoint theory gives `to_fp` several + // overloads. All but one take a rounding mode followed by a + // source operand; the exception is the bit-pattern + // reinterpretation + // ((_ to_fp eb sb) (_ BitVec eb+sb)) + // which takes a single bit-vector operand and no rounding mode, + // interpreting those bits directly as the IEEE-754 interchange + // encoding. We have just parsed the single operand into + // `rounding_mode`; if the form closes here it was actually that + // operand. + if(smt2_tokenizer.peek() == smt2_tokenizert::CLOSE) + { + next_token(); // eat the ')' + + const exprt &bv_op = rounding_mode; + if( + bv_op.type().id() != ID_unsignedbv || + to_unsignedbv_type(bv_op.type()).get_width() != spec.width()) + { + throw error() << "to_fp bit-pattern reinterpretation expects a " + "bit-vector operand of width " + << spec.width(); + } + + // Reinterpret the bits as a float. Routing through a generic + // bit-vector type makes the conversion a bit-level reinterpret + // (see boolbv_typecastt) rather than a numeric integer-to-float + // conversion. + return typecast_exprt{ + typecast_exprt{bv_op, bv_typet{spec.width()}}, spec.to_type()}; + } + auto source_op = expression(); if(next_token() != smt2_tokenizert::CLOSE) diff --git a/src/util/bitvector_types.h b/src/util/bitvector_types.h index f3db88296f1..29a3440d027 100644 --- a/src/util/bitvector_types.h +++ b/src/util/bitvector_types.h @@ -409,6 +409,12 @@ class floatbv_typet : public bitvector_typet std::size_t get_e() const { + // x86 80-bit extended precision has a hardware-fixed 15-bit + // exponent; the value-width formula below would otherwise compute + // get_width()-f-1, which over-counts because of the explicit + // integer bit (and any storage padding above the 80-bit value). + if(get_bool(ID_x86_extended)) + return 15; // subtract one for sign bit return get_width() - get_f() - 1; } diff --git a/src/util/c_types.cpp b/src/util/c_types.cpp index 2b18e322f59..542e16d47bf 100644 --- a/src/util/c_types.cpp +++ b/src/util/c_types.cpp @@ -193,30 +193,59 @@ floatbv_typet double_type() floatbv_typet long_double_type() { floatbv_typet result; - if(config.ansi_c.long_double_width==128) - result=ieee_float_spect::quadruple_precision().to_type(); - else if(config.ansi_c.long_double_width==64) - result=ieee_float_spect::double_precision().to_type(); - else if(config.ansi_c.long_double_width==80) + // x86_64 typically stores `long double` as the x87 80-bit extended + // format in 16 bytes (Linux, macOS, FreeBSD). i386 stores it in 12 + // bytes. Accurate modelling matters because programs (notably the + // Apple SDK's inline helpers) read the bytes via union-based + // bit twiddling. + // + // The decision is driven primarily by `long_double_width`: + // - 96 bits is exclusively the i386 (32-bit x86) ABI: the x87 + // 80-bit extended value in a 12-byte container. This also arises + // when running a non-x86 build with `--32` (set_ILP32 changes the + // widths but leaves `arch` at the host value, e.g. "arm64"), so + // we accept width 96 regardless of `arch`. + // - 128 bits is ambiguous: x86_64 uses the x87 80-bit extended + // value in a 16-byte container, while non-x86 targets (ppc64le, + // AArch64 Linux configured that way) use true IEEE 754 binary128. + // Here `arch` is the tie-breaker. + // - 80 bits is the bare x87 extended value with no padding. + // - 64 bits is `long double == double` (AArch64 macOS, MSVC, ...). + const auto &c = config.ansi_c; + if(c.long_double_width == 96) { - // x86 extended precision has 80 bits in total, and - // deviating from IEEE, does not use a hidden bit. - // We use the closest we have got, but the below isn't accurate. - result=ieee_float_spect(63, 15).to_type(); + result = ieee_float_spect::x86_96().to_type(); } - else if(config.ansi_c.long_double_width==96) + else if(c.long_double_width == 128 && c.arch == "x86_64") { - result=ieee_float_spect(80, 15).to_type(); - // not quite right. The extra bits beyond 80 are usually padded. + // 16-byte storage container is x86_64-specific; i386 uses the 96-bit + // container handled above, so this branch deliberately excludes it. + result = ieee_float_spect::x86_128().to_type(); } + else if(c.long_double_width == 80) + { + result = ieee_float_spect::x86_80().to_type(); + } + else if(c.long_double_width == 128) + { + // Non-x86 128-bit long double: IEEE 754 binary128 (PowerPC, AArch64 + // when configured that way, etc.). + result = ieee_float_spect::quadruple_precision().to_type(); + } + else if(c.long_double_width == 64) + result = ieee_float_spect::double_precision().to_type(); else INVARIANT(false, "width of long double"); result.set(ID_C_c_type, ID_long_double); - return result; } +bool long_double_is_x86_extended() +{ + return long_double_type().get_bool(ID_x86_extended); +} + signedbv_typet pointer_diff_type() { // The pointer-diff type varies. This is signed int on some systems, diff --git a/src/util/c_types.h b/src/util/c_types.h index a16151d17d8..1ac65ab31e8 100644 --- a/src/util/c_types.h +++ b/src/util/c_types.h @@ -498,6 +498,12 @@ unsignedbv_typet char32_t_type(); floatbv_typet float_type(); floatbv_typet double_type(); floatbv_typet long_double_type(); +/// Returns true iff `long_double_type()` returns the x86 80-bit extended +/// format (i.e., its `ID_x86_extended` flag is set). This is the canonical +/// way for callers outside `c_types.cpp` to ask "does the current target +/// use x86-extended long double" without having to re-derive the answer +/// from `(arch, long_double_width)`. +bool long_double_is_x86_extended(); unsignedbv_typet size_type(); signedbv_typet signed_size_type(); signedbv_typet pointer_diff_type(); diff --git a/src/util/ieee_float.cpp b/src/util/ieee_float.cpp index 497a9558471..36df41af4c7 100644 --- a/src/util/ieee_float.cpp +++ b/src/util/ieee_float.cpp @@ -44,17 +44,31 @@ mp_integer ieee_float_spect::max_fraction() const void ieee_float_spect::from_type(const floatbv_typet &type) { - std::size_t width=type.get_width(); - f=type.get_f(); + std::size_t type_width = type.get_width(); + f = type.get_f(); DATA_INVARIANT(f != 0, "mantissa must be at least 1 bit"); DATA_INVARIANT( - f < width, + f < type_width, "mantissa bits must be less than " "originating type width"); - e=width-f-1; - x86_extended=type.get_bool(ID_x86_extended); + x86_extended = type.get_bool(ID_x86_extended); if(x86_extended) - e=e-1; // no hidden bit + { + // x86 80-bit extended precision: the value width is f+e+2 (sign, + // exponent, explicit integer bit, fraction). The exponent width + // for x86 80-bit is always 15. + e = 15; + std::size_t value_w = f + e + 2; + DATA_INVARIANT( + type_width >= value_w, + "type width must be at least the x86-extended value width"); + storage_width_bits = (type_width > value_w) ? type_width : 0; + } + else + { + e = type_width - f - 1; + storage_width_bits = 0; + } } ieee_float_valuet ieee_float_valuet::abs() const @@ -324,6 +338,78 @@ void ieee_float_valuet::unpack(const mp_integer &i) PRECONDITION(spec.f != 0); PRECONDITION(spec.e != 0); + if(spec.x86_extended) + { + // Mirror of the pack() layout above. We read only the low + // value_width() bits and ignore any storage padding above. + mp_integer tmp = i; + mp_integer pf = power(2, spec.f); + fraction = tmp % pf; + tmp /= pf; + + mp_integer integer_bit = tmp % 2; + tmp /= 2; + + mp_integer pe = power(2, spec.e); + exponent = tmp % pe; + tmp /= pe; + + sign_flag = (tmp % 2) != 0; + + if(exponent == spec.max_exponent()) + { + if(fraction == 0 && integer_bit == 1) + { + NaN_flag = false; + infinity_flag = true; + } + else + { + // Includes "real" NaNs (integer bit = 1, fraction != 0) and the + // historical x86 pseudo-encodings (integer bit = 0). All are + // treated as NaN. + make_NaN(); + } + } + else if(exponent == 0 && integer_bit == 0) + { + if(fraction == 0) + { + // signed zero + NaN_flag = false; + infinity_flag = false; + } + else + { + // denormal + NaN_flag = false; + infinity_flag = false; + exponent = -spec.bias() + 1; + } + } + else if(integer_bit == 1) + { + // normal: re-add the integer bit so `fraction` carries the full + // (f+1)-bit significand, matching the IEEE convention used + // elsewhere in this class. + NaN_flag = false; + infinity_flag = false; + fraction += pf; + exponent -= spec.bias(); + } + else + { + // Pseudo-denormal: integer bit = 0 with non-zero exponent. Modern + // x86 hardware traps these as invalid, but for analysis purposes we + // map them to the value implied by the bit pattern (denormal-like). + NaN_flag = false; + infinity_flag = false; + exponent -= spec.bias(); + } + + return; + } + { mp_integer tmp=i; @@ -376,7 +462,52 @@ bool ieee_float_valuet::is_normal() const mp_integer ieee_float_valuet::pack() const { - mp_integer result=0; + mp_integer result = 0; + + if(spec.x86_extended) + { + // x86 80-bit extended precision layout (low to high): + // bits 0 .. f-1 : fraction (without integer bit) + // bit f : explicit integer bit (1 for normal, 0 for + // denormal/zero) + // bits f+1 .. f+e : biased exponent + // bit f+1+e : sign + // Higher bits up to spec.width()-1 are storage padding (zero). + if(NaN_flag) + { + // Quiet NaN: integer bit set, exponent all-ones, fraction MSB set. + result += power(2, spec.f); // integer bit + result += spec.max_exponent() * power(2, spec.f + 1); // exponent + result += power(2, spec.f - 1); // quiet bit + } + else if(infinity_flag) + { + result += power(2, spec.f); // integer bit + result += spec.max_exponent() * power(2, spec.f + 1); // exponent + } + else if(fraction == 0 && exponent == 0) + { + // signed zero -- nothing to do + } + else if(is_normal()) + { + // fraction holds f+1 bits including the leading 1; the leading 1 is + // the explicit integer bit. + result += fraction; + result += (exponent + spec.bias()) * power(2, spec.f + 1); + } + else // denormal + { + // Denormal: integer bit is 0, exponent is 0; fraction holds the + // sub-normal mantissa with no leading 1. + result += fraction; + } + + if(sign_flag) + result += power(2, spec.f + 1 + spec.e); + + return result; + } // sign bit if(sign_flag) diff --git a/src/util/ieee_float.h b/src/util/ieee_float.h index 4965742eb4e..c5b1d37a4aa 100644 --- a/src/util/ieee_float.h +++ b/src/util/ieee_float.h @@ -29,6 +29,13 @@ class ieee_float_spect // integer bit. bool x86_extended; + // Storage size in bits, when different from the natural value width. + // x86 80-bit extended precision is stored in 12-byte (96-bit) or + // 16-byte (128-bit) memory locations on i386 and x86_64 respectively; + // the upper bits are padding. A value of 0 means "no padding, storage + // size equals value_width()". + std::size_t storage_width_bits; + mp_integer bias() const; explicit ieee_float_spect(const floatbv_typet &type) @@ -38,20 +45,32 @@ class ieee_float_spect void from_type(const floatbv_typet &type); - ieee_float_spect():f(0), e(0), x86_extended(false) + ieee_float_spect() : f(0), e(0), x86_extended(false), storage_width_bits(0) { } - ieee_float_spect(std::size_t _f, std::size_t _e): - f(_f), e(_e), x86_extended(false) + ieee_float_spect(std::size_t _f, std::size_t _e) + : f(_f), e(_e), x86_extended(false), storage_width_bits(0) { } - std::size_t width() const + /// Number of bits actually used to encode a value: sign + exponent + + /// fraction (+ explicit integer bit for x86 extended). + std::size_t value_width() const { // Add one for the sign bit. // Add one if x86 explicit integer bit is used. - return f+e+1+(x86_extended?1:0); + return f + e + 1 + (x86_extended ? 1 : 0); + } + + /// Storage width in bits. For most formats this equals value_width(), + /// but x86 80-bit extended precision is stored in larger containers + /// (96 bits on i386, 128 bits on x86_64) with the high bits as padding. + std::size_t width() const + { + if(storage_width_bits != 0) + return storage_width_bits; + return value_width(); } mp_integer max_exponent() const; @@ -87,23 +106,38 @@ class ieee_float_spect static ieee_float_spect x86_80() { - // Intel, not IEEE + // Intel x87 80-bit extended precision, used as the value with no + // additional padding (e.g. for solver-internal arithmetic). ieee_float_spect result(63, 15); - result.x86_extended=true; + result.x86_extended = true; return result; } static ieee_float_spect x86_96() { - // Intel, not IEEE + // x86 80-bit extended precision stored in a 96-bit (12-byte) + // container, as used for `long double` on i386. + ieee_float_spect result(63, 15); + result.x86_extended = true; + result.storage_width_bits = 96; + return result; + } + + static ieee_float_spect x86_128() + { + // x86 80-bit extended precision stored in a 128-bit (16-byte) + // container, as used for `long double` on x86_64 (Linux, macOS, + // FreeBSD, etc.). ieee_float_spect result(63, 15); - result.x86_extended=true; + result.x86_extended = true; + result.storage_width_bits = 128; return result; } bool operator==(const ieee_float_spect &other) const { - return f==other.f && e==other.e && x86_extended==other.x86_extended; + return f == other.f && e == other.e && x86_extended == other.x86_extended && + storage_width_bits == other.storage_width_bits; } bool operator!=(const ieee_float_spect &other) const diff --git a/unit/util/ieee_float.cpp b/unit/util/ieee_float.cpp index 5d5ddafbd14..9c7139f87ec 100644 --- a/unit/util/ieee_float.cpp +++ b/unit/util/ieee_float.cpp @@ -129,3 +129,179 @@ TEST_CASE("round_to_integral", "[unit][util][ieee_float]") REQUIRE(round_to_integral(from_double(0x1.0p+52), away) == 0x1.0p+52); REQUIRE(round_to_integral(from_double(dmax), away) == dmax); } + +TEST_CASE( + "ieee_float_spect: x86 extended specs", + "[core][util][ieee_float][x86_extended]") +{ + // The 80-bit value contains 1 sign + 15 exponent + 1 explicit integer + // bit + 63 fraction bits. Storage may add padding to reach 96 or 128 + // bits. + + const auto x86_80 = ieee_float_spect::x86_80(); + REQUIRE(x86_80.f == 63); + REQUIRE(x86_80.e == 15); + REQUIRE(x86_80.x86_extended); + REQUIRE(x86_80.value_width() == 80); + REQUIRE(x86_80.width() == 80); + + const auto x86_96 = ieee_float_spect::x86_96(); + REQUIRE(x86_96.f == 63); + REQUIRE(x86_96.e == 15); + REQUIRE(x86_96.x86_extended); + REQUIRE(x86_96.value_width() == 80); + REQUIRE(x86_96.width() == 96); + + const auto x86_128 = ieee_float_spect::x86_128(); + REQUIRE(x86_128.f == 63); + REQUIRE(x86_128.e == 15); + REQUIRE(x86_128.x86_extended); + REQUIRE(x86_128.value_width() == 80); + REQUIRE(x86_128.width() == 128); + + // The bias is 2^(e-1)-1 = 16383. + REQUIRE(x86_128.bias() == 16383); +} + +TEST_CASE( + "ieee_float_valuet: pack/unpack on x86 80-bit extended (all storage widths)", + "[core][util][ieee_float][x86_extended]") +{ + // Reference encodings captured on macOS-15 x86_64 hardware (Apple + // clang 17.0.0). The 80-bit value pattern is the same across all + // three storage widths; only the leading zero padding differs. + // + // 1.0L -> 00 00 00 00 00 00 00 80 ff 3f [+ 0..6 zero pad bytes] + // -1.0L -> 00 00 00 00 00 00 00 80 ff bf [+ 0..6 zero pad bytes] + // 2.0L -> 00 00 00 00 00 00 00 80 00 40 [+ 0..6 zero pad bytes] + // 0.5L -> 00 00 00 00 00 00 00 80 fe 3f [+ 0..6 zero pad bytes] + // 0.0L -> 00 00 00 00 00 00 00 00 00 00 [+ 0..6 zero pad bytes] + // + // Storage padding is zero on Linux/macOS, so the mp_integer + // representation is just the bottom 80 bits of the value, regardless + // of whether the storage container is 80, 96 or 128 bits. + + auto from_hex = [](const char *hex) -> mp_integer + { return string2integer(hex, 16); }; + + // The bottom 80 bits of each reference value, written as a hex + // big-endian word. These are reused for x86_80, x86_96 and x86_128 + // because the storage padding above the value bits is zero. + const mp_integer one_80 = from_hex("3FFF8000000000000000"); + const mp_integer neg_one_80 = from_hex("BFFF8000000000000000"); + const mp_integer two_80 = from_hex("40008000000000000000"); + const mp_integer half_80 = from_hex("3FFE8000000000000000"); + const mp_integer zero_80 = 0; + + // Smallest positive denormal: J=0, exp=0, frac=1. Hardware encodes + // this as 80 bits with the low fraction bit set and the J-bit clear, + // i.e. mp_integer value 1. Verifies the canonical denormal pattern + // from `float_utilst::pack` / `float_bvt::pack` (J=0 for exp=0). + const mp_integer smallest_denormal_80 = 1; + + // Positive infinity: J=1, exp=all-1s, frac=0. + const mp_integer pos_inf_80 = from_hex("7FFF8000000000000000"); + + // A representable quiet NaN: J=1, exp=all-1s, frac with at least one + // bit set (here, the fraction MSB so the NaN is "quiet" on hardware + // that interprets that bit as the quiet flag). + const mp_integer quiet_nan_80 = from_hex("7FFFC000000000000000"); + + for(const auto &spec : + {ieee_float_spect::x86_80(), + ieee_float_spect::x86_96(), + ieee_float_spect::x86_128()}) + { + auto make = [&spec](mp_integer i) -> ieee_floatt + { + ieee_floatt v{spec, ieee_floatt::ROUND_TO_EVEN}; + v.from_integer(i); + return v; + }; + + // Hardware-reference round-trip for representable integer values. + { + const auto v = make(1); + REQUIRE(v.pack() == one_80); + + ieee_float_valuet u{spec}; + u.unpack(one_80); + REQUIRE(u == v); + } + + { + const auto v = make(-1); + REQUIRE(v.pack() == neg_one_80); + + ieee_float_valuet u{spec}; + u.unpack(neg_one_80); + REQUIRE(u == v); + } + + { + const auto v = make(2); + REQUIRE(v.pack() == two_80); + + ieee_float_valuet u{spec}; + u.unpack(two_80); + REQUIRE(u == v); + } + + { + const auto v = make(0); + REQUIRE(v.pack() == zero_80); + + ieee_float_valuet u{spec}; + u.unpack(zero_80); + REQUIRE(u == v); + } + + // Round-trip a range of integers through pack/unpack. + for(int i : {-1024, -3, -2, -1, 0, 1, 2, 3, 7, 1024}) + { + const auto v = make(i); + ieee_float_valuet u{spec}; + u.unpack(v.pack()); + REQUIRE(u == v); + } + + // Non-integer encodings: NaN, infinity, smallest denormal. These + // exercise the special-case branches in `unpack`/`pack` that are + // not reached by the integer round-trip above. + + // Positive infinity: unpack -> classify as infinity -> pack + // re-emits the same canonical pattern. + { + ieee_float_valuet u{spec}; + u.unpack(pos_inf_80); + REQUIRE(u.is_infinity()); + REQUIRE(!u.get_sign()); + REQUIRE(u.pack() == pos_inf_80); + } + + // A quiet NaN. Pack re-emits a canonical NaN pattern, which is + // not necessarily bitwise equal to the input (the exact NaN + // payload is implementation-defined), so we only check that the + // round-trip preserves the NaN classification. + { + ieee_float_valuet u{spec}; + u.unpack(quiet_nan_80); + REQUIRE(u.is_NaN()); + + ieee_float_valuet u2{spec}; + u2.unpack(u.pack()); + REQUIRE(u2.is_NaN()); + } + + // Smallest positive denormal: J=0, exp=0, frac=1. This catches + // the pseudo-denormal bug in the pack path: an x86 denormal must + // be encoded with J=0, and `unpack` of J=0+exp=0 must take the + // denormal branch (exponent -bias+1) rather than the + // pseudo-denormal/normal branch (-bias). + { + ieee_float_valuet u{spec}; + u.unpack(smallest_denormal_80); + REQUIRE(u.pack() == smallest_denormal_80); + } + } +}