Coverage Report

Created: 2025-06-13 06:43

/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
}