random: Fix γ-section implementation for Randomizer::getFloat() (#12402)

The reference implementation of the "Drawing Random Floating-Point Numbers from
an Interval" paper contains two mistakes that will result in erroneous values
being returned under certain circumstances:

- For large values of `g` the multiplication of `k * g` might overflow to
  infinity.
- The value of `ceilint()` might exceed 2^53, possibly leading to a rounding
  error when promoting `k` to double within the multiplication of `k * g`.

This commit updates the implementation based on Prof. Goualard suggestions
after reaching out to him. It will correctly handle inputs larger than 2^-1020
in absolute values. This limitation will be documented and those inputs
possibly be rejected in a follow-up commit depending on performance concerns.

NEWS

@@ -13,6 +13,10 @@ PHP                                                                        NEWS
   . Fixed bug GH-8143 (Crashes in zend_accel_inheritance_cache_find since
     upgrading to 8.1.3 due to corrupt on-disk file cache). (turchanov)
 
+- Random:
+  . Fix Randomizer::getFloat() returning incorrect results under
+    certain circumstances. (timwolla)
+
 - SOAP:
   . Fixed bug GH-12392 (Segmentation fault on SoapClient::__getTypes).
     (nielsdos)

ext/random/gammasection.c

@@ -13,6 +13,9 @@
    | Authors: Tim Düsterhus <timwolla@php.net>                            |
    |                                                                      |
    | Based on code from: Frédéric Goualard                                |
+   | Corrected to handled appropriately random integers larger than 2^53  |
+   | converted to double precision floats, and to avoid overflows for     |
+   | large gammas.                                                        |
    +----------------------------------------------------------------------+
 */
 
@@ -46,6 +49,12 @@ static double gamma_max(double x, double y)
 	return (fabs(x) > fabs(y)) ? gamma_high(x) : gamma_low(y);
 }
 
