Fix GH-12143: Optimize round

SakiTakamachi · SakiTakamachi · commit 78970ef6b2f4 · 2024-02-03T22:23:02.000+09:00
Fixed an error in the result due to "pre-rounding" of the round function. "Pre-rounding" has been abolished and the method of comparing numbers has been changed. Closes GH-12268.
diff --git a/NEWS b/NEWS
@@ -153,6 +153,7 @@ Standard:
     the precision is not lost. (Marc Bennewitz)
   . Add support for 4 new rounding modes to the round() function. (Jorg Sowa)
   . debug_zval_dump() now indicates whether an array is packed. (Max Semenik)
+  . Fix GH-12143 (Optimize round). (SakiTakamachi)
 
 XML:
   . Added XML_OPTION_PARSE_HUGE parser option. (nielsdos)
diff --git a/UPGRADING b/UPGRADING
@@ -343,6 +343,11 @@ PDO_SQLITE:
 
     RFC: https://wiki.php.net/rfc/new_rounding_modes_to_round_function
   . debug_zval_dump() now indicates whether an array is packed.
+  . Fixed a bug caused by "pre-rounding" of the round() function. Previously, using
+    "pre-rounding" to treat a value like 0.285 (actually 0.28499999999999998) as a
+    decimal number and round it to 0.29. However, "pre-rounding" incorrectly rounds
+    certain numbers, so this fix removes "pre-rounding" and changes the way numbers
+    are compared, so that the values ​​are correctly rounded as decimal numbers.
 
 ========================================
 6. New Functions
diff --git a/ext/standard/math.c b/ext/standard/math.c
@@ -27,53 +27,10 @@
 #include <float.h>
 #include <math.h>
 #include <stdlib.h>
+#include <fenv.h>
 
 #include "basic_functions.h"
 
