From cd13e7f6889dc8b8580d446820a1b36c24f4dbfa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tim=20D=C3=BCsterhus?= Date: Tue, 10 Oct 2023 16:34:25 +0200 Subject: [PATCH] =?UTF-8?q?random:=20Fix=20=CE=B3-section=20implementation?= =?UTF-8?q?=20for=20Randomizer::getFloat()?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 | 4 ++ ext/random/gammasection.c | 65 ++++++++++++++++--- .../methods/getFloat_extreme_range.phpt | 41 ++++++++++++ .../methods/getFloat_opposite_signs.phpt | 52 +++++++++++++++ 4 files changed, 154 insertions(+), 8 deletions(-) create mode 100644 ext/random/tests/03_randomizer/methods/getFloat_extreme_range.phpt create mode 100644 ext/random/tests/03_randomizer/methods/getFloat_opposite_signs.phpt diff --git a/NEWS b/NEWS index 01984977e236f..f2a819d33963f 100644 --- a/NEWS +++ b/NEWS @@ -35,6 +35,10 @@ PHP NEWS . Fixed bug GH-12380 (JIT+private array property access inside closure accesses private property in child class). (nielsdos) +- Random: + . Fix Randomizer::getFloat() returning incorrect results under + certain circumstances. (timwolla) + - SimpleXML: . Apply iterator fixes only on master. (nielsdos) diff --git a/ext/random/gammasection.c b/ext/random/gammasection.c index aa4531fba22f7..79bf63c866992 100644 --- a/ext/random/gammasection.c +++ b/ext/random/gammasection.c @@ -13,6 +13,9 @@ | Authors: Tim Düsterhus | | | | 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; } } diff --git a/ext/random/tests/03_randomizer/methods/getFloat_extreme_range.phpt b/ext/random/tests/03_randomizer/methods/getFloat_extreme_range.phpt new file mode 100644 index 0000000000000..f71960510c823 --- /dev/null +++ b/ext/random/tests/03_randomizer/methods/getFloat_extreme_range.phpt @@ -0,0 +1,41 @@ +--TEST-- +Random: Randomizer: getFloat(): Extreme ranges are handled correctly +--FILE-- +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 diff --git a/ext/random/tests/03_randomizer/methods/getFloat_opposite_signs.phpt b/ext/random/tests/03_randomizer/methods/getFloat_opposite_signs.phpt new file mode 100644 index 0000000000000..f9d2fae77804e --- /dev/null +++ b/ext/random/tests/03_randomizer/methods/getFloat_opposite_signs.phpt @@ -0,0 +1,52 @@ +--TEST-- +Random: Randomizer: getFloat(): Opposite signs are handled correctly +--FILE-- +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)