Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2018 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 :
28 : namespace PLMD {
29 :
30 : const double Random::fact=5.9604644775390625e-8; /* 1 / 2^24 */
31 : const double Random::EPS=3.0e-16;
32 : const double Random::AM=1.0/IM;
33 : const double Random::RNMX=(1.0-EPS); // 1.0-EPS;
34 841 : const std::string Random::noname="noname";
35 :
36 1461 : Random::Random(const std::string & name):
37 : switchGaussian(false),
38 : saveGaussian(0.0),
39 : // this is required because if a Random object is created during
40 : // initialization than Random::noname could still be initialized.
41 : // In practice: without this it is not possible to declare
42 : // a static Random object without enforcing the order of the
43 : // static constructors.
44 1461 : name(&name!=&noname?name:"noname")
45 : {
46 1461 : iy=0;
47 1461 : iv[0]=0;
48 1461 : setSeed(0);
49 1461 : }
50 :
51 1472 : void Random::setSeed(int idum_) {
52 1472 : idum=idum_;
53 1472 : incPrec=false;
54 1472 : }
55 :
56 196335 : double Random::RandU01 ()
57 : {
58 196335 : if (incPrec)
59 0 : return U01d();
60 : else
61 196335 : return U01();
62 : }
63 :
64 :
65 0 : double Random::U01d ()
66 : {
67 : double u;
68 0 : u = U01();
69 0 : u += U01() * fact;
70 0 : return (u < 1.0) ? u : (u - 1.0);
71 : }
72 :
73 :
74 3932722 : double Random::U01() {
75 : int j,k;
76 : double temp;
77 3932722 : if (idum <= 0 || !iy) {
78 265 : if (-idum < 1) idum=1;
79 7 : else idum = -idum;
80 10865 : for (j=NTAB+7; j>=0; j--) {
81 10600 : k=idum/IQ;
82 10600 : idum=IA*(idum-k*IQ)-IR*k;
83 10600 : if (idum < 0) idum += IM;
84 10600 : if (j < NTAB) iv[j] = idum;
85 : }
86 265 : iy=iv[0];
87 : }
88 3932722 : k=idum/IQ;
89 3932722 : idum=IA*(idum-k*IQ)-IR*k;
90 3932722 : if (idum < 0) idum += IM;
91 3932722 : j=iy/NDIV;
92 3932722 : iy=iv[j];
93 3932722 : iv[j] = idum;
94 3932722 : if ((temp=AM*iy) > RNMX) return RNMX;
95 3932722 : else return temp;
96 : }
97 :
98 0 : void Random::WriteStateFull(std::ostream & out)const {
99 0 : out<<name<<std::endl;
100 0 : out<<idum<<" "<<iy;
101 0 : for(int i=0; i<NTAB; i++) {
102 0 : out<<" "<<iv[i];
103 : };
104 0 : out<<" "<<switchGaussian;
105 0 : out<<" "<<saveGaussian;
106 0 : out<<std::endl;
107 0 : }
108 :
109 0 : void Random::ReadStateFull (std::istream & in) {
110 0 : getline(in,name);
111 0 : in>>idum>>iy;
112 0 : for (int i = 0; i < NTAB; i++) in>>iv[i];
113 0 : in>>switchGaussian;
114 0 : in>>saveGaussian;
115 0 : }
116 :
117 0 : void Random::toString(std::string & str)const {
118 0 : std::ostringstream ostr;
119 0 : ostr<<idum<<"|"<<iy;
120 0 : for(int i=0; i<NTAB; i++) {
121 0 : ostr<<"|"<<iv[i];
122 : };
123 0 : str=ostr.str();
124 0 : }
125 :
126 0 : void Random::fromString(const std::string & str) {
127 0 : std::string s=str;
128 0 : for(unsigned i=0; i<s.length(); i++) if(s[i]=='|') s[i]=' ';
129 0 : std::istringstream istr(s.c_str());
130 0 : istr>>idum>>iy;
131 0 : for (int i = 0; i < NTAB; i++) istr>>iv[i];
132 0 : }
133 :
134 : // This allows to have the same stream of random numbers
135 : // with different compilers:
136 : #ifdef __INTEL_COMPILER
137 : #pragma intel optimization_level 0
138 : #endif
139 :
140 154287 : double Random::Gaussian() {
141 : double v1,v2,rsq;
142 154287 : if(switchGaussian) {
143 77143 : switchGaussian=false;
144 77143 : return saveGaussian;
145 : }
146 : while(true) {
147 98167 : v1=2.0*RandU01()-1.0;
148 98167 : v2=2.0*RandU01()-1.0;
149 98167 : rsq=v1*v1+v2*v2;
150 98167 : if(rsq<1.0 && rsq>0.0) break;
151 : }
152 77144 : double fac=sqrt(-2.*std::log(rsq)/rsq);
153 77144 : saveGaussian=v1*fac;
154 77144 : switchGaussian=true;
155 98167 : return v2*fac;
156 : }
157 :
158 2523 : }
|