-/* {{{ php_intlog10abs
-   Returns floor(log10(fabs(val))), uses fast binary search */
-static inline int php_intlog10abs(double value) {
-	value = fabs(value);
-
-	if (value < 1e-8 || value > 1e22) {
-		return (int)floor(log10(value));
-	} else {
-		/* Do a binary search with 5 steps */
-		int result = 15;
-		static const double values[] = {
-				1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2,
-				1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13,
-				1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22};
-
-		if (value < values[result]) {
-			result -= 8;
-		} else {
-			result += 8;
-		}
-		if (value < values[result]) {
-			result -= 4;
-		} else {
-			result += 4;
-		}
-		if (value < values[result]) {
-			result -= 2;
-		} else {
-			result += 2;
-		}
-		if (value < values[result]) {
-			result -= 1;
-		} else {
-			result += 1;
-		}
-		if (value < values[result]) {
-			result -= 1;
-		}
-		result -= 8;
-		return result;
-	}
-}
-/* }}} */
-
 /* {{{ php_intpow10
        Returns pow(10.0, (double)power), uses fast lookup table for exact powers */
 static inline double php_intpow10(int power) {
@@ -90,22 +47,30 @@ static inline double php_intpow10(int power) {
 }
 /* }}} */
 
-/* {{{ php_round_helper
-       Actually performs the rounding of a value to integer in a certain mode */
-static inline double php_round_helper(double value, int mode) {
-	double integral, fractional;
+static zend_always_inline double php_round_get_basic_edge_case(double integral, double exponent, int places)
+{
+	return (places > 0)
+		? fabs((integral + copysign(0.5, integral)) / exponent)
+		: fabs((integral + copysign(0.5, integral)) * exponent);
+}
 
-	/* Split the input value into the integral and fractional part.
-	 *
-	 * Both parts will have the same sign as the input value. We take
-	 * the absolute value of the fractional part (which will not result
-	 * in branches in the assembly) to make the following cases simpler.
-	 */
-	fractional = fabs(modf(value, &integral));
+static zend_always_inline double php_round_get_zero_edge_case(double integral, double exponent, int places)
+{
+	return (places > 0)
+		? fabs((integral) / exponent)
+		: fabs((integral) * exponent);
+}
+
+/* {{{ php_round_helper
+	   Actually performs the rounding of a value to integer in a certain mode */
+static inline double php_round_helper(double integral, double value, double exponent, int places, int mode) {
+	double value_abs = fabs(value);
+	double edge_case;
 
 	switch (mode) {
 		case PHP_ROUND_HALF_UP:
-			if (fractional >= 0.5) {
+			edge_case = php_round_get_basic_edge_case(integral, exponent, places);
+			if (value_abs >= edge_case) {
 				/* We must increase the magnitude of the integral part
 				 * (rounding up / towards infinity). copysign(1.0, integral)
 				 * will either result in 1.0 or -1.0 depending on the sign
@@ -120,21 +85,24 @@ static inline double php_round_helper(double value, int mode) {
 			return integral;
 
 		case PHP_ROUND_HALF_DOWN:
-			if (fractional > 0.5) {
+			edge_case = php_round_get_basic_edge_case(integral, exponent, places);
+			if (value_abs > edge_case) {
 				return integral + copysign(1.0, integral);
 			}
 
 			return integral;
 
 		case PHP_ROUND_CEILING:
-			if (value > 0.0 && fractional > 0.0) {
+			edge_case = php_round_get_zero_edge_case(integral, exponent, places);
+			if (value > 0.0 && value_abs > edge_case) {
 				return integral + 1.0;
 			}
 
 			return integral;
 
 		case PHP_ROUND_FLOOR:
-			if (value < 0.0 && fractional > 0.0) {
+			edge_case = php_round_get_zero_edge_case(integral, exponent, places);
+			if (value < 0.0 && value_abs > edge_case) {
 				return integral - 1.0;
 			}
 
@@ -144,18 +112,18 @@ static inline double php_round_helper(double value, int mode) {
 			return integral;
 
 		case PHP_ROUND_AWAY_FROM_ZERO:
-			if (fractional > 0.0) {
+			edge_case = php_round_get_zero_edge_case(integral, exponent, places);
+			if (value_abs > edge_case) {
 				return integral + copysign(1.0, integral);
 			}
 
 			return integral;
 
 		case PHP_ROUND_HALF_EVEN:
-			if (fractional > 0.5) {
+			edge_case = php_round_get_basic_edge_case(integral, exponent, places);
+			if (value_abs > edge_case) {
 				return integral + copysign(1.0, integral);
-			}
-
-			if (UNEXPECTED(fractional == 0.5)) {
+			} else if (UNEXPECTED(value_abs == edge_case)) {
 				bool even = !fmod(integral, 2.0);
 
 				/* If the integral part is not even we can make it even
@@ -169,11 +137,10 @@ static inline double php_round_helper(double value, int mode) {
 			return integral;
 
 		case PHP_ROUND_HALF_ODD:
-			if (fractional > 0.5) {
+			edge_case = php_round_get_basic_edge_case(integral, exponent, places);
+			if (value_abs > edge_case) {
 				return integral + copysign(1.0, integral);
-			}
-
-			if (UNEXPECTED(fractional == 0.5)) {
+			} else if (UNEXPECTED(value_abs == edge_case)) {
 				bool even = !fmod(integral, 2.0);
 
 				if (even) {
@@ -196,63 +163,55 @@ static inline double php_round_helper(double value, int mode) {
  * mode. For the specifics of the algorithm, see http://wiki.php.net/rfc/rounding
  */
 PHPAPI double _php_math_round(double value, int places, int mode) {
-	double f1, f2;
+	double exponent;
 	double tmp_value;
-	int precision_places;
+	int cpu_round_mode;
 
 	if (!zend_finite(value) || value == 0.0) {
 		return value;
 	}
 
 	places = places < INT_MIN+1 ? INT_MIN+1 : places;
-	precision_places = 14 - php_intlog10abs(value);
 
-	f1 = php_intpow10(abs(places));
+	exponent = php_intpow10(abs(places));
 
-	/* If the decimal precision guaranteed by FP arithmetic is higher than
-	   the requested places BUT is small enough to make sure a non-zero value
-	   is returned, pre-round the result to the precision */
-	if (precision_places > places && precision_places - 15 < places) {
-		int64_t use_precision = precision_places < INT_MIN+1 ? INT_MIN+1 : precision_places;
-
-		f2 = php_intpow10(abs((int)use_precision));
-		if (use_precision >= 0) {
-			tmp_value = value * f2;
-		} else {
-			tmp_value = value / f2;
-		}
-		/* preround the result (tmp_value will always be something * 1e14,
-		   thus never larger than 1e15 here) */
-		tmp_value = php_round_helper(tmp_value, mode);
-
-		use_precision = places - precision_places;
-		use_precision = use_precision < INT_MIN+1 ? INT_MIN+1 : use_precision;
-		/* now correctly move the decimal point */
-		f2 = php_intpow10(abs((int)use_precision));
-		/* because places < precision_places */
-		tmp_value = tmp_value / f2;
+	/**
+	 * When extracting the integer part, the result may be incorrect as a decimal
+	 * number due to floating point errors.
+	 * e.g.
+	 * 0.285 * 10000000000 => 2849999999.9999995
+	 * floor(0.285 * 10000000000) => 2849999999
+	 *
+	 * Therefore, change the CPU rounding mode to away from 0 only from
+	 * fegetround to fesetround.
+	 * e.g.
+	 * 0.285 * 10000000000 => 2850000000.0
+	 * floor(0.285 * 10000000000) => 2850000000
+	 */
+	cpu_round_mode = fegetround();
+	if (value >= 0.0) {
+		fesetround(FE_UPWARD);
+		tmp_value = floor(places > 0  ? value * exponent : value / exponent);
 	} else {
-		/* adjust the value */
-		if (places >= 0) {
-			tmp_value = value * f1;
-		} else {
-			tmp_value = value / f1;
-		}
-		/* This value is beyond our precision, so rounding it is pointless */
-		if (fabs(tmp_value) >= 1e15) {
-			return value;
-		}
+		fesetround(FE_DOWNWARD);
+		tmp_value = ceil(places > 0  ? value * exponent : value / exponent);
+	}
+	fesetround(cpu_round_mode);
+
+	/* This value is beyond our precision, so rounding it is pointless */
+	if (fabs(tmp_value) >= 1e15) {
+		return value;
 	}
 
 	/* round the temp value */
-	tmp_value = php_round_helper(tmp_value, mode);
+	tmp_value = php_round_helper(tmp_value, value, exponent, places, mode);
 
 	/* see if it makes sense to use simple division to round the value */
 	if (abs(places) < 23) {
 		if (places > 0) {
-			tmp_value = tmp_value / f1;
+			tmp_value = tmp_value / exponent;
 		} else {
-			tmp_value = tmp_value * f1;
+			tmp_value = tmp_value * exponent;
 		}
 	} else {
 		/* Simple division can't be used since that will cause wrong results.
@@ -272,7 +231,6 @@ PHPAPI double _php_math_round(double value, int places, int mode) {
 			tmp_value = value;
 		}
 	}
-
 	return tmp_value;
 }
 /* }}} */
diff --git a/ext/standard/tests/math/bug24142.phpt b/ext/standard/tests/math/bug24142.phpt
@@ -2,19 +2,65 @@
 Bug #24142 (round() problems)
 --FILE--
 <?php
-$v = 0.005;
-for ($i = 1; $i < 10; $i++) {
-    echo "round({$v}, 2) -> ".round($v, 2)."\n";
-    $v += 0.01;
-}
+echo "round(0.005, 2)\n";
+var_dump(round(0.005, 2));
+echo "\n";
+
+echo "round(0.015, 2)\n";
+var_dump(round(0.015, 2));
+echo "\n";
+
+echo "round(0.025, 2)\n";
+var_dump(round(0.025, 2));
+echo "\n";
+
+echo "round(0.035, 2)\n";
+var_dump(round(0.035, 2));
+echo "\n";
+
+echo "round(0.045, 2)\n";
+var_dump(round(0.045, 2));
+echo "\n";
+
+echo "round(0.055, 2)\n";
+var_dump(round(0.055, 2));
+echo "\n";
+
+echo "round(0.065, 2)\n";
+var_dump(round(0.065, 2));
+echo "\n";
+
+echo "round(0.075, 2)\n";
+var_dump(round(0.075, 2));
+echo "\n";
+
+echo "round(0.085, 2)\n";
+var_dump(round(0.085, 2));
 ?>
 --EXPECT--
-round(0.005, 2) -> 0.01
-round(0.015, 2) -> 0.02
-round(0.025, 2) -> 0.03
-round(0.035, 2) -> 0.04
-round(0.045, 2) -> 0.05
-round(0.055, 2) -> 0.06
-round(0.065, 2) -> 0.07
-round(0.075, 2) -> 0.08
-round(0.085, 2) -> 0.09
+round(0.005, 2)
+float(0.01)
+
+round(0.015, 2)
+float(0.02)
+
+round(0.025, 2)
+float(0.03)
+
+round(0.035, 2)
+float(0.04)
+
+round(0.045, 2)
+float(0.05)
+
+round(0.055, 2)
+float(0.06)
+
+round(0.065, 2)
+float(0.07)
+
+round(0.075, 2)
+float(0.08)
+
+round(0.085, 2)
+float(0.09)
diff --git a/ext/standard/tests/math/round_gh12143_optimize_round.phpt b/ext/standard/tests/math/round_gh12143_optimize_round.phpt
