git.gag.com Git - debian/gcpegg/blob - lecuyer.c

   1 /*
   2
   3     L'Ecuyer's two-sequence pseudorandom generator with a 32 cell
   4     Bays-Durham shuffle on the back end.
   5
   6     When updating the individual generators, an intermediate value as
   7     large as 64 bits is required.  Many modern C compilers provide a
   8     "long long" 64-bit integer type.  If this type is available,
   9     compile this file with LONGLONG defined and it will be used,
  10     generally improving execution speed.  If the compiler does not
  11     implement "long long", compile without LONGLONG and Schrage's
  12     algorithm will be used to perform the computation using only
  13     32 bit arithmetic.
  14
  15     For additional details, see L'Ecuyer's paper:
  16
  17         L'Ecuyer, P. "Communications of the ACM", Vol. 31, p. 742 (1988).
  18
  19     and that of Bays and Durham:
  20
  21         Bays, Carter, and S.D. Durham.  "ACM Transactions on
  22             Mathematical Software", Vol. 2, p. 59 (1976).
  23
  24     Schrage's multiple-precision modular algorithm is described in:
  25
  26         Schrage, P.  "ACM Transactions on Mathematical Software",
  27             Vol. 5, p. 132 (1979).
  28
  29 */
  30
  31 #include <stdio.h>
  32 #include <assert.h>
  33
  34 #define LONGLONG
  35
  36 #ifdef LONGLONG
  37 #define int64 long long
  38 #else
  39 #define int64 long
  40 #endif
  41 #define int32 long
  42
  43 /* L'Ecuyer's recommended multiplier and modulus for the two
  44    multiplicative congruential generators.  Even though the
  45    values fit in 32 bits, we declare them as int64 so that the
  46    arithmetic in calculating the next value will be automatically
  47    done in int64 without need for casting when LONGLONG is defined.
  48    In a non-LONGLONG build int64 is defined as long, so these values
  49    participate correctly in the Schrage algorithm computation. */
  50
  51 #define mul1 ((int64) 40014)
  52 #define mod1 ((int64) 2147483563)
  53 #define mul2 ((int64) 40692)
  54 #define mod2 ((int64) 2147483399)
  55
  56 #define shuffleSize 32                /* Shuffle table size */
  57 #define warmup 19                     /* Number of initial warmup results to "burn" */
  58
  59 static int32 gen1, gen2, state;
  60 static int32 shuffle[shuffleSize];    /* Bays-Durham shuffle table */
  61
  62 /*  Update generator which, using either long long arithmetic or
  63     Schrage's algorithm depending on whether LONGLONG is defined.
  64     The definition for Schrage's algorithm relies, for efficiency,
  65     on the C compiler performing compile-time constant arithmetic
  66     for the quotient and remainder of the modulus and multiplier.
  67     If you're porting this to a language which lacks that feature,
  68     you'll want to predefine the quotient and remainder for each
  69     generator and use the explicit values.  */
  70
  71 #ifdef LONGLONG
  72 #define updgen(which)   gen##which = (int32) ((gen##which * mul##which) % mod##which)
  73 #else
  74 #define updgen(which) {                                                               \
  75         int32 t = gen##which / (mod##which / mul##which);                             \
  76                                                                                       \
  77         gen##which = mul##which * (gen##which - (t * (mod##which / mul##which))) -    \
  78                      t * ((mod##which % mul##which));                                 \
  79         if (gen##which < 0) {                                                         \
  80             gen##which += mod##which;                                                 \
  81         }                                                                             \
  82     }
  83 #endif
  84
  85 /*  LEsetSeed  --  Set seed for generator.  Subsequent values will be based
  86                    on the given nonzero seed.  */
  87
  88 void LEsetSeed(int32 seed)
  89 {
  90     int32 i;
  91
  92     assert(seed != 0);
  93
  94     gen1 = gen2 = (int32) (seed & 0x7FFFFFFFL);
  95
  96     /* "Warm up" the generator for a number of rounds to eliminate
  97        any residual inflence of the seed. */
  98
  99     for (i = 0; i < warmup; i++) {
 100         updgen(1);
 101     }
 102
 103     /* Fill the shuffle table with values.  */
 104
 105     for (i = 0; i < shuffleSize; i++) {
 106         updgen(1);
 107         shuffle[(shuffleSize - 1) - i] = gen1;
 108     }
 109     state = shuffle[0];
 110 }
 111
 112 /*  LEnextByte  --  Get next byte from generator.  */
 113
 114 unsigned char LEnextByte(void)
 115 {
 116     int i;
 117
 118     updgen(1);
 119     updgen(2);
 120
 121     /* Extract shuffle table index from most significant part
 122        of the previous result. */
 123
 124     i = state / (1 + (((int32) mod1) - 1) / shuffleSize);
 125
 126     /* New state is sum of generators modulo one of their moduli.  */
 127
 128     state = (int32) (((shuffle[i]) + ((unsigned int32) gen2)) % mod1);
 129
 130     /* Replace value in shuffle table with generator 1 result.  */
 131
 132     shuffle[i] = gen1;
 133
 134     return (unsigned char) (state / (1 + (((int32) mod1) - 1) / 256));
 135 }