aboutsummaryrefslogtreecommitdiff
path: root/math/test/rtest/random.c
blob: 56123966b8c48f8acbeb1501d1e56d5b1d5e3e2c (plain)
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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
/*
 * random.c - random number generator for producing mathlib test cases
 *
 * Copyright (c) 1998-2019, Arm Limited.
 * SPDX-License-Identifier: MIT
 */

#include "types.h"
#include "random.h"

static uint32 seedbuf[55];
static int seedptr;

void seed_random(uint32 seed) {
    int i;

    seedptr = 0;
    for (i = 0; i < 55; i++) {
        seed = seed % 44488 * 48271 - seed / 44488 * 3399;
        seedbuf[i] = seed - 1;
    }
}

uint32 base_random(void) {
    seedptr %= 55;
    seedbuf[seedptr] += seedbuf[(seedptr+31)%55];
    return seedbuf[seedptr++];
}

uint32 random32(void) {
    uint32 a, b, b1, b2;
    a = base_random();
    b = base_random();
    for (b1 = 0x80000000, b2 = 1; b1 > b2; b1 >>= 1, b2 <<= 1) {
        uint32 b3 = b1 | b2;
        if ((b & b3) != 0 && (b & b3) != b3)
            b ^= b3;
    }
    return a ^ b;
}

/*
 * random_upto: generate a uniformly randomised number in the range
 * 0,...,limit-1. (Precondition: limit > 0.)
 *
 * random_upto_biased: generate a number in the same range, but with
 * the probability skewed towards the high end by means of taking the
 * maximum of 8*bias+1 samples from the uniform distribution on the
 * same range. (I don't know why bias is given in that curious way -
 * historical reasons, I expect.)
 *
 * For speed, I separate the implementation of random_upto into the
 * two stages of (a) generate a bitmask which reduces a 32-bit random
 * number to within a factor of two of the right range, (b) repeatedly
 * generate numbers in that range until one is small enough. Splitting
 * it up like that means that random_upto_biased can do (a) only once
 * even when it does (b) lots of times.
 */

static uint32 random_upto_makemask(uint32 limit) {
    uint32 mask = 0xFFFFFFFF;
    int i;
    for (i = 16; i > 0; i >>= 1)
        if ((limit & (mask >> i)) == limit)
            mask >>= i;
    return mask;
}

static uint32 random_upto_internal(uint32 limit, uint32 mask) {
    uint32 ret;
    do {
        ret = random32() & mask;
    } while (ret > limit);
    return ret;
}

uint32 random_upto(uint32 limit) {
    uint32 mask = random_upto_makemask(limit);
    return random_upto_internal(limit, mask);
}

uint32 random_upto_biased(uint32 limit, int bias) {
    uint32 mask = random_upto_makemask(limit);

    uint32 ret = random_upto_internal(limit, mask);
    while (bias--) {
        uint32 tmp;
        tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
        tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
        tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
        tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
        tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
        tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
        tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
        tmp = random_upto_internal(limit, mask); if (tmp < ret) ret = tmp;
    }

    return ret;
}