The GLIBC random number generator
Written Jan 4, 2007.
The GNU C library's random() function provides pseudorandom
numbers via a linear additive feedback method. A description of the
exact algorithm used is hard to find, so I have documented it here.
The random(3) man page states, misleadingly, that "the
random() function uses a nonlinear additive feedback
random number generator". This is not actually true, as the feedback
is in fact linear modulo 2^{32}. The only slight
nonlinearity is introduced during the seeding stage, due to the fact
that the seeding calculation is done modulo 2^{31}  1 and not
modulo 2^{31} or 2^{32}. This means that, although
each generated random number depends linearly on previous numbers in
the sequence, the random numbers as a whole do not depend linearly on
the seed.
Description of the algorithm
For reference, and because it is easier to read than source code, here
is a mathematical description of the exact algorithm used by the GLIBC
pseudorandom number generator. Note, in the following description,
that 2147483647 = 2^{31}  1 and 4294967296 =
2^{32}. All quantities are mathematical integers.
For a given seed s, an initialization vector
r_{0}...r_{33} is calculated according to the
following scheme:
 (1)
r_{0} = s
 (2)
r_{i} = (16807 *
(signed int) r_{i1}) mod 2147483647 (for i = 1...30)
 (3)
r_{i} =
r_{i31} (for i = 31...33)
Note that the multiplication by 16807 is done in a large enough signed
integer type so that there is no overflow before the modulo
operation. Also note that r_{i1} is converted
to a signed 32bit value before the multiplication, but the only time
this value can be negative is in the case of i=1, if s
≥ 2^{31}. The modulo operation is mathematical, i.e., the result is
taken between 0 and 2147483646, even if r_{i1} is
negative.
Then a sequence of pseudorandom numbers r_{34}... is
generated by a linear feedback loop as follows:
 (4)
r_{i} = (r_{i3} +
r_{i31}) mod 4294967296 (for i ≥ 34)
r_{0}...r_{343} are thrown away. The
ith output o_{i} of rand() is
 (5)
o_{i} = r_{i+344} >> 1
Note that this is a 31bit number; the least significant bit of
r_{i+344} is thrown away.
Linearity
Despite the fact that the least significant bit is thrown away, we
still have almostlinearity of the output sequence. As a consequence
of equations (4) and (5), we have
 (6)
o_{i} =
o_{i31} + o_{i3} mod 2^{31}
or
o_{i} =
o_{i31} + o_{i3} + 1 mod
2^{31},
for all i ≥ 31.
Indeed, you can check this for yourself: the first 60 numbers returned
by random() for seed 1 (the default seed), are:
0:  1804289383 
1:  846930886 
2:  1681692777 
3:  1714636915 
4:  1957747793 
5:  424238335 
6:  719885386 
7:  1649760492 
8:  596516649 
9:  1189641421 
10:  1025202362 
11:  1350490027 

12:  783368690 
13:  1102520059 
14:  2044897763 
15:  1967513926 
16:  1365180540 
17:  1540383426 
18:  304089172 
19:  1303455736 
20:  35005211 
21:  521595368 
22:  294702567 
23:  1726956429 

24:  336465782 
25:  861021530 
26:  278722862 
27:  233665123 
28:  2145174067 
29:  468703135 
30:  1101513929 
31:  1801979802 
32:  1315634022 
33:  635723058 
34:  1369133069 
35:  1125898167 

36:  1059961393 
37:  2089018456 
38:  628175011 
39:  1656478042 
40:  1131176229 
41:  1653377373 
42:  859484421 
43:  1914544919 
44:  608413784 
45:  756898537 
46:  1734575198 
47:  1973594324 

48:  149798315 
49:  2038664370 
50:  1129566413 
51:  184803526 
52:  412776091 
53:  1424268980 
54:  1911759956 
55:  749241873 
56:  137806862 
57:  42999170 
58:  982906996 
59:  135497281 

And indeed,
 o_{31} = 1801979802 = o_{0} +
o_{28} = 1804289383 + 2145174067 (mod 2^{31})
 o_{32} = 1315634022 = o_{1} +
o_{29} + 1= 846930886 + 468703135 + 1 (mod 2^{31})
 etc.
Simplified code
The following simple C program generates exactly the same sequence of
pseudorandom numbers as the GLIBC random() function, for any
given seed.
#include <stdio.h>
#define MAX 1000
#define seed 1
main() {
int r[MAX];
int i;
r[0] = seed;
for (i=1; i<31; i++) {
r[i] = (16807LL * r[i1]) % 2147483647;
if (r[i] < 0) {
r[i] += 2147483647;
}
}
for (i=31; i<34; i++) {
r[i] = r[i31];
}
for (i=34; i<344; i++) {
r[i] = r[i31] + r[i3];
}
for (i=344; i<MAX; i++) {
r[i] = r[i31] + r[i3];
printf("%d\n", ((unsigned int)r[i]) >> 1);
}
}
Note: the GLIBC implementation allows other degrees besides 31 in the
linear feedback register, via the initstate() function. An
analogous description applies in these situations.
Back to Homepage:
Peter Selinger /
Department of Mathematics and Statistics /
Dalhousie University
selinger@mathstat.dal.ca
/ PGP key