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)); } /* }}} */ 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)