1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
|
/*
+----------------------------------------------------------------------+
| Copyright (c) The PHP Group |
+----------------------------------------------------------------------+
| This source file is subject to version 3.01 of the PHP license, |
| that is bundled with this package in the file LICENSE, and is |
| available through the world-wide-web at the following url: |
| https://www.php.net/license/3_01.txt |
| If you did not receive a copy of the PHP license and are unable to |
| obtain it through the world-wide-web, please send a note to |
| license@php.net so we can mail you a copy immediately. |
+----------------------------------------------------------------------+
| Authors: Tim Düsterhus <timwolla@php.net> |
| |
| 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. |
+----------------------------------------------------------------------+
*/
#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include "php.h"
#include "php_random.h"
#include <math.h>
/* This file implements the γ-section algorithm as published in:
*
* Drawing Random Floating-Point Numbers from an Interval. Frédéric
* Goualard, ACM Trans. Model. Comput. Simul., 32:3, 2022.
* https://doi.org/10.1145/3503512
*/
static double gamma_low(double x)
{
return x - nextafter(x, -DBL_MAX);
}
static double gamma_high(double x)
{
return nextafter(x, DBL_MAX) - x;
}
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;
double e;
if (fabs(a) <= fabs(b)) {
e = -a / g - (s - b / g);
} else {
e = b / g - (s + a / g);
}
double si = ceil(s);
return (s != si) ? (uint64_t)si : (uint64_t)si + (e > 0);
}
PHPAPI double php_random_gammasection_closed_open(php_random_algo_with_state engine, double min, double max)
{
double g = gamma_max(min, max);
uint64_t hi = ceilint(min, max, g);
if (UNEXPECTED(max <= min || hi < 1)) {
return NAN;
}
uint64_t k = 1 + php_random_range64(engine, hi - 1); /* [1, hi] */
if (fabs(min) <= fabs(max)) {
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 {
double k_hi, k_lo;
splitint64(k - 1, &k_hi, &k_lo);
return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g;
}
}
PHPAPI double php_random_gammasection_closed_closed(php_random_algo_with_state engine, double min, double max)
{
double g = gamma_max(min, max);
uint64_t hi = ceilint(min, max, g);
if (UNEXPECTED(max < min)) {
return NAN;
}
uint64_t k = php_random_range64(engine, hi); /* [0, hi] */
if (fabs(min) <= fabs(max)) {
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 {
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;
}
}
}
PHPAPI double php_random_gammasection_open_closed(php_random_algo_with_state engine, double min, double max)
{
double g = gamma_max(min, max);
uint64_t hi = ceilint(min, max, g);
if (UNEXPECTED(max <= min || hi < 1)) {
return NAN;
}
uint64_t k = php_random_range64(engine, hi - 1); /* [0, hi - 1] */
if (fabs(min) <= fabs(max)) {
double k_hi, k_lo;
splitint64(k, &k_hi, &k_lo);
return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g;
} else {
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;
}
}
}
PHPAPI double php_random_gammasection_open_open(php_random_algo_with_state engine, double min, double max)
{
double g = gamma_max(min, max);
uint64_t hi = ceilint(min, max, g);
if (UNEXPECTED(max <= min || hi < 2)) {
return NAN;
}
uint64_t k = 1 + php_random_range64(engine, hi - 2); /* [1, hi - 1] */
if (fabs(min) <= fabs(max)) {
double k_hi, k_lo;
splitint64(k, &k_hi, &k_lo);
return 0x1p+2 * (max * 0x1p-2 - k_hi * g) - k_lo * g;
} else {
double k_hi, k_lo;
splitint64(k, &k_hi, &k_lo);
return 0x1p+2 * (min * 0x1p-2 + k_hi * g) + k_lo * g;
}
}
|