COPASI API  4.16.103
Cr250.cpp
Go to the documentation of this file.
1 /* Begin CVS Header
2  $Source: /Volumes/Home/Users/shoops/cvs/copasi_dev/copasi/randomGenerator/Cr250.cpp,v $
3  $Revision: 1.9 $
4  $Name: $
5  $Author: shoops $
6  $Date: 2006/04/27 01:31:00 $
7  End CVS Header */
8 
9 // Copyright 2005 by Pedro Mendes, Virginia Tech Intellectual
10 // Properties, Inc. and EML Research, gGmbH.
11 // All rights reserved.
12 
13 /******************************************************************************
14  * Module: Cr250.cpp
15  * Description: implements R250 random number generator, from S.
16  * Kirkpatrick and E. Stoll, Journal of Computational Physics, 40,
17  * p. 517 (1981).
18  * Written by: W. L. Maier
19  *
20  * METHOD...
21  * 16 parallel copies of a linear shift register with
22  * period 2^250 - 1. FAR longer period than the usual
23  * linear congruent generator, and commonly faster as
24  * well. (For details see the above paper, and the
25  * article in DDJ referenced below.)
26  *
27  * HISTORY...
28  * Sep 02: Converted to c++ class by Stefan Hoops
29  * Sep 92: Number returned by dr250() is in range [0,1) instead
30  * of [0,1], so for example a random angle in the
31  * interval [0, 2*PI) can be calculated
32  * conveniently. (J. R. Van Zandt <jrv@mbunix.mitre.org>)
33  * Aug 92: Initialization is optional. Default condition is
34  * equivalent to initializing with the seed 12345,
35  * so that the sequence of random numbers begins:
36  * 1173, 53403, 52352, 35341... (9996 numbers
37  * skipped) ...57769, 14511, 46930, 11942, 7978,
38  * 56163, 46506, 45768, 21162, 43113... Using ^=
39  * operator rather than two separate statements.
40  * Initializing with own linear congruent
41  * pseudorandom number generator for portability.
42  * Function prototypes moved to a header file.
43  * Implemented r250n() to generate numbers
44  * uniformly distributed in a specific range
45  * [0,n), because the common expedient rand()%n is
46  * incorrect. (J. R. Van Zandt <jrv@mbunix.mitre.org>)
47  * May 91: Published by W. L. Maier, "A Fast Pseudo Random Number
48  * Generator", Dr. Dobb's Journal #176.
49  ******************************************************************************/
50 
51 #include "copasi.h"
52 #include "CRandom.h"
53 
54 Cr250::Cr250(unsigned C_INT32 seed):
55  CRandom(),
56  mIndex(0)
57 {
58  setModulus(65535);
59  initialize(seed);
60 }
62 
63 void Cr250::initialize(unsigned C_INT32 seed)
64 {
65  /*--------------------------------------------------------------------------*/
66  int j, k;
67  unsigned int mask;
68  unsigned int msb;
69  /*--------------------------------------------------------------------------*/
70 
71  mIndex = 0;
72  mSeed = seed;
73 
74  /* Fill the r250 buffer with 15-bit values */
75  for (j = 0; j < 250; j++)
76  mBuffer[j] = myrand();
77 
78  /* Set some of the MS bits to 1 */
79  for (j = 0; j < 250; j++)
80  if (myrand() > 16384)
81  mBuffer[j] |= 0x8000;
82 
83  msb = 0x8000; /* To turn on the diagonal bit */
84  mask = 0xffff; /* To turn off the leftmost bits */
85 
86  for (j = 0; j < 16; j++)
87  {
88  k = 11 * j + 3; /* Select a word to operate on */
89  mBuffer[k] &= mask; /* Turn off bits left of the diagonal */
90  mBuffer[k] |= msb; /* Turn on the diagonal bit */
91  mask >>= 1;
92  msb >>= 1;
93  }
94  return;
95 }
96 
98 {return r250();}
99 
101 {return r250();}
102 
104 {return r250() * mModulusInv;}
105 
107 {return dr250();}
108 
110 {return (r250() + .5) * mModulusInv1;}
111 
112 unsigned C_INT32 Cr250::r250(void)
113 {
114  register C_INT16 j;
115 
116  if (mIndex > 146)
117  j = mIndex - 147; /* Wrap pointer around */
118  else
119  j = mIndex + 103;
120 
121  mNumberU = mBuffer[mIndex] ^= mBuffer[j];
122 
123  if (mIndex > 248) /* Increment pointer for next time */
124  mIndex = 0;
125  else
126  mIndex++;
127 
128  return mNumberU;
129 }
130 
131 unsigned C_INT32 Cr250::r250n(const unsigned C_INT16 & max)
132 {
133  register unsigned C_INT16 limit;
134 
135  limit = (65535U / max) * max;
136 
137  do
138  {
139  r250();
140  r250(); // Why a second call?
141  }
142  while (mNumberU >= limit);
143 
144  return mNumberU % max;
145 }
146 
148 {
149  mNumberU = r250();
150  return mNumberU / 65536.; /* Return a number in [0.0 to 1.0) */
151 }
152 
154 {
155  mSeed = mSeed * 0x015a4e35L + 1;
156  return (mSeed >> 16)&0x7fff;
157 }
C_FLOAT64 getRandomCO()
Definition: Cr250.cpp:106
unsigned C_INT32 mSeed
Definition: Cr250.h:32
void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
Definition: Cr250.cpp:63
unsigned C_INT32 r250n(const unsigned C_INT16 &max)
Definition: Cr250.cpp:131
C_INT32 getRandomS()
Definition: Cr250.cpp:100
C_FLOAT64 getRandomCC()
Definition: Cr250.cpp:103
C_FLOAT64 getRandomOO()
Definition: Cr250.cpp:109
#define C_INT32
Definition: copasi.h:90
~Cr250()
Definition: Cr250.cpp:61
C_FLOAT64 dr250(void)
Definition: Cr250.cpp:147
unsigned C_INT32 r250(void)
Definition: Cr250.cpp:112
Cr250(unsigned C_INT32 seed)
Definition: Cr250.cpp:54
C_FLOAT64 mModulusInv1
Definition: CRandom.h:108
C_FLOAT64 mModulusInv
Definition: CRandom.h:103
void setModulus(const unsigned C_INT32 &modulus)
Definition: CRandom.cpp:129
#define C_INT16
Definition: copasi.h:91
C_INT32 mIndex
Definition: Cr250.h:30
unsigned C_INT16 mBuffer[250]
Definition: Cr250.h:34
unsigned C_INT16 myrand()
Definition: Cr250.cpp:153
#define C_FLOAT64
Definition: copasi.h:92
unsigned C_INT32 mNumberU
Definition: CRandom.h:78
unsigned C_INT32 getRandomU()
Definition: Cr250.cpp:97
#define max(a, b)
Definition: f2c.h:176