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
|
/* See LICENSE file for copyright and license details. */
#include "common.h"
#ifndef TEST
#include <math.h>
long double
libsimple_random_float(uintmax_t (*rng)(size_t bits, void *user), void *user, long double min, long double postmax) /* TODO add test */
{
long double range, offset, mantissa, value;
uintmax_t r;
int exp;
size_t bits;
if (min < postmax) {
range = postmax - min;
offset = min;
} else if (min > postmax) {
range = min - postmax;
offset = postmax;
} else {
return min;
}
mantissa = frexpl(range, &exp);
do {
bits = sizeof(uintmax_t) * CHAR_BIT;
r = (*rng)(bits, user);
while (bits && (r & 1)) {
exp -= 1;
bits -= 1;
r >>= 1;
}
} while (!bits);
bits = sizeof(uintmax_t) * CHAR_BIT;
r = (*rng)(bits, user);
value = (long double)r;
value = fmal(value, 0.5, 0.5);
value = ldexpl(value, -(int)bits);
value = ldexpl(value, exp);
value *= mantissa;
return value + offset;
}
#else
#include "test.h"
int
main(void)
{
return 0;
}
#endif
|
