source: vis_dev/cusp-1.1/src/util/random.c @ 103

Last change on this file since 103 was 12, checked in by cecile, 13 years ago

cusp added
File size: 7.5 KB
RevLine 
[12]1/**CFile***********************************************************************
2
3  FileName    [random.c]
4
5  PackageName [util]
6
7  Synopsis    [Our own portable random number generator.]
8
9  Author      [Fabio Somenzi <Fabio@Colorado.EDU>]
10
11  Copyright   [Copyright (c) 1994-1996 The Regents of the Univ. of California.
12  All rights reserved.
13
14  Permission is hereby granted, without written agreement and without license
15  or royalty fees, to use, copy, modify, and distribute this software and its
16  documentation for any purpose, provided that the above copyright notice and
17  the following two paragraphs appear in all copies of this software.
18
19  IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY FOR
20  DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT
21  OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF
22  CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
23
24  THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
25  INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
26  FITNESS FOR A PARTICULAR PURPOSE.  THE SOFTWARE PROVIDED HEREUNDER IS ON AN
27  \"AS IS\" BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATION TO PROVIDE
28  MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.]
29
30******************************************************************************/
31
32#include "util.h"
33
34static char rcsid[] UNUSED = "$Id: random.c,v 1.1.1.1 2008-11-14 20:40:10 hhkim Exp $";
35
36/*---------------------------------------------------------------------------*/
37/* Constant declarations                                                     */
38/*---------------------------------------------------------------------------*/
39
40/* Random generator constants. */
41#define MODULUS1 2147483563
42#define LEQA1 40014
43#define LEQQ1 53668
44#define LEQR1 12211
45#define MODULUS2 2147483399
46#define LEQA2 40692
47#define LEQQ2 52774
48#define LEQR2 3791
49#define STAB_SIZE 64
50#define STAB_DIV (1 + (MODULUS1 - 1) / STAB_SIZE)
51
52/*---------------------------------------------------------------------------*/
53/* Type declarations                                                         */
54/*---------------------------------------------------------------------------*/
55
56
57/*---------------------------------------------------------------------------*/
58/* Structure declarations                                                    */
59/*---------------------------------------------------------------------------*/
60
61
62/*---------------------------------------------------------------------------*/
63/* Variable declarations                                                     */
64/*---------------------------------------------------------------------------*/
65static long utilRand = 0;
66static long utilRand2;
67static long shuffleSelect;
68static long shuffleTable[STAB_SIZE];
69
70/*---------------------------------------------------------------------------*/
71/* Macro declarations                                                        */
72/*---------------------------------------------------------------------------*/
73
74
75/**AutomaticStart*************************************************************/
76
77/*---------------------------------------------------------------------------*/
78/* Static function prototypes                                                */
79/*---------------------------------------------------------------------------*/
80
81/**AutomaticEnd***************************************************************/
82
83
84/*---------------------------------------------------------------------------*/
85/* Definition of exported functions                                          */
86/*---------------------------------------------------------------------------*/
87
88
89/**Function********************************************************************
90
91  Synopsis    [Initializer for the portable random number generator.]
92
93  Description [Initializer for the portable number generator based on
94  ran2 in "Numerical Recipes in C." The input is the seed for the
95  generator. If it is negative, its absolute value is taken as seed.
96  If it is 0, then 1 is taken as seed. The initialized sets up the two
97  recurrences used to generate a long-period stream, and sets up the
98  shuffle table.]
99
100  SideEffects [None]
101
102  SeeAlso     [util_random]
103
104******************************************************************************/
105void
106util_srandom(long seed)
107{
108    int i;
109
110    if (seed < 0)       utilRand = -seed;
111    else if (seed == 0) utilRand = 1;
112    else                utilRand = seed;
113    utilRand2 = utilRand;
114    /* Load the shuffle table (after 11 warm-ups). */
115    for (i = 0; i < STAB_SIZE + 11; i++) {
116        long int w;
117        w = utilRand / LEQQ1;
118        utilRand = LEQA1 * (utilRand - w * LEQQ1) - w * LEQR1;
119        utilRand += (utilRand < 0) * MODULUS1;
120        shuffleTable[i % STAB_SIZE] = utilRand;
121    }
122    shuffleSelect = shuffleTable[1 % STAB_SIZE];
123} /* end of util_srandom */
124
125/**Function********************************************************************
126
127  Synopsis    [Portable random number generator.]
128
129  Description [Portable number generator based on ran2 from "Numerical
130  Recipes in C." It is a long period (> 2 * 10^18) random number generator
131  of L'Ecuyer with Bays-Durham shuffle. Returns a long integer uniformly
132  distributed between 0 and 2147483561 (inclusive of the endpoint values).
133  The random generator can be explicitly initialized by calling
134  util_srandom. If no explicit initialization is performed, then the
135  seed 1 is assumed.]
136
137  SideEffects []
138
139  SeeAlso     [util_srandom]
140
141******************************************************************************/
142long
143util_random(void)
144{
145    int i;      /* index in the shuffle table */
146    long int w; /* work variable */
147
148    /* utilRand == 0 if the geneartor has not been initialized yet. */
149    if (utilRand == 0) util_srandom(1);
150
151    /* Compute utilRand = (utilRand * LEQA1) % MODULUS1 avoiding
152    ** overflows by Schrage's method.
153    */
154    w          = utilRand / LEQQ1;
155    utilRand   = LEQA1 * (utilRand - w * LEQQ1) - w * LEQR1;
156    utilRand  += (utilRand < 0) * MODULUS1;
157
158    /* Compute utilRand2 = (utilRand2 * LEQA2) % MODULUS2 avoiding
159    ** overflows by Schrage's method.
160    */
161    w          = utilRand2 / LEQQ2;
162    utilRand2  = LEQA2 * (utilRand2 - w * LEQQ2) - w * LEQR2;
163    utilRand2 += (utilRand2 < 0) * MODULUS2;
164
165    /* utilRand is shuffled with the Bays-Durham algorithm.
166    ** shuffleSelect and utilRand2 are combined to generate the output.
167    */
168
169    /* Pick one element from the shuffle table; "i" will be in the range
170    ** from 0 to STAB_SIZE-1.
171    */
172    i = shuffleSelect / STAB_DIV;
173    /* Mix the element of the shuffle table with the current iterate of
174    ** the second sub-generator, and replace the chosen element of the
175    ** shuffle table with the current iterate of the first sub-generator.
176    */
177    shuffleSelect   = shuffleTable[i] - utilRand2;
178    shuffleTable[i] = utilRand;
179    shuffleSelect  += (shuffleSelect < 1) * (MODULUS1 - 1);
180    /* Since shuffleSelect != 0, and we want to be able to return 0,
181    ** here we subtract 1 before returning.
182    */
183    return(shuffleSelect - 1);
184
185} /* end of util_random */
186
187
188/*---------------------------------------------------------------------------*/
189/* Definition of internal functions                                          */
190/*---------------------------------------------------------------------------*/
191
192
193/*---------------------------------------------------------------------------*/
194/* Definition of static functions                                            */
195/*---------------------------------------------------------------------------*/
196
197
198
Note: See TracBrowser for help on using the repository browser.