+static void splitint64(uint64_t v, double *vhi, double *vlo)
+{
+	*vhi = v >> 2;
+	*vlo = v & UINT64_C(0x3);
+}
+
 static uint64_t ceilint(double a, double b, double g)
 {
 	double s = b / g - a / g;
@@ -74,9 +83,19 @@ PHPAPI double php_random_gammasection_closed_open(const php_random_algo *algo, p
 	uint64_t k = 1 + php_random_range64(algo, status, hi - 1); /* [1, hi] */
 
 	if (fabs(min) <= fabs(max)) {
-		return k == hi ? min : max - k * g;
+		if (k == hi) {
+			return min;
+		} else {
+			double k_hi, k_lo;
+			splitint64(k, &k_hi, &k_lo);
+
+			return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g;
+		}
 	} else {
-		return min + (k - 1) * g;
+		double k_hi, k_lo;
+		splitint64(k - 1, &k_hi, &k_lo);
+
+		return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g;
 	}
 }
 
@@ -92,9 +111,23 @@ PHPAPI double php_random_gammasection_closed_closed(const php_random_algo *algo,
 	uint64_t k = php_random_range64(algo, status, hi); /* [0, hi] */
 
 	if (fabs(min) <= fabs(max)) {
-		return k == hi ? min : max - k * g;
+		if (k == hi) {
+			return min;
+		} else {
+			double k_hi, k_lo;
+			splitint64(k, &k_hi, &k_lo);
+
+			return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g;
+		}
 	} else {
-		return k == hi ? max : min + k * g;
+		if (k == hi) {
+			return max;
+		} else {
+			double k_hi, k_lo;
+			splitint64(k, &k_hi, &k_lo);
+
+			return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g;
+		}
 	}
 }
 
@@ -110,9 +143,19 @@ PHPAPI double php_random_gammasection_open_closed(const php_random_algo *algo, p
 	uint64_t k = php_random_range64(algo, status, hi - 1); /* [0, hi - 1] */
 
 	if (fabs(min) <= fabs(max)) {
-		return max - k * g;
+		double k_hi, k_lo;
+		splitint64(k, &k_hi, &k_lo);
+
+		return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g;
 	} else {
-		return k == (hi - 1) ? max : min + (k + 1) * g;
+		if (k == (hi - 1)) {
+			return max;
+		} else {
+			double k_hi, k_lo;
+			splitint64(k + 1, &k_hi, &k_lo);
+
+			return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g;
+		}
 	}
 }
 
@@ -128,8 +171,14 @@ PHPAPI double php_random_gammasection_open_open(const php_random_algo *algo, php
 	uint64_t k = 1 + php_random_range64(algo, status, hi - 2); /* [1, hi - 1] */
 
 	if (fabs(min) <= fabs(max)) {
-		return max - k * g;
+		double k_hi, k_lo;
+		splitint64(k, &k_hi, &k_lo);
+
+		return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g;
 	} else {
-		return min + k * g;
+		double k_hi, k_lo;
+		splitint64(k, &k_hi, &k_lo);
+
+		return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g;
 	}
 }

ext/random/tests/03_randomizer/methods/getFloat_extreme_range.phpt

@@ -0,0 +1,41 @@
+--TEST--
+Random: Randomizer: getFloat(): Extreme ranges are handled correctly
+--FILE--
+<?php
+
+use Random\Engine;
+use Random\Randomizer;
+
+final class TestEngine implements Engine
+{
+    private array $values = [
+        "\x64\xe8\x58\x79\x3e\xf6\x38\x00",
+        "\x65\xe8\x58\x79\x3e\xf6\x38\x00",
+        "\x00\x00\x00\x00\x00\x00\x00\x00",
+    ];
+
+    public function generate(): string
+    {
+        return \array_shift($this->values);
+    }
+}
+
+$r = new Randomizer(new TestEngine());
+
+$min = -1.6e308;
+$max = 1.6e308;
+printf("%.17g\n", $min);
+printf("%.17g\n\n", $max);
+
+printf("%.17g\n", $r->getFloat($min, $max));
+printf("%.17g\n", $r->getFloat($min, $max));
+printf("%.17g\n", $r->getFloat($min, $max));
+
+?>
+--EXPECT--
+-1.6e+308
+1.6e+308
+
+-1.5999999999999998e+308
+-1.6e+308
+1.5999999999999998e+308

ext/random/tests/03_randomizer/methods/getFloat_opposite_signs.phpt

@@ -0,0 +1,52 @@
+--TEST--
+Random: Randomizer: getFloat(): Opposite signs are handled correctly
+--FILE--
+<?php
+
+use Random\Engine;
+use Random\Randomizer;
+
+final class TestEngine implements Engine
+{
+    private array $values = [
+        "\x63\xe8\x58\x79\x3e\xf6\x38\x00",
+        "\x64\xe8\x58\x79\x3e\xf6\x38\x00",
+        "\x65\xe8\x58\x79\x3e\xf6\x38\x00",
+        "\x66\xe8\x58\x79\x3e\xf6\x38\x00",
+        "\x67\xe8\x58\x79\x3e\xf6\x38\x00",
+    ];
+
+    public function generate(): string
+    {
+        return \array_shift($this->values);
+    }
+}
+
+$r = new Randomizer(new TestEngine());
+
+$min = -1;
+$max = 1;
+printf("%.17g\n", $min);
+printf("%.17g\n\n", $max);
+
+$prev = null;
+for ($i = 0; $i < 5; $i++) {
+    $float = $r->getFloat($min, $max);
+    printf("%.17f", $float);
+    if ($prev !== null) {
+        printf(" (%+.17g)", ($float - $prev));
+    }
+    printf("\n");
+
+    $prev = $float;
+}
+?>
+--EXPECT--
+-1
+1
+
+-0.78005908680576086
+-0.78005908680576097 (-1.1102230246251565e-16)
+-0.78005908680576108 (-1.1102230246251565e-16)
+-0.78005908680576119 (-1.1102230246251565e-16)
+-0.78005908680576130 (-1.1102230246251565e-16)