CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // ---------------------------------------------------------------------- 00002 // 00003 // testBug58950 -- test problem with RanecuEngine on 64bit machines 00004 // 00005 // R. Weller 11/11/09 initial test from Vanderbilt 00006 // L. Garren 12/1/09 rewritten for test suite 00007 // 00008 // ---------------------------------------------------------------------- 00009 #include <iostream> 00010 #include <stdexcept> 00011 #include <cmath> 00012 #include <stdlib.h> 00013 #include <limits> 00014 #include <complex> 00015 #include "CLHEP/Random/RanecuEngine.h" 00016 #include "CLHEP/Random/Random.h" 00017 #include "pretend.h" 00018 00019 bool printCheck( int & i, double & r, std::ofstream & os ) 00020 { 00021 os << i << " " << r << std::endl; 00022 if (r < 0 || r > 1.0 ) { 00023 std::cout << "Error: bad random number " << r << std::endl; 00024 return false; 00025 } 00026 return true; 00027 } 00028 00029 int main() { 00030 00031 std::ofstream output("testBug58950.cout"); 00032 00033 output << std::endl << "short " << sizeof(short) << std::endl; 00034 output << "int " << sizeof(int) << std::endl; 00035 output << "unsigned int " << sizeof(unsigned int) << std::endl; 00036 output << "long " << sizeof(long) << std::endl; 00037 output << "float " << sizeof(float) << std::endl; 00038 output << "double " << sizeof(double) << std::endl; 00039 output << "long double " << sizeof(long double) << std::endl << std::endl; 00040 00041 CLHEP::RanecuEngine *eng = new CLHEP::RanecuEngine; 00042 CLHEP::HepRandom::setTheEngine(eng); 00043 CLHEP::HepRandom *g; 00044 g=CLHEP::HepRandom::getTheGenerator(); 00045 00046 long rvals[2]; 00047 try { 00048 std::ifstream in("/dev/urandom", std::ios::in | std::ios::binary); 00049 if(in.is_open()) { 00050 in.read((char *)(&rvals), 2*sizeof(long)); 00051 in.close(); 00052 if(in.fail()) { 00053 throw std::runtime_error("File read error"); 00054 } 00055 } else throw std::runtime_error("File open error"); 00056 } catch(std::runtime_error e) { 00057 std::ostringstream dStr; 00058 dStr << "Error: " << e.what() 00059 << " processing seed from file \"" << "/dev/urandom" << "\"."; 00060 throw std::runtime_error(dStr.str().c_str()); 00061 } 00062 00063 int nNumbers = 20; 00064 int badcount = 0; 00065 00066 long seeds[3]; 00067 const long *pseeds; 00068 //*********************************************************************** 00069 // Seeds are expected to be positive. Therefore, if either seed is 00070 // negative then prior to 2.0.4.5 the generator set initial conditions 00071 // and generated the same sequence of numbers no matter what the seeds were. 00072 seeds[0]=rvals[0]; 00073 seeds[1]=rvals[1]; 00074 seeds[2]=0; 00075 if( rvals[0] > 0 ) seeds[0] = -rvals[0]; 00076 00077 double negseq[20] = { 0.154707, 0.587114, 0.702059, 0.566, 0.988325, 00078 0.525921, 0.191554, 0.269338, 0.234277, 0.358997, 00079 0.549936, 0.296877, 0.162243, 0.227732, 0.528862, 00080 0.631571, 0.176462, 0.247858, 0.170025, 0.284483 }; 00081 double eps = 1.0E-6; 00082 00083 output << std::endl << "********************" << std::endl; 00084 output << "This is the case that may or may not fail." << std::endl; 00085 output << "However, if it has values in (0,1), they are a " << std::endl 00086 << "deterministic sequence beginning with 0.154707." << std::endl; 00087 output << "seeds[0] = " << seeds[0] << "\n" 00088 << "seeds[1] = " << seeds[1] << std::endl << std::endl; 00089 00090 g->setTheSeeds(seeds); 00091 int rseq = 0; 00092 for (int i=0; i < nNumbers; ++i) { 00093 double r = g->flat(); 00094 if( ! printCheck(i,r,output) ) ++badcount; 00095 // before the change, the random number sequence was reliably the same 00096 if( std::fabs(r-negseq[i]) < eps ) { 00097 std::cout << " reproducing sequence " << i << " " 00098 << r << " " << negseq[i] << std::endl; 00099 ++rseq; 00100 } 00101 } 00102 if( rseq == 20 ) ++badcount; 00103 pseeds=g->getTheSeeds(); 00104 output << "Final seeds[0] = " << pseeds[0] << "\n" 00105 << "Final seeds[1] = " << pseeds[1] << std::endl << std::endl; 00106 00107 //*********************************************************************** 00108 // Prior to the 2.0.4.5 bug fix, 64bit seeds resulted in incorrect randoms 00109 seeds[0]=labs(rvals[0]); 00110 seeds[1]=labs(rvals[1]); 00111 seeds[2]=0; 00112 00113 output << std::endl << "********************" << std::endl; 00114 output << "This is the case that always fails." << std::endl; 00115 output << "seeds[0] = " << seeds[0] << "\n" 00116 << "seeds[1] = " << seeds[1] << std::endl << std::endl; 00117 00118 g->setTheSeeds(seeds); 00119 for (int i=0; i < nNumbers; ++i) { 00120 double r = g->flat(); 00121 if( ! printCheck(i,r,output) ) ++badcount; 00122 } 00123 pseeds=g->getTheSeeds(); 00124 output << "Final seeds[0] = " << pseeds[0] << "\n" 00125 << "Final seeds[1] = " << pseeds[1] << std::endl << std::endl; 00126 00127 //*********************************************************************** 00128 // recover and reuse seeds 00129 seeds[0]=labs(rvals[0]); 00130 seeds[1]=labs(rvals[1]); 00131 seeds[2]=0; 00132 00133 output << std::endl << "********************" << std::endl; 00134 output << "Check rolling back a random number seed." << std::endl; 00135 output << "seeds[0] = " << seeds[0] << "\n" 00136 << "seeds[1] = " << seeds[1] << std::endl << std::endl; 00137 std::vector<double> v; 00138 g->setTheSeeds(seeds); 00139 00140 for (int i=0; i < nNumbers; ++i) { 00141 double r = g->flat(); 00142 if( ! printCheck(i,r,output) ) ++badcount; 00143 } 00144 pseeds=g->getTheSeeds(); 00145 seeds[0] = pseeds[0]; 00146 seeds[1] = pseeds[1]; 00147 output << " pseeds[0] = " << pseeds[0] << "\n" 00148 << "pseeds[1] = " << pseeds[1] << std::endl; 00149 for (int i=0; i < nNumbers; ++i) { 00150 double r = g->flat(); 00151 v.push_back(r); 00152 } 00153 g->setTheSeeds(seeds); 00154 for (int i=0; i < nNumbers; ++i) { 00155 double r = g->flat(); 00156 // if(v[i] != r ) { 00157 if(std::abs(v[i] - r) >= std::numeric_limits<double>::epsilon()) { 00158 ++badcount; 00159 std::cerr << " rollback fails: i, v[i], r " 00160 << i << " " << v[i] << " " << r << std::endl; 00161 } 00162 } 00163 output << std::endl; 00164 00165 //*********************************************************************** 00166 // 4-byte positive integers generate valid sequences, which remain within bounds. 00167 seeds[0]= labs(static_cast<int>(rvals[0])); 00168 seeds[1]= labs(static_cast<int>(rvals[1])); 00169 seeds[2]=0; 00170 00171 output << std::endl << "********************" << std::endl; 00172 output << "This is the case that works." << std::endl; 00173 output << std::endl << "seeds[0] = " << seeds[0] << "\n" 00174 << "seeds[1] = " << seeds[1] << "\n" 00175 << "seeds[2] = " << seeds[2] << std::endl << std::endl; 00176 00177 g->setTheSeeds(seeds); 00178 for (int i=0; i < nNumbers; ++i) { 00179 double r = g->flat(); 00180 if( ! printCheck(i,r,output) ) ++badcount; 00181 } 00182 pseeds=g->getTheSeeds(); 00183 output << "Final seeds[0] = " << pseeds[0] << "\n" 00184 << "Final seeds[1] = " << pseeds[1] << std::endl << std::endl; 00185 00186 //*********************************************************************** 00187 // Before the fix, a bad 64bit sequence would eventually rectify itself. 00188 // This starts with seeds that would have failed before the 64bit corrections 00189 // were applied and loops until both seeds are positive 32-bit integers. 00190 // This looping should no longer occur. 00191 seeds[0]=labs(rvals[0]); 00192 seeds[1]=labs(rvals[1]); 00193 seeds[2]=0; 00194 00195 output << std::endl << "********************" << std::endl; 00196 output << "This case loops until valid short seeds occur." << std::endl; 00197 output << "seeds[0] = " << seeds[0] << "\n" 00198 << "seeds[1] = " << seeds[1] << std::endl << std::endl; 00199 00200 g->setTheSeeds(seeds); 00201 // Loop as long as the values are bad. 00202 double r; 00203 unsigned int low = ~0; 00204 unsigned long mask = (~0) << 31; 00205 unsigned long skipcount = 0; 00206 output << "low = " << low << " mask = " << mask << std::endl; 00207 do { 00208 r = g->flat(); 00209 pretend_to_use( r ); 00210 pseeds = g->getTheSeeds(); 00211 ++skipcount; 00212 } while((pseeds[0]&mask) || (pseeds[1]&mask)); 00213 if ( skipcount > 1 ) ++badcount; 00214 00215 output << std::endl << "Loop terminates on two short seeds." << std::endl; 00216 output << "Skipcount = " << skipcount << std::endl; 00217 output << "pseeds[0]&mask = " << (pseeds[0]&mask) << std::endl; 00218 output << "pseeds[1]&mask = " << (pseeds[1]&mask) << std::endl; 00219 output << "Final seeds[0] = " << pseeds[0] << "\n" 00220 << "Final seeds[1] = " << pseeds[1] << std::endl << std::endl; 00221 00222 output << "This should be a valid sequence." << std::endl; 00223 for (int i=0; i < nNumbers; ++i) { 00224 double r1 = g->flat(); 00225 if( ! printCheck(i,r1,output) ) ++badcount; 00226 } 00227 pseeds=g->getTheSeeds(); 00228 output << "seeds[0] = " << pseeds[0] << "\n" 00229 << "seeds[1] = " << pseeds[1] << std::endl << std::endl; 00230 00231 if( badcount > 0 ) std::cout << "Error count is " << badcount << std::endl; 00232 return badcount; 00233 }