1 /* 2 This file is part of MADNESS. 3 4 Copyright (C) 2007,2010 Oak Ridge National Laboratory 5 6 This program is free software; you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation; either version 2 of the License, or 9 (at your option) any later version. 10 11 This program is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with this program; if not, write to the Free Software 18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 19 20 For more information please contact: 21 22 Robert J. Harrison 23 Oak Ridge National Laboratory 24 One Bethel Valley Road 25 P.O. Box 2008, MS-6367 26 27 email: harrisonrj@ornl.gov 28 tel: 865-241-3937 29 fax: 865-572-0680 30 31 $Id$ 32 */ 33 #include <madness/misc/ran.h> 34 #include <iostream> 35 #include <algorithm> 36 #include <cstdio> 37 #include <limits.h> 38 using std::fopen; 39 using std::fwrite; 40 41 42 namespace madness { 43 44 generate()45 void Random::generate() { 46 // Assume we have the lock when we come in here 47 48 double * MADNESS_RESTRICT ur = const_cast<double*>(u); 49 double * MADNESS_RESTRICT us = const_cast<double*>(u)+r-s; 50 for (int i=0; i<s; ++i) { 51 double t = ur[i] + beta*us[i]; 52 ur[i] = t - int(t); 53 } 54 55 ur = const_cast<double*>(u)+s; 56 us = const_cast<double*>(u); 57 int rs = r-s; 58 for (int i=0; i<rs; ++i) { 59 double t = ur[i] + beta*us[i]; 60 ur[i] = t - int(t); 61 } 62 63 cur = 0; 64 } 65 simple()66 unsigned int Random::simple() { 67 return simple_state = 1103515245U*simple_state + 12345U; 68 } 69 70 Random(unsigned int seed)71 Random::Random(unsigned int seed) : r(1279), s(861), beta(7.0), cur(0), u(new double [r]) { 72 // If you switch r don't forget to change size in RandomState 73 74 // a) not clear if beta != 1 is an improvement. 75 // b) must ensure s >= r/2. 76 // c) r=19937, s=10095 or r=1279, s=861 seem worse than 4423/3004 ? 77 // Random(unsigned int seed = 5461) : r(4423), s(3004), beta(7.0), cur(0) { 78 79 setstate(seed); 80 } 81 setstate(unsigned int seed)82 void Random::setstate(unsigned int seed) { 83 ScopedMutex<Mutex> safe(this); 84 // Initialize startup generator 85 if ((seed&1) == 0) seed += 1; 86 simple_state = seed; 87 for (int i=0; i<10*r; ++i) simple(); 88 89 // Initialize stream with 48 bit values by first generating 90 // roughly 52 bit values and then truncating to exactly 48 bits 91 double two21 = 2097152.0; 92 double two52 = 4503599627370496.0; 93 double two24 = 16777216.0; 94 double rtwo24 = 5.9604644775390625e-08; 95 for (int i=0; i<r; ++i) u[i] = double(simple()); 96 for (int i=0; i<r; ++i) u[i] += double(simple())*two21; 97 for (int i=0; i<r; ++i) u[i] /= two52; 98 // Next line breaks on Cray X1 CC 5.3 ... sigh 99 //for (int i=0; i<r; ++i) u[i] -= int(u[i]); 100 for (int i=0; i<r; ++i) {int tmp=int(u[i]); u[i] -= double(tmp);} 101 for (int i=0; i<r; ++i) { 102 int high = int(two24*u[i]); 103 int lo = int(two24*(two24*u[i]-high)); 104 u[i] = rtwo24*(high + rtwo24*lo); 105 } 106 107 // Verify that we have only set 48 bits and that at least 108 // some of the 48th bits are set. 109 double n48 = 0; 110 for (int i=0; i<r; ++i) { 111 double rem = u[i]*two24; 112 rem -= int(rem); 113 double rem48 = rem*two24; 114 rem48 -= int(rem48); 115 if (rem48 != 0) throw "Random: bad bits?"; 116 117 double rem47 = rem*two24*0.5; 118 rem47 -= int(rem47); 119 if (rem47 != 0) ++n48; 120 } 121 if (n48 == 0) throw "Random: bad 48'th bit?"; 122 123 // Warm up 124 for (int i=0; i<2000; ++i) generate(); 125 } 126 ~Random()127 Random::~Random() {delete [] u;} 128 129 130 // // Might be faster than the one-byte-per-word version 131 // void getbytes(int n, unsigned char * MADNESS_RESTRICT v) { 132 // const double two24 = 16777216.0; 133 // int n6 = n%6; 134 // n -= n6; 135 // while (n>0) { 136 // if (cur >= r) generate(); 137 // int ndo = std::min(n/6,r-cur); 138 // for (int i=0; i<ndo; ++i) { 139 // unsigned int high = int(two24*u[i+cur]); 140 // unsigned int lo = int(two24*(two24*u[i+cur]-high)); 141 // *v++ = (high>>16)&0xff; 142 // *v++ = (high>>8 )&0xff; 143 // *v++ = (high )&0xff; 144 // *v++ = (lo>>16 )&0xff; 145 // *v++ = (lo>>8 )&0xff; 146 // *v++ = (lo )&0xff; 147 // }; 148 // n -= ndo*6; 149 // cur += ndo; 150 // } 151 // for (int i=0; i<n6; ++i) *v++ = (unsigned char) (256*get()); 152 // }; 153 getbytes(int n,unsigned char * MADNESS_RESTRICT v)154 void Random::getbytes(int n, unsigned char * MADNESS_RESTRICT v) { 155 ScopedMutex<Mutex> safe(this); 156 while (n) { 157 if (cur >= r) generate(); 158 int ndo = std::min(n,r-cur); 159 const double* ucur = const_cast<const double*>(u) + cur; 160 for (int i=0; i<ndo; ++i) v[i] = (unsigned char) (256.0*ucur[i]); 161 n -= ndo; 162 v += ndo; 163 cur += ndo; 164 } 165 } 166 getstate() const167 RandomState Random::getstate() const { 168 ScopedMutex<Mutex> safe(this); 169 RandomState s; 170 s.cur = cur; 171 for (int i=0; i<r; ++i) s.u[i] = u[i]; 172 return s; 173 } 174 setstate(const RandomState & s)175 void Random::setstate(const RandomState &s) { 176 ScopedMutex<Mutex> safe(this); 177 cur = s.cur; 178 for (int i=0; i<r; ++i) u[i] = s.u[i]; 179 } 180 181 test()182 void Random::test() { 183 // Crude tests just to find gross errors ... the file written 184 // at the end is used as input to diehard. 185 Random r; 186 double sum1 = 0; 187 double sum2 = 0; 188 double x12 = 0; 189 double two24 = 16777216.0; 190 191 double n = 100000000.0; 192 double nn = n; 193 const int nbuf = 64; 194 nn /= nbuf; 195 n = nn*nbuf; 196 197 double buf[nbuf]; 198 199 while (nn--) { 200 r.getv(nbuf,buf); 201 for (int i=0; i<nbuf; ++i) { 202 double d1 = buf[i]; 203 sum1 += d1; 204 double d2 = d1 * two24; 205 d2 -= int(d2); 206 sum2 += d2; 207 x12 += (d1-0.5)*(d2-0.5); 208 } 209 } 210 211 std::cout << "high " << sum1/n << std::endl; 212 std::cout << "lo " << sum2/n << std::endl; 213 std::cout << "hi-lo " << x12/n << std::endl; 214 215 const int nb = 100000000; 216 unsigned char *b = new unsigned char[nb+1]; 217 218 b[nb-1] = 0; 219 b[nb] = 99; 220 r.getbytes(nb,b); 221 222 std::cout << int(b[nb-1]) << std::endl; 223 std::cout << int(b[nb]) << std::endl; 224 225 FILE *f = fopen("stream","w+"); 226 if (!f) {std::cout << "fopen?\n"; std::exit(1);} 227 fwrite(b, 1, nb, f); 228 fclose(f); 229 } 230 231 Random default_random_generator; 232 233 RandomValue()234 template <> double RandomValue<double> () { 235 return default_random_generator.get(); 236 } 237 RandomValue()238 template <> float RandomValue<float> () { 239 return float(default_random_generator.get()); 240 } 241 RandomValue()242 template <> double_complex RandomValue<double_complex> () { 243 return double_complex(RandomValue<double>(),RandomValue<double>()); 244 } 245 RandomValue()246 template <> float_complex RandomValue<float_complex> () { 247 return float_complex(RandomValue<float>(),RandomValue<float>()); 248 } 249 RandomValue()250 template <> int RandomValue<int> () { 251 return int(2147483648.0 * RandomValue<double>()); 252 } 253 RandomValue()254 template <> long RandomValue<long> () { 255 return long(2147483648.0 * RandomValue<double>()); 256 } 257 RandomVector(int n,double * t)258 template <> void RandomVector<double>(int n, double* t) { 259 default_random_generator.getv(n, t); 260 } 261 RandomVector(int n,float * t)262 template <> void RandomVector<float>(int n, float* t) { 263 default_random_generator.getv(n, t); 264 } 265 RandomVector(int n,double_complex * t)266 template <> void RandomVector<double_complex>(int n, double_complex* t) { 267 default_random_generator.getv(2*n, (double*)(t)); 268 } 269 RandomVector(int n,float_complex * t)270 template <> void RandomVector<float_complex>(int n, float_complex* t) { 271 default_random_generator.getv(2*n, (float*)(t)); 272 } 273 } 274 275 // using namespace madness; 276 277 // int main() { 278 // Random::test(); 279 // return 0; 280 // } 281