Coverage Report

Created: 2026-06-02 06:39

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