randmt.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #include "bu.h"
00035 #include "bn.h"
00036
00037
00038 #define N 624
00039 #define M 397
00040 #define MATRIX_A 0x9908b0df
00041 #define UPPER_MASK 0x80000000
00042 #define LOWER_MASK 0x7fffffff
00043
00044
00045
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;
00058 uint32_t mt[N];
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
00077 void
00078 bn_randmt_state_seed(struct _internal_state_s *is, uint32_t seed)
00079 {
00080
00081
00082
00083
00084
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
00098 if (is->mti >= N) {
00099 int kk;
00100
00101
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
00131
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
00163
00164
00165
00166
00167
00168
00169