diff options
author | Jens Arnold <amiconn@rockbox.org> | 2006-01-20 01:21:26 +0000 |
---|---|---|
committer | Jens Arnold <amiconn@rockbox.org> | 2006-01-20 01:21:26 +0000 |
commit | c05cd1676f323f1346099f436aaa0212fd18e178 (patch) | |
tree | 26583afe9c2e529fc3c96c1836dad3e10eba1d89 | |
parent | cd7779c229f89f732869c6b7228833ae859aca88 (diff) | |
download | rockbox-c05cd1676f323f1346099f436aaa0212fd18e178.tar.gz rockbox-c05cd1676f323f1346099f436aaa0212fd18e178.zip |
Better random number generator: Mersenne twiser with improved initialisation. BSD-style license.
git-svn-id: svn://svn.rockbox.org/rockbox/trunk@8396 a1c6a512-1295-4272-9138-f99709370657
-rw-r--r-- | docs/LICENSES | 35 | ||||
-rw-r--r-- | firmware/common/random.c | 234 |
2 files changed, 132 insertions, 137 deletions
diff --git a/docs/LICENSES b/docs/LICENSES index 652a276eed..47ddcb19c2 100644 --- a/docs/LICENSES +++ b/docs/LICENSES | |||
@@ -34,3 +34,38 @@ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY | |||
34 | OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF | 34 | OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF |
35 | SUCH DAMAGE. | 35 | SUCH DAMAGE. |
36 | @(#)gmon.c 5.3 (Berkeley) 5/22/91 | 36 | @(#)gmon.c 5.3 (Berkeley) 5/22/91 |
37 | |||
38 | ************************************************************************* | ||
39 | In: random.c | ||
40 | From: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html | ||
41 | mt19937ar-cok.c | ||
42 | ************************************************************************* | ||
43 | Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, | ||
44 | All rights reserved. | ||
45 | |||
46 | Redistribution and use in source and binary forms, with or without | ||
47 | modification, are permitted provided that the following conditions | ||
48 | are met: | ||
49 | |||
50 | 1. Redistributions of source code must retain the above copyright | ||
51 | notice, this list of conditions and the following disclaimer. | ||
52 | |||
53 | 2. Redistributions in binary form must reproduce the above copyright | ||
54 | notice, this list of conditions and the following disclaimer in the | ||
55 | documentation and/or other materials provided with the distribution. | ||
56 | |||
57 | 3. The names of its contributors may not be used to endorse or promote | ||
58 | products derived from this software without specific prior written | ||
59 | permission. | ||
60 | |||
61 | THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS | ||
62 | "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT | ||
63 | LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR | ||
64 | A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR | ||
65 | CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, | ||
66 | EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, | ||
67 | PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR | ||
68 | PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF | ||
69 | LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING | ||
70 | NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS | ||
71 | SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | ||
diff --git a/firmware/common/random.c b/firmware/common/random.c index d23c29ee3d..f3efe89351 100644 --- a/firmware/common/random.c +++ b/firmware/common/random.c | |||
@@ -1,159 +1,119 @@ | |||
1 | /* | 1 | /* |
2 | * This is the ``Mersenne Twister'' random number generator MT19937, which | 2 | A C-program for MT19937, with initialization improved 2002/2/10. |
3 | * generates pseudorandom integers uniformly distributed in 0..(2^32 - 1) | 3 | Coded by Takuji Nishimura and Makoto Matsumoto. |
4 | * starting from any odd seed in 0..(2^32 - 1). This version is a recode | 4 | This is a faster version by taking Shawn Cokus's optimization, |
5 | * by Shawn Cokus (Cokus@math.washington.edu) on March 8, 1998 of a version by | 5 | Matthe Bellew's simplification. |
6 | * Takuji Nishimura (who had suggestions from Topher Cooper and Marc Rieffel in | ||
7 | * July-August 1997). | ||
8 | * | ||
9 | * Effectiveness of the recoding (on Goedel2.math.washington.edu, a DEC Alpha | ||
10 | * running OSF/1) using GCC -O3 as a compiler: before recoding: 51.6 sec. to | ||
11 | * generate 300 million random numbers; after recoding: 24.0 sec. for the same | ||
12 | * (i.e., 46.5% of original time), so speed is now about 12.5 million random | ||
13 | * number generations per second on this machine. | ||
14 | * | ||
15 | * According to the URL <http://www.math.keio.ac.jp/~matumoto/emt.html> | ||
16 | * (and paraphrasing a bit in places), the Mersenne Twister is ``designed | ||
17 | * with consideration of the flaws of various existing generators,'' has | ||
18 | * a period of 2^19937 - 1, gives a sequence that is 623-dimensionally | ||
19 | * equidistributed, and ``has passed many stringent tests, including the | ||
20 | * die-hard test of G. Marsaglia and the load test of P. Hellekalek and | ||
21 | * S. Wegenkittl.'' It is efficient in memory usage (typically using 2506 | ||
22 | * to 5012 bytes of static data, depending on data type sizes, and the code | ||
23 | * is quite short as well). It generates random numbers in batches of 624 | ||
24 | * at a time, so the caching and pipelining of modern systems is exploited. | ||
25 | * It is also divide- and mod-free. | ||
26 | * | ||
27 | * This library is free software; you can redistribute it and/or modify it | ||
28 | * under the terms of the GNU Library General Public License as published by | ||
29 | * the Free Software Foundation (either version 2 of the License or, at your | ||
30 | * option, any later version). This library is distributed in the hope that | ||
31 | * it will be useful, but WITHOUT ANY WARRANTY, without even the implied | ||
32 | * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See | ||
33 | * the GNU Library General Public License for more details. You should have | ||
34 | * received a copy of the GNU Library General Public License along with this | ||
35 | * library; if not, write to the Free Software Foundation, Inc., 59 Temple | ||
36 | * Place, Suite 330, Boston, MA 02111-1307, USA. | ||
37 | * | ||
38 | * The code as Shawn received it included the following notice: | ||
39 | * | ||
40 | * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. When | ||
41 | * you use this, send an e-mail to <matumoto@math.keio.ac.jp> with | ||
42 | * an appropriate reference to your work. | ||
43 | * | ||
44 | * It would be nice to CC: <Cokus@math.washington.edu> when you write. | ||
45 | * | ||
46 | */ | ||
47 | 6 | ||
7 | Before using, initialize the state by using srand(seed). | ||
8 | |||
9 | Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, | ||
10 | All rights reserved. | ||
11 | |||
12 | Redistribution and use in source and binary forms, with or without | ||
13 | modification, are permitted provided that the following conditions | ||
14 | are met: | ||
15 | |||
16 | 1. Redistributions of source code must retain the above copyright | ||
17 | notice, this list of conditions and the following disclaimer. | ||
18 | |||
19 | 2. Redistributions in binary form must reproduce the above copyright | ||
20 | notice, this list of conditions and the following disclaimer in the | ||
21 | documentation and/or other materials provided with the distribution. | ||
22 | |||
23 | 3. The names of its contributors may not be used to endorse or promote | ||
24 | products derived from this software without specific prior written | ||
25 | permission. | ||
26 | |||
27 | THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS | ||
28 | "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT | ||
29 | LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR | ||
30 | A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR | ||
31 | CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, | ||
32 | EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, | ||
33 | PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR | ||
34 | PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF | ||
35 | LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING | ||
36 | NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS | ||
37 | SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | ||
38 | |||
39 | Any feedback is very welcome. | ||
40 | http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html | ||
41 | email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) | ||
42 | */ | ||
43 | |||
44 | /* | ||
45 | Adapted to Rockbox by Jens Arnold | ||
46 | */ | ||
47 | |||
48 | #include <stdlib.h> | 48 | #include <stdlib.h> |
49 | 49 | ||
50 | #define N (624) /* length of state vector */ | 50 | /* Period parameters */ |
51 | #define M (397) /* a period parameter */ | 51 | #define N 624 |
52 | #define K (0x9908B0DFUL) /* a magic constant */ | 52 | #define M 397 |
53 | #define hiBit(u) ((u) & 0x80000000UL) /* mask all but highest bit of u */ | 53 | #define MATRIX_A 0x9908b0dfUL /* constant vector a */ |
54 | #define loBit(u) ((u) & 0x00000001UL) /* mask all but lowest bit of u */ | 54 | #define UMASK 0x80000000UL /* most significant w-r bits */ |
55 | #define loBits(u) ((u) & 0x7FFFFFFFUL) /* mask the highest bit of u */ | 55 | #define LMASK 0x7fffffffUL /* least significant r bits */ |
56 | #define mixBits(u, v) (hiBit(u)|loBits(v)) /* move highest bit of u to | 56 | #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) |
57 | highest bit of v */ | 57 | #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) |
58 | 58 | ||
59 | static unsigned long state[N+1]; /* state vector + 1 to not violate ANSI C */ | 59 | static unsigned long state[N]; /* the array for the state vector */ |
60 | static unsigned long *next; /* next random value is computed from here */ | 60 | static int left = 0; |
61 | static int left = -1; /* can *next++ this many times before reloading */ | 61 | static unsigned long *next; |
62 | 62 | ||
63 | /* initializes state[N] with a seed */ | ||
63 | void srand(unsigned int seed) | 64 | void srand(unsigned int seed) |
64 | { | 65 | { |
65 | /* | 66 | unsigned long x = seed & 0xffffffffUL; |
66 | * We initialize state[0..(N-1)] via the generator | 67 | unsigned long *s = state; |
67 | * | ||
68 | * x_new = (69069 * x_old) mod 2^32 | ||
69 | * | ||
70 | * from Line 15 of Table 1, p. 106, Sec. 3.3.4 of Knuth's | ||
71 | * _The Art of Computer Programming_, Volume 2, 3rd ed. | ||
72 | * | ||
73 | * Notes (SJC): I do not know what the initial state requirements | ||
74 | * of the Mersenne Twister are, but it seems this seeding generator | ||
75 | * could be better. It achieves the maximum period for its modulus | ||
76 | * (2^30) iff x_initial is odd (p. 20-21, Sec. 3.2.1.2, Knuth); if | ||
77 | * x_initial can be even, you have sequences like 0, 0, 0, ...; | ||
78 | * 2^31, 2^31, 2^31, ...; 2^30, 2^30, 2^30, ...; 2^29, 2^29 + 2^31, | ||
79 | * 2^29, 2^29 + 2^31, ..., etc. so I force seed to be odd below. | ||
80 | * | ||
81 | * Even if x_initial is odd, if x_initial is 1 mod 4 then | ||
82 | * | ||
83 | * the lowest bit of x is always 1, | ||
84 | * the next-to-lowest bit of x is always 0, | ||
85 | * the 2nd-from-lowest bit of x alternates ... 0 1 0 1 0 1 0 1 ... , | ||
86 | * the 3rd-from-lowest bit of x 4-cycles ... 0 1 1 0 0 1 1 0 ... , | ||
87 | * the 4th-from-lowest bit of x has the 8-cycle ... 0 0 0 1 1 1 1 0 ... , | ||
88 | * ... | ||
89 | * | ||
90 | * and if x_initial is 3 mod 4 then | ||
91 | * | ||
92 | * the lowest bit of x is always 1, | ||
93 | * the next-to-lowest bit of x is always 1, | ||
94 | * the 2nd-from-lowest bit of x alternates ... 0 1 0 1 0 1 0 1 ... , | ||
95 | * the 3rd-from-lowest bit of x 4-cycles ... 0 0 1 1 0 0 1 1 ... , | ||
96 | * the 4th-from-lowest bit of x has the 8-cycle ... 0 0 1 1 1 1 0 0 ... , | ||
97 | * ... | ||
98 | * | ||
99 | * The generator's potency (min. s>=0 with (69069-1)^s = 0 mod 2^32) is | ||
100 | * 16, which seems to be alright by p. 25, Sec. 3.2.1.3 of Knuth. It | ||
101 | * also does well in the dimension 2..5 spectral tests, but it could be | ||
102 | * better in dimension 6 (Line 15, Table 1, p. 106, Sec. 3.3.4, Knuth). | ||
103 | * | ||
104 | * Note that the random number user does not see the values generated | ||
105 | * here directly since reloadMT() will always munge them first, so maybe | ||
106 | * none of all of this matters. In fact, the seed values made here could | ||
107 | * even be extra-special desirable if the Mersenne Twister theory says | ||
108 | * so-- that's why the only change I made is to restrict to odd seeds. | ||
109 | */ | ||
110 | |||
111 | unsigned long x = (seed | 1UL) & 0xFFFFFFFFUL, *s = state; | ||
112 | int j; | 68 | int j; |
113 | 69 | ||
114 | for(left=0, *s++=x, j=N; --j; | 70 | for (*s++ = x, j = 1; j < N; j++) { |
115 | *s++ = (x*=69069UL) & 0xFFFFFFFFUL); | 71 | x = (1812433253UL * (x ^ (x >> 30)) + j) & 0xffffffffUL; |
72 | *s++ = x; | ||
73 | /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ | ||
74 | /* In the previous versions, MSBs of the seed affect */ | ||
75 | /* only MSBs of the array state[]. */ | ||
76 | /* 2002/01/09 modified by Makoto Matsumoto */ | ||
77 | } | ||
78 | left = 1; | ||
116 | } | 79 | } |
117 | 80 | ||
118 | static int rand_reload(void) | 81 | static void next_state(void) |
119 | { | 82 | { |
120 | unsigned long *p0=state, *p2=state+2, *pM=state+M, s0, s1; | 83 | unsigned long *p = state; |
121 | int j; | 84 | int j; |
122 | 85 | ||
123 | if(left < -1) | 86 | /* if srand() has not been called, */ |
124 | srand(4357UL); | 87 | /* a default initial seed is used */ |
125 | 88 | if (left < 0) | |
126 | left=N-1, next=state+1; | 89 | srand(5489UL); |
127 | 90 | ||
128 | for(s0=state[0], s1=state[1], j=N-M+1; --j; s0=s1, s1=*p2++) | 91 | left = N; |
129 | *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0UL); | 92 | next = state; |
130 | 93 | ||
131 | for(pM=state, j=M; --j; s0=s1, s1=*p2++) | 94 | for (j = N - M + 1; --j; p++) |
132 | *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0UL); | 95 | *p = p[M] ^ TWIST(p[0], p[1]); |
133 | 96 | ||
134 | s1=state[0], *p0 = *pM ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0UL); | 97 | for (j = M; --j; p++) |
135 | s1 ^= (s1 >> 11); | 98 | *p = p[M-N] ^ TWIST(p[0], p[1]); |
136 | s1 ^= (s1 << 7) & 0x9D2C5680UL; | 99 | |
137 | s1 ^= (s1 << 15) & 0xEFC60000UL; | 100 | *p = p[M-N] ^ TWIST(p[0], state[0]); |
138 | return (long)s1 ^ (s1 >> 18); | ||
139 | } | 101 | } |
140 | 102 | ||
103 | /* generates a random number on [0,RAND_MAX]-interval */ | ||
141 | int rand(void) | 104 | int rand(void) |
142 | { | 105 | { |
143 | unsigned long y; | 106 | unsigned long y; |
144 | 107 | ||
145 | if(--left < 0) { | 108 | if (--left <= 0) |
146 | y = rand_reload(); | 109 | next_state(); |
147 | } | 110 | y = *next++; |
148 | else { | 111 | |
149 | y = *next++; | 112 | /* Tempering */ |
150 | y ^= (y >> 11); | 113 | y ^= (y >> 11); |
151 | y ^= (y << 7) & 0x9D2C5680UL; | 114 | y ^= (y << 7) & 0x9d2c5680UL; |
152 | y ^= (y << 15) & 0xEFC60000UL; | 115 | y ^= (y << 15) & 0xefc60000UL; |
153 | y ^= (y >> 18); | 116 | y ^= (y >> 18); |
154 | } | ||
155 | 117 | ||
156 | return ((unsigned int)y) >> 1; | 118 | return ((unsigned int)y) >> 1; |
157 | /* 31-bit limit by Björn Stenberg*/ | ||
158 | /* 16-bit architectures compatibility by Jean-Philippe Bernardy */ | ||
159 | } | 119 | } |