SphinxBase  0.6
src/libsphinxbase/util/genrand.c
00001 /* -*- c-basic-offset: 4; indent-tabs-mode: nil -*- */
00002 /* ====================================================================
00003  * Copyright (c) 1999-2004 Carnegie Mellon University.  All rights
00004  * reserved.
00005  *
00006  * Redistribution and use in source and binary forms, with or without
00007  * modification, are permitted provided that the following conditions
00008  * are met:
00009  *
00010  * 1. Redistributions of source code must retain the above copyright
00011  *    notice, this list of conditions and the following disclaimer. 
00012  *
00013  * 2. Redistributions in binary form must reproduce the above copyright
00014  *    notice, this list of conditions and the following disclaimer in
00015  *    the documentation and/or other materials provided with the
00016  *    distribution.
00017  *
00018  * This work was supported in part by funding from the Defense Advanced 
00019  * Research Projects Agency and the National Science Foundation of the 
00020  * United States of America, and the CMU Sphinx Speech Consortium.
00021  *
00022  * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
00023  * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
00024  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00025  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
00026  * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00027  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
00028  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
00029  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
00030  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
00031  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
00032  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00033  *
00034  * ====================================================================
00035  *
00036  */
00037 /* 
00038    A C-program for MT19937, with initialization improved 2002/1/26.
00039    Coded by Takuji Nishimura and Makoto Matsumoto.
00040 
00041    Before using, initialize the state by using init_genrand(seed)  
00042    or init_by_array(init_key, key_length).
00043 
00044    Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
00045    All rights reserved.                          
00046 
00047    Redistribution and use in source and binary forms, with or without
00048    modification, are permitted provided that the following conditions
00049    are met:
00050 
00051      1. Redistributions of source code must retain the above copyright
00052         notice, this list of conditions and the following disclaimer.
00053 
00054      2. Redistributions in binary form must reproduce the above copyright
00055 `        notice, this list of conditions and the following disclaimer in the
00056         documentation and/or other materials provided with the distribution.
00057 
00058      3. The names of its contributors may not be used to endorse or promote 
00059         products derived from this software without specific prior written 
00060         permission.
00061 
00062    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00063    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00064    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
00065    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
00066    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00067    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00068    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00069    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00070    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00071    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00072    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00073 
00074 
00075    Any feedback is very welcome.
00076    http://www.math.keio.ac.jp/matumoto/emt.html
00077    email: matumoto@math.keio.ac.jp
00078 */
00079 
00080 
00081 /*
00082  * randgen.c   : a portable random generator 
00083  * 
00084  *
00085  * **********************************************
00086  * CMU ARPA Speech Project
00087  *
00088  * Copyright (c) 1999 Carnegie Mellon University.
00089  * ALL RIGHTS RESERVED.
00090  * **********************************************
00091  * 
00092  * HISTORY
00093  * $Log: genrand.c,v $
00094  * Revision 1.2  2005/06/22 03:01:50  arthchan2003
00095  * Added  keyword
00096  *
00097  * Revision 1.3  2005/03/30 01:22:48  archan
00098  * Fixed mistakes in last updates. Add
00099  *
00100  * 
00101  * 18-Nov-04 ARCHAN (archan@cs.cmu.edu) at Carnegie Mellon University
00102  *              First incorporated from the Mersenne Twister Random
00103  *              Number Generator package. It was chosen because it is
00104  *              in BSD-license and its performance is quite
00105  *              reasonable. Of course if you look at the inventors's
00106  *              page.  This random generator can actually gives
00107  *              19937-bits period.  This is already far from we need. 
00108  *              This will possibly good enough for the next 10 years. 
00109  *
00110  *              I also downgrade the code a little bit to avoid Sphinx's
00111  *              developers misused it. 
00112  */
00113 
00114 
00115 
00116 #include <stdio.h>
00117 
00118 #include "sphinxbase/genrand.h"
00119 
00120 /* Period parameters */
00121 #define N 624
00122 #define M 397
00123 #define MATRIX_A 0x9908b0dfUL   /* constant vector a */
00124 #define UPPER_MASK 0x80000000UL /* most significant w-r bits */
00125 #define LOWER_MASK 0x7fffffffUL /* least significant r bits */
00126 
00127 void init_genrand(unsigned long s);
00128 
00129 void
00130 genrand_seed(unsigned long s)
00131 {
00132     init_genrand(s);
00133 }
00134 
00135 
00136 static unsigned long mt[N];     /* the array for the state vector  */
00137 static int mti = N + 1;         /* mti==N+1 means mt[N] is not initialized */
00138 
00139 /* initializes mt[N] with a seed */
00140 void
00141 init_genrand(unsigned long s)
00142 {
00143     mt[0] = s & 0xffffffffUL;
00144     for (mti = 1; mti < N; mti++) {
00145         mt[mti] =
00146             (1812433253UL * (mt[mti - 1] ^ (mt[mti - 1] >> 30)) + mti);
00147         /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
00148         /* In the previous versions, MSBs of the seed affect   */
00149         /* only MSBs of the array mt[].                        */
00150         /* 2002/01/09 modified by Makoto Matsumoto             */
00151         mt[mti] &= 0xffffffffUL;
00152         /* for >32 bit machines */
00153     }
00154 }
00155 
00156 /* generates a random number on [0,0xffffffff]-interval */
00157 unsigned long
00158 genrand_int32(void)
00159 {
00160     unsigned long y;
00161     static unsigned long mag01[2] = { 0x0UL, MATRIX_A };
00162     /* mag01[x] = x * MATRIX_A  for x=0,1 */
00163 
00164     if (mti >= N) {             /* generate N words at one time */
00165         int kk;
00166 
00167         if (mti == N + 1)       /* if init_genrand() has not been called, */
00168             init_genrand(5489UL);       /* a default initial seed is used */
00169 
00170         for (kk = 0; kk < N - M; kk++) {
00171             y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
00172             mt[kk] = mt[kk + M] ^ (y >> 1) ^ mag01[y & 0x1UL];
00173         }
00174         for (; kk < N - 1; kk++) {
00175             y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
00176             mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
00177         }
00178         y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
00179         mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ mag01[y & 0x1UL];
00180 
00181         mti = 0;
00182     }
00183 
00184     y = mt[mti++];
00185 
00186     /* Tempering */
00187     y ^= (y >> 11);
00188     y ^= (y << 7) & 0x9d2c5680UL;
00189     y ^= (y << 15) & 0xefc60000UL;
00190     y ^= (y >> 18);
00191 
00192     return y;
00193 }
00194 
00195 /* generates a random number on [0,0x7fffffff]-interval */
00196 long
00197 genrand_int31(void)
00198 {
00199     return (long) (genrand_int32() >> 1);
00200 }
00201 
00202 /* generates a random number on [0,1]-real-interval */
00203 double
00204 genrand_real1(void)
00205 {
00206     return genrand_int32() * (1.0 / 4294967295.0);
00207     /* divided by 2^32-1 */
00208 }
00209 
00210 /* generates a random number on [0,1)-real-interval */
00211 double
00212 genrand_real2(void)
00213 {
00214     return genrand_int32() * (1.0 / 4294967296.0);
00215     /* divided by 2^32 */
00216 }
00217 
00218 /* generates a random number on (0,1)-real-interval */
00219 double
00220 genrand_real3(void)
00221 {
00222     return (((double) genrand_int32()) + 0.5) * (1.0 / 4294967296.0);
00223     /* divided by 2^32 */
00224 }
00225 
00226 /* generates a random number on [0,1) with 53-bit resolution*/
00227 double
00228 genrand_res53(void)
00229 {
00230     unsigned long a = genrand_int32() >> 5, b = genrand_int32() >> 6;
00231     return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
00232 }
00233 
00234 /* These real versions are due to Isaku Wada, 2002/01/09 added */