source: trunk/optimizer/src/pgapack/pgapack/source/random.c @ 816

Last change on this file since 816 was 816, checked in by liveletlive, 15 years ago

Committing the newer version of PGAPack.

File size: 12.7 KB
Line 
1/*
2 * 
3 *  *********************************************************************
4 *  (C) COPYRIGHT 1995 UNIVERSITY OF CHICAGO
5 *  *********************************************************************
6 * 
7 *  This software was authored by
8 * 
9 *  D. Levine
10 *  Mathematics and Computer Science Division Argonne National Laboratory
11 *  Argonne IL 60439
12 *  levine@mcs.anl.gov
13 *  (708) 252-6735
14 *  (708) 252-5986 (FAX)
15 * 
16 *  with programming assistance of participants in Argonne National
17 *  Laboratory's SERS program.
18 * 
19 *  This program contains material protectable under copyright laws of the
20 *  United States.  Permission is hereby granted to use it, reproduce it,
21 *  to translate it into another language, and to redistribute it to
22 *  others at no charge except a fee for transferring a copy, provided
23 *  that you conspicuously and appropriately publish on each copy the
24 *  University of Chicago's copyright notice, and the disclaimer of
25 *  warranty and Government license included below.  Further, permission
26 *  is hereby granted, subject to the same provisions, to modify a copy or
27 *  copies or any portion of it, and to distribute to others at no charge
28 *  materials containing or derived from the material.
29 * 
30 *  The developers of the software ask that you acknowledge its use in any
31 *  document referencing work based on the  program, such as published
32 *  research.  Also, they ask that you supply to Argonne National
33 *  Laboratory a copy of any published research referencing work based on
34 *  the software.
35 * 
36 *  Any entity desiring permission for further use must contact:
37 * 
38 *  J. Gleeson
39 *  Industrial Technology Development Center Argonne National Laboratory
40 *  Argonne IL 60439
41 *  gleesonj@smtplink.eid.anl.gov
42 *  (708) 252-6055
43 * 
44 *  ********************************************************************
45 *  DISCLAIMER
46 * 
47 *  THIS PROGRAM WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED BY AN AGENCY
48 *  OF THE UNITED STATES GOVERNMENT.  NEITHER THE UNIVERSITY OF CHICAGO,
49 *  THE UNITED STATES GOVERNMENT NOR ANY OF THEIR EMPLOYEES MAKE ANY
50 *  WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR
51 *  RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR USEFULNESS OF ANY
52 *  INFORMATION OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT
53 *  INFRINGE PRIVATELY OWNED RIGHTS.
54 * 
55 *  **********************************************************************
56 *  GOVERNMENT LICENSE
57 * 
58 *  The Government is granted for itself and others acting on its behalf a
59 *  paid-up, non-exclusive, irrevocable worldwide license in this computer
60 *  software to reproduce, prepare derivative works, and perform publicly
61 *  and display publicly.
62 */
63
64/*****************************************************************************
65*     FILE: random.c: This file contains routines to generate randomness.
66*
67*     Authors: David M. Levine, Philip L. Hallstrom, David M. Noelle,
68*              Brian P. Walenz
69*****************************************************************************/
70
71#include "pgapack.h"
72
73/*U****************************************************************************
74   PGARandomFlip - flip a biased coin and return PGA_TRUE if the coin is
75   a "winner."  Otherwise, return PGA_FALSE.
76
77   Category: Utility
78
79   Inputs:
80      ctx - context variable
81      p   - biased probability (.5 is a fair coin)
82
83   Outputs:
84      PGA_TRUE or PGA_FALSE
85
86   Example:
87      To return PGA_TRUE approximately seventy percent of the time, use
88
89      PGAContext *ctx;
90      int p;
91      :
92      PGARandomFlip(ctx, 0.7)
93
94****************************************************************************U*/
95int PGARandomFlip ( PGAContext *ctx, double p )
96{
97    PGADebugEntered("PGARandomFlip");
98
99    PGADebugExited("PGARandomFlip");
100
101    return( (PGARandom01(ctx, 0) < p) ? PGA_TRUE : PGA_FALSE);
102}
103
104
105/*U****************************************************************************
106   PGARandomInterval - returns a uniform random number on the specified
107   interval
108
109   Category: Utility
110
111   Inputs:
112      ctx - context variable
113      start - starting (integer) value of the interval
114      end   - ending   (integer) value of the interval
115
116   Outputs:
117      A uniformly distributed random number in the interval [start, end].
118
119   Example:
120      Generate a value uniformly random from the interval [0,99]
121
122      PGAContext *ctx;
123      :
124      PGARandomInterval(ctx, 0, 99);
125
126****************************************************************************U*/
127int PGARandomInterval( PGAContext *ctx, int start, int end)
128{
129    PGADebugEntered("PGARandomInterval");
130   
131    PGADebugExited("PGARandomInterval");
132
133    return( (int)floor(PGARandom01(ctx, 0) * (double)(end-start+1) ) + start );
134   
135/*
136   The original call...
137
138   return(ceil( (double)(end-start+1) * PGARandom01(ctx, 0) )+start-1);
139*/
140}
141
142
143/*****************************************************************************
144*  This is a C language implementation of the universal random number        *
145*  generator proposed by G. Marsaglia and A. Zaman and translated from       *
146*  F. James' version.                                                        *
147*                                                                            *
148*  F. James                                                                  *
149*  A review of pseudorandom number generators                                *
150*  Computer Physics Communication                                            *
151*  60 (1990) 329-344                                                         *
152*                                                                            *
153*  G. Marsaglia, A. Zaman, W. Tseng                                          *
154*  Stat Prob. Letter                                                         *
155*  9 (1990) 35.                                                              *
156*                                                                            *
157*  G. Marsaglia, A. Zaman                                                    *
158*  FSU-SCRI-87-50                                                            *
159*                                                                            *
160*  This algorithm is a combination of a lagged Fibonacci and arithmetic      *
161*  sequence (F. James) generator with period of 2^144.  It provides 32-bit   *
162*  floating point numbers in the range from zero to one.  It is claimed to   *
163*  be portable and provides bit-identical results on all machines with at    *
164*  least 24-bit mantissas.                                                   *
165*                                                                            *
166*  PGARandom01 should be initialized with a 32-bit integer seed such that    *
167*  0 <= seed <= 900,000,000.  Each of these 900,000,000 values gives rise    *
168*  to an independent sequence of ~ 10^30.                                    *
169*                                                                            *
170*  warning on use of static storage class on thread shared memory machines   *
171*****************************************************************************/
172/*U****************************************************************************
173   PGARandom01 - generates a uniform random number on the interval [0,1)
174   If the second argument is 0 it returns the next random number in the
175   sequence.  Otherwise, the second argument is used as a new seed for the
176   population
177
178   Category: Utility
179
180   Inputs:
181      ctx     - context variable
182      newseed - either 0 to get the next random number, or nonzero
183                to reseed
184   Outputs:
185      A random number on the interval [0,1)
186
187   Example:
188      To get the next random number use
189
190      PGAContext *ctx;
191      double r;
192      :
193      r = PGARandom01(ctx,0);
194
195****************************************************************************U*/
196double PGARandom01( PGAContext *ctx, int newseed )
197{
198
199    /* initialization variables */
200    int ij, kl, i, j, k, l, m, ii, jj;
201    float s, t;
202
203    /* random number variables */
204    static int seed=1;      /* default seed if none specified */
205    static int i96, j96;
206    static float u[97], uni, c, cd, cm;
207
208
209    PGADebugEntered("PGARandom01");
210
211    /* initialization */
212/*     printf("i96 = %d\tj96 = %d\n", i96, j96); */
213
214    if ( newseed != 0 ) {
215
216        seed = newseed % 900000000;
217        ij   = seed / 30082;
218        kl   = seed - 30082 * ij;
219        i    = ( (ij/177) % 177 ) + 2;
220        j    = (  ij      % 177 ) + 2;
221        k    = ( (kl/169) % 178 ) + 1;
222        l    = (  kl      % 169 );
223
224        for ( ii=0; ii<97; ii++ ) {
225
226            s = 0.0;
227            t = 0.5;
228
229            for ( jj=0; jj<24; jj++ ) {
230
231                m = ( ((i*j) % 179) * k ) % 179;
232                i = j;
233                j = k;
234                k = m;
235                l = ( (53*l) + 1 ) % 169;
236                if ( ( (l*m) % 64 ) >= 32 )
237                    s += t;
238                t *= .5;
239            }
240
241            u[ii] = s;
242        }
243
244        c   = 362436.  /16777216.;
245        cd  = 7654321. /16777216.;
246        cm  = 16777213./16777216.;
247        i96 = 96;
248        j96 = 32;
249    }
250
251    /* random number generation */
252    uni = u[i96] - u[j96];
253    if ( uni < 0. ) uni += 1.0;
254    u[i96] = uni;
255    i96--;
256    if ( i96 < 0  ) i96 = 96;
257    j96--;
258    if ( j96 < 0  ) j96 = 96;
259    c   -= cd;
260    if ( c   < 0. ) c += cm;
261    uni -= c;
262    if ( uni < 0. ) uni += 1.0;
263
264    PGADebugExited("PGARandom01");
265    return( (double) uni);
266}
267
268/*U****************************************************************************
269   PGARandomUniform - returns a uniform random number on the interval
270   [start,end]
271
272   Category: Utility
273
274   Inputs:
275      ctx - context variable
276      start - starting (double) value of the interval
277      end   - ending   (double) value of the interval
278
279   Outputs:
280      A random number on the interval [start,end]
281
282   Example:
283      Generate a uniform random number on the interval [-0.5, 1.5]
284
285      PGAContext *ctx;
286      double r;
287      :
288      r = PGARandomUniform(ctx, -0.5, 1.5);
289
290****************************************************************************U*/
291double PGARandomUniform( PGAContext *ctx, double start, double end)
292{
293    double val, r;
294
295    PGADebugEntered("PGARandomUniform");
296
297    r = PGARandom01(ctx, 0);
298    val = (end-start) * r + start;
299
300    PGADebugExited("PGARandomUniform");
301
302    return(val);
303}
304
305
306/*U****************************************************************************
307   PGARandomGaussian - returns an approximation to a Gaussian random number
308
309   Category: Utility
310
311   Inputs:
312       mean  - the mean of the Gaussian distribution
313       sigma - the standard deviation of the Gaussian distribution
314
315   Outputs:
316      A random number selected from a Gaussian distribution with given
317      mean and standard deviation
318
319   Example:
320      To generate a Gaussian random number with mean 0.0 and standard
321      deviation 1.0 use
322
323      PGAContext *ctx;
324      :
325      r = PGARandomGaussian(ctx, 0.0, 1.0);
326
327****************************************************************************U*/
328double PGARandomGaussian( PGAContext *ctx, double mean, double sigma)
329{
330    int i;
331    double sum = 0.;
332
333    PGADebugEntered("PGARandomGaussian");
334
335    for (i=11;i>=0; i--)
336        sum += PGARandom01(ctx, 0);
337
338    PGADebugExited("PGARandomGaussian");
339
340    return ( (sum-6.0) * sigma + mean );
341}
342
343/*U***************************************************************************
344   PGAGetRandomSeed - returns the integer to seed random numbers with
345
346   Category: Utility
347
348   Inputs:
349      ctx - context variable
350
351   Outputs:
352      The seed for the random number generator
353
354   Example:
355      PGAContext *ctx;
356      int seed;
357      :
358      seed = PGAGetRandomSeed(ctx);
359
360***************************************************************************U*/
361int PGAGetRandomSeed(PGAContext *ctx)
362{
363    PGADebugEntered("PGAGetRandomSeed");
364
365    PGADebugExited("PGAGetRandomSeed");
366
367    return(ctx->init.RandomSeed);
368}
369
370/*U****************************************************************************
371   PGASetRandomSeed - set a seed for the random number generator.  The
372   default is to use a random seed.  Specifying a seed exlicitly allows
373   for reproducibility of runs.
374
375   Category: Utility
376
377   Inputs:
378      ctx  - context variable
379      seed - seed  for the random number generator
380
381   Outputs:
382      None
383
384   Example:
385      PGAContext *ctx;
386      :
387      PGASetRandomSeed(ctx,1);
388
389****************************************************************************U*/
390void PGASetRandomSeed(PGAContext *ctx, int seed)
391{
392#define MAX_PROCESSORS 2048
393
394    PGADebugEntered("PGASetRandomSeed");
395    PGAFailIfSetUp("PGASetRandomSeed");
396
397    if ((seed < 1) || (seed + MAX_PROCESSORS > 900000000))
398        PGAError ( ctx, "PGASetRandomSeed: Invalid value of seed:",
399                  PGA_FATAL, PGA_INT, (void *) &seed);
400    else
401        ctx->init.RandomSeed = seed;
402   
403    PGADebugExited("PGASetRandomSeed");
404}
Note: See TracBrowser for help on using the repository browser.