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