devs/foxbox
1
0
Fork 0
forked from len0rd/rockbox
Code Pull requests Activity

Better random number generator: Mersenne twiser with improved initialisation. BSD-style license.

Browse source 
git-svn-id: svn://svn.rockbox.org/rockbox/trunk@8396 a1c6a512-1295-4272-9138-f99709370657
This commit is contained in:
Jens Arnold 2006-01-20 01:21:26 +00:00
parent cd7779c229
commit c05cd1676f
2 changed files with 132 additions and 137 deletions
Download patch file Download diff file Expand all files Collapse all files

35
docs/LICENSES
View file

@ -34,3 +34,38 @@ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH DAMAGE. SUCH DAMAGE.
@(#)gmon.c 5.3 (Berkeley) 5/22/91 @(#)gmon.c 5.3 (Berkeley) 5/22/91
*************************************************************************
In: random.c
From: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
mt19937ar-cok.c
*************************************************************************
Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. The names of its contributors may not be used to endorse or promote
products derived from this software without specific prior written
permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

222
firmware/common/random.c
View file

@ -1,159 +1,119 @@
/* /*
* This is the ``Mersenne Twister'' random number generator MT19937, which A C-program for MT19937, with initialization improved 2002/2/10.
* generates pseudorandom integers uniformly distributed in 0..(2^32 - 1) Coded by Takuji Nishimura and Makoto Matsumoto.
* starting from any odd seed in 0..(2^32 - 1). This version is a recode This is a faster version by taking Shawn Cokus's optimization,
* by Shawn Cokus (Cokus@math.washington.edu) on March 8, 1998 of a version by Matthe Bellew's simplification.
* Takuji Nishimura (who had suggestions from Topher Cooper and Marc Rieffel in
* July-August 1997). Before using, initialize the state by using srand(seed).
*
* Effectiveness of the recoding (on Goedel2.math.washington.edu, a DEC Alpha Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
* running OSF/1) using GCC -O3 as a compiler: before recoding: 51.6 sec. to All rights reserved.
* generate 300 million random numbers; after recoding: 24.0 sec. for the same
* (i.e., 46.5% of original time), so speed is now about 12.5 million random Redistribution and use in source and binary forms, with or without
* number generations per second on this machine. modification, are permitted provided that the following conditions
* are met:
* According to the URL <http://www.math.keio.ac.jp/~matumoto/emt.html>
* (and paraphrasing a bit in places), the Mersenne Twister is ``designed 1. Redistributions of source code must retain the above copyright
* with consideration of the flaws of various existing generators,'' has notice, this list of conditions and the following disclaimer.
* a period of 2^19937 - 1, gives a sequence that is 623-dimensionally
* equidistributed, and ``has passed many stringent tests, including the 2. Redistributions in binary form must reproduce the above copyright
* die-hard test of G. Marsaglia and the load test of P. Hellekalek and notice, this list of conditions and the following disclaimer in the
* S. Wegenkittl.'' It is efficient in memory usage (typically using 2506 documentation and/or other materials provided with the distribution.
* to 5012 bytes of static data, depending on data type sizes, and the code
* is quite short as well). It generates random numbers in batches of 624 3. The names of its contributors may not be used to endorse or promote
* at a time, so the caching and pipelining of modern systems is exploited. products derived from this software without specific prior written
* It is also divide- and mod-free. permission.
*
* This library is free software; you can redistribute it and/or modify it THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
* under the terms of the GNU Library General Public License as published by "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
* the Free Software Foundation (either version 2 of the License or, at your LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
* option, any later version). This library is distributed in the hope that A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
* it will be useful, but WITHOUT ANY WARRANTY, without even the implied CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
* the GNU Library General Public License for more details. You should have PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
* received a copy of the GNU Library General Public License along with this PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
* library; if not, write to the Free Software Foundation, Inc., 59 Temple LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
* Place, Suite 330, Boston, MA 02111-1307, USA. NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
* The code as Shawn received it included the following notice:
* Any feedback is very welcome.
* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. When http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
* you use this, send an e-mail to <matumoto@math.keio.ac.jp> with email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
* an appropriate reference to your work. */
*
* It would be nice to CC: <Cokus@math.washington.edu> when you write. /*
* Adapted to Rockbox by Jens Arnold
*/ */
#include <stdlib.h> #include <stdlib.h>
#define N (624) /* length of state vector */ /* Period parameters */
#define M (397) /* a period parameter */ #define N 624
#define K (0x9908B0DFUL) /* a magic constant */ #define M 397
#define hiBit(u) ((u) & 0x80000000UL) /* mask all but highest bit of u */ #define MATRIX_A 0x9908b0dfUL /* constant vector a */
#define loBit(u) ((u) & 0x00000001UL) /* mask all but lowest bit of u */ #define UMASK 0x80000000UL /* most significant w-r bits */
#define loBits(u) ((u) & 0x7FFFFFFFUL) /* mask the highest bit of u */ #define LMASK 0x7fffffffUL /* least significant r bits */
#define mixBits(u, v) (hiBit(u)|loBits(v)) /* move highest bit of u to #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
highest bit of v */ #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
static unsigned long state[N+1]; /* state vector + 1 to not violate ANSI C */ static unsigned long state[N]; /* the array for the state vector */
static unsigned long *next; /* next random value is computed from here */ static int left = 0;
static int left = -1; /* can *next++ this many times before reloading */ static unsigned long *next;
/* initializes state[N] with a seed */
void srand(unsigned int seed) void srand(unsigned int seed)
{ {
/* unsigned long x = seed & 0xffffffffUL;
* We initialize state[0..(N-1)] via the generator unsigned long *s = state;
*
* x_new = (69069 * x_old) mod 2^32
*
* from Line 15 of Table 1, p. 106, Sec. 3.3.4 of Knuth's
* _The Art of Computer Programming_, Volume 2, 3rd ed.
*
* Notes (SJC): I do not know what the initial state requirements
* of the Mersenne Twister are, but it seems this seeding generator
* could be better. It achieves the maximum period for its modulus
* (2^30) iff x_initial is odd (p. 20-21, Sec. 3.2.1.2, Knuth); if
* x_initial can be even, you have sequences like 0, 0, 0, ...;
* 2^31, 2^31, 2^31, ...; 2^30, 2^30, 2^30, ...; 2^29, 2^29 + 2^31,
* 2^29, 2^29 + 2^31, ..., etc. so I force seed to be odd below.
*
* Even if x_initial is odd, if x_initial is 1 mod 4 then
*
* the lowest bit of x is always 1,
* the next-to-lowest bit of x is always 0,
* the 2nd-from-lowest bit of x alternates ... 0 1 0 1 0 1 0 1 ... ,
* the 3rd-from-lowest bit of x 4-cycles ... 0 1 1 0 0 1 1 0 ... ,
* the 4th-from-lowest bit of x has the 8-cycle ... 0 0 0 1 1 1 1 0 ... ,
* ...
*
* and if x_initial is 3 mod 4 then
*
* the lowest bit of x is always 1,
* the next-to-lowest bit of x is always 1,
* the 2nd-from-lowest bit of x alternates ... 0 1 0 1 0 1 0 1 ... ,
* the 3rd-from-lowest bit of x 4-cycles ... 0 0 1 1 0 0 1 1 ... ,
* the 4th-from-lowest bit of x has the 8-cycle ... 0 0 1 1 1 1 0 0 ... ,
* ...
*
* The generator's potency (min. s>=0 with (69069-1)^s = 0 mod 2^32) is
* 16, which seems to be alright by p. 25, Sec. 3.2.1.3 of Knuth. It
* also does well in the dimension 2..5 spectral tests, but it could be
* better in dimension 6 (Line 15, Table 1, p. 106, Sec. 3.3.4, Knuth).
*
* Note that the random number user does not see the values generated
* here directly since reloadMT() will always munge them first, so maybe
* none of all of this matters. In fact, the seed values made here could
* even be extra-special desirable if the Mersenne Twister theory says
* so-- that's why the only change I made is to restrict to odd seeds.
*/
unsigned long x = (seed | 1UL) & 0xFFFFFFFFUL, *s = state;
int j; int j;
for(left=0, *s++=x, j=N; --j; for (*s++ = x, j = 1; j < N; j++) {
*s++ = (x*=69069UL) & 0xFFFFFFFFUL); x = (1812433253UL * (x ^ (x >> 30)) + j) & 0xffffffffUL;
*s++ = x;
/* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
/* In the previous versions, MSBs of the seed affect */
/* only MSBs of the array state[]. */
/* 2002/01/09 modified by Makoto Matsumoto */
}
left = 1;
} }
static int rand_reload(void) static void next_state(void)
{ {
unsigned long *p0=state, *p2=state+2, *pM=state+M, s0, s1; unsigned long *p = state;
int j; int j;
if(left < -1) /* if srand() has not been called, */
srand(4357UL); /* a default initial seed is used */
if (left < 0)
srand(5489UL);
left=N-1, next=state+1; left = N;
next = state;
for(s0=state[0], s1=state[1], j=N-M+1; --j; s0=s1, s1=*p2++) for (j = N - M + 1; --j; p++)
*p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0UL); *p = p[M] ^ TWIST(p[0], p[1]);
for(pM=state, j=M; --j; s0=s1, s1=*p2++) for (j = M; --j; p++)
*p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0UL); *p = p[M-N] ^ TWIST(p[0], p[1]);
s1=state[0], *p0 = *pM ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0UL); *p = p[M-N] ^ TWIST(p[0], state[0]);
s1 ^= (s1 >> 11);
s1 ^= (s1 << 7) & 0x9D2C5680UL;
s1 ^= (s1 << 15) & 0xEFC60000UL;
return (long)s1 ^ (s1 >> 18);
} }
/* generates a random number on [0,RAND_MAX]-interval */
int rand(void) int rand(void)
{ {
unsigned long y; unsigned long y;
if(--left < 0) { if (--left <= 0)
y = rand_reload(); next_state();
} y = *next++;
else {
y = *next++; /* Tempering */
y ^= (y >> 11); y ^= (y >> 11);
y ^= (y << 7) & 0x9D2C5680UL; y ^= (y << 7) & 0x9d2c5680UL;
y ^= (y << 15) & 0xEFC60000UL; y ^= (y << 15) & 0xefc60000UL;
y ^= (y >> 18); y ^= (y >> 18);
}
return ((unsigned int)y) >> 1; return ((unsigned int)y) >> 1;
/* 31-bit limit by Björn Stenberg*/
/* 16-bit architectures compatibility by Jean-Philippe Bernardy */
} }