BRL-CAD
randmt.c
Go to the documentation of this file.
1 /* R A N D M T . C
2  * BRL-CAD
3  *
4  * A C-program for MT19937: Real number version
5  * genrand() generates one pseudorandom real number (double)
6  * which is uniformly distributed on [0, 1]-interval, for each
7  * call. sgenrand(seed) set initial values to the working area
8  * of 624 words. Before genrand(), sgenrand(seed) must be
9  * called once. (seed is any 32-bit integer except for 0).
10  * Integer generator is obtained by modifying two lines.
11  * Coded by Takuji Nishimura, considering the suggestions by
12  * Topher Cooper and Marc Rieffel in July-Aug. 1997.
13  *
14  * This library is free software; you can redistribute it and/or
15  * modify it under the terms of the GNU Lesser General Public
16  * License as published by the Free Software Foundation; either
17  * version 2 of the License, or (at your option) any later version.
18  *
19  * This library is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
22  * See the GNU Lesser General Public License for more details.
23  * You should have received a copy of the GNU Library General
24  * Public License along with this library; if not, write to the
25  * Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
26  * 02111-1307 USA
27  *
28  * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
29  * Any feedback is very welcome. For any question, comments,
30  * see http://www.math.keio.ac.jp/matumoto/emt.html or email
31  * matumoto@math.keio.ac.jp
32  */
33 
34 #include "bu/log.h"
35 #include "bu/malloc.h"
36 #include "bn/randmt.h"
37 
38 /* Period parameters */
39 #define N 624
40 #define M 397
41 #define MATRIX_A 0x9908b0df /* constant vector a */
42 #define UPPER_MASK 0x80000000 /* most significant w-r bits */
43 #define LOWER_MASK 0x7fffffff /* least significant r bits */
44 
45 
46 /* Tempering parameters */
47 #define TEMPERING_MASK_B 0x9d2c5680
48 #define TEMPERING_MASK_C 0xefc60000
49 #define TEMPERING_SHIFT_U(y) (y >> 11)
50 #define TEMPERING_SHIFT_S(y) (y << 7)
51 #define TEMPERING_SHIFT_T(y) (y << 15)
52 #define TEMPERING_SHIFT_L(y) (y >> 18)
53 
54 #define MERSENNE_MAGIC 0x4D54524E
55 
57  uint32_t magic;
58  int mti; /* state index */
59  uint32_t mt[N]; /* state vector */
60 };
61 
62 static struct _internal_state_s global_state = { MERSENNE_MAGIC, N+1, {0} };
63 
64 void *
66 {
67  struct _internal_state_s *is;
68 
69  BU_ALLOC(is, struct _internal_state_s);
70  is->magic = MERSENNE_MAGIC;
71  is->mti = N+1;
72  return (void *)is;
73 }
74 
75 /* initializing the array with a NONZERO seed */
76 void
77 bn_randmt_state_seed(struct _internal_state_s *is, uint32_t seed)
78 {
79  /*
80  * setting initial seeds to mt[N] using
81  * the generator Line 25 of Table 1 in
82  * [KNUTH 1981, The Art of Computer Programming
83  * Vol. 2 (2nd Ed.), pp102]
84  */
85  is->mt[0] = seed & 0xffffffff;
86  for (is->mti = 1; is->mti<N; is->mti++)
87  is->mt[is->mti] = (69069 * is->mt[is->mti-1]) & 0xffffffff;
88 }
89 
90 double
92 {
93  uint32_t y;
94  static uint32_t mag01[2]={0x0, MATRIX_A};
95 
96  /* generate N words at one time */
97  if (is->mti >= N) {
98  int kk;
99 
100  /* if sgenrand() has not been called, a default initial seed is used */
101  if (is->mti == N+1)
102  bn_randmt_seed(4357);
103 
104  for (kk = 0; kk < N-M; kk++) {
105  y = (is->mt[kk]&UPPER_MASK)|(is->mt[kk+1]&LOWER_MASK);
106  is->mt[kk] = is->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
107  }
108 
109  for (; kk < N-1; kk++) {
110  y = (is->mt[kk]&UPPER_MASK)|(is->mt[kk+1]&LOWER_MASK);
111  is->mt[kk] = is->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
112  }
113 
114  y = (is->mt[N-1]&UPPER_MASK)|(is->mt[0]&LOWER_MASK);
115  is->mt[N-1] = is->mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
116 
117  is->mti = 0;
118  }
119 
120  y = is->mt[is->mti++];
121  y ^= TEMPERING_SHIFT_U(y);
124  y ^= TEMPERING_SHIFT_L(y);
125 
126  return (double)y / (uint32_t)0xffffffff;
127 }
128 
129 /* straight binhex would result in a 8+4992 character encoding. We should be able
130  * to compress it quite a bit with better encoding? maybe uuencode? */
131 void
133 {
134  bu_bomb("Not implemented yet.\n");
135 }
136 
137 void
139 {
140  bu_bomb("Not implemented yet.\n");
141 }
142 
143 void
144 bn_randmt_seed(unsigned long seed)
145 {
146  bn_randmt_state_seed(&global_state, (uint32_t)seed);
147 }
148 
149 double bn_randmt(void)
150 {
151  return bn_randmt_state(&global_state);
152 }
153 
154 /*
155  * Local Variables:
156  * mode: C
157  * tab-width: 8
158  * indent-tabs-mode: t
159  * c-file-style: "stroustrup"
160  * End:
161  * ex: shiftwidth=4 tabstop=8
162  */
#define MERSENNE_MAGIC
Definition: randmt.c:54
#define LOWER_MASK
Definition: randmt.c:43
#define TEMPERING_MASK_B
Definition: randmt.c:47
#define MATRIX_A
Definition: randmt.c:41
if lu s
Definition: nmg_mod.c:3860
#define TEMPERING_SHIFT_L(y)
Definition: randmt.c:52
#define N
Definition: randmt.c:39
void bn_randmt_state_serialize(struct _internal_state_s *is, struct bu_vls *s)
Definition: randmt.c:132
#define TEMPERING_SHIFT_U(y)
Definition: randmt.c:49
#define UPPER_MASK
Definition: randmt.c:42
double bn_randmt_state(struct _internal_state_s *is)
Definition: randmt.c:91
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
#define M
Definition: randmt.c:40
#define TEMPERING_MASK_C
Definition: randmt.c:48
uint32_t mt[N]
Definition: randmt.c:59
#define UNUSED(parameter)
Definition: common.h:239
#define TEMPERING_SHIFT_T(y)
Definition: randmt.c:51
void bn_randmt_seed(unsigned long seed)
Definition: randmt.c:144
void bn_randmt_state_seed(struct _internal_state_s *is, uint32_t seed)
Definition: randmt.c:77
void * bn_randmt_state_create(void)
Definition: randmt.c:65
#define TEMPERING_SHIFT_S(y)
Definition: randmt.c:50
uint32_t magic
Definition: randmt.c:57
Definition: vls.h:56
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
void bn_randmt_state_deserialize(struct _internal_state_s *is, struct bu_vls *s)
Definition: randmt.c:138
double bn_randmt(void)
Mersenne Twister random number generation as defined by MT19937.
Definition: randmt.c:149