Fix pow(), ** and fpow() losing precision on 32-bit x86

Zend/tests/pow_fpu_precision.phpt

@@ -0,0 +1,73 @@
+--TEST--
+The ** operator must not lose precision when the FPU is clamped to double precision
+--INI--
+serialize_precision=-1
+--FILE--
+<?php
+/* 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 pow()
+ * 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.
+ *
+ * This file covers the ** operator (ZEND_POW). The same cases are covered for
+ * pow() in ext/standard/tests/math/pow_fpu_precision.phpt.
+ *
+ * Every expected value is exactly the correctly rounded double of its decimal
+ * literal, so a build that keeps full precision prints the short literal back
+ * and compares identical to it. */
+
+$ten = 10;
+$five = 5;
+$two = 2;
+
+echo "-- runtime --\n";
+var_dump($ten ** -2);
+var_dump($ten ** -3);
+var_dump($ten ** -4);
+var_dump($five ** -2);
+var_dump($five ** -3);
+var_dump($five ** -4);
+echo "-- compile-time constant folding --\n";
+var_dump(10 ** -2);
+var_dump(10 ** -3);
+var_dump(10 ** -4);
+var_dump(5 ** -2);
+var_dump(5 ** -3);
+var_dump(5 ** -4);
+echo "-- exact powers of two are unaffected --\n";
+var_dump($two ** -1);
+var_dump($two ** -10);
+echo "-- identical to the decimal literal --\n";
+var_dump($ten ** -2 === 0.01);
+var_dump($ten ** -3 === 0.001);
+var_dump($ten ** -4 === 0.0001);
+var_dump($five ** -2 === 0.04);
+var_dump($five ** -3 === 0.008);
+var_dump($five ** -4 === 0.0016);
+?>
+--EXPECT--
+-- runtime --
+float(0.01)
+float(0.001)
+float(0.0001)
+float(0.04)
+float(0.008)
+float(0.0016)
+-- compile-time constant folding --
+float(0.01)
+float(0.001)
+float(0.0001)
+float(0.04)
+float(0.008)
+float(0.0016)
+-- exact powers of two are unaffected --
+float(0.5)
+float(0.0009765625)
+-- identical to the decimal literal --
+bool(true)
+bool(true)
+bool(true)
+bool(true)
+bool(true)
+bool(true)

Zend/zend_float.h

@@ -20,6 +20,8 @@
 
 #include "zend_portability.h"
 
+#include <math.h>
+
 BEGIN_EXTERN_C()
 
 /*
@@ -413,4 +415,22 @@ 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 pow() relies on the FPU running at extended
+ * precision internally, so while that clamp is in place its result is off by
+ * one or more ULP compared to platforms whose FPU has no precision control.
+ * Restore extended precision for the call and truncate the result back to
+ * double. Compiles down to a plain pow() everywhere XPFPA_HAVE_CW is 0, which
+ * includes x86-64 and arm64. */
+static zend_always_inline double zend_pow(double base, double exponent)
+{
+#if XPFPA_HAVE_CW
+	XPFPA_DECLARE
+	XPFPA_SWITCH_DOUBLE_EXTENDED();
+	XPFPA_RETURN_DOUBLE(pow(base, exponent));
+#else
+	return pow(base, exponent);
+#endif
+}
+
 #endif

Zend/zend_operators.c

@@ -1347,7 +1347,7 @@ static double safe_pow(double base, double exponent)
 		zend_power_base_0_exponent_lt_0_error();
 	}
 
-	return pow(base, exponent);
+	return zend_pow(base, exponent);
 }
 
 static zend_result ZEND_FASTCALL pow_function_base(zval *result, zval *op1, zval *op2) /* {{{ */

ext/standard/math.c

@@ -1485,7 +1485,7 @@ PHP_FUNCTION(fpow)
 		Z_PARAM_DOUBLE(exponent)
 	ZEND_PARSE_PARAMETERS_END();
 
-	RETURN_DOUBLE(pow(base, exponent));
+	RETURN_DOUBLE(zend_pow(base, exponent));
 }
 /* }}} */
 

ext/standard/tests/math/fpow_fpu_precision.phpt

