1*0928368fSKristof Beyls /*
2*0928368fSKristof Beyls * random.c - random number generator for producing mathlib test cases
3*0928368fSKristof Beyls *
4*0928368fSKristof Beyls * Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
5*0928368fSKristof Beyls * See https://llvm.org/LICENSE.txt for license information.
6*0928368fSKristof Beyls * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
7*0928368fSKristof Beyls */
8*0928368fSKristof Beyls
9*0928368fSKristof Beyls #include "types.h"
10*0928368fSKristof Beyls #include "random.h"
11*0928368fSKristof Beyls
12*0928368fSKristof Beyls static uint32 seedbuf[55];
13*0928368fSKristof Beyls static int seedptr;
14*0928368fSKristof Beyls
seed_random(uint32 seed)15*0928368fSKristof Beyls void seed_random(uint32 seed) {
16*0928368fSKristof Beyls int i;
17*0928368fSKristof Beyls
18*0928368fSKristof Beyls seedptr = 0;
19*0928368fSKristof Beyls for (i = 0; i < 55; i++) {
20*0928368fSKristof Beyls seed = seed % 44488 * 48271 - seed / 44488 * 3399;
21*0928368fSKristof Beyls seedbuf[i] = seed - 1;
22*0928368fSKristof Beyls }
23*0928368fSKristof Beyls }
24*0928368fSKristof Beyls
base_random(void)25*0928368fSKristof Beyls uint32 base_random(void) {
26*0928368fSKristof Beyls seedptr %= 55;
27*0928368fSKristof Beyls seedbuf[seedptr] += seedbuf[(seedptr+31)%55];
28*0928368fSKristof Beyls return seedbuf[seedptr++];
29*0928368fSKristof Beyls }
30*0928368fSKristof Beyls
random32(void)31*0928368fSKristof Beyls uint32 random32(void) {
32*0928368fSKristof Beyls uint32 a, b, b1, b2;
33*0928368fSKristof Beyls a = base_random();
34*0928368fSKristof Beyls b = base_random();
35*0928368fSKristof Beyls for (b1 = 0x80000000, b2 = 1; b1 > b2; b1 >>= 1, b2 <<= 1) {
36*0928368fSKristof Beyls uint32 b3 = b1 | b2;
37*0928368fSKristof Beyls if ((b & b3) != 0 && (b & b3) != b3)
38*0928368fSKristof Beyls b ^= b3;
39*0928368fSKristof Beyls }
40*0928368fSKristof Beyls return a ^ b;
41*0928368fSKristof Beyls }
42*0928368fSKristof Beyls
43*0928368fSKristof Beyls /*
44*0928368fSKristof Beyls * random_upto: generate a uniformly randomised number in the range
45*0928368fSKristof Beyls * 0,...,limit-1. (Precondition: limit > 0.)
46*0928368fSKristof Beyls *
47*0928368fSKristof Beyls * random_upto_biased: generate a number in the same range, but with
48*0928368fSKristof Beyls * the probability skewed towards the high end by means of taking the
49*0928368fSKristof Beyls * maximum of 8*bias+1 samples from the uniform distribution on the
50*0928368fSKristof Beyls * same range. (I don't know why bias is given in that curious way -
51*0928368fSKristof Beyls * historical reasons, I expect.)
52*0928368fSKristof Beyls *
53*0928368fSKristof Beyls * For speed, I separate the implementation of random_upto into the
54*0928368fSKristof Beyls * two stages of (a) generate a bitmask which reduces a 32-bit random
55*0928368fSKristof Beyls * number to within a factor of two of the right range, (b) repeatedly
56*0928368fSKristof Beyls * generate numbers in that range until one is small enough. Splitting
57*0928368fSKristof Beyls * it up like that means that random_upto_biased can do (a) only once
58*0928368fSKristof Beyls * even when it does (b) lots of times.
59*0928368fSKristof Beyls */
60*0928368fSKristof Beyls
random_upto_makemask(uint32 limit)61*0928368fSKristof Beyls static uint32 random_upto_makemask(uint32 limit) {
62*0928368fSKristof Beyls uint32 mask = 0xFFFFFFFF;
63*0928368fSKristof Beyls int i;
64*0928368fSKristof Beyls for (i = 16; i > 0; i >>= 1)
65*0928368fSKristof Beyls if ((limit & (mask >> i)) == limit)
66*0928368fSKristof Beyls mask >>= i;
67*0928368fSKristof Beyls return mask;
68*0928368fSKristof Beyls }
69*0928368fSKristof Beyls
random_upto_internal(uint32 limit,uint32 mask)70*0928368fSKristof Beyls static uint32 random_upto_internal(uint32 limit, uint32 mask) {
71*0928368fSKristof Beyls uint32 ret;
72*0928368fSKristof Beyls do {
73*0928368fSKristof Beyls ret = random32() & mask;
74*0928368fSKristof Beyls } while (ret > limit);
75*0928368fSKristof Beyls return ret;
76*0928368fSKristof Beyls }
77*0928368fSKristof Beyls
random_upto(uint32 limit)78*0928368fSKristof Beyls uint32 random_upto(uint32 limit) {
79*0928368fSKristof Beyls uint32 mask = random_upto_makemask(limit);
80*0928368fSKristof Beyls return random_upto_internal(limit, mask);
81*0928368fSKristof Beyls }
82*0928368fSKristof Beyls
random_upto_biased(uint32 limit,int bias)83*0928368fSKristof Beyls uint32 random_upto_biased(uint32 limit, int bias) {
84*0928368fSKristof Beyls uint32 mask = random_upto_makemask(limit);
85*0928368fSKristof Beyls
86*0928368fSKristof Beyls uint32 ret = random_upto_internal(limit, mask);
87*0928368fSKristof Beyls while (bias--) {
88*0928368fSKristof Beyls uint32 tmp;
89*0928368fSKristof Beyls tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
90*0928368fSKristof Beyls tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
91*0928368fSKristof Beyls tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
92*0928368fSKristof Beyls tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
93*0928368fSKristof Beyls tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
94*0928368fSKristof Beyls tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
95*0928368fSKristof Beyls tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
96*0928368fSKristof Beyls tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
97*0928368fSKristof Beyls }
98*0928368fSKristof Beyls
99*0928368fSKristof Beyls return ret;
100*0928368fSKristof Beyls }
101