1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2012-2021 The plumed team
3 (see the PEOPLE file at the root of the distribution for a list of names)
4
5 See http://www.plumed.org for more information.
6
7 This file is part of plumed, version 2.
8
9 plumed is free software: you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
13
14 plumed is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 #include "Random.h"
23 #include <cmath>
24 #include <cstdlib>
25 #include <sstream>
26 #include <iostream>
27 #include <iterator>
28 #include <functional>
29
30 namespace PLMD {
31
32 const double Random::fact=5.9604644775390625e-8; /* 1 / 2^24 */
33 const double Random::EPS=3.0e-16;
34 const double Random::AM=1.0/IM;
35 const double Random::RNMX=(1.0-EPS); // 1.0-EPS;
36 const std::string Random::noname="noname";
37
Random(const std::string & name)38 Random::Random(const std::string & name):
39 switchGaussian(false),
40 saveGaussian(0.0),
41 // this is required because if a Random object is created during
42 // initialization than Random::noname could still be initialized.
43 // In practice: without this it is not possible to declare
44 // a static Random object without enforcing the order of the
45 // static constructors.
46 name(&name!=&noname?name:"noname")
47 {
48 iy=0;
49 for(unsigned i=0; i<NTAB; i++) iv[i]=0;
50 setSeed(0);
51 }
52
setSeed(int idum_)53 void Random::setSeed(int idum_) {
54 if(idum_>0) idum_=-idum_;
55 idum=idum_;
56 incPrec=false;
57 }
58
RandU01()59 double Random::RandU01 ()
60 {
61 if (incPrec)
62 return U01d();
63 else
64 return U01();
65 }
66
67
U01d()68 double Random::U01d ()
69 {
70 double u;
71 u = U01();
72 u += U01() * fact;
73 return (u < 1.0) ? u : (u - 1.0);
74 }
75
U01()76 double Random::U01() {
77 int j,k;
78 double temp;
79 if (idum <= 0 || !iy) {
80 if (-idum < 1) idum=1;
81 else idum = -idum;
82 for (j=NTAB+7; j>=0; j--) {
83 k=idum/IQ;
84 idum=IA*(idum-k*IQ)-IR*k;
85 if (idum < 0) idum += IM;
86 if (j < NTAB) iv[j] = idum;
87 }
88 iy=iv[0];
89 }
90 k=idum/IQ;
91 idum=IA*(idum-k*IQ)-IR*k;
92 if (idum < 0) idum += IM;
93 j=iy/NDIV;
94 iy=iv[j];
95 iv[j] = idum;
96 if ((temp=AM*iy) > RNMX) return RNMX;
97 else return temp;
98 }
99
WriteStateFull(std::ostream & out) const100 void Random::WriteStateFull(std::ostream & out)const {
101 out<<name<<std::endl;
102 out<<idum<<" "<<iy;
103 for(int i=0; i<NTAB; i++) {
104 out<<" "<<iv[i];
105 };
106 out<<" "<<switchGaussian;
107 out<<" "<<saveGaussian;
108 out<<std::endl;
109 }
110
ReadStateFull(std::istream & in)111 void Random::ReadStateFull (std::istream & in) {
112 getline(in,name);
113 in>>idum>>iy;
114 for (int i = 0; i < NTAB; i++) in>>iv[i];
115 in>>switchGaussian;
116 in>>saveGaussian;
117 }
118
toString(std::string & str) const119 void Random::toString(std::string & str)const {
120 std::ostringstream ostr;
121 ostr<<idum<<"|"<<iy;
122 for(int i=0; i<NTAB; i++) {
123 ostr<<"|"<<iv[i];
124 };
125 str=ostr.str();
126 }
127
fromString(const std::string & str)128 void Random::fromString(const std::string & str) {
129 std::string s=str;
130 for(unsigned i=0; i<s.length(); i++) if(s[i]=='|') s[i]=' ';
131 std::istringstream istr(s.c_str());
132 istr>>idum>>iy;
133 for (int i = 0; i < NTAB; i++) istr>>iv[i];
134 }
135
136 // This allows to have the same stream of random numbers
137 // with different compilers:
138 #ifdef __INTEL_COMPILER
139 #pragma intel optimization_level 0
140 #endif
141
Gaussian()142 double Random::Gaussian() {
143 double v1,v2,rsq;
144 if(switchGaussian) {
145 switchGaussian=false;
146 return saveGaussian;
147 }
148 while(true) {
149 v1=2.0*RandU01()-1.0;
150 v2=2.0*RandU01()-1.0;
151 rsq=v1*v1+v2*v2;
152 if(rsq<1.0 && rsq>0.0) break;
153 }
154 double fac=sqrt(-2.*std::log(rsq)/rsq);
155 saveGaussian=v1*fac;
156 switchGaussian=true;
157 return v2*fac;
158 }
159
Shuffle(std::vector<unsigned> & vec)160 void Random::Shuffle(std::vector<unsigned>& vec) {
161 std::iterator_traits<std::vector<unsigned>::iterator >::difference_type i, n;
162 n = vec.end() - vec.begin();
163 for(i=n-1; i>0; --i) {
164 std::swap(vec[i], vec[(int)round(RandU01() * IM) % i]);
165 }
166 }
167
168
169 }
170