BRL-CAD
msr.c
Go to the documentation of this file.
1 /* M S R . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1994-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 
21 /** @addtogroup msr */
22 /** @{ */
23 /** @file libbn/msr.c
24  *
25  * @brief
26  * Minimal Standard RANdom number generator
27  *
28  * @par From:
29  * Stephen K. Park and Keith W. Miller
30  * @n "Random number generators: good ones are hard to find"
31  * @n CACM vol 31 no 10, Oct 88
32  *
33  */
34 
35 #include "common.h"
36 
37 #include <stdio.h>
38 #include <math.h>
39 
40 #include "bu/malloc.h"
41 #include "bu/log.h"
42 #include "vmath.h"
43 #include "bn/msr.h"
44 
45 /**
46  * Note: BN_MSR_MAXTBL must be an even number, preferably a power of two.
47  */
48 #define BN_MSR_MAXTBL 4096 /* Size of random number tables. */
49 
50 
51 #define A 16807
52 #define M 2147483647
53 #define DM 2147483647.0
54 #define Q 127773 /* Q = M / A */
55 #define R 2836 /* R = M % A */
56 
57 struct bn_unif *
58 bn_unif_init(long int setseed, int method)
59 {
60  struct bn_unif *p;
61  BU_ALLOC(p, struct bn_unif);
62  p->msr_longs = (long *) bu_malloc(BN_MSR_MAXTBL*sizeof(long), "msr long table");
63  p->msr_doubles=(double *) bu_malloc(BN_MSR_MAXTBL*sizeof(double), "msr double table");
64  p->msr_seed = 1;
65  p->msr_long_ptr = 0;
66  p->msr_double_ptr = 0;
67 
68  if (method != 0)
69  bu_bomb("Method not yet supported in bn_unif_init()");
70 
71  if (setseed&0x7fffffff) p->msr_seed=setseed&0x7fffffff;
72  p->magic = BN_UNIF_MAGIC;
73  return p;
74 }
75 
76 
77 long
79 {
80  register long test, work_seed;
81  register int i;
82 
83  if (!p)
84  return 1;
85 
86  /*
87  * Gauss and uniform structures have the same format for the
88  * first part (gauss is an extension of uniform) so that a gauss
89  * structure can be passed to a uniform routine. This means
90  * that we only maintain one structure for gaussian distributions
91  * rather than two. It also means that the user can pull
92  * uniform numbers from a gauss structure when the user wants.
93  */
94  if (p->magic != BN_UNIF_MAGIC && p->magic != BN_GAUSS_MAGIC) {
95  BN_CK_UNIF(p);
96  }
97 
98  work_seed = p->msr_seed;
99 
100  if ( p->msr_longs) {
101  for (i=0; i < BN_MSR_MAXTBL; i++) {
102  test = A*(work_seed % Q) - R*(work_seed / Q);
103  p->msr_longs[i] = work_seed = (test < 0) ?
104  test+M : test;
105  }
107  }
108  test = A*(work_seed % Q) - R*(work_seed / Q);
109  p->msr_seed = (test < 0) ? test+M : test;
110  return p->msr_seed;
111 }
112 
113 
114 double
116 {
117  register long test, work_seed;
118  register int i;
119 
120  if (!p)
121  return 0.0;
122 
123  /*
124  * Gauss and uniform structures have the same format for the
125  * first part (gauss is an extension of uniform) so that a gauss
126  * structure can be passed to a uniform routine. This means
127  * that we only maintain one structure for gaussian distributions
128  * rather than two. It also means that the user can pull
129  * uniform numbers from a gauss structure when the user wants.
130  */
131  if (p->magic != BN_UNIF_MAGIC && p->magic != BN_GAUSS_MAGIC) {
132  BN_CK_UNIF(p);
133  }
134 
135  work_seed = p->msr_seed;
136 
137  if (p->msr_doubles) {
138  for (i=0; i < BN_MSR_MAXTBL; i++) {
139  test = A*(work_seed % Q) - R*(work_seed / Q);
140  work_seed = (test < 0) ? test+M : test;
141  p->msr_doubles[i] = ( work_seed - M/2) * 1.0/DM;
142  }
144  }
145  test = A*(work_seed % Q) - R*(work_seed / Q);
146  p->msr_seed = (test < 0) ? test+M : test;
147 
148  return ((p->msr_seed - M/2) * 1.0/DM);
149 }
150 
151 /* bn_unif_free free random number table
152  *
153  */
154 void
156 {
157  if (!p)
158  return;
159 
160  bu_free(p->msr_doubles, "msr double table");
161  bu_free(p->msr_longs, "msr long table");
162  p->magic = 0;
163  bu_free(p, "bn_unif");
164 }
165 
166 
167 
168 struct bn_gauss *
169 bn_gauss_init(long int setseed, int method)
170 {
171  struct bn_gauss *p;
172 
173  if (method != 0)
174  bu_bomb("Method not yet supported in bn_unif_init()");
175 
176  BU_ALLOC(p, struct bn_gauss);
177  p->msr_gausses=(double *) bu_malloc(BN_MSR_MAXTBL*sizeof(double), "msr gauss table");
178  p->msr_gauss_doubles=(double *) bu_malloc(BN_MSR_MAXTBL*sizeof(double), "msr gauss doubles");
179  p->msr_gauss_seed = 1;
180  p->msr_gauss_ptr = 0;
181  p->msr_gauss_dbl_ptr = 0;
182 
183  if (setseed&0x7fffffff) p->msr_gauss_seed=setseed&0x7fffffff;
184  p->magic = BN_GAUSS_MAGIC;
185  return p;
186 }
187 
188 
189 double
191 {
192  register int i;
193  /* register */ double v1, v2, r, fac;
194 
195  if (!p)
196  return 0.0;
197 
198  BN_CK_GAUSS(p);
199 
200  if (p->msr_gausses) {
201  for (i=0; i< BN_MSR_MAXTBL-1; ) {
202  BN_UNIF_CIRCLE((struct bn_unif *)p, v1, v2, r);
203  if (r<0.00001) continue;
204  fac = sqrt(-2.0*log(r)/r);
205  p->msr_gausses[i++] = v1*fac;
206  p->msr_gausses[i++] = v2*fac;
207  }
209  }
210 
211  do {
212  BN_UNIF_CIRCLE((struct bn_unif *)p, v1, v2, r);
213  } while (r < 0.00001);
214  fac = sqrt(-2.0*log(r)/r);
215  return v1*fac;
216 }
217 
218 
219 void
221 {
222  bu_free(p->msr_gauss_doubles, "msr gauss doubles");
223  bu_free(p->msr_gausses, "msr gauss table");
224  bu_free(p, "bn_msr_gauss");
225 }
226 
227 
228 #undef A
229 #undef M
230 #undef DM
231 #undef Q
232 #undef R
233 
234 /** @} */
235 /*
236  * Local Variables:
237  * mode: C
238  * tab-width: 8
239  * indent-tabs-mode: t
240  * c-file-style: "stroustrup"
241  * End:
242  * ex: shiftwidth=4 tabstop=8
243  */
Minimal Standard RANdom number generator.
Definition: msr.h:53
int msr_long_ptr
Definition: msr.h:58
struct bn_unif * bn_unif_init(long int setseed, int method)
Definition: msr.c:58
int msr_double_ptr
Definition: msr.h:56
#define BN_UNIF_CIRCLE(_p, _x, _y, _r)
Definition: msr.h:250
#define DM
Definition: msr.c:53
Header file for the BRL-CAD common definitions.
void bn_unif_free(struct bn_unif *p)
Definition: msr.c:155
long bn_unif_long_fill(struct bn_unif *p)
Definition: msr.c:78
void bn_gauss_free(struct bn_gauss *p)
Definition: msr.c:220
Definition: msr.h:71
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
double * msr_gauss_doubles
Definition: msr.h:75
long msr_seed
Definition: msr.h:55
int msr_gauss_ptr
Definition: msr.h:76
#define BN_MSR_MAXTBL
Definition: msr.c:48
double * msr_gausses
Definition: msr.h:77
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
#define BN_CK_UNIF(_p)
Definition: msr.h:62
#define BN_UNIF_MAGIC
Definition: magic.h:75
double * msr_doubles
Definition: msr.h:57
#define BN_CK_GAUSS(_p)
Definition: msr.h:63
int msr_gauss_dbl_ptr
Definition: msr.h:74
uint32_t magic
Definition: msr.h:72
uint32_t magic
Definition: msr.h:54
long * msr_longs
Definition: msr.h:59
#define BN_GAUSS_MAGIC
Definition: magic.h:69
double bn_unif_double_fill(struct bn_unif *p)
Definition: msr.c:115
#define R
Definition: msr.c:55
#define A
Definition: msr.c:51
double bn_gauss_fill(struct bn_gauss *p)
Definition: msr.c:190
struct bn_gauss * bn_gauss_init(long int setseed, int method)
Definition: msr.c:169
#define Q
Definition: msr.c:54
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
long msr_gauss_seed
Definition: msr.h:73
#define M
Definition: msr.c:52