about summary refs log tree commit diff
path: root/lib/util/rand.c
blob: 6a0a2cdb660b273e4df20aac1b275a44d98ce0c9 (plain) (blame)
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
/*

Pseudo-random number generator for Netpbm

The interface provided herein should be flexible enough for anybody
who wishes to use some other random number generator.

---

If you desire to implement a different generator, or writing an original
one, first take a look at the random number generator section of the
GNU Scientific Library package (GSL).

GNU Scientific Library
https://www.gnu.org/software/gsl/

GSL Random Number Generators
https://wnww.gnu.org/software/gsl/doc/html/rng.html

*/

#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>
#include <strings.h>
#include <time.h>
#include <float.h>
#include <math.h>

#include "netpbm/pm_c_util.h"
#include "netpbm/mallocvar.h"
#include "netpbm/pm.h"
#include "netpbm/rand.h"

/*-----------------------------------------------------------------------------
                              Use
-------------------------------------------------------------------------------
  Typical usage:

      #include "rand.h"

      ...

      myfunction( ... , unsigned int const seed , ... ) {

          struct randSt;

          ...

          pm_randinit(&randSt);
          pm_srand(&randSt, seed);  // pm_srand2() is often more useful

          ...

          pm_rand(&randSt);

          ...

          pm_randterm(&randSt);

      }
-----------------------------------------------------------------------------*/



/*-----------------------------------------------------------------------------
                            Design note
-------------------------------------------------------------------------------

  Netpbm code contains multiple random number generators.  Stock Netpbm always
  uses an internal pseudo-random number generator that implements the Mersenne
  Twister method and does not rely on any randomness facility of the operating
  system, but it is easy to compile an alternative version that uses others.

  The Mersenne Twister method was new to Netpbm in Netpbm 10.94 (March 2021).
  Before that, Netpbm used standard OS-provided facilities.

  Programs that use random numbers have existed in Netpbm since PBMPlus days.
  The system rand() function was used in instances randomness was required;
  exceptions were rare and all of them appear to be errors on the part of the
  original author.

  Although the rand() function is available in every system on which Netpbm
  runs, differences exist in the underlying algorithm, so that Netpbm programs
  produce different output on different systems even when the user specifies
  the same random number seed.

  This was not considered a problem in the early days.  Deterministic
  operation was not a feature users requested and it was impossible regardless
  of the random number generation method on most programs because they did not
  allow a user to specify a seed for the generator.

  This state of affairs changed as Netpbm got firmly established as a
  base-level system package.  Security became critical for many users.  A
  crucial component of quality control is automated regression tests (="make
  check").  Unpredictable behavior gets in the way of testing.  One by one
  programs were given the -randomseed (or -seed) option to ensure reproducible
  results.  Often this was done as new test cases were written.  However,
  inconsistent output caused by system-level differences in rand()
  implementation remained a major obstacle.

  In 2020 the decision was made to replace all calls to rand() in the Netpbm
  source code with an internal random number generator.  We decided to use the
  Mersenne Twister, which is concise, enjoys a fine reputation and is
  available under liberal conditions (see below.)
-----------------------------------------------------------------------------*/


void
pm_srand(struct pm_randSt * const randStP,
         unsigned int       const seed) {
/*----------------------------------------------------------------------------
  Initialize (or "seed") the random number generation sequence with value
  'seed'.
-----------------------------------------------------------------------------*/
    pm_randinit(randStP);

    randStP->vtable.srand(randStP, seed);

    randStP->seed = seed;
}



void
pm_srand2(struct pm_randSt * const randStP,
          bool               const seedValid,
          unsigned int       const seed) {
/*----------------------------------------------------------------------------
  Seed the random number generator.  If 'seedValid' is true, use 'seed"..
  Otherwise, use pm_randseed().

  For historical reasons pm_randseed() is defined in libpm.c rather than
  this source file.
-----------------------------------------------------------------------------*/
    pm_srand(randStP, seedValid ? seed : pm_randseed() );

}



unsigned long int
pm_rand(struct pm_randSt * const randStP) {
/*----------------------------------------------------------------------------
  An integer random number in the interval [0, randStP->max].
-----------------------------------------------------------------------------*/
    return randStP->vtable.rand(randStP);
}



