13
13
| Authors: Tim Düsterhus <timwolla@php.net> |
14
14
| |
15
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. |
16
19
+----------------------------------------------------------------------+
17
20
*/
18
21
@@ -46,6 +49,12 @@ static double gamma_max(double x, double y)
46
49
return (fabs (x ) > fabs (y )) ? gamma_high (x ) : gamma_low (y );
47
50
}
48
51
52
+ static void splitint64 (uint64_t v , double * vhi , double * vlo )
53
+ {
54
+ * vhi = v >> 2 ;
55
+ * vlo = v & UINT64_C (0x3 );
56
+ }
57
+
49
58
static uint64_t ceilint (double a , double b , double g )
50
59
{
51
60
double s = b / g - a / g ;
@@ -74,9 +83,19 @@ PHPAPI double php_random_gammasection_closed_open(const php_random_algo *algo, p
74
83
uint64_t k = 1 + php_random_range64 (algo , status , hi - 1 ); /* [1, hi] */
75
84
76
85
if (fabs (min ) <= fabs (max )) {
77
- return k == hi ? min : max - k * g ;
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
+ }
78
94
} else {
79
- return min + (k - 1 ) * g ;
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 ;
80
99
}
81
100
}
82
101
@@ -92,9 +111,23 @@ PHPAPI double php_random_gammasection_closed_closed(const php_random_algo *algo,
92
111
uint64_t k = php_random_range64 (algo , status , hi ); /* [0, hi] */
93
112
94
113
if (fabs (min ) <= fabs (max )) {
95
- return k == hi ? min : max - k * g ;
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
+ }
96
122
} else {
97
- return k == hi ? max : min + k * g ;
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
+ }
98
131
}
99
132
}
100
133
@@ -110,9 +143,19 @@ PHPAPI double php_random_gammasection_open_closed(const php_random_algo *algo, p
110
143
uint64_t k = php_random_range64 (algo , status , hi - 1 ); /* [0, hi - 1] */
111
144
112
145
if (fabs (min ) <= fabs (max )) {
113
- return max - k * g ;
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 ;
114
150
} else {
115
- return k == (hi - 1 ) ? max : min + (k + 1 ) * g ;
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
+ }
116
159
}
117
160
}
118
161
@@ -128,8 +171,14 @@ PHPAPI double php_random_gammasection_open_open(const php_random_algo *algo, php
128
171
uint64_t k = 1 + php_random_range64 (algo , status , hi - 2 ); /* [1, hi - 1] */
129
172
130
173
if (fabs (min ) <= fabs (max )) {
131
- return max - k * g ;
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 ;
132
178
} else {
133
- return min + k * g ;
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 ;
134
183
}
135
184
}
0 commit comments