/src/php-src/ext/random/gammasection.c

Source (jump to first uncovered line)
/*
   +----------------------------------------------------------------------+
   | 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;
  }
}

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	}

Coverage Report

Created: 2025-06-13 06:43