File: lecuyer.c

package info (click to toggle)
gcpegg 5.1-9
  • links: PTS
  • area: main
  • in suites: etch, etch-m68k
  • size: 328 kB
  • ctags: 445
  • sloc: ansic: 3,706; makefile: 95; sh: 27; csh: 21
file content (135 lines) | stat: -rw-r--r-- 4,524 bytes parent folder | download | duplicates (6)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
/*

    L'Ecuyer's two-sequence pseudorandom generator with a 32 cell
    Bays-Durham shuffle on the back end.

    When updating the individual generators, an intermediate value as
    large as 64 bits is required.  Many modern C compilers provide a
    "long long" 64-bit integer type.  If this type is available,
    compile this file with LONGLONG defined and it will be used,
    generally improving execution speed.  If the compiler does not
    implement "long long", compile without LONGLONG and Schrage's
    algorithm will be used to perform the computation using only
    32 bit arithmetic.

    For additional details, see L'Ecuyer's paper:

        L'Ecuyer, P. "Communications of the ACM", Vol. 31, p. 742 (1988).

    and that of Bays and Durham:

        Bays, Carter, and S.D. Durham.  "ACM Transactions on
            Mathematical Software", Vol. 2, p. 59 (1976).

    Schrage's multiple-precision modular algorithm is described in:

        Schrage, P.  "ACM Transactions on Mathematical Software",
	    Vol. 5, p. 132 (1979).

*/

#include <stdio.h>
#include <assert.h>

#define LONGLONG

#ifdef LONGLONG
#define int64 long long
#else
#define int64 long
#endif
#define int32 long

/* L'Ecuyer's recommended multiplier and modulus for the two
   multiplicative congruential generators.  Even though the
   values fit in 32 bits, we declare them as int64 so that the
   arithmetic in calculating the next value will be automatically
   done in int64 without need for casting when LONGLONG is defined.
   In a non-LONGLONG build int64 is defined as long, so these values
   participate correctly in the Schrage algorithm computation. */

#define mul1 ((int64) 40014)
#define mod1 ((int64) 2147483563)
#define mul2 ((int64) 40692)
#define mod2 ((int64) 2147483399)

#define shuffleSize 32		      /* Shuffle table size */
#define warmup 19                     /* Number of initial warmup results to "burn" */

static int32 gen1, gen2, state;
static int32 shuffle[shuffleSize];    /* Bays-Durham shuffle table */

/*  Update generator which, using either long long arithmetic or
    Schrage's algorithm depending on whether LONGLONG is defined.
    The definition for Schrage's algorithm relies, for efficiency,
    on the C compiler performing compile-time constant arithmetic
    for the quotient and remainder of the modulus and multiplier.
    If you're porting this to a language which lacks that feature,
    you'll want to predefine the quotient and remainder for each
    generator and use the explicit values.  */

#ifdef LONGLONG
#define updgen(which)	gen##which = (int32) ((gen##which * mul##which) % mod##which)
#else
#define updgen(which) {                                                               \
        int32 t = gen##which / (mod##which / mul##which);                             \
                                                                                      \
        gen##which = mul##which * (gen##which - (t * (mod##which / mul##which))) -    \
                     t * ((mod##which % mul##which));                                 \
        if (gen##which < 0) {                                                         \
            gen##which += mod##which;                                                 \
        }                                                                             \
    }
#endif

/*  LEsetSeed  --  Set seed for generator.  Subsequent values will be based
		   on the given nonzero seed.  */

void LEsetSeed(int32 seed)
{
    int32 i;

    assert(seed != 0);

    gen1 = gen2 = (int32) (seed & 0x7FFFFFFFL);

    /* "Warm up" the generator for a number of rounds to eliminate
       any residual inflence of the seed. */

    for (i = 0; i < warmup; i++) {
	updgen(1);
    }

    /* Fill the shuffle table with values.  */

    for (i = 0; i < shuffleSize; i++) {
	updgen(1);
	shuffle[(shuffleSize - 1) - i] = gen1;
    }
    state = shuffle[0];
}

/*  LEnextByte	--  Get next byte from generator.  */

unsigned char LEnextByte(void)
{
    int i;

    updgen(1);
    updgen(2);

    /* Extract shuffle table index from most significant part
       of the previous result. */

    i = state / (1 + (((int32) mod1) - 1) / shuffleSize);

    /* New state is sum of generators modulo one of their moduli.  */

    state = (int32) (((shuffle[i]) + ((unsigned int32) gen2)) % mod1);

    /* Replace value in shuffle table with generator 1 result.	*/

    shuffle[i] = gen1;

    return (unsigned char) (state / (1 + (((int32) mod1) - 1) / 256));
}