/src/php-src/ext/random/gammasection.c
Line | Count | Source (jump to first uncovered line) |
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 | | |
37 | | static double gamma_low(double x) |
38 | 0 | { |
39 | 0 | return x - nextafter(x, -DBL_MAX); |
40 | 0 | } |
41 | | |
42 | | static double gamma_high(double x) |
43 | 0 | { |
44 | 0 | return nextafter(x, DBL_MAX) - x; |
45 | 0 | } |
46 | | |
47 | | static double gamma_max(double x, double y) |
48 | 0 | { |
49 | 0 | return (fabs(x) > fabs(y)) ? gamma_high(x) : gamma_low(y); |
50 | 0 | } |
51 | | |
52 | | static void splitint64(uint64_t v, double *vhi, double *vlo) |
53 | 0 | { |
54 | 0 | *vhi = v >> 2; |
55 | 0 | *vlo = v & UINT64_C(0x3); |
56 | 0 | } |
57 | | |
58 | | static uint64_t ceilint(double a, double b, double g) |
59 | 0 | { |
60 | 0 | double s = b / g - a / g; |
61 | 0 | double e; |
62 | |
|
63 | 0 | if (fabs(a) <= fabs(b)) { |
64 | 0 | e = -a / g - (s - b / g); |
65 | 0 | } else { |
66 | 0 | e = b / g - (s + a / g); |
67 | 0 | } |
68 | |
|
69 | 0 | double si = ceil(s); |
70 | |
|
71 | 0 | return (s != si) ? (uint64_t)si : (uint64_t)si + (e > 0); |
72 | 0 | } |
73 | | |
74 | | PHPAPI double php_random_gammasection_closed_open(php_random_algo_with_state engine, double min, double max) |
75 | 0 | { |
76 | 0 | double g = gamma_max(min, max); |
77 | 0 | uint64_t hi = ceilint(min, max, g); |
78 | |
|
79 | 0 | if (UNEXPECTED(max <= min || hi < 1)) { |
80 | 0 | return NAN; |
81 | 0 | } |
82 | | |
83 | 0 | uint64_t k = 1 + php_random_range64(engine, hi - 1); /* [1, hi] */ |
84 | |
|
85 | 0 | if (fabs(min) <= fabs(max)) { |
86 | 0 | if (k == hi) { |
87 | 0 | return min; |
88 | 0 | } else { |
89 | 0 | double k_hi, k_lo; |
90 | 0 | splitint64(k, &k_hi, &k_lo); |
91 | |
|
92 | 0 | return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g; |
93 | 0 | } |
94 | 0 | } else { |
95 | 0 | double k_hi, k_lo; |
96 | 0 | splitint64(k - 1, &k_hi, &k_lo); |
97 | |
|
98 | 0 | return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g; |
99 | 0 | } |
100 | 0 | } |
101 | | |
102 | | PHPAPI double php_random_gammasection_closed_closed(php_random_algo_with_state engine, double min, double max) |
103 | 0 | { |
104 | 0 | double g = gamma_max(min, max); |
105 | 0 | uint64_t hi = ceilint(min, max, g); |
106 | |
|
107 | 0 | if (UNEXPECTED(max < min)) { |
108 | 0 | return NAN; |
109 | 0 | } |
110 | | |
111 | 0 | uint64_t k = php_random_range64(engine, hi); /* [0, hi] */ |
112 | |
|
113 | 0 | if (fabs(min) <= fabs(max)) { |
114 | 0 | if (k == hi) { |
115 | 0 | return min; |
116 | 0 | } else { |
117 | 0 | double k_hi, k_lo; |
118 | 0 | splitint64(k, &k_hi, &k_lo); |
119 | |
|
120 | 0 | return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g; |
121 | 0 | } |
122 | 0 | } else { |
123 | 0 | if (k == hi) { |
124 | 0 | return max; |
125 | 0 | } else { |
126 | 0 | double k_hi, k_lo; |
127 | 0 | splitint64(k, &k_hi, &k_lo); |
128 | |
|
129 | 0 | return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g; |
130 | 0 | } |
131 | 0 | } |
132 | 0 | } |
133 | | |
134 | | PHPAPI double php_random_gammasection_open_closed(php_random_algo_with_state engine, double min, double max) |
135 | 0 | { |
136 | 0 | double g = gamma_max(min, max); |
137 | 0 | uint64_t hi = ceilint(min, max, g); |
138 | |
|
139 | 0 | if (UNEXPECTED(max <= min || hi < 1)) { |
140 | 0 | return NAN; |
141 | 0 | } |
142 | | |
143 | 0 | uint64_t k = php_random_range64(engine, hi - 1); /* [0, hi - 1] */ |
144 | |
|
145 | 0 | if (fabs(min) <= fabs(max)) { |
146 | 0 | double k_hi, k_lo; |
147 | 0 | splitint64(k, &k_hi, &k_lo); |
148 | |
|
149 | 0 | return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g; |
150 | 0 | } else { |
151 | 0 | if (k == (hi - 1)) { |
152 | 0 | return max; |
153 | 0 | } else { |
154 | 0 | double k_hi, k_lo; |
155 | 0 | splitint64(k + 1, &k_hi, &k_lo); |
156 | |
|
157 | 0 | return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g; |
158 | 0 | } |
159 | 0 | } |
160 | 0 | } |
161 | | |
162 | | PHPAPI double php_random_gammasection_open_open(php_random_algo_with_state engine, double min, double max) |
163 | 0 | { |
164 | 0 | double g = gamma_max(min, max); |
165 | 0 | uint64_t hi = ceilint(min, max, g); |
166 | |
|
167 | 0 | if (UNEXPECTED(max <= min || hi < 2)) { |
168 | 0 | return NAN; |
169 | 0 | } |
170 | | |
171 | 0 | uint64_t k = 1 + php_random_range64(engine, hi - 2); /* [1, hi - 1] */ |
172 | |
|
173 | 0 | if (fabs(min) <= fabs(max)) { |
174 | 0 | double k_hi, k_lo; |
175 | 0 | splitint64(k, &k_hi, &k_lo); |
176 | |
|
177 | 0 | return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g; |
178 | 0 | } else { |
179 | 0 | double k_hi, k_lo; |
180 | 0 | splitint64(k, &k_hi, &k_lo); |
181 | |
|
182 | 0 | return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g; |
183 | 0 | } |
184 | 0 | } |