summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--docs/LICENSES35
-rw-r--r--firmware/common/random.c234
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
34OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 34OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
35SUCH DAMAGE. 35SUCH DAMAGE.
36@(#)gmon.c 5.3 (Berkeley) 5/22/91 36@(#)gmon.c 5.3 (Berkeley) 5/22/91
37
38*************************************************************************
39In: random.c
40From: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
41 mt19937ar-cok.c
42*************************************************************************
43Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
44All rights reserved.
45
46Redistribution and use in source and binary forms, with or without
47modification, are permitted provided that the following conditions
48are 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
61THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
62"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
63LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
64A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
65CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
66EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
67PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
68PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
69LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
70NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
71SOFTWARE, 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
59static unsigned long state[N+1]; /* state vector + 1 to not violate ANSI C */ 59static unsigned long state[N]; /* the array for the state vector */
60static unsigned long *next; /* next random value is computed from here */ 60static int left = 0;
61static int left = -1; /* can *next++ this many times before reloading */ 61static unsigned long *next;
62 62
63/* initializes state[N] with a seed */
63void srand(unsigned int seed) 64void 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
118static int rand_reload(void) 81static 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 */
141int rand(void) 104int 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}