ARTS 2.5.0 (git: 9ee3ac6c)
rng.cc
Go to the documentation of this file.
1/* Copyright (C) 2003-2012 Cory Davis <cory@met.ed.ac.uk>
2
3 This program is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 USA. */
17
31#include "rng.h"
32#include <algorithm>
33#include <climits>
34#include <cstddef>
35#include <cstdio>
36#include <cstdlib>
37#include <cstring>
38#include <iostream>
39#include <vector>
40
41#include "arts.h"
42#include "messages.h"
43
48
53
54void Rng::seed(unsigned long int n, const Verbosity &verbosity) {
55 // Static pool of previously used seeds.
56 static vector<unsigned long int> seeds;
57
58#pragma omp critical(Rng_seed)
59 {
60 unsigned long int n_orig = n;
61 //cout << "Requested seed: " << n;
62 while (find(seeds.begin(), seeds.end(), n) != seeds.end()) {
63 if (n >= ULONG_MAX - 1)
64 n = 0;
65 else
66 n++;
67
68 // If all possible seeds were already used, we empty the pool and
69 // start over.
70 if (n == n_orig) {
72 out0
73 << "Rng Warning: Couldn't find an unused seed. Clearing seed pool.\n";
74 seeds.clear();
75 break;
76 }
77 }
78 seeds.push_back(n);
79 seed_no = n;
80 //cout << " Got seed: " << seed_no << endl;
81 }
82
84}
85
89void Rng::force_seed(unsigned long int n) {
90 seed_no = n;
92}
93
97double Rng::draw() { return gsl_rng_uniform(r); }
98
102unsigned long int Rng::showseed() const { return seed_no; }
103
104/*
105 * rng/mt.c
106 *
107 * This program is free software; you can redistribute it and/or
108 * modify it under the terms of the GNU General Public License as
109 * published by the Free Software Foundation; either version 2 of the
110 * License, or (at your option) any later version.
111 *
112 * This program is distributed in the hope that it will be useful, but
113 * WITHOUT ANY WARRANTY; without even the implied warranty of
114 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
115 * General Public License for more details. You should have received
116 * a copy of the GNU General Public License along with this program;
117 * if not, write to the Free Foundation, Inc., 59 Temple Place, Suite
118 * 330, Boston, MA 02111-1307 USA.
119 *
120 * Original implementation was copyright (C) 1997 Makoto Matsumoto and
121 * Takuji Nishimura. Coded by Takuji Nishimura, considering the
122 * suggestions by Topher Cooper and Marc Rieffel in July-Aug. 1997, "A
123 * C-program for MT19937: Integer version (1998/4/6)".
124 *
125 * This implementation copyright (C) 1998 Brian Gough. I reorganized
126 * the code to use the module framework of GSL. The license on this
127 * implementation was changed from LGPL to GPL, following paragraph 3
128 * of the LGPL, version 2.
129 *
130 * Update:
131 *
132 * The seeding procedure has been updated to match the 10/99 release
133 * of MT19937.
134 *
135 * Update:
136 *
137 * The seeding procedure has been updated again to match the 2002
138 * release of MT19937.
139 *
140 * The original code included the comment: "When you use this, send an
141 * email to: matumoto@math.keio.ac.jp with an appropriate reference to
142 * your work".
143 *
144 * Makoto Matsumoto has a web page with more information about the
145 * generator, http://www.math.keio.ac.jp/~matumoto/emt.html.
146 *
147 * The paper below has details of the algorithm.
148 *
149 * From: Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A
150 * 623-dimensionally equidistributerd uniform pseudorandom number
151 *generator". ACM Transactions on Modeling and Computer Simulation,
152 * Vol. 8, No. 1 (Jan. 1998), Pages 3-30.
153 *
154 * You can obtain the paper directly from Makoto Matsumoto's web page.
155 *
156 * The period of this generator is 2^{19937} - 1.
157 *
158 */
159
160static inline unsigned long int mt_get(void *vstate);
161static double mt_get_double(void *vstate);
162static void mt_set(void *state, unsigned long int s);
163
164#define N 624 /* Period parameters */
165#define M 397
166
167/* most significant w-r bits */
168static const unsigned long UPPER_MASK = 0x80000000UL;
169
170/* least significant r bits */
171static const unsigned long LOWER_MASK = 0x7fffffffUL;
172
173typedef struct {
174 unsigned long mt[N];
175 int mti;
176} mt_state_t;
177
178static inline unsigned long mt_get(void *vstate) {
179 mt_state_t *state = (mt_state_t *)vstate;
180
181 unsigned long k;
182 unsigned long int *const mt = state->mt;
183
184#define MAGIC(y) (((y)&0x1) ? 0x9908b0dfUL : 0)
185
186 if (state->mti >= N) { /* generate N words at one time */
187 int kk;
188
189 for (kk = 0; kk < N - M; kk++) {
190 unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
191 mt[kk] = mt[kk + M] ^ (y >> 1) ^ MAGIC(y);
192 }
193 for (; kk < N - 1; kk++) {
194 unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
195 mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ MAGIC(y);
196 }
197
198 {
199 unsigned long y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
200 mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ MAGIC(y);
201 }
202
203 state->mti = 0;
204 }
205
206 /* Tempering */
207
208 k = mt[state->mti];
209 k ^= (k >> 11);
210 k ^= (k << 7) & 0x9d2c5680UL;
211 k ^= (k << 15) & 0xefc60000UL;
212 k ^= (k >> 18);
213
214 state->mti++;
215
216 return k;
217}
218
219static double mt_get_double(void *vstate) {
220 return (double)mt_get(vstate) / 4294967296.0;
221}
222
223static void mt_set(void *vstate, unsigned long int s) {
224 mt_state_t *state = (mt_state_t *)vstate;
225 int i;
226
227 if (s == 0) s = 4357; /* the default seed is 4357 */
228
229 state->mt[0] = s & 0xffffffffUL;
230
231 for (i = 1; i < N; i++) {
232 /* See Knuth's "Art of Computer Programming" Vol. 2, 3rd
233 Ed. p.106 for multiplier. */
234
235 state->mt[i] =
236 (1812433253UL * (state->mt[i - 1] ^ (state->mt[i - 1] >> 30)) + i);
237
238 state->mt[i] &= 0xffffffffUL;
239 }
240
241 state->mti = i;
242}
243
244static const gsl_rng_type mt_type = {"mt19937", /* name */
245 0xffffffffUL, /* RAND_MAX */
246 0, /* RAND_MIN */
247 sizeof(mt_state_t),
248 &mt_set,
249 &mt_get,
250 &mt_get_double};
251
253unsigned long int gsl_rng_default_seed = 0;
254/* rng/types.c
255 *
256 * Copyright (C) 2001 Brian Gough
257 *
258 * This program is free software; you can redistribute it and/or modify
259 * it under the terms of the GNU General Public License as published by
260 * the Free Software Foundation; either version 2 of the License, or (at
261 * your option) any later version.
262 *
263 * This program is distributed in the hope that it will be useful, but
264 * WITHOUT ANY WARRANTY; without even the implied warranty of
265 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
266 * General Public License for more details.
267 *
268 * You should have received a copy of the GNU General Public License
269 * along with this program; if not, write to the Free Software
270 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
271 */
272
273#define N1 100
274
276
277#define ADD(t) \
278 { \
279 if (i == N1) abort(); \
280 gsl_rng_generator_types[i] = (t); \
281 i++; \
282 };
283
285 int i = 0;
286
288 ADD(0);
289
291}
292
293/* err/stream.c
294 *
295 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
296 *
297 * This program is free software; you can redistribute it and/or modify
298 * it under the terms of the GNU General Public License as published by
299 * the Free Software Foundation; either version 2 of the License, or (at
300 * your option) any later version.
301 *
302 * This program is distributed in the hope that it will be useful, but
303 * WITHOUT ANY WARRANTY; without even the implied warranty of
304 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
305 * General Public License for more details.
306 *
307 * You should have received a copy of the GNU General Public License
308 * along with this program; if not, write to the Free Software
309 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
310 */
311
312FILE *gsl_stream = NULL;
314
315void gsl_stream_printf(const char *label,
316 const char *file,
317 int line,
318 const char *reason) {
319 if (gsl_stream == NULL) {
320 gsl_stream = stderr;
321 }
322 if (gsl_stream_handler) {
323 (*gsl_stream_handler)(label, file, line, reason);
324 return;
325 }
326 fprintf(gsl_stream, "gsl: %s:%d: %s: %s\n", file, line, label, reason);
327}
328
330 gsl_stream_handler_t *new_handler) {
331 gsl_stream_handler_t *previous_handler = gsl_stream_handler;
332 gsl_stream_handler = new_handler;
333 return previous_handler;
334}
335
336FILE *gsl_set_stream(FILE *new_stream) {
337 FILE *previous_stream;
338 if (gsl_stream == NULL) {
339 gsl_stream = stderr;
340 }
341 previous_stream = gsl_stream;
342 gsl_stream = new_stream;
343 return previous_stream;
344}
345
346/* err/error.c
347 *
348 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
349 *
350 * This program is free software; you can redistribute it and/or modify
351 * it under the terms of the GNU General Public License as published by
352 * the Free Software Foundation; either version 2 of the License, or (at
353 * your option) any later version.
354 *
355 * This program is distributed in the hope that it will be useful, but
356 * WITHOUT ANY WARRANTY; without even the implied warranty of
357 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
358 * General Public License for more details.
359 *
360 * You should have received a copy of the GNU General Public License
361 * along with this program; if not, write to the Free Software
362 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
363 */
364
366
367static void no_error_handler(const char *reason,
368 const char *file,
369 int line,
370 int gsl_errno);
371
372void gsl_error(const char *reason, const char *file, int line, int gsl_errno) {
373 if (gsl_error_handler) {
374 (*gsl_error_handler)(reason, file, line, gsl_errno);
375 return;
376 }
377
378 gsl_stream_printf("ERROR", file, line, reason);
379
380 fprintf(stderr, "Default GSL error handler invoked.\n");
381 abort();
382}
383
385 gsl_error_handler_t *previous_handler = gsl_error_handler;
386 gsl_error_handler = new_handler;
387 return previous_handler;
388}
389
391 gsl_error_handler_t *previous_handler = gsl_error_handler;
392 gsl_error_handler = no_error_handler;
393 return previous_handler;
394}
395
396static void no_error_handler(const char *reason _U_,
397 const char *file _U_,
398 int line _U_,
399 int gsl_errno _U_) {
400 /* do nothing */
401 return;
402}
403
404/* rng/rng.c
405 *
406 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
407 *
408 * This program is free software; you can redistribute it and/or modify
409 * it under the terms of the GNU General Public License as published by
410 * the Free Software Foundation; either version 2 of the License, or (at
411 * your option) any later version.
412 *
413 * This program is distributed in the hope that it will be useful, but
414 * WITHOUT ANY WARRANTY; without even the implied warranty of
415 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
416 * General Public License for more details.
417 *
418 * You should have received a copy of the GNU General Public License
419 * along with this program; if not, write to the Free Software
420 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
421 */
422
424 gsl_rng *r = (gsl_rng *)malloc(sizeof(gsl_rng));
425
426 if (r == 0) {
427 GSL_ERROR_VAL("failed to allocate space for rng struct", GSL_ENOMEM, 0);
428 };
429
430 r->state = malloc(T->size);
431
432 if (r->state == 0) {
433 free(r); /* exception in constructor, avoid memory leak */
434
435 GSL_ERROR_VAL("failed to allocate space for rng state", GSL_ENOMEM, 0);
436 };
437
438 r->type = T;
439
440 gsl_rng_set(r, gsl_rng_default_seed); /* seed the generator */
441
442 return r;
443}
444
445int gsl_rng_memcpy(gsl_rng *dest, const gsl_rng *src) {
446 if (dest->type != src->type) {
447 GSL_ERROR("generators must be of the same type", GSL_EINVAL);
448 }
449
450 memcpy(dest->state, src->state, src->type->size);
451
452 return GSL_SUCCESS;
453}
454
456 gsl_rng *r = (gsl_rng *)malloc(sizeof(gsl_rng));
457
458 if (r == 0) {
459 GSL_ERROR_VAL("failed to allocate space for rng struct", GSL_ENOMEM, 0);
460 };
461
462 r->state = malloc(q->type->size);
463
464 if (r->state == 0) {
465 free(r); /* exception in constructor, avoid memory leak */
466
467 GSL_ERROR_VAL("failed to allocate space for rng state", GSL_ENOMEM, 0);
468 };
469
470 r->type = q->type;
471
472 memcpy(r->state, q->state, q->type->size);
473
474 return r;
475}
476
477void gsl_rng_set(const gsl_rng *r, unsigned long int seed) {
478 (r->type->set)(r->state, seed);
479}
480
481#ifndef HIDE_INLINE_STATIC
482unsigned long int gsl_rng_get(const gsl_rng *r) {
483 return (r->type->get)(r->state);
484}
485
486double gsl_rng_uniform(const gsl_rng *r) {
487 return (r->type->get_double)(r->state);
488}
489
490double gsl_rng_uniform_pos(const gsl_rng *r) {
491 double x;
492 do {
493 x = (r->type->get_double)(r->state);
494 } while (x == 0);
495
496 return x;
497}
498
499unsigned long int gsl_rng_uniform_int(const gsl_rng *r, unsigned long int n) {
500 unsigned long int offset = r->type->min;
501 unsigned long int range = r->type->max - offset;
502 unsigned long int scale = range / n;
503 unsigned long int k;
504
505 if (n > range) {
506 GSL_ERROR_VAL("n exceeds maximum value of generator", GSL_EINVAL, 0);
507 }
508
509 do {
510 k = (((r->type->get)(r->state)) - offset) / scale;
511 } while (k >= n);
512
513 return k;
514}
515#endif
516
517unsigned long int gsl_rng_max(const gsl_rng *r) { return r->type->max; }
518
519unsigned long int gsl_rng_min(const gsl_rng *r) { return r->type->min; }
520
521const char *gsl_rng_name(const gsl_rng *r) { return r->type->name; }
522
523size_t gsl_rng_size(const gsl_rng *r) { return r->type->size; }
524
525void *gsl_rng_state(const gsl_rng *r) { return r->state; }
526
528 size_t i;
529 unsigned char *p = (unsigned char *)(r->state);
530 const size_t n = r->type->size;
531
532 for (i = 0; i < n; i++) {
533 /* FIXME: we're assuming that a char is 8 bits */
534 printf("%.2x", *(p + i));
535 }
536}
537
539 free(r->state);
540 free(r);
541}
The global header file for ARTS.
double draw()
Draws a double from the uniform distribution [0,1)
Definition: rng.cc:97
Rng()
Constructor creates instance of gsl_rng of type gsl_rng_mt19937.
Definition: rng.cc:47
unsigned long int seed_no
Definition: rng.h:557
void seed(unsigned long int n, const Verbosity &verbosity)
Seeds the Rng with the integer argument.
Definition: rng.cc:54
unsigned long int showseed() const
Returns the seed number.
Definition: rng.cc:102
~Rng()
Destructor frees memory allocated to gsl_rng.
Definition: rng.cc:52
gsl_rng * r
Definition: rng.h:555
void force_seed(unsigned long int n)
Seeds the Rng with the integer argument.
Definition: rng.cc:89
#define _U_
Definition: config.h:180
#define q
Declarations having to do with the four output streams.
#define CREATE_OUT0
Definition: messages.h:204
#define N
Definition: rng.cc:164
void gsl_stream_printf(const char *label, const char *file, int line, const char *reason)
Definition: rng.cc:315
void gsl_rng_print_state(const gsl_rng *r)
Definition: rng.cc:527
unsigned long int gsl_rng_min(const gsl_rng *r)
Definition: rng.cc:519
double gsl_rng_uniform_pos(const gsl_rng *r)
Definition: rng.cc:490
#define MAGIC(y)
gsl_error_handler_t * gsl_set_error_handler(gsl_error_handler_t *new_handler)
Definition: rng.cc:384
gsl_rng * gsl_rng_alloc(const gsl_rng_type *T)
Definition: rng.cc:423
void gsl_rng_set(const gsl_rng *r, unsigned long int seed)
Definition: rng.cc:477
#define ADD(t)
Definition: rng.cc:277
void gsl_error(const char *reason, const char *file, int line, int gsl_errno)
Definition: rng.cc:372
unsigned long int gsl_rng_max(const gsl_rng *r)
Definition: rng.cc:517
#define M
Definition: rng.cc:165
gsl_error_handler_t * gsl_set_error_handler_off(void)
Definition: rng.cc:390
unsigned long int gsl_rng_uniform_int(const gsl_rng *r, unsigned long int n)
Definition: rng.cc:499
void * gsl_rng_state(const gsl_rng *r)
Definition: rng.cc:525
const gsl_rng_type * gsl_rng_mt19937
Definition: rng.cc:252
FILE * gsl_stream
Definition: rng.cc:312
unsigned long int gsl_rng_get(const gsl_rng *r)
Definition: rng.cc:482
gsl_stream_handler_t * gsl_set_stream_handler(gsl_stream_handler_t *new_handler)
Definition: rng.cc:329
gsl_rng * gsl_rng_clone(const gsl_rng *q)
Definition: rng.cc:455
unsigned long int gsl_rng_default_seed
Definition: rng.cc:253
const gsl_rng_type ** gsl_rng_types_setup(void)
Definition: rng.cc:284
gsl_error_handler_t * gsl_error_handler
Definition: rng.cc:365
#define N1
Definition: rng.cc:273
const gsl_rng_type * gsl_rng_generator_types[N1]
Definition: rng.cc:275
FILE * gsl_set_stream(FILE *new_stream)
Definition: rng.cc:336
size_t gsl_rng_size(const gsl_rng *r)
Definition: rng.cc:523
int gsl_rng_memcpy(gsl_rng *dest, const gsl_rng *src)
Definition: rng.cc:445
gsl_stream_handler_t * gsl_stream_handler
Definition: rng.cc:313
double gsl_rng_uniform(const gsl_rng *r)
Definition: rng.cc:486
const char * gsl_rng_name(const gsl_rng *r)
Definition: rng.cc:521
void gsl_rng_free(gsl_rng *r)
Definition: rng.cc:538
member functions of the Rng class and gsl_rng code
void gsl_stream_handler_t(const char *label, const char *file, int line, const char *reason)
Definition: rng.h:366
@ GSL_ENOMEM
Definition: rng.h:323
@ GSL_SUCCESS
Definition: rng.h:313
@ GSL_EINVAL
Definition: rng.h:319
void gsl_error_handler_t(const char *reason, const char *file, int line, int gsl_errno)
Definition: rng.h:361
#define GSL_ERROR(reason, gsl_errno)
Definition: rng.h:381
#define GSL_ERROR_VAL(reason, gsl_errno, value)
Definition: rng.h:389
unsigned long int max
Definition: rng.h:120
size_t size
Definition: rng.h:122
unsigned long int(* get)(void *state)
Definition: rng.h:124
void(* set)(void *state, unsigned long int seed)
Definition: rng.h:123
double(* get_double)(void *state)
Definition: rng.h:125
const char * name
Definition: rng.h:119
unsigned long int min
Definition: rng.h:121
Definition: rng.h:128
void * state
Definition: rng.h:130
const gsl_rng_type * type
Definition: rng.h:129
unsigned long mt[N]
Definition: rng.cc:174
int mti
Definition: rng.cc:175