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 | Corrected to handled appropriately random integers larger than 2^53 |
17 | converted to double precision floats, and to avoid overflows for |
18 | large gammas. |
19 +----------------------------------------------------------------------+
20 */
21
22 #ifdef HAVE_CONFIG_H
23 # include "config.h"
24 #endif
25
26 #include "php.h"
27 #include "php_random.h"
28 #include <math.h>
29
30 /* This file implements the γ-section algorithm as published in:
31 *
32 * Drawing Random Floating-Point Numbers from an Interval. Frédéric
33 * Goualard, ACM Trans. Model. Comput. Simul., 32:3, 2022.
34 * https://doi.org/10.1145/3503512
35 */
36
gamma_low(double x)37 static double gamma_low(double x)
38 {
39 return x - nextafter(x, -DBL_MAX);
40 }
41
gamma_high(double x)42 static double gamma_high(double x)
43 {
44 return nextafter(x, DBL_MAX) - x;
45 }
46
gamma_max(double x,double y)47 static double gamma_max(double x, double y)
48 {
49 return (fabs(x) > fabs(y)) ? gamma_high(x) : gamma_low(y);
50 }
51
splitint64(uint64_t v,double * vhi,double * vlo)52 static void splitint64(uint64_t v, double *vhi, double *vlo)
53 {
54 *vhi = v >> 2;
55 *vlo = v & UINT64_C(0x3);
56 }
57
ceilint(double a,double b,double g)58 static uint64_t ceilint(double a, double b, double g)
59 {
60 double s = b / g - a / g;
61 double e;
62
63 if (fabs(a) <= fabs(b)) {
64 e = -a / g - (s - b / g);
65 } else {
66 e = b / g - (s + a / g);
67 }
68
69 double si = ceil(s);
70
71 return (s != si) ? (uint64_t)si : (uint64_t)si + (e > 0);
72 }
73
php_random_gammasection_closed_open(php_random_algo_with_state engine,double min,double max)74 PHPAPI double php_random_gammasection_closed_open(php_random_algo_with_state engine, double min, double max)
75 {
76 double g = gamma_max(min, max);
77 uint64_t hi = ceilint(min, max, g);
78
79 if (UNEXPECTED(max <= min || hi < 1)) {
80 return NAN;
81 }
82
83 uint64_t k = 1 + php_random_range64(engine, hi - 1); /* [1, hi] */
84
85 if (fabs(min) <= fabs(max)) {
86 if (k == hi) {
87 return min;
88 } else {
89 double k_hi, k_lo;
90 splitint64(k, &k_hi, &k_lo);
91
92 return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g;
93 }
94 } else {
95 double k_hi, k_lo;
96 splitint64(k - 1, &k_hi, &k_lo);
97
98 return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g;
99 }
100 }
101
php_random_gammasection_closed_closed(php_random_algo_with_state engine,double min,double max)102 PHPAPI double php_random_gammasection_closed_closed(php_random_algo_with_state engine, double min, double max)
103 {
104 double g = gamma_max(min, max);
105 uint64_t hi = ceilint(min, max, g);
106
107 if (UNEXPECTED(max < min)) {
108 return NAN;
109 }
110
111 uint64_t k = php_random_range64(engine, hi); /* [0, hi] */
112
113 if (fabs(min) <= fabs(max)) {
114 if (k == hi) {
115 return min;
116 } else {
117 double k_hi, k_lo;
118 splitint64(k, &k_hi, &k_lo);
119
120 return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g;
121 }
122 } else {
123 if (k == hi) {
124 return max;
125 } else {
126 double k_hi, k_lo;
127 splitint64(k, &k_hi, &k_lo);
128
129 return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g;
130 }
131 }
132 }
133
php_random_gammasection_open_closed(php_random_algo_with_state engine,double min,double max)134 PHPAPI double php_random_gammasection_open_closed(php_random_algo_with_state engine, double min, double max)
135 {
136 double g = gamma_max(min, max);
137 uint64_t hi = ceilint(min, max, g);
138
139 if (UNEXPECTED(max <= min || hi < 1)) {
140 return NAN;
141 }
142
143 uint64_t k = php_random_range64(engine, hi - 1); /* [0, hi - 1] */
144
145 if (fabs(min) <= fabs(max)) {
146 double k_hi, k_lo;
147 splitint64(k, &k_hi, &k_lo);
148
149 return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g;
150 } else {
151 if (k == (hi - 1)) {
152 return max;
153 } else {
154 double k_hi, k_lo;
155 splitint64(k + 1, &k_hi, &k_lo);
156
157 return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g;
158 }
159 }
160 }
161
php_random_gammasection_open_open(php_random_algo_with_state engine,double min,double max)162 PHPAPI double php_random_gammasection_open_open(php_random_algo_with_state engine, double min, double max)
163 {
164 double g = gamma_max(min, max);
165 uint64_t hi = ceilint(min, max, g);
166
167 if (UNEXPECTED(max <= min || hi < 2)) {
168 return NAN;
169 }
170
171 uint64_t k = 1 + php_random_range64(engine, hi - 2); /* [1, hi - 1] */
172
173 if (fabs(min) <= fabs(max)) {
174 double k_hi, k_lo;
175 splitint64(k, &k_hi, &k_lo);
176
177 return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g;
178 } else {
179 double k_hi, k_lo;
180 splitint64(k, &k_hi, &k_lo);
181
182 return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g;
183 }
184 }
185
