From c05cd1676f323f1346099f436aaa0212fd18e178 Mon Sep 17 00:00:00 2001 From: Jens Arnold Date: Fri, 20 Jan 2006 01:21:26 +0000 Subject: 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 --- docs/LICENSES | 35 +++++++ 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 OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @(#)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. 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 @@ /* - * This is the ``Mersenne Twister'' random number generator MT19937, which - * generates pseudorandom integers uniformly distributed in 0..(2^32 - 1) - * starting from any odd seed in 0..(2^32 - 1). This version is a recode - * by Shawn Cokus (Cokus@math.washington.edu) on March 8, 1998 of a version by - * Takuji Nishimura (who had suggestions from Topher Cooper and Marc Rieffel in - * July-August 1997). - * - * Effectiveness of the recoding (on Goedel2.math.washington.edu, a DEC Alpha - * running OSF/1) using GCC -O3 as a compiler: before recoding: 51.6 sec. to - * 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 - * number generations per second on this machine. - * - * According to the URL - * (and paraphrasing a bit in places), the Mersenne Twister is ``designed - * with consideration of the flaws of various existing generators,'' has - * a period of 2^19937 - 1, gives a sequence that is 623-dimensionally - * equidistributed, and ``has passed many stringent tests, including the - * die-hard test of G. Marsaglia and the load test of P. Hellekalek and - * S. Wegenkittl.'' It is efficient in memory usage (typically using 2506 - * 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 - * at a time, so the caching and pipelining of modern systems is exploited. - * It is also divide- and mod-free. - * - * This library is free software; you can redistribute it and/or modify it - * under the terms of the GNU Library General Public License as published by - * the Free Software Foundation (either version 2 of the License or, at your - * option, any later version). This library is distributed in the hope that - * it will be useful, but WITHOUT ANY WARRANTY, without even the implied - * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See - * the GNU Library General Public License for more details. You should have - * received a copy of the GNU Library General Public License along with this - * library; if not, write to the Free Software Foundation, Inc., 59 Temple - * Place, Suite 330, Boston, MA 02111-1307, USA. - * - * The code as Shawn received it included the following notice: - * - * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. When - * you use this, send an e-mail to with - * an appropriate reference to your work. - * - * It would be nice to CC: when you write. - * - */ + A C-program for MT19937, with initialization improved 2002/2/10. + Coded by Takuji Nishimura and Makoto Matsumoto. + This is a faster version by taking Shawn Cokus's optimization, + Matthe Bellew's simplification. + Before using, initialize the state by using srand(seed). + + 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. + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) +*/ + +/* + Adapted to Rockbox by Jens Arnold +*/ + #include -#define N (624) /* length of state vector */ -#define M (397) /* a period parameter */ -#define K (0x9908B0DFUL) /* a magic constant */ -#define hiBit(u) ((u) & 0x80000000UL) /* mask all but highest bit of u */ -#define loBit(u) ((u) & 0x00000001UL) /* mask all but lowest bit of u */ -#define loBits(u) ((u) & 0x7FFFFFFFUL) /* mask the highest bit of u */ -#define mixBits(u, v) (hiBit(u)|loBits(v)) /* move highest bit of u to - highest bit of v */ +/* Period parameters */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0dfUL /* constant vector a */ +#define UMASK 0x80000000UL /* most significant w-r bits */ +#define LMASK 0x7fffffffUL /* least significant r bits */ +#define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) +#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 *next; /* next random value is computed from here */ -static int left = -1; /* can *next++ this many times before reloading */ +static unsigned long state[N]; /* the array for the state vector */ +static int left = 0; +static unsigned long *next; +/* initializes state[N] with a seed */ void srand(unsigned int seed) { - /* - * We initialize state[0..(N-1)] via the generator - * - * 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; + unsigned long x = seed & 0xffffffffUL; + unsigned long *s = state; int j; - for(left=0, *s++=x, j=N; --j; - *s++ = (x*=69069UL) & 0xFFFFFFFFUL); + for (*s++ = x, j = 1; j < N; j++) { + 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; - - if(left < -1) - srand(4357UL); - - left=N-1, next=state+1; - - for(s0=state[0], s1=state[1], j=N-M+1; --j; s0=s1, s1=*p2++) - *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0UL); - - for(pM=state, j=M; --j; s0=s1, s1=*p2++) - *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0UL); - - s1=state[0], *p0 = *pM ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0UL); - s1 ^= (s1 >> 11); - s1 ^= (s1 << 7) & 0x9D2C5680UL; - s1 ^= (s1 << 15) & 0xEFC60000UL; - return (long)s1 ^ (s1 >> 18); + + /* if srand() has not been called, */ + /* a default initial seed is used */ + if (left < 0) + srand(5489UL); + + left = N; + next = state; + + for (j = N - M + 1; --j; p++) + *p = p[M] ^ TWIST(p[0], p[1]); + + for (j = M; --j; p++) + *p = p[M-N] ^ TWIST(p[0], p[1]); + + *p = p[M-N] ^ TWIST(p[0], state[0]); } +/* generates a random number on [0,RAND_MAX]-interval */ int rand(void) { unsigned long y; - - if(--left < 0) { - y = rand_reload(); - } - else { - y = *next++; - y ^= (y >> 11); - y ^= (y << 7) & 0x9D2C5680UL; - y ^= (y << 15) & 0xEFC60000UL; - y ^= (y >> 18); - } + + if (--left <= 0) + next_state(); + y = *next++; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); return ((unsigned int)y) >> 1; - /* 31-bit limit by Björn Stenberg*/ - /* 16-bit architectures compatibility by Jean-Philippe Bernardy */ } -- cgit v1.2.3