@@ -0,0 +1,115 @@
1+/*
2+ +----------------------------------------------------------------------+
3+ | Copyright (c) The PHP Group |
4+ +----------------------------------------------------------------------+
5+ | This source file is subject to version 3.01 of the PHP license, |
6+ | that is bundled with this package in the file LICENSE, and is |
7+ | available through the world-wide-web at the following url: |
8+ | https://www.php.net/license/3_01.txt |
9+ | If you did not receive a copy of the PHP license and are unable to |
10+ | obtain it through the world-wide-web, please send a note to |
11+ | license@php.net so we can mail you a copy immediately. |
12+ +----------------------------------------------------------------------+
13+ | Authors: Tim Düsterhus <timwolla@php.net> |
14+ | |
15+ | Based on code from: Frédéric Goualard |
16+ +----------------------------------------------------------------------+
17+*/
18+19+#ifdef HAVE_CONFIG_H
20+# include "config.h"
21+#endif
22+23+#include "php.h"
24+#include "php_random.h"
25+#include <math.h>
26+27+/* This file implements the γ-section algorithm as published in:
28+ *
29+ * Drawing Random Floating-Point Numbers from an Interval. Frédéric
30+ * Goualard, ACM Trans. Model. Comput. Simul., 32:3, 2022.
31+ * https://doi.org/10.1145/3503512
32+ */
33+34+static double gamma_low(double x)
35+{
36+return x - nextafter(x, -DBL_MAX);
37+}
38+39+static double gamma_high(double x)
40+{
41+return nextafter(x, DBL_MAX) - x;
42+}
43+44+static double gamma_max(double x, double y)
45+{
46+return (fabs(x) > fabs(y)) ? gamma_high(x) : gamma_low(y);
47+}
48+49+static uint64_t ceilint(double a, double b, double g)
50+{
51+double s = b / g - a / g;
52+double e;
53+54+if (fabs(a) <= fabs(b)) {
55+e = -a / g - (s - b / g);
56+ } else {
57+e = b / g - (s + a / g);
58+ }
59+60+double si = ceil(s);
61+62+return (s != si) ? (uint64_t)si : (uint64_t)si + (e > 0);
63+}
64+65+PHPAPI double php_random_gammasection_closed_open(const php_random_algo *algo, php_random_status *status, double min, double max)
66+{
67+double g = gamma_max(min, max);
68+uint64_t hi = ceilint(min, max, g);
69+uint64_t k = 1 + php_random_range64(algo, status, hi - 1); /* [1, hi] */
70+71+if (fabs(min) <= fabs(max)) {
72+return k == hi ? min : max - k * g;
73+ } else {
74+return min + (k - 1) * g;
75+ }
76+}
77+78+PHPAPI double php_random_gammasection_closed_closed(const php_random_algo *algo, php_random_status *status, double min, double max)
79+{
80+double g = gamma_max(min, max);
81+uint64_t hi = ceilint(min, max, g);
82+uint64_t k = php_random_range64(algo, status, hi); /* [0, hi] */
83+84+if (fabs(min) <= fabs(max)) {
85+return k == hi ? min : max - k * g;
86+ } else {
87+return k == hi ? max : min + k * g;
88+ }
89+}
90+91+PHPAPI double php_random_gammasection_open_closed(const php_random_algo *algo, php_random_status *status, double min, double max)
92+{
93+double g = gamma_max(min, max);
94+uint64_t hi = ceilint(min, max, g);
95+uint64_t k = php_random_range64(algo, status, hi - 1); /* [0, hi - 1] */
96+97+if (fabs(min) <= fabs(max)) {
98+return max - k * g;
99+ } else {
100+return k == (hi - 1) ? max : min + (k + 1) * g;
101+ }
102+}
103+104+PHPAPI double php_random_gammasection_open_open(const php_random_algo *algo, php_random_status *status, double min, double max)
105+{
106+double g = gamma_max(min, max);
107+uint64_t hi = ceilint(min, max, g);
108+uint64_t k = 1 + php_random_range64(algo, status, hi - 2); /* [1, hi - 1] */
109+110+if (fabs(min) <= fabs(max)) {
111+return max - k * g;
112+ } else {
113+return min + k * g;
114+ }
115+}