Carve booleans library integration

==================================

Merging Carve library integration project into the trunk.

This commit switches Boolean modifier to another library which handles
mesh boolean operations in much stable and faster way, resolving old
well-known limitations of intern boolop library.

Carve is integrating as alternative interface for boolop library and
which makes it totally transparent for blender sources to switch between
old-fashioned boolop and new Carve backends.

Detailed changes in this commit:

- Integrated needed subset of Carve library sources into extern/
  Added script for re-bundling it (currently works only if repo
  was cloned by git-svn).
- Added BOP_CarveInterface for boolop library which can be used by
  Boolean modifier.
- Carve backend is enabled by default, can be disabled by WITH_BF_CARVE
  SCons option and WITH_CARVE CMake option.
- If Boost library is found in build environment it'll be used for
  unordered collections. If Boost isn't found, it'll fallback to TR1
  implementation for GCC compilers. Boost is obligatory if MSVC is used.

Tested on Linux 64bit and Windows 7 64bit.

NOTE: behavior of flat objects was changed. E.g. Plane-Sphere now gives
      plane with circle hole, not plane with semisphere. Don't think
      it's really issue because it's not actually defined behavior in
      such situations and both of ways might be useful. Since it's
      only known "regression" think it's OK to deal with it.

Details are there http://wiki.blender.org/index.php/User:Nazg-gul/CarveBooleans

Special thanks to:

- Ken Hughes: author of original carve integration patch.
- Campbell Barton: help in project development, review tests.
- Tobias Sargeant: author of Carve library, help in resolving some
                   merge stoppers, bug fixing.

2 files changed, 785 insertions, 0 deletions

diff --git a/extern/carve/patches/files/config.h b/extern/carve/patches/files/config.h
new file mode 100644
index 00000000000..fdae2d2843f
--- /dev/null
+++ b/extern/carve/patches/files/config.h
@@ -0,0 +1,12 @@
+#define CARVE_VERSION "2.0.0a"
+
+#undef CARVE_DEBUG
+#undef CARVE_DEBUG_WRITE_PLY_DATA
+
+#if defined(__GNUC__)
+#  if !defined(HAVE_BOOST_UNORDERED_COLLECTIONS)
+#    define HAVE_TR1_UNORDERED_COLLECTIONS
+#  endif
+
+#  define HAVE_STDINT_H
+#endif
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