CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

testBug58950.cc
Go to the documentation of this file.
1 // ----------------------------------------------------------------------
2 //
3 // testBug58950 -- test problem with RanecuEngine on 64bit machines
4 //
5 // R. Weller 11/11/09 initial test from Vanderbilt
6 // L. Garren 12/1/09 rewritten for test suite
7 //
8 // ----------------------------------------------------------------------
9 #include <iostream>
10 #include <stdexcept>
11 #include <cmath>
12 #include <stdlib.h>
13 #include "CLHEP/Random/RanecuEngine.h"
14 #include "CLHEP/Random/Random.h"
15 #include "pretend.h"
16 
17 // Absolutely Safe Equals Without Registers Screwing Us Up
18 bool equals01(const std::vector<double> &ab) {
19  return ab[1]==ab[0];
20 }
21 bool equals(double a, double b) {
22  std::vector<double> ab(2);
23  ab[0]=a; ab[1]=b;
24  return (equals01(ab));
25 }
26 
27 bool printCheck( int & i, double & r, std::ofstream & os )
28 {
29  os << i << " " << r << std::endl;
30  if (r < 0 || r > 1.0 ) {
31  std::cout << "Error: bad random number " << r << std::endl;
32  return false;
33  }
34  return true;
35 }
36 
37 int main() {
38 
39  std::ofstream output("testBug58950.cout");
40 
41  output << std::endl << "short " << sizeof(short) << std::endl;
42  output << "int " << sizeof(int) << std::endl;
43  output << "unsigned int " << sizeof(unsigned int) << std::endl;
44  output << "long " << sizeof(long) << std::endl;
45  output << "float " << sizeof(float) << std::endl;
46  output << "double " << sizeof(double) << std::endl;
47  output << "long double " << sizeof(long double) << std::endl << std::endl;
48 
53 
54  long rvals[2];
55  try {
56  std::ifstream in("/dev/urandom", std::ios::in | std::ios::binary);
57  if(in.is_open()) {
58  in.read((char *)(&rvals), 2*sizeof(long));
59  in.close();
60  if(in.fail()) {
61  throw std::runtime_error("File read error");
62  }
63  } else throw std::runtime_error("File open error");
64  } catch(std::runtime_error e) {
65  std::ostringstream dStr;
66  dStr << "Error: " << e.what()
67  << " processing seed from file \"" << "/dev/urandom" << "\".";
68  throw std::runtime_error(dStr.str().c_str());
69  }
70 
71  int nNumbers = 20;
72  int badcount = 0;
73 
74  long seeds[3];
75  const long *pseeds;
76  //***********************************************************************
77  // Seeds are expected to be positive. Therefore, if either seed is
78  // negative then prior to 2.0.4.5 the generator set initial conditions
79  // and generated the same sequence of numbers no matter what the seeds were.
80  seeds[0]=rvals[0];
81  seeds[1]=rvals[1];
82  seeds[2]=0;
83  if( rvals[0] > 0 ) seeds[0] = -rvals[0];
84 
85  double negseq[20] = { 0.154707, 0.587114, 0.702059, 0.566, 0.988325,
86  0.525921, 0.191554, 0.269338, 0.234277, 0.358997,
87  0.549936, 0.296877, 0.162243, 0.227732, 0.528862,
88  0.631571, 0.176462, 0.247858, 0.170025, 0.284483 };
89  double eps = 1.0E-6;
90 
91  output << std::endl << "********************" << std::endl;
92  output << "This is the case that may or may not fail." << std::endl;
93  output << "However, if it has values in (0,1), they are a " << std::endl
94  << "deterministic sequence beginning with 0.154707." << std::endl;
95  output << "seeds[0] = " << seeds[0] << "\n"
96  << "seeds[1] = " << seeds[1] << std::endl << std::endl;
97 
98  g->setTheSeeds(seeds);
99  int rseq = 0;
100  for (int i=0; i < nNumbers; ++i) {
101  double r = g->flat();
102  if( ! printCheck(i,r,output) ) ++badcount;
103  // before the change, the random number sequence was reliably the same
104  if( std::fabs(r-negseq[i]) < eps ) {
105  std::cout << " reproducing sequence " << i << " "
106  << r << " " << negseq[i] << std::endl;
107  ++rseq;
108  }
109  }
110  if( rseq == 20 ) ++badcount;
111  pseeds=g->getTheSeeds();
112  output << "Final seeds[0] = " << pseeds[0] << "\n"
113  << "Final seeds[1] = " << pseeds[1] << std::endl << std::endl;
114 
115  //***********************************************************************
116  // Prior to the 2.0.4.5 bug fix, 64bit seeds resulted in incorrect randoms
117  seeds[0]=labs(rvals[0]);
118  seeds[1]=labs(rvals[1]);
119  seeds[2]=0;
120 
121  output << std::endl << "********************" << std::endl;
122  output << "This is the case that always fails." << std::endl;
123  output << "seeds[0] = " << seeds[0] << "\n"
124  << "seeds[1] = " << seeds[1] << std::endl << std::endl;
125 
126  g->setTheSeeds(seeds);
127  for (int i=0; i < nNumbers; ++i) {
128  double r = g->flat();
129  if( ! printCheck(i,r,output) ) ++badcount;
130  }
131  pseeds=g->getTheSeeds();
132  output << "Final seeds[0] = " << pseeds[0] << "\n"
133  << "Final seeds[1] = " << pseeds[1] << std::endl << std::endl;
134 
135  //***********************************************************************
136  // recover and reuse seeds
137  seeds[0]=labs(rvals[0]);
138  seeds[1]=labs(rvals[1]);
139  seeds[2]=0;
140 
141  output << std::endl << "********************" << std::endl;
142  output << "Check rolling back a random number seed." << std::endl;
143  output << "seeds[0] = " << seeds[0] << "\n"
144  << "seeds[1] = " << seeds[1] << std::endl << std::endl;
145  std::vector<double> v;
146  g->setTheSeeds(seeds);
147 
148  for (int i=0; i < nNumbers; ++i) {
149  double r = g->flat();
150  if( ! printCheck(i,r,output) ) ++badcount;
151  }
152  pseeds=g->getTheSeeds();
153  seeds[0] = pseeds[0];
154  seeds[1] = pseeds[1];
155  output << " pseeds[0] = " << pseeds[0] << "\n"
156  << "pseeds[1] = " << pseeds[1] << std::endl;
157  for (int i=0; i < nNumbers; ++i) {
158  double r = g->flat();
159  v.push_back(r);
160  }
161  g->setTheSeeds(seeds);
162  for (int i=0; i < nNumbers; ++i) {
163  double r = g->flat();
164  if(! equals(v[i],r) ) {
165  ++badcount;
166  std::cerr << " rollback fails: i, v[i], r "
167  << i << " " << v[i] << " " << r << std::endl;
168  }
169  }
170  output << std::endl;
171 
172  //***********************************************************************
173  // 4-byte positive integers generate valid sequences, which remain within bounds.
174  seeds[0]= labs(static_cast<int>(rvals[0]));
175  seeds[1]= labs(static_cast<int>(rvals[1]));
176  seeds[2]=0;
177 
178  output << std::endl << "********************" << std::endl;
179  output << "This is the case that works." << std::endl;
180  output << std::endl << "seeds[0] = " << seeds[0] << "\n"
181  << "seeds[1] = " << seeds[1] << "\n"
182  << "seeds[2] = " << seeds[2] << std::endl << std::endl;
183 
184  g->setTheSeeds(seeds);
185  for (int i=0; i < nNumbers; ++i) {
186  double r = g->flat();
187  if( ! printCheck(i,r,output) ) ++badcount;
188  }
189  pseeds=g->getTheSeeds();
190  output << "Final seeds[0] = " << pseeds[0] << "\n"
191  << "Final seeds[1] = " << pseeds[1] << std::endl << std::endl;
192 
193  //***********************************************************************
194  // Before the fix, a bad 64bit sequence would eventually rectify itself.
195  // This starts with seeds that would have failed before the 64bit corrections
196  // were applied and loops until both seeds are positive 32-bit integers.
197  // This looping should no longer occur.
198  seeds[0]=labs(rvals[0]);
199  seeds[1]=labs(rvals[1]);
200  seeds[2]=0;
201 
202  output << std::endl << "********************" << std::endl;
203  output << "This case loops until valid short seeds occur." << std::endl;
204  output << "seeds[0] = " << seeds[0] << "\n"
205  << "seeds[1] = " << seeds[1] << std::endl << std::endl;
206 
207  g->setTheSeeds(seeds);
208  // Loop as long as the values are bad.
209  double r;
210  unsigned int low = ~0;
211  unsigned long mask = (~0) << 31;
212  unsigned long skipcount = 0;
213  output << "low = " << low << " mask = " << mask << std::endl;
214  do {
215  r = g->flat();
216  pretend_to_use( r );
217  pseeds = g->getTheSeeds();
218  ++skipcount;
219  } while((pseeds[0]&mask) || (pseeds[1]&mask));
220  if ( skipcount > 1 ) ++badcount;
221 
222  output << std::endl << "Loop terminates on two short seeds." << std::endl;
223  output << "Skipcount = " << skipcount << std::endl;
224  output << "pseeds[0]&mask = " << (pseeds[0]&mask) << std::endl;
225  output << "pseeds[1]&mask = " << (pseeds[1]&mask) << std::endl;
226  output << "Final seeds[0] = " << pseeds[0] << "\n"
227  << "Final seeds[1] = " << pseeds[1] << std::endl << std::endl;
228 
229  output << "This should be a valid sequence." << std::endl;
230  for (int i=0; i < nNumbers; ++i) {
231  double r1 = g->flat();
232  if( ! printCheck(i,r1,output) ) ++badcount;
233  }
234  pseeds=g->getTheSeeds();
235  output << "seeds[0] = " << pseeds[0] << "\n"
236  << "seeds[1] = " << pseeds[1] << std::endl << std::endl;
237 
238  if( badcount > 0 ) std::cout << "Error count is " << badcount << std::endl;
239  return badcount;
240 }
double
#define double(obj)
Definition: excDblThrow.cc:32
a
@ a
Definition: testCategories.cc:125
g
int g(shared_ptr< X >)
Definition: testSharedPtrConvertible.cc:46
equals
bool equals(double a, double b)
Definition: testBug58950.cc:21
b
@ b
Definition: testCategories.cc:125
CLHEP::RanecuEngine
Definition: Matrix/CLHEP/Random/RanecuEngine.h:48
CLHEP::HepRandom::setTheEngine
static void setTheEngine(HepRandomEngine *theNewEngine)
Definition: Random.cc:171
equals01
bool equals01(const std::vector< double > &ab)
Definition: testBug58950.cc:18
CLHEP::HepRandom::getTheGenerator
static HepRandom * getTheGenerator()
Definition: Random.cc:161
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
pretend_to_use
void pretend_to_use(T const &)
Definition: pretend.h:15
output
std::ofstream output("ranRestoreTest.cout")
seeds
Technical Maintenance Note for CLHEP Random Consequences of seeding JamesRandom with positive seed values greater than In the source code JamesRandom The usual way of seeding a generator is via the default which makes use of the table of seeds(with some trickery to ensure that the values won 't repeat after the table rows are exhausted). The trickery preserves the fact that sees are never negative(because the table values are never negative
printCheck
bool printCheck(int &i, double &r, std::ofstream &os)
Definition: testBug58950.cc:27
i
long i
Definition: JamesRandomSeeding.txt:27
pretend.h
main
int main()
Definition: testBug58950.cc:37
in
it has advantages For I leave the ZMthrows in
Definition: keyMergeIssues.doc:62
CLHEP::HepRandom
Definition: Matrix/CLHEP/Random/Random.h:50