randmt.c

Go to the documentation of this file.
00001 /*                          R A N D . C
00002  * BRL-CAD
00003  *
00004  * A C-program for MT19937: Real number version
00005  *   genrand() generates one pseudorandom real number (double)
00006  * which is uniformly distributed on [0, 1]-interval, for each
00007  * call. sgenrand(seed) set initial values to the working area
00008  * of 624 words. Before genrand(), sgenrand(seed) must be
00009  * called once. (seed is any 32-bit integer except for 0).
00010  * Integer generator is obtained by modifying two lines.
00011  *   Coded by Takuji Nishimura, considering the suggestions by
00012  * Topher Cooper and Marc Rieffel in July-Aug. 1997.
00013  *
00014  * This library is free software; you can redistribute it and/or
00015  * modify it under the terms of the GNU Lesser General Public
00016  * License as published by the Free Software Foundation; either
00017  * version 2 of the License, or (at your option) any later version.
00018  *
00019  * This library is distributed in the hope that it will be useful,
00020  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00021  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
00022  * See the GNU Lesser General Public License for more details.
00023  * You should have received a copy of the GNU Library General
00024  * Public License along with this library; if not, write to the
00025  * Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00026  * 02111-1307  USA
00027  *
00028  * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
00029  * Any feedback is very welcome. For any question, comments,
00030  * see http://www.math.keio.ac.jp/matumoto/emt.html or email
00031  * matumoto@math.keio.ac.jp
00032  */
00033 
00034 #include "bu.h"
00035 #include "bn.h"
00036 
00037 /* Period parameters */
00038 #define N 624
00039 #define M 397
00040 #define MATRIX_A 0x9908b0df   /* constant vector a */
00041 #define UPPER_MASK 0x80000000 /* most significant w-r bits */
00042 #define LOWER_MASK 0x7fffffff /* least significant r bits */
00043 
00044 
00045 /* Tempering parameters */
00046 #define TEMPERING_MASK_B 0x9d2c5680
00047 #define TEMPERING_MASK_C 0xefc60000
00048 #define TEMPERING_SHIFT_U(y)  (y >> 11)
00049 #define TEMPERING_SHIFT_S(y)  (y << 7)
00050 #define TEMPERING_SHIFT_T(y)  (y << 15)
00051 #define TEMPERING_SHIFT_L(y)  (y >> 18)
00052 
00053 #define MERSENNE_MAGIC 0x4D54524E
00054 
00055 static struct _internal_state_s {
00056     uint32_t magic;
00057     int mti;            /* state index */
00058     uint32_t mt[N];     /* state vector */
00059 } global_state_static = { MERSENNE_MAGIC, N+1, {0} };
00060 
00061 static struct _internal_state_s *global_state = &global_state_static;
00062 
00063 void *
00064 bn_randmt_state_create()
00065 {
00066     struct _internal_state_s *is;
00067 
00068     is = bu_malloc(sizeof(struct _internal_state_s), "Mersenne Twister state");
00069     if(is == NULL)
00070         return NULL;
00071     is->magic = MERSENNE_MAGIC;
00072     is->mti = N+1;
00073     return (void *)is;
00074 }
00075 
00076 /* initializing the array with a NONZERO seed */
00077 void
00078 bn_randmt_state_seed(struct _internal_state_s *is, uint32_t seed)
00079 {
00080     /*
00081      * setting initial seeds to mt[N] using
00082      * the generator Line 25 of Table 1 in
00083      * [KNUTH 1981, The Art of Computer Programming
00084      *    Vol. 2 (2nd Ed.), pp102]
00085      */
00086     is->mt[0] = seed & 0xffffffff;
00087     for (is->mti = 1; is->mti<N; is->mti++)
00088         is->mt[is->mti] = (69069 * is->mt[is->mti-1]) & 0xffffffff;
00089 }
00090 
00091 double
00092 bn_randmt_state(struct _internal_state_s *is)
00093 {
00094     uint32_t y;
00095     static uint32_t mag01[2]={0x0, MATRIX_A};
00096 
00097     /* generate N words at one time */
00098     if (is->mti >= N) {
00099         int kk;
00100 
00101         /* if sgenrand() has not been called, a default initial seed is used */
00102         if (is->mti == N+1)
00103             bn_randmt_seed(4357);
00104 
00105         for (kk = 0; kk < N-M; kk++) {
00106             y = (is->mt[kk]&UPPER_MASK)|(is->mt[kk+1]&LOWER_MASK);
00107             is->mt[kk] = is->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
00108         }
00109 
00110         for (; kk < N-1; kk++) {
00111             y = (is->mt[kk]&UPPER_MASK)|(is->mt[kk+1]&LOWER_MASK);
00112             is->mt[kk] = is->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
00113         }
00114 
00115         y = (is->mt[N-1]&UPPER_MASK)|(is->mt[0]&LOWER_MASK);
00116         is->mt[N-1] = is->mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
00117 
00118         is->mti = 0;
00119     }
00120 
00121     y = is->mt[is->mti++];
00122     y ^= TEMPERING_SHIFT_U(y);
00123     y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
00124     y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
00125     y ^= TEMPERING_SHIFT_L(y);
00126 
00127     return (double)y / (uint32_t)0xffffffff;
00128 }
00129 
00130 /* straight binhex would result in a 8+4992 character encoding. We should be able
00131  * to compress it quite a bit with better encoding? maybe uuencode? */
00132 void
00133 bn_randmt_state_serialize(struct _internal_state_s *UNUSED(is), struct bu_vls *UNUSED(s))
00134 {
00135     bu_bomb("Not implemented yet.\n");
00136 }
00137 
00138 void
00139 bn_randmt_state_deserialize(struct _internal_state_s *UNUSED(is), struct bu_vls *UNUSED(s))
00140 {
00141     bu_bomb("Not implemented yet.\n");
00142 }
00143 
00144 void
00145 bn_rand_mt_state_set_global(struct _internal_state_s *is)
00146 {
00147     global_state = is;
00148 }
00149 
00150 void
00151 bn_randmt_seed(unsigned long seed)
00152 {
00153     bn_randmt_state_seed(global_state, (uint32_t)seed);
00154 }
00155 
00156 double bn_randmt()
00157 {
00158     return bn_randmt_state(global_state);
00159 }
00160 
00161 /*
00162  * Local Variables:
00163  * mode: C
00164  * tab-width: 8
00165  * indent-tabs-mode: t
00166  * c-file-style: "stroustrup"
00167  * End:
00168  * ex: shiftwidth=4 tabstop=8
00169  */
Generated on Tue Dec 11 13:14:28 2012 for LIBBN by  doxygen 1.6.3