msr.c

Go to the documentation of this file.
00001 /*                           M S R . C
00002  * BRL-CAD
00003  *
00004  * Copyright (c) 1994-2012 United States Government as represented by
00005  * the U.S. Army Research Laboratory.
00006  *
00007  * This library is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU Lesser General Public License
00009  * version 2.1 as published by the Free Software Foundation.
00010  *
00011  * This library is distributed in the hope that it will be useful, but
00012  * WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014  * Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public
00017  * License along with this file; see the file named COPYING for more
00018  * information.
00019  */
00020 /** @addtogroup msr */
00021 /** @{ */
00022 /** @file libbn/msr.c
00023  *
00024  * @brief
00025  * Minimal Standard RANdom number generator
00026  *
00027  * @par From:
00028  *      Stephen K. Park and Keith W. Miller
00029  * @n   "Random number generators: good ones are hard to find"
00030  * @n   CACM vol 31 no 10, Oct 88
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  * Note: BN_MSR_MAXTBL must be an even number, preferably a power of two.
00044  */
00045 #define BN_MSR_MAXTBL   4096    /**< Size of random number tables. */
00046 
00047 
00048 /*      bn_unif_init    Initialize a random number structure.
00049  *
00050  * @par Entry
00051  *      setseed seed to use
00052  *      method  method to use to generate numbers;
00053  *
00054  * @par Exit
00055  *      returns a pointer to a bn_unif structure.
00056  *      returns 0 on error.
00057  *
00058  * @par Uses
00059  *      None.
00060  *
00061  * @par Calls
00062  *      bu_malloc
00063  *
00064  * @par Method @code
00065  *      malloc up a structure with pointers to the numbers
00066  *      get space for the integer table.
00067  *      get space for the floating table.
00068  *      set pointer counters
00069  *      set seed if one was given and setseed != 1
00070  @endcode
00071  *
00072  */
00073 #define A       16807
00074 #define M       2147483647
00075 #define DM      2147483647.0
00076 #define Q       127773          /**< Q = M / A */
00077 #define R       2836            /**< R = M % A */
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 /**     bn_unif_long_fill       fill a random number table.
00098  *
00099  * Use the msrad algorithm to fill a random number table
00100  * with values from 1 to 2^31-1.  These numbers can (and are) extracted from
00101  * the random number table via high speed macros and bn_unif_long_fill called
00102  * when the table is exauseted.
00103  *
00104  * @par Entry
00105  *      p       pointer to a bn_unif structure.
00106  *
00107  * @par Exit
00108  *      if (!p) returns 1 else returns a value between 1 and 2^31-1
00109  *
00110  * @par Calls
00111  *      None.  msran is inlined for speed reasons.
00112  *
00113  * @par Uses
00114  *      None.
00115  *
00116  * @par Method @code
00117  *      if (!p) return 1;
00118  *      if p->msr_longs != NULL
00119  *              msr_longs is reloaded with random numbers;
00120  *              msr_long_ptr is set to BN_MSR_MAXTBL
00121  *      endif
00122  *      msr_seed is updated.
00123  @endcode
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      * Gauss and uniform structures have the same format for the
00133      * first part (gauss is an extention of uniform) so that a gauss
00134      * structure can be passed to a uniform routine.  This means
00135      * that we only maintain one structure for gaussian distributions
00136      * rather than two.  It also means that the user can pull
00137      * uniform numbers from a guass structure when the user wants.
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 /**     bn_unif_double_fill     fill a random number table.
00160  *
00161  * Use the msrad algorithm to fill a random number table
00162  * with values from -0.5 to 0.5.  These numbers can (and are) extracted from
00163  * the random number table via high speed macros and bn_unif_double_fill
00164  * called when the table is exauseted.
00165  *
00166  * @par Entry
00167  *      p       pointer to a bn_unif structure.
00168  *
00169  * @par Exit
00170  *      if (!p) returns 0.0 else returns a value between -0.5 and 0.5
00171  *
00172  * @par Calls
00173  *      None.  msran is inlined for speed reasons.
00174  *
00175  * @par Uses
00176  *      None.
00177  *
00178  * @par Method @code
00179  *      if (!p) return (0.0)
00180  *      if p->msr_longs != NULL
00181  *              msr_longs is reloaded with random numbers;
00182  *              msr_long_ptr is set to BN_MSR_MAXTBL
00183  *      endif
00184  *      msr_seed is updated.
00185  @endcode
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      * Gauss and uniform structures have the same format for the
00195      * first part (gauss is an extention of uniform) so that a gauss
00196      * structure can be passed to a uniform routine.  This means
00197      * that we only maintain one structure for gaussian distributions
00198      * rather than two.  It also means that the user can pull
00199      * uniform numbers from a guass structure when the user wants.
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 /*      bn_unif_free    free random number table
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 /*      bn_gauss_init   Initialize a random number struct for gaussian
00236  *      numbers.
00237  *
00238  * @par Entry
00239  *      setseed         Seed to use.
00240  *      method          method to use to generate numbers (not used)
00241  *
00242  * @par Exit
00243  *      Returns a pointer toa bn_msr_guass structure.
00244  *      returns 0 on error.
00245  *
00246  * @par Calls
00247  *      bu_malloc
00248  *
00249  * @par Uses
00250  *      None.
00251  *
00252  * @par Method @code
00253  *      malloc up a structure
00254  *      get table space
00255  *      set seed and pointer.
00256  *      if setseed != 0 then seed = setseed
00257  @endcode
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 /*      bn_gauss_fill   fill a random number table.
00280  *
00281  * Use the msrad algorithm to fill a random number table.
00282  * hese numbers can (and are) extracted from
00283  * the random number table via high speed macros and bn_msr_guass_fill
00284  * called when the table is exauseted.
00285  *
00286  * @par Entry
00287  *      p       pointer to a bn_msr_guass structure.
00288  *
00289  * @par Exit
00290  *      if (!p) returns 0.0 else returns a value with a mean of 0 and
00291  *          a variance of 1.0.
00292  *
00293  * @par Calls
00294  *      BN_UNIF_CIRCLE to get to uniform random number whos radius is
00295  *      <= 1.0. I.e. sqrt(v1*v1 + v2*v2) <= 1.0
00296  *      BN_UNIF_CIRCLE is a macro which can call bn_unif_double_fill.
00297  *
00298  * @par Uses
00299  *      None.
00300  *
00301  * @par Method @code
00302  *      if (!p) return (0.0)
00303  *      if p->msr_longs != NULL
00304  *              msr_longs is reloaded with random numbers;
00305  *              msr_long_ptr is set to BN_MSR_MAXTBL
00306  *      endif
00307  *      msr_seed is updated.
00308  @endcode
00309 */
00310 double
00311 bn_gauss_fill(struct bn_gauss *p)
00312 {
00313     register int i;
00314     /* register */ 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 /*      bn_gauss_free   free random number table
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  * Local Variables:
00356  * mode: C
00357  * tab-width: 8
00358  * indent-tabs-mode: t
00359  * c-file-style: "stroustrup"
00360  * End:
00361  * ex: shiftwidth=4 tabstop=8
00362  */
Generated on Tue Dec 11 13:14:27 2012 for LIBBN by  doxygen 1.6.3