From f158d8663ee2e56e6a8ced3524c8cdab90b64833 Mon Sep 17 00:00:00 2001 From: Marc Bennewitz Date: Fri, 4 Sep 2026 16:36:29 +0200 Subject: [PATCH 1/2] Add test for expm1() precision on a double-precision-clamped FPU On i386 the x87 FPU is the platform default and zend_init_fpu() clamps it to double precision (53-bit mantissa) for the whole request. libm's expm1() computes on the x87 unit and relies on the extended range to round correctly to double, so while that clamp is in place expm1() is off by one or more ULP. This matters because it makes expm1() silently architecture-dependent: the same script prints different floats on i386 than on x86-64 or arm64, comparisons against precomputed constants fail, and the error propagates into anything built on expm1(). Every expected value is the correctly rounded double, verified against a high-precision reference rather than against another platform, and written as the shortest decimal that round-trips back to it. All arguments are exact binary fractions. The first group is inputs that are already correct everywhere, so a failure points at the affected inputs rather than at the whole function. This test fails on i386 and passes wherever XPFPA_HAVE_CW is 0, which includes x86-64 and arm64. The following commit fixes i386. --- .../tests/math/expm1_fpu_precision.phpt | 48 +++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 ext/standard/tests/math/expm1_fpu_precision.phpt diff --git a/ext/standard/tests/math/expm1_fpu_precision.phpt b/ext/standard/tests/math/expm1_fpu_precision.phpt new file mode 100644 index 000000000000..f2f664672c1a --- /dev/null +++ b/ext/standard/tests/math/expm1_fpu_precision.phpt @@ -0,0 +1,48 @@ +--TEST-- +expm1(): results must not lose precision when the FPU is clamped to double precision +--INI-- +serialize_precision=-1 +--FILE-- + +--EXPECT-- +-- reference points already correct on every platform -- +float(0) +float(6.38905609893065) +float(53.598150033144236) +-- inputs the FPU clamp gets wrong -- +float(19.085536923187668) +float(147.4131591025766) +float(402.4287934927351) +float(1095.6331584284585) +float(2979.9579870417283) +float(8102.083927575384) +float(22025.465794806718) +float(59873.14171519782) +float(162753.79141900392) +float(442412.3920089205) From 9cf6877e52f797bfcd014ce2bae16762c1ed3110 Mon Sep 17 00:00:00 2001 From: Marc Bennewitz Date: Fri, 4 Sep 2026 16:37:04 +0200 Subject: [PATCH 2/2] Use XPFPA for expm1() Restore extended FPU precision around the libm expm1() call and truncate the result back to double. Why this matters: zend_init_fpu() clamps the x87 FPU to a 53-bit mantissa for the whole request so that PHP's own double arithmetic behaves like IEEE-754 binary64. glibc's expm1(), however, computes on the x87 unit and relies on the extended range to round correctly to double, so the clamp costs expm1() accuracy. The loss is observable from userland and makes expm1() silently architecture-dependent, which breaks float comparisons, cached or serialized results, and any computation that accumulates the error. Measured against a high-precision reference rounded to nearest double, over 720 sampled inputs: i386 before this change 81/720 correctly rounded (11.2%) i386 after this change 714/720 correctly rounded (99.2%) arm64, unchanged 641/720 correctly rounded (89.0%) Note that arm64 is itself not a correctly rounded baseline for expm1(), so these figures are stated against the high-precision reference rather than against another platform. This is an improvement, not a guarantee: expm1(67) regresses by 1 ULP, because the clamped result there is genuinely the correctly rounded one. The accompanying test therefore asserts only inputs that this change fixes and that are also correct on x86-64 and arm64. zend_expm1() in zend_float.h compiles down to a plain expm1() wherever XPFPA_HAVE_CW is 0, which includes x86-64 and arm64, so no other platform is affected. --- Zend/zend_float.h | 19 +++++++++++++++++++ ext/standard/math.c | 2 +- 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/Zend/zend_float.h b/Zend/zend_float.h index 12aa02005c2a..c03f26ad52ea 100644 --- a/Zend/zend_float.h +++ b/Zend/zend_float.h @@ -20,6 +20,8 @@ #include "zend_portability.h" +#include + BEGIN_EXTERN_C() /* @@ -413,4 +415,21 @@ END_EXTERN_C() #endif /* FPU CONTROL */ +/* zend_init_fpu() clamps the x87 FPU to double precision (53-bit mantissa) for + * the whole request. libm's expm1() computes on the x87 unit and relies on the + * extended range to round correctly to double, so while that clamp is in place + * its result is off by one or more ULP. Restore extended precision for the call + * and truncate the result back to double. Compiles down to a plain expm1() + * everywhere XPFPA_HAVE_CW is 0, which includes x86-64 and arm64. */ +static zend_always_inline double zend_expm1(double num) +{ +#if XPFPA_HAVE_CW + XPFPA_DECLARE + XPFPA_SWITCH_DOUBLE_EXTENDED(); + XPFPA_RETURN_DOUBLE(expm1(num)); +#else + return expm1(num); +#endif +} + #endif diff --git a/ext/standard/math.c b/ext/standard/math.c index e13f153f63e2..bccf0dff533a 100644 --- a/ext/standard/math.c +++ b/ext/standard/math.c @@ -680,7 +680,7 @@ PHP_FUNCTION(expm1) Z_PARAM_DOUBLE(num) ZEND_PARSE_PARAMETERS_END(); - RETURN_DOUBLE(expm1(num)); + RETURN_DOUBLE(zend_expm1(num)); } /* }}} */