xref: /llvm-project/libc/AOR_v20.02/math/test/rtest/random.c (revision 0928368f623a0f885894f9c3ef1b740b060c0d9c)
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