@@ -0,0 +1,75 @@
+--TEST--
+fpow(): results must not lose precision when the FPU is clamped to double precision
+--INI--
+serialize_precision=-1
+--FILE--
+<?php
+/* fpow() calls libm pow() directly instead of going through the ** operator's
+ * safe_pow(), so it needs the same extended-precision handling: on i386 PHP
+ * clamps the x87 FPU to a 53-bit mantissa for the whole request, which costs
+ * pow() one or more ULP. See ext/standard/tests/math/pow_fpu_precision.phpt. */
+
+$ten = 10.0;
+$five = 5.0;
+$two = 2.0;
+
+echo "-- precision --\n";
+var_dump(fpow($ten, -2.0));
+var_dump(fpow($ten, -3.0));
+var_dump(fpow($ten, -4.0));
+var_dump(fpow($five, -2.0));
+var_dump(fpow($five, -3.0));
+var_dump(fpow($five, -4.0));
+
+echo "-- exact powers of two are unaffected --\n";
+var_dump(fpow($two, -1.0));
+var_dump(fpow($two, -10.0));
+
+echo "-- identical to the decimal literal --\n";
+var_dump(fpow($ten, -2.0) === 0.01);
+var_dump(fpow($ten, -3.0) === 0.001);
+var_dump(fpow($ten, -4.0) === 0.0001);
+var_dump(fpow($five, -2.0) === 0.04);
+var_dump(fpow($five, -3.0) === 0.008);
+var_dump(fpow($five, -4.0) === 0.0016);
+
+echo "-- matches the ** operator --\n";
+var_dump(fpow($ten, -2.0) === $ten ** -2);
+var_dump(fpow($five, -4.0) === $five ** -4);
+
+echo "-- IEEE-754 semantics are preserved --\n";
+var_dump(fpow(0.0, -1.0));
+var_dump(fpow(-0.0, -3.0));
+var_dump(fpow(1.0, NAN));
+var_dump(fpow(NAN, 0.0));
+var_dump(fpow(-1.0, INF));
+var_dump(fpow(-8.0, 1 / 3));
+?>
+--EXPECT--
+-- precision --
+float(0.01)
+float(0.001)
+float(0.0001)
+float(0.04)
+float(0.008)
+float(0.0016)
+-- exact powers of two are unaffected --
+float(0.5)
+float(0.0009765625)
+-- identical to the decimal literal --
+bool(true)
+bool(true)
+bool(true)
+bool(true)
+bool(true)
+bool(true)
+-- matches the ** operator --
+bool(true)
+bool(true)
+-- IEEE-754 semantics are preserved --
+float(INF)
+float(-INF)
+float(1)
+float(1)
+float(1)
+float(NAN)

ext/standard/tests/math/pow_fpu_precision.phpt

@@ -0,0 +1,73 @@
+--TEST--
+pow() must not lose precision when the FPU is clamped to double precision
+--INI--
+serialize_precision=-1
+--FILE--
+<?php
+/* 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 pow()
+ * 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.
+ *
+ * This file covers the pow() function. The same cases are covered for
+ * the ** operator in Zend/tests/pow_fpu_precision.phpt.
+ *
+ * Every expected value is exactly the correctly rounded double of its decimal
+ * literal, so a build that keeps full precision prints the short literal back
+ * and compares identical to it. */
+
+$ten = 10;
+$five = 5;
+$two = 2;
+
+echo "-- runtime --\n";
+var_dump(pow($ten, -2));
+var_dump(pow($ten, -3));
+var_dump(pow($ten, -4));
+var_dump(pow($five, -2));
+var_dump(pow($five, -3));
+var_dump(pow($five, -4));
+echo "-- compile-time constant folding --\n";
+var_dump(pow(10, -2));
+var_dump(pow(10, -3));
+var_dump(pow(10, -4));
+var_dump(pow(5, -2));
+var_dump(pow(5, -3));
+var_dump(pow(5, -4));
+echo "-- exact powers of two are unaffected --\n";
+var_dump(pow($two, -1));
+var_dump(pow($two, -10));
+echo "-- identical to the decimal literal --\n";
+var_dump(pow($ten, -2) === 0.01);
+var_dump(pow($ten, -3) === 0.001);
+var_dump(pow($ten, -4) === 0.0001);
+var_dump(pow($five, -2) === 0.04);
+var_dump(pow($five, -3) === 0.008);
+var_dump(pow($five, -4) === 0.0016);
+?>
+--EXPECT--
+-- runtime --
+float(0.01)
+float(0.001)
+float(0.0001)
+float(0.04)
+float(0.008)
+float(0.0016)
+-- compile-time constant folding --
+float(0.01)
+float(0.001)
+float(0.0001)
+float(0.04)
+float(0.008)
+float(0.0016)
+-- exact powers of two are unaffected --
+float(0.5)
+float(0.0009765625)
+-- identical to the decimal literal --
+bool(true)
+bool(true)
+bool(true)
+bool(true)
+bool(true)
+bool(true)