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

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 <http://www.math.keio.ac.jp/~matumoto/emt.html>
- * (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 <matumoto@math.keio.ac.jp> with
- *   an appropriate reference to your work.
- *
- * It would be nice to CC: <Cokus@math.washington.edu> 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 <stdlib.h>
 
-#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 */
 }