summaryrefslogtreecommitdiff
path: root/src/libc64/qrand.c
blob: dafddc864ffbf3a0ed0b856b3c5ca4baa735c80f (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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
/**
 * @file rand.c
 *
 * This file implements the primary random number generator the game relies on. The generator is pseudo-random and
 * implemented as a Linear Congruential Generator (LCG).
 *
 * A LCG computes random numbers sequentially via the relation
 *      X(n+1) = (a * X(n) + c) mod m
 * where m is the modulus, a is the multiplier and c is the increment.
 *
 * These three parameters (a,c,m) completely specify the LCG and should be chosen such that
 *  - m > 0
 *  - 0 < a < m
 *  - 0 <= c < m
 *
 * The period of the LCG (a, c, m) is the smallest period p such that X(n + p) = X(n), past n=p the sequence will repeat
 * itself in its outputs.
 * A good LCG should have the maximum possible period, which will be equal to m as there are at most m possible values
 * for X. This occurs when (Hull, T.E., & Dobell, A.R. (1962). Random Number Generators. Siam Review, 4, 230-254.):
 *  - m,c are relatively prime, that is the only integer that divides both m and c with no remainder is 1.
 *  - a - 1 is divisible by all prime factors of m.
 *  - a - 1 is divisible by 4 if m is divisible by 4.
 *
 * Ideally m is chosen to be a large power of 2 so that the modulo operation is inexpensive to compute. In this case the
 * prime factors of m = 2^k are just k copies of 2. For k > 1 m is divisible by 4, so a - 1 must be divisible by 4. 2^k
 * and c can easily be made relatively prime by making c an odd number.
 * If we let k=32 to match the size of an integer, the modulo operation is made implicit by the width of the data type
 * and becomes free to compute.
 *
 * The parameter a should be selected such that a-1 is divisible by 4 (and hence divisible by 2) and c should be any odd
 * number. The precise values should fare well against the spectral test, a measure of "how random" a particular LCG is.
 * A pair (a,c) that satisfies these requirements is (1664525, 1013904223), recommended by "Numerical Recipes in C: The
 * Art of Scientific Computing" (p. 284).
 *
 * Therefore, the LCG with parameters (1664525, 1013904223, 2^32) that is implemented in this file has a maximal period
 * of 2^32 and produces high-quality pseudo-random numbers.
 *
 * @note If sampling the LCG for a n-bit number it is important to use the upper n bits instead of the lower n bits of
 * the LCG output. The lower n bits only have a period of 2^n which may significantly worsen the quality of the
 * resulting random numbers compared to the quality of the full 32-bit result.
 *
 * @note Original name: qrand.c
 */
#include "libc64/qrand.h"
#include "z_math.h"

#define RAND_MULTIPLIER 1664525
#define RAND_INCREMENT 1013904223

/**
 * The latest generated random number, used to generate the next number in the sequence.
 *
 * @note Original name: __qrand_idum
 */
static u32 sRandInt = 1;

#if !PLATFORM_N64
/**
 * Space to store a value to be re-interpreted as a float.
 *
 * @note Orignal name: __qrand_itemp
 */
static FloatInt sRandFloat;
#endif

/**
 * Gets the next integer in the sequence of pseudo-random numbers.
 *
 * @note Original name: qrand
 */
u32 Rand_Next(void) {
#if PLATFORM_N64
    u32 next = sRandInt * RAND_MULTIPLIER + RAND_INCREMENT;

    sRandInt = next;
    return next;
#else
    return sRandInt = sRandInt * RAND_MULTIPLIER + RAND_INCREMENT;
#endif
}

/**
 * Seeds the pseudo-random number generator by providing a starting value.
 *
 * @note Original name: sqrand
 */
void Rand_Seed(u32 seed) {
    sRandInt = seed;
}

/**
 * Returns a pseudo-random floating-point number between 0.0f and 1.0f, by generating the next integer and masking it
 * to an IEEE-754 compliant floating-point number between 1.0f and 2.0f, returning the result subtract 1.0f.
 *
 * @note This technique for generating pseudo-random floats is recommended as a particularly fast but potentially
 *       non-portable generator in "Numerical Recipes in C: The Art of Scientic Computing", pp. 284-5.
 *
 * @note Original name: fqrand
 */
f32 Rand_ZeroOne(void) {
#if PLATFORM_N64
    fu v;
    f32 vf;

    // Note this samples the lower 23 bits, effectively reducing the LCG period from 2^32 to 2^23.
    // This was fixed in Gamecube versions and Majora's Mask.
    v.i = (Rand_Next() & 0x007FFFFF) | 0x3F800000;
    vf = v.f - 1.0f;
    return vf;
#else
    sRandInt = sRandInt * RAND_MULTIPLIER + RAND_INCREMENT;
    // Samples the upper 23 bits to avoid effectively reducing the LCG period.
    sRandFloat.i = (sRandInt >> 9) | 0x3F800000;
    return sRandFloat.f - 1.0f;
#endif
}

#if !PLATFORM_N64
/**
 * Returns a pseudo-random floating-point number between -0.5f and 0.5f by the same manner in which Rand_ZeroOne
 * generates its result.
 *
 * @see Rand_ZeroOne
 *
 * @note Original name: fqrand2
 */
f32 Rand_Centered(void) {
    sRandInt = sRandInt * RAND_MULTIPLIER + RAND_INCREMENT;
    sRandFloat.i = (sRandInt >> 9) | 0x3F800000;
    return sRandFloat.f - 1.5f;
}
#endif

//! All functions below are unused variants of the above four, that use a provided random number variable instead of the
//! internal `sRandInt`

/**
 * Seeds a pseudo-random number at rndNum with a provided starting value.
 *
 * @see Rand_Seed
 *
 * @note Original name: sqrand_r
 */
void Rand_Seed_Variable(u32* rndNum, u32 seed) {
    *rndNum = seed;
}

/**
 * Generates the next pseudo-random integer from the provided rndNum.
 *
 * @see Rand_Next
 *
 * @note Original name: qrand_r
 */
u32 Rand_Next_Variable(u32* rndNum) {
    return *rndNum = (*rndNum) * RAND_MULTIPLIER + RAND_INCREMENT;
}

/**
 * Generates the next pseudo-random floating-point number between 0.0f and 1.0f from the provided rndNum.
 *
 * @see Rand_ZeroOne
 *
 * @note Original name: fqrand_r
 */
f32 Rand_ZeroOne_Variable(u32* rndNum) {
#if PLATFORM_N64
    fu v;
    f32 vf;
    u32 next = Rand_Next_Variable(rndNum);

    v.i = (next & 0x007FFFFF) | 0x3F800000;
    vf = v.f - 1.0f;
    return vf;
#else
    u32 next = (*rndNum) * RAND_MULTIPLIER + RAND_INCREMENT;

    sRandFloat.i = ((*rndNum = next) >> 9) | 0x3F800000;
    return sRandFloat.f - 1.0f;
#endif
}

#if !PLATFORM_N64
/**
 * Generates the next pseudo-random floating-point number between -0.5f and 0.5f from the provided rndNum.
 *
 * @see Rand_ZeroOne, Rand_Centered
 *
 * @note Original name: fqrand2_r
 */
f32 Rand_Centered_Variable(u32* rndNum) {
    u32 next = (*rndNum) * RAND_MULTIPLIER + RAND_INCREMENT;

    sRandFloat.i = ((*rndNum = next) >> 9) | 0x3F800000;
    return sRandFloat.f - 1.5f;
}
#endif