msr.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 "common.h"
00035
00036 #include <stdio.h>
00037 #include <math.h>
00038 #include "bu.h"
00039 #include "vmath.h"
00040 #include "bn.h"
00041
00042
00043
00044
00045 #define BN_MSR_MAXTBL 4096
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073 #define A 16807
00074 #define M 2147483647
00075 #define DM 2147483647.0
00076 #define Q 127773
00077 #define R 2836
00078 struct bn_unif *
00079 bn_unif_init(long int setseed, int method)
00080 {
00081 struct bn_unif *p;
00082 p = (struct bn_unif *) bu_malloc(sizeof(struct bn_unif), "bn_unif");
00083 p->msr_longs = (long *) bu_malloc(BN_MSR_MAXTBL*sizeof(long), "msr long table");
00084 p->msr_doubles=(double *) bu_malloc(BN_MSR_MAXTBL*sizeof(double), "msr double table");
00085 p->msr_seed = 1;
00086 p->msr_long_ptr = 0;
00087 p->msr_double_ptr = 0;
00088
00089 if (method != 0)
00090 bu_bomb("Method not yet supported in bn_unif_init()");
00091
00092 if (setseed&0x7fffffff) p->msr_seed=setseed&0x7fffffff;
00093 p->magic = BN_UNIF_MAGIC;
00094 return p;
00095 }
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125 long
00126 bn_unif_long_fill(struct bn_unif *p)
00127 {
00128 register long test, work_seed;
00129 register int i;
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139 if (!p || (p->magic != BN_UNIF_MAGIC &&
00140 p->magic != BN_GAUSS_MAGIC)) {
00141 BN_CK_UNIF(p);
00142 }
00143
00144 work_seed = p->msr_seed;
00145
00146 if ( p->msr_longs) {
00147 for (i=0; i < BN_MSR_MAXTBL; i++) {
00148 test = A*(work_seed % Q) - R*(work_seed / Q);
00149 p->msr_longs[i] = work_seed = (test < 0) ?
00150 test+M : test;
00151 }
00152 p->msr_long_ptr = BN_MSR_MAXTBL;
00153 }
00154 test = A*(work_seed % Q) - R*(work_seed / Q);
00155 p->msr_seed = (test < 0) ? test+M : test;
00156 return p->msr_seed;
00157 }
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187 double
00188 bn_unif_double_fill(struct bn_unif *p)
00189 {
00190 register long test, work_seed;
00191 register int i;
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201 if (!p || (p->magic != BN_UNIF_MAGIC &&
00202 p->magic != BN_GAUSS_MAGIC)) {
00203 BN_CK_UNIF(p);
00204 }
00205
00206 work_seed = p->msr_seed;
00207
00208 if (p->msr_doubles) {
00209 for (i=0; i < BN_MSR_MAXTBL; i++) {
00210 test = A*(work_seed % Q) - R*(work_seed / Q);
00211 work_seed = (test < 0) ? test+M : test;
00212 p->msr_doubles[i] = ( work_seed - M/2) * 1.0/DM;
00213 }
00214 p->msr_double_ptr = BN_MSR_MAXTBL;
00215 }
00216 test = A*(work_seed % Q) - R*(work_seed / Q);
00217 p->msr_seed = (test < 0) ? test+M : test;
00218
00219 return((p->msr_seed - M/2) * 1.0/DM);
00220 }
00221
00222
00223
00224
00225 void
00226 bn_unif_free(struct bn_unif *p)
00227 {
00228 bu_free(p->msr_doubles, "msr double table");
00229 bu_free(p->msr_longs, "msr long table");
00230 p->magic = 0;
00231 bu_free(p, "bn_unif");
00232 }
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259 struct bn_gauss *
00260 bn_gauss_init(long int setseed, int method)
00261 {
00262 struct bn_gauss *p;
00263
00264 if (method != 0)
00265 bu_bomb("Method not yet supported in bn_unif_init()");
00266
00267 p = (struct bn_gauss *) bu_malloc(sizeof(struct bn_gauss), "bn_msr_guass");
00268 p->msr_gausses=(double *) bu_malloc(BN_MSR_MAXTBL*sizeof(double), "msr guass table");
00269 p->msr_gauss_doubles=(double *) bu_malloc(BN_MSR_MAXTBL*sizeof(double), "msr guass doubles");
00270 p->msr_gauss_seed = 1;
00271 p->msr_gauss_ptr = 0;
00272 p->msr_gauss_dbl_ptr = 0;
00273
00274 if (setseed&0x7fffffff) p->msr_gauss_seed=setseed&0x7fffffff;
00275 p->magic = BN_GAUSS_MAGIC;
00276 return p;
00277 }
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310 double
00311 bn_gauss_fill(struct bn_gauss *p)
00312 {
00313 register int i;
00314 double v1, v2, r, fac;
00315
00316 BN_CK_GAUSS(p);
00317
00318 if (p->msr_gausses) {
00319 for (i=0; i< BN_MSR_MAXTBL-1; ) {
00320 BN_UNIF_CIRCLE((struct bn_unif *)p, v1, v2, r);
00321 if (r<0.00001) continue;
00322 fac = sqrt(-2.0*log(r)/r);
00323 p->msr_gausses[i++] = v1*fac;
00324 p->msr_gausses[i++] = v2*fac;
00325 }
00326 p->msr_gauss_ptr = BN_MSR_MAXTBL;
00327 }
00328
00329 do {
00330 BN_UNIF_CIRCLE((struct bn_unif *)p, v1, v2, r);
00331 } while (r < 0.00001);
00332 fac = sqrt(-2.0*log(r)/r);
00333 return v1*fac;
00334 }
00335
00336
00337
00338 void
00339 bn_gauss_free(struct bn_gauss *p)
00340 {
00341 bu_free(p->msr_gauss_doubles, "msr guass doubles");
00342 bu_free(p->msr_gausses, "msr guass table");
00343 bu_free(p, "bn_msr_guass");
00344 }
00345
00346
00347 #undef A
00348 #undef M
00349 #undef DM
00350 #undef Q
00351 #undef R
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362