double
pm_drand(struct pm_randSt * const randStP) {
/*----------------------------------------------------------------------------
  A floating point random number in the interval [0, 1).

  Although the return value is declared as double, the actual value will have
  no more precision than a single call to pm_rand() provides.  This is 32 bits
  for Mersenne Twister.
-----------------------------------------------------------------------------*/
    return (double) pm_rand(randStP) / randStP->max;
}



void
pm_gaussrand2(struct pm_randSt * const randStP,
              double *           const r1P,
              double *           const r2P) {
/*----------------------------------------------------------------------------
  Generate two Gaussian (or normally) distributed random numbers *r1P and
  *r2P.

  Mean = 0, Standard deviation = 1.

  This is called the Box-Muller method.

  For details of this algorithm and other methods for producing
  Gaussian random numbers see:

  http://www.doc.ic.ac.uk/~wl/papers/07/csur07dt.pdf
-----------------------------------------------------------------------------*/
    double u1, u2;

    u1 = pm_drand(randStP);
    u2 = pm_drand(randStP);

    if (u1 < DBL_EPSILON)
        u1 = DBL_EPSILON;

    *r1P = sqrt(-2.0 * log(u1)) * cos(2.0 * M_PI * u2);
    *r2P = sqrt(-2.0 * log(u1)) * sin(2.0 * M_PI * u2);
}



double
pm_gaussrand(struct pm_randSt * const randStP) {
/*----------------------------------------------------------------------------
  A Gaussian (or normally) distributed random number.

  Mean = 0, Standard deviation = 1.

  If a randStP->gaussCache has a value, return that value.  Otherwise call
  pm_gaussrand2; return one generated value, remember the other.
-----------------------------------------------------------------------------*/
    double retval;

    if (!randStP->gaussCacheValid) {
        pm_gaussrand2(randStP, &retval, &randStP->gaussCache);
        randStP->gaussCacheValid = true;
    } else {
        retval = randStP->gaussCache;
        randStP->gaussCacheValid = false;
    }

    return retval;
}



uint32_t
pm_rand32(struct pm_randSt * const randStP) {
/*-----------------------------------------------------------------------------
  Generate a 32-bit random number.

  This is a provision for users who select a non-default random number
  generator which returns less than 32 bits per call.  Many system generators
  are known to return 31 bits (max = 2147483647 or 0x7FFFFFFF).

  This does not work with generators that return less than 11 bits per call.
  The least we know of is the archaic RANDU, which generates 15 bits (max =
  32767 or 0x7FFF).
-----------------------------------------------------------------------------*/
    unsigned int const randMax = randStP->max;

    uint32_t retval;

    if (randMax >= 0xFFFFFFFF)
        retval = pm_rand(randStP);
    else {
        uint32_t scale;

        retval = 0;  /* Initial value */

        for (scale = 0xFFFFFFFF; scale > 0; scale /= (randMax +1))
            retval *= (randMax + 1) + pm_rand(randStP);
    }

    return retval;;
}



void
pm_randinit(struct pm_randSt * const randStP) {
/*----------------------------------------------------------------------------
  Initialize the random number generator.
-----------------------------------------------------------------------------*/
    switch (PM_RANDOM_NUMBER_GENERATOR) {
    case PM_RAND_SYS_RAND:
        randStP->vtable = pm_randsysrand_vtable;
        break;
    case PM_RAND_SYS_RANDOM:
        randStP->vtable = pm_randsysrandom_vtable;
        break;
    case PM_RAND_MERSENNETWISTER:
        randStP->vtable = pm_randmersenne_vtable;
        break;
    default:
        pm_error("INTERNAL ERROR: Invalid value of "
                 "PM_RANDOM_NUMBER_GENERATOR (random number generator "
                 "engine type): %u", PM_RANDOM_NUMBER_GENERATOR);
    }

    randStP->vtable.init(randStP);

    randStP->gaussCacheValid = false;
}



void
pm_randterm(struct pm_randSt * const randStP) {
/*----------------------------------------------------------------------------
  Tear down the random number generator.
-----------------------------------------------------------------------------*/
     if (randStP->stateP)
         free(randStP->stateP);
}