diff options
Diffstat (limited to 'extern/carve/patches/files/random.hpp')
| -rw-r--r-- | extern/carve/patches/files/random.hpp | 773 |
1 files changed, 773 insertions, 0 deletions
diff --git a/extern/carve/patches/files/random.hpp b/extern/carve/patches/files/random.hpp new file mode 100644 index 00000000000..15866bb496e --- /dev/null +++ b/extern/carve/patches/files/random.hpp @@ -0,0 +1,773 @@ +#pragma once + +#include <iostream> +#include <vector> +#include <limits> +#include <stdexcept> +#include <cmath> +#include <algorithm> + +#if !defined(_MSC_VER) +#include <stdint.h> +#endif + +namespace boost { + +// type_traits could help here, but I don't want to depend on type_traits. +template<class T> +struct ptr_helper +{ + typedef T value_type; + typedef T& reference_type; + typedef const T& rvalue_type; + static reference_type ref(T& r) { return r; } + static const T& ref(const T& r) { return r; } +}; + +template<class T> +struct ptr_helper<T&> +{ + typedef T value_type; + typedef T& reference_type; + typedef T& rvalue_type; + static reference_type ref(T& r) { return r; } + static const T& ref(const T& r) { return r; } +}; + +template<class T> +struct ptr_helper<T*> +{ + typedef T value_type; + typedef T& reference_type; + typedef T* rvalue_type; + static reference_type ref(T * p) { return *p; } + static const T& ref(const T * p) { return *p; } +}; + +template<class UniformRandomNumberGenerator> +class pass_through_engine +{ +private: + typedef ptr_helper<UniformRandomNumberGenerator> helper_type; + +public: + typedef typename helper_type::value_type base_type; + typedef typename base_type::result_type result_type; + + explicit pass_through_engine(UniformRandomNumberGenerator rng) + // make argument an rvalue to avoid matching Generator& constructor + : _rng(static_cast<typename helper_type::rvalue_type>(rng)) + { } + + result_type min () const { return (base().min)(); } + result_type max () const { return (base().max)(); } + base_type& base() { return helper_type::ref(_rng); } + const base_type& base() const { return helper_type::ref(_rng); } + + result_type operator()() { return base()(); } + +private: + UniformRandomNumberGenerator _rng; +}; + +template<class RealType> +class new_uniform_01 +{ +public: + typedef RealType input_type; + typedef RealType result_type; + // compiler-generated copy ctor and copy assignment are fine + result_type min () const { return result_type(0); } + result_type max () const { return result_type(1); } + void reset() { } + + template<class Engine> + result_type operator()(Engine& eng) { + for (;;) { + typedef typename Engine::result_type base_result; + result_type factor = result_type(1) / + (result_type((eng.max)()-(eng.min)()) + + result_type(std::numeric_limits<base_result>::is_integer ? 1 : 0)); + result_type result = result_type(eng() - (eng.min)()) * factor; + if (result < result_type(1)) + return result; + } + } + + template<class CharT, class Traits> + friend std::basic_ostream<CharT,Traits>& + operator<<(std::basic_ostream<CharT,Traits>& os, const new_uniform_01&) + { + return os; + } + + template<class CharT, class Traits> + friend std::basic_istream<CharT,Traits>& + operator>>(std::basic_istream<CharT,Traits>& is, new_uniform_01&) + { + return is; + } +}; + +template<class UniformRandomNumberGenerator, class RealType> +class backward_compatible_uniform_01 +{ + typedef ptr_helper<UniformRandomNumberGenerator> traits; + typedef pass_through_engine<UniformRandomNumberGenerator> internal_engine_type; +public: + typedef UniformRandomNumberGenerator base_type; + typedef RealType result_type; + + static const bool has_fixed_range = false; + + explicit backward_compatible_uniform_01(typename traits::rvalue_type rng) + : _rng(rng), + _factor(result_type(1) / + (result_type((_rng.max)()-(_rng.min)()) + + result_type(std::numeric_limits<base_result>::is_integer ? 1 : 0))) + { + } + // compiler-generated copy ctor and copy assignment are fine + + result_type min () const { return result_type(0); } + result_type max () const { return result_type(1); } + typename traits::value_type& base() { return _rng.base(); } + const typename traits::value_type& base() const { return _rng.base(); } + void reset() { } + + result_type operator()() { + for (;;) { + result_type result = result_type(_rng() - (_rng.min)()) * _factor; + if (result < result_type(1)) + return result; + } + } + + template<class CharT, class Traits> + friend std::basic_ostream<CharT,Traits>& + operator<<(std::basic_ostream<CharT,Traits>& os, const backward_compatible_uniform_01& u) + { + os << u._rng; + return os; + } + + template<class CharT, class Traits> + friend std::basic_istream<CharT,Traits>& + operator>>(std::basic_istream<CharT,Traits>& is, backward_compatible_uniform_01& u) + { + is >> u._rng; + return is; + } + +private: + typedef typename internal_engine_type::result_type base_result; + internal_engine_type _rng; + result_type _factor; +}; + +// A definition is required even for integral static constants +template<class UniformRandomNumberGenerator, class RealType> +const bool backward_compatible_uniform_01<UniformRandomNumberGenerator, RealType>::has_fixed_range; + +template<class UniformRandomNumberGenerator> +struct select_uniform_01 +{ + template<class RealType> + struct apply + { + typedef backward_compatible_uniform_01<UniformRandomNumberGenerator, RealType> type; + }; +}; + +template<> +struct select_uniform_01<float> +{ + template<class RealType> + struct apply + { + typedef new_uniform_01<float> type; + }; +}; + +template<> +struct select_uniform_01<double> +{ + template<class RealType> + struct apply + { + typedef new_uniform_01<double> type; + }; +}; + +template<> +struct select_uniform_01<long double> +{ + template<class RealType> + struct apply + { + typedef new_uniform_01<long double> type; + }; +}; + +// Because it is so commonly used: uniform distribution on the real [0..1) +// range. This allows for specializations to avoid a costly int -> float +// conversion plus float multiplication +template<class UniformRandomNumberGenerator = double, class RealType = double> +class uniform_01 + : public select_uniform_01<UniformRandomNumberGenerator>::template apply<RealType>::type +{ + typedef typename select_uniform_01<UniformRandomNumberGenerator>::template apply<RealType>::type impl_type; + typedef ptr_helper<UniformRandomNumberGenerator> traits; +public: + + uniform_01() {} + + explicit uniform_01(typename traits::rvalue_type rng) + : impl_type(rng) + { + } + + template<class CharT, class Traits> + friend std::basic_ostream<CharT,Traits>& + operator<<(std::basic_ostream<CharT,Traits>& os, const uniform_01& u) + { + os << static_cast<const impl_type&>(u); + return os; + } + + template<class CharT, class Traits> + friend std::basic_istream<CharT,Traits>& + operator>>(std::basic_istream<CharT,Traits>& is, uniform_01& u) + { + is >> static_cast<impl_type&>(u); + return is; + } +}; + +template<class UniformRandomNumberGenerator, class IntType = unsigned long> +class uniform_int_float +{ +public: + typedef UniformRandomNumberGenerator base_type; + typedef IntType result_type; + + uniform_int_float(base_type rng, IntType min_arg = 0, IntType max_arg = 0xffffffff) + : _rng(rng), _min(min_arg), _max(max_arg) + { + init(); + } + + result_type min () const { return _min; } + result_type max () const { return _max; } + base_type& base() { return _rng.base(); } + const base_type& base() const { return _rng.base(); } + + result_type operator()() + { + return static_cast<IntType>(_rng() * _range) + _min; + } + + template<class CharT, class Traits> + friend std::basic_ostream<CharT,Traits>& + operator<<(std::basic_ostream<CharT,Traits>& os, const uniform_int_float& ud) + { + os << ud._min << " " << ud._max; + return os; + } + + template<class CharT, class Traits> + friend std::basic_istream<CharT,Traits>& + operator>>(std::basic_istream<CharT,Traits>& is, uniform_int_float& ud) + { + is >> std::ws >> ud._min >> std::ws >> ud._max; + ud.init(); + return is; + } + +private: + void init() + { + _range = static_cast<base_result>(_max-_min)+1; + } + + typedef typename base_type::result_type base_result; + uniform_01<base_type> _rng; + result_type _min, _max; + base_result _range; +}; + + +template<class UniformRandomNumberGenerator, class CharT, class Traits> +std::basic_ostream<CharT,Traits>& +operator<<( + std::basic_ostream<CharT,Traits>& os + , const pass_through_engine<UniformRandomNumberGenerator>& ud + ) +{ + return os << ud.base(); +} + +template<class UniformRandomNumberGenerator, class CharT, class Traits> +std::basic_istream<CharT,Traits>& +operator>>( + std::basic_istream<CharT,Traits>& is + , const pass_through_engine<UniformRandomNumberGenerator>& ud + ) +{ + return is >> ud.base(); +} + + + +template<class RealType = double> +class normal_distribution +{ +public: + typedef RealType input_type; + typedef RealType result_type; + + explicit normal_distribution(const result_type& mean_arg = result_type(0), + const result_type& sigma_arg = result_type(1)) + : _mean(mean_arg), _sigma(sigma_arg), _valid(false) + { + //assert(_sigma >= result_type(0)); + } + + // compiler-generated copy constructor is NOT fine, need to purge cache + normal_distribution(const normal_distribution& other) + : _mean(other._mean), _sigma(other._sigma), _valid(false) + { + } + + // compiler-generated copy ctor and assignment operator are fine + + RealType mean() const { return _mean; } + RealType sigma() const { return _sigma; } + + void reset() { _valid = false; } + + template<class Engine> + result_type operator()(Engine& eng) + { +#ifndef BOOST_NO_STDC_NAMESPACE + // allow for Koenig lookup + using std::sqrt; using std::log; using std::sin; using std::cos; +#endif + if(!_valid) { + _r1 = eng(); + _r2 = eng(); + _cached_rho = sqrt(-result_type(2) * log(result_type(1)-_r2)); + _valid = true; + } else { + _valid = false; + } + // Can we have a boost::mathconst please? + const result_type pi = result_type(3.14159265358979323846); + + return _cached_rho * (_valid ? + cos(result_type(2)*pi*_r1) : + sin(result_type(2)*pi*_r1)) + * _sigma + _mean; + } + +#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS + template<class CharT, class Traits> + friend std::basic_ostream<CharT,Traits>& + operator<<(std::basic_ostream<CharT,Traits>& os, const normal_distribution& nd) + { + os << nd._mean << " " << nd._sigma << " " + << nd._valid << " " << nd._cached_rho << " " << nd._r1; + return os; + } + + template<class CharT, class Traits> + friend std::basic_istream<CharT,Traits>& + operator>>(std::basic_istream<CharT,Traits>& is, normal_distribution& nd) + { + is >> std::ws >> nd._mean >> std::ws >> nd._sigma + >> std::ws >> nd._valid >> std::ws >> nd._cached_rho + >> std::ws >> nd._r1; + return is; + } +#endif +private: + result_type _mean, _sigma; + result_type _r1, _r2, _cached_rho; + bool _valid; +}; + +// http://www.math.keio.ac.jp/matumoto/emt.html +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +class mersenne_twister +{ +public: + typedef UIntType result_type; + static const int word_size = w; + static const int state_size = n; + static const int shift_size = m; + static const int mask_bits = r; + static const UIntType parameter_a = a; + static const int output_u = u; + static const int output_s = s; + static const UIntType output_b = b; + static const int output_t = t; + static const UIntType output_c = c; + static const int output_l = l; + + static const bool has_fixed_range = false; + + mersenne_twister() { seed(); } + + explicit mersenne_twister(const UIntType& value) + { seed(value); } + template<class It> mersenne_twister(It& first, It last) { seed(first,last); } + + template<class Generator> \ + explicit mersenne_twister(Generator& gen) + { seed(gen); } + + // compiler-generated copy ctor and assignment operator are fine + + void seed() { seed(UIntType(5489)); } + + void seed(const UIntType& value) + { + // New seeding algorithm from + // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html + // In the previous versions, MSBs of the seed affected only MSBs of the + // state x[]. + const UIntType mask = ~0u; + x[0] = value & mask; + for (i = 1; i < n; i++) { + // See Knuth "The Art of Computer Programming" Vol. 2, 3rd ed., page 106 + x[i] = (1812433253UL * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask; + } + } + + // For GCC, moving this function out-of-line prevents inlining, which may + // reduce overall object code size. However, MSVC does not grok + // out-of-line definitions of member function templates. + template<class Generator> \ + void seed(Generator& gen) + { +/*#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS + BOOST_STATIC_ASSERT(!std::numeric_limits<result_type>::is_signed); +#endif*/ + // I could have used std::generate_n, but it takes "gen" by value + for(int j = 0; j < n; j++) + x[j] = gen(); + i = n; + } + + template<class It> + void seed(It& first, It last) + { + int j; + for(j = 0; j < n && first != last; ++j, ++first) + x[j] = *first; + i = n; + if(first == last && j < n) + throw std::invalid_argument("mersenne_twister::seed"); + } + + result_type min () const { return 0; } + result_type max () const + { + // avoid "left shift count >= with of type" warning + result_type res = 0; + for(int j = 0; j < w; ++j) + res |= (1u << j); + return res; + } + + result_type operator()(); + static bool validation(result_type v) { return val == v; } + + template<class CharT, class Traits> + friend std::basic_ostream<CharT,Traits>& + operator<<(std::basic_ostream<CharT,Traits>& os, const mersenne_twister& mt) + { + for(int j = 0; j < mt.state_size; ++j) + os << mt.compute(j) << " "; + return os; + } + + template<class CharT, class Traits> + friend std::basic_istream<CharT,Traits>& + operator>>(std::basic_istream<CharT,Traits>& is, mersenne_twister& mt) + { + for(int j = 0; j < mt.state_size; ++j) + is >> mt.x[j] >> std::ws; + // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template + // value parameter "n" available from the class template scope, so use + // the static constant with the same value + mt.i = mt.state_size; + return is; + } + + friend bool operator==(const mersenne_twister& x, const mersenne_twister& y) + { + for(int j = 0; j < state_size; ++j) + if(x.compute(j) != y.compute(j)) + return false; + return true; + } + + friend bool operator!=(const mersenne_twister& x, const mersenne_twister& y) + { return !(x == y); } + +private: + // returns x(i-n+index), where index is in 0..n-1 + UIntType compute(unsigned int index) const + { + // equivalent to (i-n+index) % 2n, but doesn't produce negative numbers + return x[ (i + n + index) % (2*n) ]; + } + void twist(int block); + + // state representation: next output is o(x(i)) + // x[0] ... x[k] x[k+1] ... x[n-1] x[n] ... x[2*n-1] represents + // x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) x(i-k-n) ... x[i(i-k-1)] + // The goal is to always have x(i-n) ... x(i-1) available for + // operator== and save/restore. + + UIntType x[2*n]; + int i; +}; + +// A definition is required even for integral static constants + +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +const bool mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::has_fixed_range; +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::state_size; +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::shift_size; +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::mask_bits; +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +const UIntType mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::parameter_a; +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_u; +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_s; +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +const UIntType mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_b; +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_t; +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +const UIntType mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_c; +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_l; + +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +void mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::twist(int block) +{ + const UIntType upper_mask = (~0u) << r; + const UIntType lower_mask = ~upper_mask; + + if(block == 0) { + for(int j = n; j < 2*n; j++) { + UIntType y = (x[j-n] & upper_mask) | (x[j-(n-1)] & lower_mask); + x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0); + } + } else if (block == 1) { + // split loop to avoid costly modulo operations + { // extra scope for MSVC brokenness w.r.t. for scope + for(int j = 0; j < n-m; j++) { + UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask); + x[j] = x[j+n+m] ^ (y >> 1) ^ (y&1 ? a : 0); + } + } + + for(int j = n-m; j < n-1; j++) { + UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask); + x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0); + } + // last iteration + UIntType y = (x[2*n-1] & upper_mask) | (x[0] & lower_mask); + x[n-1] = x[m-1] ^ (y >> 1) ^ (y&1 ? a : 0); + i = 0; + } +} + +template<class UIntType, int w, int n, int m, int r, UIntType a, int u, + int s, UIntType b, int t, UIntType c, int l, UIntType val> +inline typename mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::result_type +mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::operator()() +{ + if(i == n) + twist(0); + else if(i >= 2*n) + twist(1); + // Step 4 + UIntType z = x[i]; + ++i; + z ^= (z >> u); + z ^= ((z << s) & b); + z ^= ((z << t) & c); + z ^= (z >> l); + return z; +} + +typedef mersenne_twister<uint32_t,32,351,175,19,0xccab8ee7,11, + 7,0x31b6ab00,15,0xffe50000,17, 0xa37d3c92> mt11213b; + +// validation by experiment from mt19937.c +typedef mersenne_twister<uint32_t,32,624,397,31,0x9908b0df,11, + 7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937; + + +template<class RealType = double, class Cont = std::vector<RealType> > +class uniform_on_sphere +{ +public: + typedef RealType input_type; + typedef Cont result_type; + + explicit uniform_on_sphere(int dim = 2) : _container(dim), _dim(dim) { } + + // compiler-generated copy ctor and assignment operator are fine + + void reset() { _normal.reset(); } + + template<class Engine> + const result_type & operator()(Engine& eng) + { + RealType sqsum = 0; + for(typename Cont::iterator it = _container.begin(); + it != _container.end(); + ++it) { + RealType val = _normal(eng); + *it = val; + sqsum += val * val; + } + using std::sqrt; + // for all i: result[i] /= sqrt(sqsum) + std::transform(_container.begin(), _container.end(), _container.begin(), + std::bind2nd(std::divides<RealType>(), sqrt(sqsum))); + return _container; + } + + template<class CharT, class Traits> + friend std::basic_ostream<CharT,Traits>& + operator<<(std::basic_ostream<CharT,Traits>& os, const uniform_on_sphere& sd) + { + os << sd._dim; + return os; + } + + template<class CharT, class Traits> + friend std::basic_istream<CharT,Traits>& + operator>>(std::basic_istream<CharT,Traits>& is, uniform_on_sphere& sd) + { + is >> std::ws >> sd._dim; + sd._container.resize(sd._dim); + return is; + } + +private: + normal_distribution<RealType> _normal; + result_type _container; + int _dim; +}; + + + +template<bool have_int, bool want_int> +struct engine_helper; + + +template<> +struct engine_helper<true, true> +{ + template<class Engine, class DistInputType> + struct impl + { + typedef pass_through_engine<Engine> type; + }; +}; + +template<> +struct engine_helper<false, false> +{ + template<class Engine, class DistInputType> + struct impl + { + typedef uniform_01<Engine, DistInputType> type; + }; +}; + +template<> +struct engine_helper<true, false> +{ + template<class Engine, class DistInputType> + struct impl + { + typedef uniform_01<Engine, DistInputType> type; + }; +}; + +template<> +struct engine_helper<false, true> +{ + template<class Engine, class DistInputType> + struct impl + { + typedef uniform_int_float<Engine, unsigned long> type; + }; +}; + +template<class Engine, class Distribution> +class variate_generator +{ +private: + typedef pass_through_engine<Engine> decorated_engine; + +public: + typedef typename decorated_engine::base_type engine_value_type; + typedef Engine engine_type; + typedef Distribution distribution_type; + typedef typename Distribution::result_type result_type; + + variate_generator(Engine e, Distribution d) + : _eng(decorated_engine(e)), _dist(d) { } + + result_type operator()() { return _dist(_eng); } + template<class T> + result_type operator()(T value) { return _dist(_eng, value); } + + engine_value_type& engine() { return _eng.base().base(); } + const engine_value_type& engine() const { return _eng.base().base(); } + + distribution_type& distribution() { return _dist; } + const distribution_type& distribution() const { return _dist; } + + result_type min () const { return (distribution().min)(); } + result_type max () const { return (distribution().max)(); } + +private: + enum { + have_int = std::numeric_limits<typename decorated_engine::result_type>::is_integer, + want_int = std::numeric_limits<typename Distribution::input_type>::is_integer + }; + typedef typename engine_helper<have_int, want_int>::template impl<decorated_engine, typename Distribution::input_type>::type internal_engine_type; + + internal_engine_type _eng; + distribution_type _dist; +}; + +} // namespace boost |