CORSIKA  @c8_version@
The framework to simulate particle cascades for astroparticle physics
Engine.hpp
1 /*
2 Copyright 2010-2011, D. E. Shaw Research.
3 All rights reserved.
4 
5 Redistribution and use in source and binary forms, with or without
6 modification, are permitted provided that the following conditions are
7 met:
8 
9 * Redistributions of source code must retain the above copyright
10  notice, this list of conditions, and the following disclaimer.
11 
12 * Redistributions in binary form must reproduce the above copyright
13  notice, this list of conditions, and the following disclaimer in the
14  documentation and/or other materials provided with the distribution.
15 
16 * Neither the name of D. E. Shaw Research nor the names of its
17  contributors may be used to endorse or promote products derived from
18  this software without specific prior written permission.
19 
20 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32 #ifndef __Engine_dot_hpp_
33 #define __Engine_dot_hpp_
34 
35 #include "../features/compilerfeatures.h"
36 #include "../array.h"
37 #include <limits>
38 #include <stdexcept>
39 #include <sstream>
40 #include <algorithm>
41 #include <vector>
42 #if RANDOM_ITERATOR_R123_USE_CXX11_TYPE_TRAITS
43 #include <type_traits>
44 #endif
45 
46 namespace random_iterator_r123 {
67  template <typename CBRNG>
68  struct Engine {
69  typedef CBRNG cbrng_type;
70  typedef typename CBRNG::ctr_type ctr_type;
71  typedef typename CBRNG::key_type key_type;
72  typedef typename CBRNG::ukey_type ukey_type;
73  typedef typename ctr_type::value_type result_type;
74 
75  protected:
76  cbrng_type b;
77  key_type key;
78  ctr_type c;
79  ctr_type v;
80 
81  void fix_invariant() {
82  if (v.back() != 0) {
83  result_type vv = v.back();
84  v = b(c, key);
85  v.back() = vv;
86  }
87  }
88 
89  public:
90  explicit Engine()
91  : b()
92  , c() {
93  ukey_type x = {{}};
94  v.back() = 0;
95  key = x;
96  }
97  explicit Engine(result_type r)
98  : b()
99  , c() {
100  ukey_type x = {{typename ukey_type::value_type(r)}};
101  v.back() = 0;
102  key = x;
103  }
104  // 26.5.3 says that the SeedSeq templates shouldn't particpate in
105  // overload resolution unless the type qualifies as a SeedSeq.
106  // How that is determined is unspecified, except that "as a
107  // minimum a type shall not qualify as a SeedSeq if it is
108  // implicitly convertible to a result_type."
109  //
110  // First, we make sure that even the non-const copy constructor
111  // works as expected. In addition, if we've got C++11
112  // type_traits, we use enable_if and is_convertible to implement
113  // the convertible-to-result_type restriction. Otherwise, the
114  // template is unconditional and will match in some surpirsing
115  // and undesirable situations.
116  Engine(Engine& e)
117  : b(e.b)
118  , key(e.key)
119  , c(e.c) {
120  v.back() = e.v.back();
121  fix_invariant();
122  }
123  Engine(const Engine& e)
124  : b(e.b)
125  , key(e.key)
126  , c(e.c) {
127  v.back() = e.v.back();
128  fix_invariant();
129  }
130 
131  template <typename SeedSeq>
132  explicit Engine(SeedSeq& s
133 #if RANDOM_ITERATOR_R123_USE_CXX11_TYPE_TRAITS
134  ,
135  typename std::enable_if<
136  !std::is_convertible<SeedSeq, result_type>::value>::type* = 0
137 #endif
138  )
139  : b()
140  , c() {
141  ukey_type ukey = ukey_type::seed(s);
142  key = ukey;
143  v.back() = 0;
144  }
145  void seed(result_type r) { *this = Engine(r); }
146  template <typename SeedSeq>
147  void seed(SeedSeq& s
148 #if RANDOM_ITERATOR_R123_USE_CXX11_TYPE_TRAITS
149  ,
150  typename std::enable_if<
151  !std::is_convertible<SeedSeq, result_type>::value>::type* = 0
152 #endif
153  ) {
154  *this = Engine(s);
155  }
156  void seed() { *this = Engine(); }
157  friend bool operator==(const Engine& lhs, const Engine& rhs) {
158  return lhs.c == rhs.c && lhs.v.back() == rhs.v.back() && lhs.key == rhs.key;
159  }
160  friend bool operator!=(const Engine& lhs, const Engine& rhs) {
161  return lhs.c != rhs.c || lhs.v.back() != rhs.v.back() || lhs.key != rhs.key;
162  }
163 
164  friend std::ostream& operator<<(std::ostream& os, const Engine& be) {
165  return os << be.c << " " << be.key << " " << be.v.back();
166  }
167 
168  friend std::istream& operator>>(std::istream& is, Engine& be) {
169  is >> be.c >> be.key >> be.v.back();
170  be.fix_invariant();
171  return is;
172  }
173 
174  // The <random> shipped with MacOS Xcode 4.5.2 imposes a
175  // non-standard requirement that URNGs also have static data
176  // members: _Min and _Max. Later versions of libc++ impose the
177  // requirement only when constexpr isn't supported. Although the
178  // Xcode 4.5.2 requirement is clearly non-standard, it is unlikely
179  // to be fixed and it is very easy work around. We certainly
180  // don't want to go to great lengths to accommodate every buggy
181  // library we come across, but in this particular case, the effort
182  // is low and the benefit is high, so it's worth doing. Thanks to
183  // Yan Zhou for pointing this out to us. See similar code in
184  // ../MicroURNG.hpp
185  const static result_type _Min = 0;
186  const static result_type _Max = ~((result_type)0);
187 
188  static RANDOM_ITERATOR_R123_CONSTEXPR result_type min
189  RANDOM_ITERATOR_R123_NO_MACRO_SUBST() {
190  return _Min;
191  }
192  static RANDOM_ITERATOR_R123_CONSTEXPR result_type max
193  RANDOM_ITERATOR_R123_NO_MACRO_SUBST() {
194  return _Max;
195  }
196 
197  result_type operator()() {
198  if (c.size() == 1) // short-circuit the scalar case. Compilers aren't mind-readers.
199  return b(c.incr(), key)[0];
200  result_type& elem = v.back();
201  if (elem == 0) {
202  v = b(c.incr(), key);
203  result_type ret = v.back();
204  elem = c.size() - 1;
205  return ret;
206  }
207  return v[--elem];
208  }
209 
210  void discard(RANDOM_ITERATOR_R123_ULONG_LONG skip) {
211  // don't forget: elem counts down
212  size_t nelem = c.size();
213  size_t sub = skip % nelem;
214  result_type& elem = v.back();
215  skip /= nelem;
216  if (elem < sub) {
217  elem += nelem;
218  skip++;
219  }
220  elem -= sub;
221  c.incr(skip);
222  fix_invariant();
223  }
224 
225  //--------------------------
226  // Some bonus methods, not required for a Random Number
227  // Engine
228 
229  // Constructors and seed() method for ukey_type seem useful
230  // We need const and non-const to supersede the SeedSeq template.
231  explicit Engine(const ukey_type& uk)
232  : key(uk)
233  , c() {
234  v.back() = 0;
235  }
236  explicit Engine(ukey_type& uk)
237  : key(uk)
238  , c() {
239  v.back() = 0;
240  }
241  void seed(const ukey_type& uk) { *this = Engine(uk); }
242  void seed(ukey_type& uk) { *this = Engine(uk); }
243 
244 #if RANDOM_ITERATOR_R123_USE_CXX11_TYPE_TRAITS
245  template <typename DUMMY = void>
246  explicit Engine(const key_type& k,
247  typename std::enable_if<!std::is_same<ukey_type, key_type>::value,
248  DUMMY>::type* = 0)
249  : key(k)
250  , c() {
251  v.back() = 0;
252  }
253 
254  template <typename DUMMY = void>
255  void seed(const key_type& k,
256  typename std::enable_if<!std::is_same<ukey_type, key_type>::value,
257  DUMMY>::type* = 0) {
258  *this = Engine(k);
259  }
260 #endif
261 
262  // Forward the e(counter) to the CBRNG we are templated
263  // on, using the current value of the key.
264  ctr_type operator()(const ctr_type& c) const { return b(c, key); }
265 
266  key_type getkey() const { return key; }
267 
268  // N.B. setkey(k) is different from seed(k) because seed(k) zeros
269  // the counter (per the C++11 requirements for an Engine), whereas
270  // setkey does not.
271  void setkey(const key_type& k) {
272  key = k;
273  fix_invariant();
274  }
275 
276  // Maybe the caller want's to know the details of
277  // the internal state, e.g., so it can call a different
278  // bijection with the same counter.
279  std::pair<ctr_type, result_type> getcounter() const {
280  return std::make_pair(c, v.back());
281  }
282 
283  // And the inverse.
284  void setcounter(const ctr_type& _c, result_type _elem) {
285  static const size_t nelem = c.size();
286  if (_elem >= nelem)
287  throw std::range_error("Engine::setcounter called with elem out of range");
288  c = _c;
289  v.back() = _elem;
290  fix_invariant();
291  }
292 
293  void setcounter(const std::pair<ctr_type, result_type>& ce) {
294  setcounter(ce.first, ce.second);
295  }
296  };
297 } // namespace random_iterator_r123
298 
299 #endif
If G satisfies the requirements of a CBRNG, and has a ctr_type whose value_type is an unsigned integr...
Definition: Engine.hpp:68