MCKL
Monte Carlo Kernel Library
random/internal/common.hpp
Go to the documentation of this file.
1 //============================================================================
2 // MCKL/include/mckl/random/internal/common.hpp
3 //----------------------------------------------------------------------------
4 // MCKL: Monte Carlo Kernel Library
5 //----------------------------------------------------------------------------
6 // Copyright (c) 2013-2018, Yan Zhou
7 // All rights reserved.
8 //
9 // Redistribution and use in source and binary forms, with or without
10 // modification, are permitted provided that the following conditions are met:
11 //
12 // Redistributions of source code must retain the above copyright notice,
13 // this list of conditions and the following disclaimer.
14 //
15 // Redistributions in binary form must reproduce the above copyright notice,
16 // this list of conditions and the following disclaimer in the documentation
17 // and/or other materials provided with the distribution.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS AS IS
20 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 // POSSIBILITY OF SUCH DAMAGE.
30 //============================================================================
31 
32 #ifndef MCKL_RANDOM_INTERNAL_COMMON_HPP
33 #define MCKL_RANDOM_INTERNAL_COMMON_HPP
34 
35 #include <mckl/internal/common.hpp>
36 
37 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_ASSERT_REAL_TYPE(Name) \
38  static_assert(std::is_floating_point<RealType>::value, \
39  "**" #Name \
40  "Distribution** used with RealType other than floating point types");
41 
42 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_ASSERT_BLAS_TYPE(Name) \
43  static_assert(::mckl::internal::is_blas_floating_point<RealType>::value, \
44  "**" #Name \
45  "Distribution** used with RealType other than float or double");
46 
47 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_ASSERT_INT_TYPE(Name, MinBits) \
48  static_assert(std::is_integral<IntType>::value, \
49  "**" #Name \
50  "Distribution** used with IntType other than integer types"); \
51  static_assert(CHAR_BIT * sizeof(IntType) >= MinBits, \
52  "**" #Name "Distribution** used with IntType smaller than " #MinBits \
53  " bits");
54 
55 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_ASSERT_UINT_TYPE(Name, MinBits) \
56  static_assert(std::is_unsigned<UIntType>::value, \
57  "**" #Name \
58  "Distribution** used with UIntType other than unsigned integer " \
59  "types"); \
60  static_assert(std::numeric_limits<UIntType>::digits >= MinBits, \
61  "**" #Name "Distribution** used with IntType smaller than " #MinBits \
62  " bits");
63 
64 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_ASSERT_REAL_TYPE(Name) \
65  static_assert(std::is_floating_point<RealType>::value, \
66  "**" #Name \
67  "Distribution** used with RealType other than floating point types");
68 
69 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_BATCH_0(Name, name, T) \
70  template <typename T, typename RNGType> \
71  inline void name##_distribution(RNGType &rng, std::size_t n, T *r) \
72  { \
73  const std::size_t K = BufferSize<T>::value; \
74  const std::size_t M = n / K; \
75  const std::size_t L = n % K; \
76  for (std::size_t i = 0; i != M; ++i, r += K) \
77  ::mckl::internal::name##_distribution_impl<K>(rng, K, r); \
78  ::mckl::internal::name##_distribution_impl<K>(rng, L, r); \
79  } \
80  \
81  template <typename T, typename RNGType> \
82  inline void name##_distribution(RNGType &rng, std::size_t N, T *r, \
83  const typename Name##Distribution<T>::param_type &) \
84  { \
85  name##_distribution(rng, N, r); \
86  }
87 
88 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_BATCH_1(Name, name, T, T1, p1) \
89  template <typename T, typename RNGType> \
90  inline void name##_distribution(RNGType &rng, std::size_t N, T *r, T1 p1) \
91  { \
92  const std::size_t K = BufferSize<T>::value; \
93  const std::size_t M = N / K; \
94  const std::size_t L = N % K; \
95  for (std::size_t i = 0; i != M; ++i, r += K) \
96  ::mckl::internal::name##_distribution_impl<K>(rng, K, r, p1); \
97  ::mckl::internal::name##_distribution_impl<K>(rng, L, r, p1); \
98  } \
99  \
100  template <typename T, typename RNGType> \
101  inline void name##_distribution(RNGType &rng, std::size_t N, T *r, \
102  const typename Name##Distribution<T>::param_type &param) \
103  { \
104  name##_distribution(rng, N, r, param.p1()); \
105  }
106 
107 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_BATCH_2( \
108  Name, name, T, T1, p1, T2, p2) \
109  template <typename T, typename RNGType> \
110  inline void name##_distribution( \
111  RNGType &rng, std::size_t N, T *r, T1 p1, T2 p2) \
112  { \
113  const std::size_t K = BufferSize<T>::value; \
114  const std::size_t M = N / K; \
115  const std::size_t L = N % K; \
116  for (std::size_t i = 0; i != M; ++i, r += K) \
117  ::mckl::internal::name##_distribution_impl<K>(rng, K, r, p1, p2); \
118  ::mckl::internal::name##_distribution_impl<K>(rng, L, r, p1, p2); \
119  } \
120  \
121  template <typename T, typename RNGType> \
122  inline void name##_distribution(RNGType &rng, std::size_t N, T *r, \
123  const typename Name##Distribution<T>::param_type &param) \
124  { \
125  name##_distribution(rng, N, r, param.p1(), param.p2()); \
126  }
127 
128 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_BATCH_4( \
129  Name, name, T, T1, p1, T2, p2, T3, p3, T4, p4) \
130  template <typename T, typename RNGType> \
131  inline void name##_distribution( \
132  RNGType &rng, std::size_t N, T *r, T1 p1, T2 p2, T3 p3, T4 p4) \
133  { \
134  const std::size_t K = BufferSize<T>::value; \
135  const std::size_t M = N / K; \
136  const std::size_t L = N % K; \
137  for (std::size_t i = 0; i != M; ++i, r += K) { \
138  ::mckl::internal::name##_distribution_impl<K>( \
139  rng, K, r, p1, p2, p3, p4); \
140  } \
141  ::mckl::internal::name##_distribution_impl<K>( \
142  rng, L, r, p1, p2, p3, p4); \
143  } \
144  \
145  template <typename T, typename RNGType> \
146  inline void name##_distribution(RNGType &rng, std::size_t N, T *r, \
147  const typename Name##Distribution<T>::param_type &param) \
148  { \
149  name##_distribution( \
150  rng, N, r, param.p1(), param.p2(), param.p3(), param.p4()); \
151  }
152 
153 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_PARAM_TYPE_0(Name, T) \
154  public: \
155  class param_type \
156  { \
157  public: \
158  using result_type = T; \
159  using distribution_type = Name##Distribution<T>; \
160  \
161  friend bool operator==(const param_type &, const param_type &) \
162  { \
163  return true; \
164  } \
165  \
166  friend bool operator!=(const param_type &, const param_type &) \
167  { \
168  return false; \
169  } \
170  \
171  template <typename CharT, typename Traits> \
172  friend std::basic_ostream<CharT, Traits> &operator<<( \
173  std::basic_ostream<CharT, Traits> &os, const param_type &) \
174  { \
175  return os; \
176  } \
177  \
178  template <typename CharT, typename Traits> \
179  friend std::basic_istream<CharT, Traits> &operator>>( \
180  std::basic_istream<CharT, Traits> &is, param_type &) \
181  { \
182  return is; \
183  } \
184  \
185  private: \
186  friend distribution_type; \
187  }; // class param_type
188 
189 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_PARAM_TYPE_1( \
190  Name, name, T, T1, p1, v1) \
191  public: \
192  class param_type \
193  { \
194  public: \
195  using result_type = T; \
196  using distribution_type = Name##Distribution<T>; \
197  \
198  explicit param_type(T1 p1 = v1) : p1##_(p1) \
199  { \
200  ::mckl::runtime_assert( \
201  ::mckl::internal::name##_distribution_check_param(p1), \
202  "**" #Name \
203  "Distribution** constructed with invalid arguments"); \
204  } \
205  \
206  T1 p1() const { return p1##_; } \
207  \
208  friend bool operator==( \
209  const param_type &param1, const param_type &param2) \
210  { \
211  MCKL_PUSH_CLANG_WARNING("-Wfloat-equal") \
212  MCKL_PUSH_INTEL_WARNING(1572) /* floating-point comparison */ \
213  if (param1.p1##_ != param2.p1##_) \
214  return false; \
215  return true; \
216  MCKL_POP_CLANG_WARNING \
217  MCKL_POP_INTEL_WARNING \
218  } \
219  \
220  friend bool operator!=( \
221  const param_type &param1, const param_type &param2) \
222  { \
223  return !(param1 == param2); \
224  } \
225  \
226  template <typename CharT, typename Traits> \
227  friend std::basic_ostream<CharT, Traits> &operator<<( \
228  std::basic_ostream<CharT, Traits> &os, const param_type &param) \
229  { \
230  if (!os) \
231  return os; \
232  \
233  os << param.p1##_; \
234  \
235  return os; \
236  } \
237  \
238  template <typename CharT, typename Traits> \
239  friend std::basic_istream<CharT, Traits> &operator>>( \
240  std::basic_istream<CharT, Traits> &is, param_type &param) \
241  { \
242  if (!is) \
243  return is; \
244  \
245  T1 p1 = 0; \
246  is >> std::ws >> p1; \
247  \
248  if (is) { \
249  if (::mckl::internal::name##_distribution_check_param(p1)) { \
250  param.p1##_ = p1; \
251  } else { \
252  is.setstate(std::ios_base::failbit); \
253  } \
254  } \
255  \
256  return is; \
257  } \
258  \
259  private: \
260  T1 p1##_; \
261  \
262  friend distribution_type; \
263  }; // class param_type
264 
265 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_PARAM_TYPE_2( \
266  Name, name, T, T1, p1, v1, T2, p2, v2) \
267  public: \
268  class param_type \
269  { \
270  public: \
271  using result_type = T; \
272  using distribution_type = Name##Distribution<T>; \
273  \
274  explicit param_type(T1 p1 = v1, T2 p2 = v2) : p1##_(p1), p2##_(p2) \
275  { \
276  ::mckl::runtime_assert( \
277  ::mckl::internal::name##_distribution_check_param(p1, p2), \
278  "**" #Name \
279  "Distribution** constructed with invalid arguments"); \
280  } \
281  \
282  T1 p1() const { return p1##_; } \
283  T2 p2() const { return p2##_; } \
284  \
285  friend bool operator==( \
286  const param_type &param1, const param_type &param2) \
287  { \
288  MCKL_PUSH_CLANG_WARNING("-Wfloat-equal") \
289  MCKL_PUSH_INTEL_WARNING(1572) /* floating-point comparison */ \
290  if (param1.p1##_ != param2.p1##_) \
291  return false; \
292  if (param1.p2##_ != param2.p2##_) \
293  return false; \
294  return true; \
295  MCKL_POP_CLANG_WARNING \
296  MCKL_POP_INTEL_WARNING \
297  } \
298  \
299  friend bool operator!=( \
300  const param_type &param1, const param_type &param2) \
301  { \
302  return !(param1 == param2); \
303  } \
304  \
305  template <typename CharT, typename Traits> \
306  friend std::basic_ostream<CharT, Traits> &operator<<( \
307  std::basic_ostream<CharT, Traits> &os, const param_type &param) \
308  { \
309  if (!os) \
310  return os; \
311  \
312  os << param.p1##_ << ' '; \
313  os << param.p2##_; \
314  \
315  return os; \
316  } \
317  \
318  template <typename CharT, typename Traits> \
319  friend std::basic_istream<CharT, Traits> &operator>>( \
320  std::basic_istream<CharT, Traits> &is, param_type &param) \
321  { \
322  if (!is) \
323  return is; \
324  \
325  T1 p1 = 0; \
326  T2 p2 = 0; \
327  is >> std::ws >> p1; \
328  is >> std::ws >> p2; \
329  \
330  if (is) { \
331  if (::mckl::internal::name##_distribution_check_param( \
332  p1, p2)) { \
333  param.p1##_ = p1; \
334  param.p2##_ = p2; \
335  } else { \
336  is.setstate(std::ios_base::failbit); \
337  } \
338  } \
339  \
340  return is; \
341  } \
342  \
343  private: \
344  T1 p1##_; \
345  T2 p2##_; \
346  \
347  friend distribution_type; \
348  }; // class param_type
349 
350 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_PARAM_TYPE_4( \
351  Name, name, T, T1, p1, v1, T2, p2, v2, T3, p3, v3, T4, p4, v4) \
352  public: \
353  class param_type \
354  { \
355  public: \
356  using result_type = T; \
357  using distribution_type = Name##Distribution<T>; \
358  \
359  explicit param_type(T1 p1 = v1, T2 p2 = v2, T3 p3 = v3, T4 p4 = v4) \
360  : p1##_(p1), p2##_(p2), p3##_(p3), p4##_(p4) \
361  { \
362  ::mckl::runtime_assert( \
363  ::mckl::internal::name##_distribution_check_param( \
364  p1, p2, p3, p4), \
365  "**" #Name \
366  "Distribution** constructed with invalid arguments"); \
367  } \
368  \
369  T1 p1() const { return p1##_; } \
370  T2 p2() const { return p2##_; } \
371  T3 p3() const { return p3##_; } \
372  T4 p4() const { return p4##_; } \
373  \
374  friend bool operator==( \
375  const param_type &param1, const param_type &param2) \
376  { \
377  MCKL_PUSH_CLANG_WARNING("-Wfloat-equal") \
378  MCKL_PUSH_INTEL_WARNING(1572) /* floating-point comparison */ \
379  if (param1.p1##_ != param2.p1##_) \
380  return false; \
381  if (param1.p2##_ != param2.p2##_) \
382  return false; \
383  if (param1.p3##_ != param2.p3##_) \
384  return false; \
385  if (param1.p4##_ != param2.p4##_) \
386  return false; \
387  return true; \
388  MCKL_POP_CLANG_WARNING \
389  MCKL_POP_INTEL_WARNING \
390  } \
391  \
392  friend bool operator!=( \
393  const param_type &param1, const param_type &param2) \
394  { \
395  return !(param1 == param2); \
396  } \
397  \
398  template <typename CharT, typename Traits> \
399  friend std::basic_ostream<CharT, Traits> &operator<<( \
400  std::basic_ostream<CharT, Traits> &os, const param_type &param) \
401  { \
402  if (!os) \
403  return os; \
404  \
405  os << param.p1##_ << ' '; \
406  os << param.p2##_ << ' '; \
407  os << param.p3##_ << ' '; \
408  os << param.p4##_; \
409  \
410  return os; \
411  } \
412  \
413  template <typename CharT, typename Traits> \
414  friend std::basic_istream<CharT, Traits> &operator>>( \
415  std::basic_istream<CharT, Traits> &is, param_type &param) \
416  { \
417  if (!is) \
418  return is; \
419  \
420  T1 p1 = 0; \
421  T2 p2 = 0; \
422  T3 p3 = 0; \
423  T4 p4 = 0; \
424  is >> std::ws >> p1; \
425  is >> std::ws >> p2; \
426  is >> std::ws >> p3; \
427  is >> std::ws >> p4; \
428  \
429  if (is) { \
430  if (::mckl::internal::name##_distribution_check_param( \
431  p1, p2, p3, p4)) { \
432  param.p1##_ = p1; \
433  param.p2##_ = p2; \
434  param.p3##_ = p3; \
435  param.p4##_ = p4; \
436  } else { \
437  is.setstate(std::ios_base::failbit); \
438  } \
439  } \
440  \
441  return is; \
442  } \
443  \
444  private: \
445  T1 p1##_; \
446  T2 p2##_; \
447  T3 p3##_; \
448  T4 p4##_; \
449  \
450  friend distribution_type; \
451  }; // class param_type
452 
453 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_CONSTRUCTOR_0(Name, T) \
454  public: \
455  using result_type = T; \
456  using distribution_type = Name##Distribution<T>; \
457  \
458  Name##Distribution() = default; \
459  explicit Name##Distribution(const param_type &) { reset(); } \
460  explicit Name##Distribution(param_type &&) { reset(); }
461 
462 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_CONSTRUCTOR_1(Name, T, T1, p1, v1) \
463  public: \
464  using result_type = T; \
465  using distribution_type = Name##Distribution<T>; \
466  \
467  explicit Name##Distribution(T1 p1 = v1) : param_(p1) { reset(); } \
468  \
469  explicit Name##Distribution(const param_type &param) : param_(param) \
470  { \
471  reset(); \
472  } \
473  \
474  explicit Name##Distribution(param_type &&param) \
475  : param_(std::move(param)) \
476  { \
477  reset(); \
478  } \
479  \
480  T1 p1() const { return param_.p1(); }
481 
482 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_CONSTRUCTOR_2( \
483  Name, T, T1, p1, v1, T2, p2, v2) \
484  public: \
485  using result_type = T; \
486  using distribution_type = Name##Distribution<T>; \
487  \
488  explicit Name##Distribution(T1 p1 = v1, T2 p2 = v2) : param_(p1, p2) \
489  { \
490  reset(); \
491  } \
492  \
493  explicit Name##Distribution(const param_type &param) : param_(param) \
494  { \
495  reset(); \
496  } \
497  \
498  explicit Name##Distribution(param_type &&param) \
499  : param_(std::move(param)) \
500  { \
501  reset(); \
502  } \
503  \
504  T1 p1() const { return param_.p1(); } \
505  T2 p2() const { return param_.p2(); }
506 
507 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_CONSTRUCTOR_4( \
508  Name, T, T1, p1, v1, T2, p2, v2, T3, p3, v3, T4, p4, v4) \
509  public: \
510  using result_type = T; \
511  using distribution_type = Name##Distribution<T>; \
512  \
513  explicit Name##Distribution( \
514  T1 p1 = v1, T2 p2 = v2, T3 p3 = v3, T4 p4 = v4) \
515  : param_(p1, p2, p3, p4) \
516  { \
517  reset(); \
518  } \
519  \
520  explicit Name##Distribution(const param_type &param) : param_(param) \
521  { \
522  reset(); \
523  } \
524  \
525  explicit Name##Distribution(param_type &&param) \
526  : param_(std::move(param)) \
527  { \
528  reset(); \
529  } \
530  \
531  T1 p1() const { return param_.p1(); } \
532  T2 p2() const { return param_.p2(); } \
533  T3 p3() const { return param_.p3(); } \
534  T4 p4() const { return param_.p4(); }
535 
536 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_MEMBER_0 \
537  private: \
538  bool is_equal(const distribution_type &) const { return true; } \
539  \
540  template <typename CharT, typename Traits> \
541  void ostream(std::basic_ostream<CharT, Traits> &) const \
542  { \
543  } \
544  \
545  template <typename CharT, typename Traits> \
546  void istream(std::basic_istream<CharT, Traits> &) \
547  { \
548  }
549 
550 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_MEMBER_1(T1, m1) \
551  private: \
552  T1 m1; \
553  \
554  bool is_equal(const distribution_type &other) const \
555  { \
556  if (m1 != other.m1) \
557  return false; \
558  return true; \
559  } \
560  \
561  template <typename CharT, typename Traits> \
562  void ostream(std::basic_ostream<CharT, Traits> &os) const \
563  { \
564  if (!os) \
565  return; \
566  \
567  os << ' ' << m1; \
568  } \
569  \
570  template <typename CharT, typename Traits> \
571  void istream(std::basic_istream<CharT, Traits> &is) \
572  { \
573  if (!is) \
574  return; \
575  \
576  T1 tmp1; \
577  is >> std::ws >> tmp1; \
578  if (is) { \
579  m1 = std::move(tmp1); \
580  } \
581  }
582 
583 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_MEMBER_2(T1, m1, T2, m2) \
584  private: \
585  T1 m1; \
586  T2 m2; \
587  \
588  bool is_equal(const distribution_type &other) const \
589  { \
590  if (m1 != other.m1) \
591  return false; \
592  if (m2 != other.m2) \
593  return false; \
594  return true; \
595  } \
596  \
597  template <typename CharT, typename Traits> \
598  void ostream(std::basic_ostream<CharT, Traits> &os) const \
599  { \
600  if (!os) \
601  return; \
602  \
603  os << ' ' << m1; \
604  os << ' ' << m2; \
605  } \
606  \
607  template <typename CharT, typename Traits> \
608  void istream(std::basic_istream<CharT, Traits> &is) \
609  { \
610  if (!is) \
611  return; \
612  \
613  T1 tmp1; \
614  T2 tmp2; \
615  is >> std::ws >> tmp1; \
616  is >> std::ws >> tmp2; \
617  if (is) { \
618  m1 = std::move(tmp1); \
619  m2 = std::move(tmp2); \
620  } \
621  }
622 
623 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_OPERATOR(Name, name) \
624  public: \
625  const param_type &param() const { return param_; } \
626  \
627  void param(const param_type &param) \
628  { \
629  param_ = param; \
630  reset(); \
631  } \
632  \
633  void param(param_type &&param) \
634  { \
635  param_ = std::move(param); \
636  reset(); \
637  } \
638  \
639  template <typename RNGType> \
640  result_type operator()(RNGType &rng) \
641  { \
642  return operator()(rng, param_); \
643  } \
644  \
645  template <typename RNGType> \
646  result_type operator()(RNGType &rng, const param_type &param) \
647  { \
648  return generate(rng, param); \
649  } \
650  \
651  template <typename RNGType> \
652  void operator()(RNGType &rng, std::size_t n, result_type *r) \
653  { \
654  operator()(rng, n, r, param_); \
655  } \
656  \
657  template <typename RNGType> \
658  void operator()( \
659  RNGType &rng, std::size_t n, result_type *r, const param_type &param) \
660  { \
661  name##_distribution(rng, n, r, param); \
662  } \
663  \
664  friend bool operator==( \
665  const distribution_type &dist1, const distribution_type &dist2) \
666  { \
667  return dist1.param_ == dist2.param_ && dist1.is_equal(dist2); \
668  } \
669  \
670  friend bool operator!=( \
671  const distribution_type &dist1, const distribution_type &dist2) \
672  { \
673  return !(dist1 == dist2); \
674  } \
675  \
676  template <typename CharT, typename Traits> \
677  friend std::basic_ostream<CharT, Traits> &operator<<( \
678  std::basic_ostream<CharT, Traits> &os, const distribution_type &dist) \
679  { \
680  if (!os) \
681  return os; \
682  \
683  os << dist.param_; \
684  dist.ostream(os); \
685  \
686  return os; \
687  } \
688  \
689  template <typename CharT, typename Traits> \
690  friend std::basic_istream<CharT, Traits> &operator>>( \
691  std::basic_istream<CharT, Traits> &is, distribution_type &dist) \
692  { \
693  if (!is) \
694  return is; \
695  \
696  is >> std::ws >> dist.param_; \
697  if (is) \
698  dist.istream(is); \
699  \
700  return is; \
701  } \
702  \
703  private: \
704  param_type param_;
705 
706 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_RAND(Name, T) \
707  template <typename T, typename RNGType> \
708  inline void rand(RNGType &rng, Name##Distribution<T> &distribution, \
709  std::size_t N, T *r) \
710  { \
711  distribution(rng, N, r); \
712  }
713 
714 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_0(Name, name, T) \
715  MCKL_DEFINE_RANDOM_DISTRIBUTION_PARAM_TYPE_0(Name, T) \
716  MCKL_DEFINE_RANDOM_DISTRIBUTION_CONSTRUCTOR_0(Name, T) \
717  MCKL_DEFINE_RANDOM_DISTRIBUTION_OPERATOR(Name, name)
718 
719 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_1(Name, name, T, T1, p1, v1) \
720  MCKL_DEFINE_RANDOM_DISTRIBUTION_PARAM_TYPE_1(Name, name, T, T1, p1, v1) \
721  MCKL_DEFINE_RANDOM_DISTRIBUTION_CONSTRUCTOR_1(Name, T, T1, p1, v1) \
722  MCKL_DEFINE_RANDOM_DISTRIBUTION_OPERATOR(Name, name)
723 
724 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_2( \
725  Name, name, T, T1, p1, v1, T2, p2, v2) \
726  MCKL_DEFINE_RANDOM_DISTRIBUTION_PARAM_TYPE_2( \
727  Name, name, T, T1, p1, v1, T2, p2, v2) \
728  MCKL_DEFINE_RANDOM_DISTRIBUTION_CONSTRUCTOR_2( \
729  Name, T, T1, p1, v1, T2, p2, v2) \
730  MCKL_DEFINE_RANDOM_DISTRIBUTION_OPERATOR(Name, name)
731 
732 #define MCKL_DEFINE_RANDOM_DISTRIBUTION_4( \
733  Name, name, T, T1, p1, v1, T2, p2, v2, T3, p3, v3, T4, p4, v4) \
734  MCKL_DEFINE_RANDOM_DISTRIBUTION_PARAM_TYPE_4( \
735  Name, name, T, T1, p1, v1, T2, p2, v2, T3, p3, v3, T4, p4, v4) \
736  MCKL_DEFINE_RANDOM_DISTRIBUTION_CONSTRUCTOR_4( \
737  Name, T, T1, p1, v1, T2, p2, v2, T3, p3, v3, T4, p4, v4) \
738  MCKL_DEFINE_RANDOM_DISTRIBUTION_OPERATOR(Name, name)
739 
740 #define MCKL_DEFINE_RANDOM_TEST_OPERATOR(ResultType) \
741  template <typename DistributionType> \
742  ResultType operator()(DistributionType &&distribution) \
743  { \
744  ::mckl::internal::DummyRNG rng; \
745  ::mckl::internal::DummyDistribution<DistributionType> dist( \
746  distribution); \
747  \
748  return operator()(rng, dist); \
749  }
750 
751 namespace mckl {
752 
753 namespace internal {
754 
755 class DummyRNG
756 {
757 }; // class DummyRNG
758 
759 template <typename DistributionType>
760 class DummyDistribution
761 {
762  public:
763  using result_type = double;
764 
765  DummyDistribution(const DistributionType &distribution)
766  : distribution_(distribution)
767  {
768  }
769 
770  template <typename RNGType>
771  double operator()(RNGType &)
772  {
773  return distribution_();
774  }
775 
776  private:
777  DistributionType distribution_;
778 }; // class DummyDistribution
779 
780 template <std::size_t N, std::size_t D, std::size_t T>
781 inline std::size_t serial_index(const std::size_t *, std::false_type)
782 {
783  return 0;
784 }
785 
786 template <std::size_t N, std::size_t D, std::size_t T>
787 inline std::size_t serial_index(const std::size_t *s, std::true_type)
788 {
789  return Pow<std::size_t, D, T - 1 - N>::value * s[N] +
790  serial_index<N + 1, D, T>(
791  s, std::integral_constant<bool, N + 1 < T>());
792 }
793 
794 template <std::size_t D, std::size_t T, typename ResultType>
795 inline std::size_t serial_index(const ResultType *r)
796 {
797  std::array<std::size_t, T> s;
798  for (std::size_t i = 0; i != T; ++i) {
799  s[i] = std::min(static_cast<std::size_t>(r[i]), D - 1);
800  }
801 
802  return serial_index<0, D, T>(
803  s.data(), std::integral_constant<bool, 0 < T>());
804 }
805 
806 template <std::size_t N, typename ResultType>
807 inline std::size_t permutation_imax(ResultType *r)
808 {
809  ResultType v = r[0];
810  std::size_t s = 0;
811  for (std::size_t i = 1; i < N; ++i) {
812  if (v < r[i]) {
813  v = r[i];
814  s = i;
815  }
816  }
817 
818  return s;
819 }
820 
821 template <std::size_t N, typename ResultType>
822 inline void permutation_sort(ResultType *, std::size_t &, std::false_type)
823 {
824 }
825 
826 template <std::size_t N, typename ResultType>
827 inline void permutation_sort(ResultType *r, std::size_t &f, std::true_type)
828 {
829  std::size_t s = permutation_imax<N>(r);
830  f = N * f + s;
831  r[s] = r[N - 1];
832  permutation_sort<N - 1>(r, f, std::integral_constant<bool, 2 < N>());
833 }
834 
835 template <std::size_t T, typename ResultType>
836 inline std::size_t permutation_index(ResultType *r)
837 {
838  std::size_t f = 0;
839  permutation_sort<T>(r, f, std::integral_constant<bool, 1 < T>());
840 
841  return f;
842 }
843 
844 inline void group_np(double n, const Vector<double> &np_all,
845  Vector<double> &np, std::size_t &kmin, std::size_t &kmax)
846 {
847  double npmin = 10;
848 
849  std::size_t K = np_all.size();
850  if (K <= 2) {
851  np = np_all;
852  kmin = 0;
853  kmax = K - 1;
854  return;
855  }
856 
857  kmin = 0;
858  while (kmin < K) {
859  if (np_all[kmin] >= npmin) {
860  break;
861  }
862  ++kmin;
863  }
864  kmax = kmin + 1;
865  while (kmax < K && np_all[kmax] >= npmin) {
866  ++kmax;
867  }
868  --kmax;
869  if (kmin == kmax) {
870  --kmin;
871  }
872  np.clear();
873  np.reserve(kmax - kmin + 1);
874  np.push_back(0);
875  std::size_t k = 0;
876  while (k <= kmin) {
877  np[0] += np_all[k];
878  ++k;
879  }
880  while (k <= kmax) {
881  np.push_back(np_all[k]);
882  ++k;
883  }
884  np.back() += n - std::accumulate(np.begin(), np.end(), 0.0);
885 }
886 
887 template <typename IntType, IntType D, typename RealType>
888 inline IntType ftoi(RealType r)
889 {
890  std::size_t u = static_cast<IntType>(r);
891  u = u < D ? u : D - 1;
892 
893  return u;
894 }
895 
896 template <typename IntType, typename RealType>
897 inline IntType ftoi(RealType x, std::true_type)
898 {
899  constexpr RealType maxval =
900  static_cast<RealType>(std::numeric_limits<IntType>::max());
901 
902  return static_cast<IntType>(std::min(maxval, x));
903 }
904 
905 template <typename IntType, typename RealType>
906 inline IntType ftoi(RealType x, std::false_type)
907 {
908  constexpr RealType maxval =
909  static_cast<RealType>(std::numeric_limits<IntType>::max() / 2);
910 
911  return static_cast<IntType>(std::min(maxval, x));
912 }
913 
914 template <typename IntType, typename RealType>
915 inline IntType ftoi(RealType x)
916 {
917  constexpr int W = std::numeric_limits<IntType>::digits;
918  constexpr int M = std::numeric_limits<RealType>::digits;
919 
920  return ftoi<IntType>(x, std::integral_constant<bool, W <= M>());
921 }
922 
923 template <typename T, std::size_t, std::size_t,
924  int = std::numeric_limits<T>::digits>
925 class IncrementBlockSI128;
926 
927 template <typename T, std::size_t, std::size_t,
928  int = std::numeric_limits<T>::digits>
929 class IncrementBlockSI256;
930 
931 template <typename T, std::size_t, std::size_t,
932  int = std::numeric_limits<T>::digits>
933 class IncrementBlockSI512;
934 
935 } // namespace internal
936 
939 template <typename RNGType>
940 class RNGTraits
941 {
942  using result_type = typename RNGType::result_type;
943 
944  static constexpr typename RNGType::result_type R_ = RNGType::min() == 0 &&
945  RNGType::max() == std::numeric_limits<result_type>::max() ?
946  0 :
947  RNGType::max() - RNGType::min() + const_one<result_type>();
948 
949  public:
952  static constexpr bool is_full_range =
953  RNGType::min() == 0 && (R_ & (R_ - 1)) == 0;
954 
955  static constexpr int bits = internal::Log2<result_type, R_>::value;
956 }; // class RNGTraits;
957 
960 class Closed;
961 
964 class Open;
965 
968 template <typename RNGType>
969 inline typename RNGType::result_type rand(RNGType &rng)
970 {
971  return rng();
972 }
973 
976 template <typename RNGType>
977 inline void rand(RNGType &rng, std::size_t n, typename RNGType::result_type *r)
978 {
979  for (std::size_t i = 0; i != n; ++i) {
980  r[i] = rng();
981  }
982 }
983 
986 template <typename RNGType, typename DistributionType>
987 inline typename DistributionType::result_type rand(
988  RNGType &rng, DistributionType &distribution)
989 {
990  return distribution(rng);
991 }
992 
995 template <typename RNGType, typename DistributionType>
996 inline void rand(RNGType &rng, DistributionType &distribution, std::size_t n,
997  typename DistributionType::result_type *r)
998 {
999  for (std::size_t i = 0; i != n; ++i) {
1000  r[i] = distribution(rng);
1001  }
1002 }
1003 
1004 template <typename>
1005 class SeedTrait
1006 {
1007  public:
1008  using type = unsigned;
1009 }; // class SeedTrait
1010 
1011 template <typename RNGType>
1012 using SeedType = typename SeedTrait<RNGType>::type;
1013 
1014 template <typename = double>
1015 class ArcsineDistribution;
1016 
1017 template <typename = double>
1018 class BetaDistribution;
1019 
1020 template <typename = double>
1021 class CauchyDistribution;
1022 
1023 template <typename = double>
1024 class ChiSquaredDistribution;
1025 
1026 template <typename = double>
1027 class DirichletDistribution;
1028 
1029 template <typename = double>
1030 class ExponentialDistribution;
1031 
1032 template <typename = double>
1033 class ExtremeValueDistribution;
1034 
1035 template <typename = double>
1036 class FisherFDistribution;
1037 
1038 template <typename = double>
1039 class GammaDistribution;
1040 
1041 template <typename = double>
1042 class LaplaceDistribution;
1043 
1044 template <typename = double>
1045 class LevyDistribution;
1046 
1047 template <typename = double>
1048 class LogisticDistribution;
1049 
1050 template <typename = double>
1051 class LognormalDistribution;
1052 
1053 template <typename = double>
1054 class NormalDistribution;
1055 
1056 template <typename = double>
1057 class NormalMVDistribution;
1058 
1059 template <typename = double>
1060 class ParetoDistribution;
1061 
1062 template <typename = double>
1063 class RayleighDistribution;
1064 
1065 template <typename = double>
1066 class StableDistribution;
1067 
1068 template <typename = double>
1069 class StudentTDistribution;
1070 
1071 template <typename = double>
1072 class U01CanonicalDistribution;
1073 
1074 template <typename = double>
1075 class U01CCDistribution;
1076 
1077 template <typename = double>
1078 class U01CODistribution;
1079 
1080 template <typename = double>
1081 class U01OCDistribution;
1082 
1083 template <typename = double>
1084 class U01OODistribution;
1085 
1086 template <typename = unsigned>
1087 class UniformBitsDistribution;
1088 
1089 template <typename = double>
1090 class UniformRealDistribution;
1091 
1092 template <typename = double>
1093 class WeibullDistribution;
1094 
1095 template <typename = bool>
1096 class BernoulliDistribution;
1097 
1098 template <typename = int>
1099 class GeometricDistribution;
1100 
1101 template <typename = int>
1102 class DiscreteDistribution;
1103 
1104 template <typename = int>
1105 class SamplingDistribution;
1106 
1107 template <typename = int>
1108 class UniformIntDistribution;
1109 
1110 } // namespace mckl
1111 
1112 #endif // MCKL_RANDOM_INTERNAL_COMMON_HPP
mckl::NormalMVDistribution
Multivariate Normal distribution.
Definition: random/internal/common.hpp:1057
mckl::RNGTraits
Traits of RNG engines.
Definition: random/internal/common.hpp:940
mckl::U01CODistribution
Standard uniform distribution on [0, 1)
Definition: random/internal/common.hpp:1078
mckl::SamplingDistribution
Drawing a subset without replacement.
Definition: random/internal/common.hpp:1105
mckl::internal::permutation_imax
std::size_t permutation_imax(ResultType *r)
Definition: random/internal/common.hpp:807
mckl::internal::ftoi
IntType ftoi(RealType r)
Definition: random/internal/common.hpp:888
mckl::UniformIntDistribution
Uniform integer distribution.
Definition: random/internal/common.hpp:1108
mckl::internal::permutation_index
std::size_t permutation_index(ResultType *r)
Definition: random/internal/common.hpp:836
mckl::internal::group_np
void group_np(double n, const Vector< double > &np_all, Vector< double > &np, std::size_t &kmin, std::size_t &kmax)
Definition: random/internal/common.hpp:844
mckl::Vector
std::vector< T, Alloc > Vector
std::vector with Allocator as the default allocator
Definition: memory.hpp:435
mckl::U01OCDistribution
Standard uniform distribution on (0, 1].
Definition: random/internal/common.hpp:1081
mckl::internal::serial_index
std::size_t serial_index(const std::size_t *, std::false_type)
Definition: random/internal/common.hpp:781
mckl::U01OODistribution
Standard uniform distribution on (0, 1)
Definition: random/internal/common.hpp:1084
mckl
Definition: mcmc.hpp:40
mckl::rand
void rand(RNGType &rng, ArcsineDistribution< RealType > &distribution, std::size_t N, RealType *r)
Definition: arcsine_distribution.hpp:93
mckl::U01CanonicalDistribution
Standard uniform distribution on [0, 1)
Definition: random/internal/common.hpp:1072
mckl::internal::permutation_sort
void permutation_sort(ResultType *, std::size_t &, std::false_type)
Definition: random/internal/common.hpp:822
mckl::DiscreteDistribution
Draw a single sample given weights.
Definition: discrete_distribution.hpp:43
mckl::U01CCDistribution
Standard uniform distribution on [0, 1].
Definition: random/internal/common.hpp:1075
mckl::internal::DummyDistribution::DummyDistribution
DummyDistribution(const DistributionType &distribution)
Definition: random/internal/common.hpp:765
mckl::SeedType
typename SeedTrait< RNGType >::type SeedType
Definition: random/internal/common.hpp:1012
Generated by   doxygen 1.8.13