[Boost-commit] svn:boost r64219 - in sandbox/SOC/2010/quasi_random: boost/random boost/random/detail libs/random/doc libs/random/example libs/random/test

Subject: [Boost-commit] svn:boost r64219 - in sandbox/SOC/2010/quasi_random: boost/random boost/random/detail libs/random/doc libs/random/example libs/random/test
From: jvd_at_[hidden]
Date: 2010-07-21 09:27:57

Author: qrng
Date: 2010-07-21 09:27:55 EDT (Wed, 21 Jul 2010)
New Revision: 64219
URL: http://svn.boost.org/trac/boost/changeset/64219

Log:
Quasi-random number generators and tests

Added:
   sandbox/SOC/2010/quasi_random/boost/random/detail/gray_coded_qrng_base.hpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/boost/random/detail/qrng_base.hpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/boost/random/faure.hpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/boost/random/niederreiter_base2.hpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/boost/random/sobol.hpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/libs/random/doc/Jamfile.diff   (contents, props changed)
   sandbox/SOC/2010/quasi_random/libs/random/example/qrng_images.cpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/libs/random/test/discrepancy.cpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/libs/random/test/faure_validate.cpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/libs/random/test/gsl_validate.cpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/libs/random/test/qrng_typedefs.hpp   (contents, props changed)

Added: sandbox/SOC/2010/quasi_random/boost/random/detail/gray_coded_qrng_base.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/boost/random/detail/gray_coded_qrng_base.hpp	2010-07-21 09:27:55 EDT (Wed, 21 Jul 2010)
@@ -0,0 +1,134 @@
+/* boost random/detail/gray_coded_qrng_base.hpp header file
+ *
+ * Copyright Justinas Vygintas Daugmaudis 2010
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#ifndef BOOST_RANDOM_DETAIL_GRAY_CODED_QRNG_BASE_HPP
+#define BOOST_RANDOM_DETAIL_GRAY_CODED_QRNG_BASE_HPP
+
+#include <boost/random/detail/qrng_base.hpp>
+
+//!\file
+//!Describes the gray-coded quasi-random number generator base class template.
+
+namespace boost {
+namespace random {
+
+namespace detail {
+
+template<typename DerivedT, typename LatticeT>
+class gray_coded_qrng_base
+  : public qrng_base<DerivedT, const LatticeT> // immutable lattice!
+{
+private:
+
+  typedef gray_coded_qrng_base<DerivedT, LatticeT> self_t;
+  typedef qrng_base<DerivedT, const LatticeT> base_t;
+
+  friend class qrng_base<DerivedT, const LatticeT>;
+
+public:
+
+  typedef typename LatticeT::result_type result_type;
+
+  gray_coded_qrng_base()
+    : base_t(std::size_t())
+  {}
+
+  explicit gray_coded_qrng_base(std::size_t init)
+    : base_t(init)
+  {}
+
+  // default copy c-tor is fine
+
+  //!Requirements: *this is mutable.
+  //!
+  //!Effects: *this and rhs are equivalent as they both return
+  //!an identical sequence of quasi-random numbers.
+  //!\code
+  //!X u; u = rhs;
+  //!assert(u() == rhs());
+  //!\endcode
+  //!
+  //!Returns: *this.
+  //!
+  //!Throws: nothing.
+  self_t& operator=(const self_t& rhs)
+  {
+    this->curr_elem = rhs.curr_elem;
+    this->seq_count = rhs.seq_count;
+    std::copy(rhs.quasi_state, rhs.quasi_state + base_t::dimension_value, this->quasi_state);
+    // we do not copy lattice here!
+    return *this;
+  }
+
+protected:
+  void reset_state()
+  {
+    this->curr_elem = 0;
+    this->seq_count = 0;
+    std::fill(this->quasi_state, this->quasi_state + base_t::dimension_value, 0);
+  }
+
+  void seed(std::size_t init, const char *msg)
+  {
+    this->curr_elem = 0;
+    if ( init != this->seq_count )
+    {
+      base_t::derived().seed();
+
+      this->seq_count = init;
+      init ^= (init / 2);
+      for( int r = 0; init != 0; ++r, init >>= 1 )
+      {
+        if( init & 1 )
+          update_quasi(r, msg);
+      }
+    }
+  }
+
+  void compute_current_vector(const char *msg)
+  {
+    update_state(this->seq_count, msg);
+  }
+
+private:
+
+  void compute_next_vector()
+  {
+    update_state(LatticeT::next(this->seq_count),
+        "compute_next_vector");
+    this->seq_count++;
+  }
+
+  void update_state(std::size_t cnt, const char *msg)
+  {
+    // Find the position of the least-significant zero in sequence count.
+    // This is the bit that changes in the Gray-code representation as
+    // the count is advanced.
+    int r = 0;
+    for( ; cnt & 1; cnt >>= 1 ) {
+      ++r;
+    }
+    update_quasi(r, msg);
+  }
+
+  void update_quasi(int r, const char* msg)
+  {
+    if(r >= LatticeT::bit_count)
+      boost::throw_exception( std::overflow_error(msg) );
+
+    // Calculate the next state.
+    for(std::size_t i = 0; i != base_t::dimension_value; ++i)
+      this->quasi_state[i] ^= this->lattice(r, i);
+  }
+};
+
+}} // namespace detail::random
+
+} // namespace boost
+
+#endif // BOOST_RANDOM_DETAIL_GRAY_CODED_QRNG_BASE_HPP

Added: sandbox/SOC/2010/quasi_random/boost/random/detail/qrng_base.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/boost/random/detail/qrng_base.hpp	2010-07-21 09:27:55 EDT (Wed, 21 Jul 2010)
@@ -0,0 +1,169 @@
+/* boost random/detail/quasi_random_number_generator_base.hpp header file
+ *
+ * Copyright Justinas Vygintas Daugmaudis 2010
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#ifndef BOOST_RANDOM_DETAIL_QRNG_BASE_HPP
+#define BOOST_RANDOM_DETAIL_QRNG_BASE_HPP
+
+#include <istream>
+#include <ostream>
+
+#include <stdexcept>
+
+#include <boost/throw_exception.hpp>
+
+//!\file
+//!Describes the quasi-random number generator base class template.
+
+namespace boost {
+namespace random {
+
+namespace detail {
+
+template<typename DerivedT, typename LatticeT>
+class qrng_base
+{
+protected:
+
+  BOOST_STATIC_CONSTANT(std::size_t, dimension_value = LatticeT::dimension_value);
+
+public:
+
+  typedef typename LatticeT::result_type result_type;
+
+  qrng_base()
+    : lattice(std::size_t())
+  {
+    derived().seed();
+  }
+
+  explicit qrng_base(std::size_t init)
+    : lattice(init)
+  {
+    derived().seed(); // acquire initial values and then move on if necessary
+    derived().seed(init);
+  }
+
+  // default copy c-tor is fine
+
+  // default assignment operator is fine
+
+  //!Requirements: *this is mutable.
+  //!
+  //!Returns: Returns a successive element of an s-dimensional
+  //!(s = X::dimension()) vector at each invocation. When all elements are
+  //!exhausted, X::operator() begins anew with the starting element of a
+  //!subsequent s-dimensional vector.
+  //!
+  //!Throws: overflow_error.
+  result_type operator()()
+  {
+    return curr_elem != dimension_value ? load_saved(): next_state();
+  }
+
+  //!Requirements: *this is mutable.
+  //!
+  //!Effects: Advances *this state as if z consecutive
+  //!X::operator() invocations were executed.
+  //!
+  //!Throws: overflow_error.
+  void discard(std::size_t z)
+  {
+    std::size_t vec_n  = z / dimension_value;
+    std::size_t elem_n = z - vec_n * dimension_value; // z % Dimension
+    std::size_t vec_offset = vec_n + (curr_elem + elem_n) / dimension_value;
+    // Discards vec_offset consecutive s-dimensional vectors
+    discard_vector(vec_offset);
+    // Sets up the proper position of the element-to-read
+    curr_elem += (z - dimension_value * vec_offset);
+  }
+
+protected:
+  DerivedT& derived() throw()
+  {
+    return *static_cast<DerivedT * const>(this);
+  }
+
+  //!Writes a @c DerivedT to a @c std::ostream.
+  template<typename CharT, typename Traits>
+  std::basic_ostream<CharT, Traits>&
+  write_impl(std::basic_ostream<CharT, Traits>& os) const
+  {
+    os << curr_elem << " " << seq_count;
+    return os;
+  }
+
+  //!Reads a @c DerivedT from a @c std::istream.
+  template<typename CharT, typename Traits>
+  std::basic_istream<CharT, Traits>&
+  read_impl(std::basic_istream<CharT, Traits>& is)
+  {
+    std::size_t dim, init;
+    if( is >> dim >> std::ws >> init ) // initialize iff success!
+    {
+      derived().seed(init);
+      curr_elem = dim;
+    }
+    return is;
+  }
+
+  bool is_equal(const DerivedT& rhs) const
+  {
+    // compare normalized sequence position, but without the possible overflow
+    // that would go with seq_count * dimension_value + curr_elem.
+    return (seq_count + (curr_elem / dimension_value) ==
+              rhs.seq_count + (rhs.curr_elem / dimension_value)) &&
+           (curr_elem % dimension_value == rhs.curr_elem % dimension_value);
+  }
+
+private:
+
+  // Load the result from the saved state.
+  result_type load_saved()
+  {
+    return quasi_state[curr_elem++];
+  }
+
+  result_type next_state()
+  {
+    derived().compute_next_vector();
+
+    curr_elem = 0;
+    return load_saved();
+  }
+
+  // Discards z consecutive s-dimensional vectors,
+  // and preserves the position of the element-to-read
+  void discard_vector(std::size_t z)
+  {
+    std::size_t tmp = curr_elem;
+    std::size_t inc_seq_count = seq_count + z;
+    // Here we check that no overflow occurs before we
+    // begin seeding the new value
+    if( inc_seq_count >= seq_count )
+    {
+      derived().seed(inc_seq_count);
+    }
+    else
+    {
+      boost::throw_exception( std::overflow_error("discard_vector") );
+    }
+    curr_elem = tmp;
+  }
+
+protected:
+  LatticeT lattice;
+  std::size_t curr_elem;
+  std::size_t seq_count;
+  result_type quasi_state[dimension_value];
+};
+
+}} // namespace detail::random
+
+} // namespace boost
+
+#endif // BOOST_RANDOM_DETAIL_QRNG_BASE_HPP

Added: sandbox/SOC/2010/quasi_random/boost/random/faure.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/boost/random/faure.hpp	2010-07-21 09:27:55 EDT (Wed, 21 Jul 2010)
@@ -0,0 +1,389 @@
+/* boost random/faure.hpp header file
+ *
+ * Copyright Justinas Vygintas Daugmaudis 2010
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#ifndef BOOST_RANDOM_FAURE_HPP
+#define BOOST_RANDOM_FAURE_HPP
+
+#include <boost/random/detail/qrng_base.hpp>
+#include <boost/random/detail/config.hpp>
+#include <boost/random/detail/operators.hpp>
+
+#include <vector>
+#include <cmath>
+
+#include <boost/integer.hpp>
+
+#include <boost/assert.hpp>
+#include <boost/static_assert.hpp>
+
+#include <boost/mpl/vector/vector50_c.hpp>
+#include <boost/mpl/lower_bound.hpp>
+
+//!\file
+//!Describes the quasi-random number generator class template faure.
+
+namespace boost {
+namespace random {
+
+/** @cond */
+namespace detail {
+namespace fr {
+
+struct prime_list
+{
+  typedef mpl::vector50_c<int,
+      2,    3,    5,    7,   11,   13,   17,   19,   23,   29,
+     31,   37,   41,   43,   47,   53,   59,   61,   67,   71,
+     73,   79,   83,   89,   97,  101,  103,  107,  109,  113,
+    127,  131,  137,  139,  149,  151,  157,  163,  167,  173,
+    179,  181,  191,  193,  197,  199,  211,  223,  227,  229> type;
+};
+
+
+// Returns the smallest prime greater than or equal to N.
+template<int N>
+struct prime_ge
+{
+  typedef prime_list::type prime_list_t;
+  typedef typename mpl::lower_bound<prime_list_t, mpl::int_<N> >::type iter_t;
+
+  BOOST_STATIC_CONSTANT(int, value = mpl::deref<iter_t>::type::value);
+};
+
+// Returns the integer part of the logarithm base Base of arg.
+// In erroneous situations, e.g., integer_log<Base>(0) the function
+// returns 0 and does not report the error. This is the intended
+// behavior.
+template<std::size_t Base>
+inline std::size_t integer_log(std::size_t arg)
+{
+  BOOST_STATIC_ASSERT( 2 <= Base );
+
+  std::size_t ilog = 0;
+  while( Base <= arg )
+  {
+    arg /= Base; ++ilog;
+  }
+  return ilog;
+}
+
+inline std::size_t integer_pow(std::size_t arg, std::size_t p)
+{
+  std::size_t ipow = 1;
+  while( p-- )
+    ipow *= arg;
+  return ipow;
+}
+
+// Computes a table of binomial coefficients modulo qs.
+template<typename RealType, std::size_t Dimension>
+struct binomial_coefficients
+{
+#if defined(BOOST_NO_STATIC_ASSERT)
+  BOOST_STATIC_ASSERT( Dimension <= 229 );
+#else
+  static_assert(Dimension <= 229, "The Faure quasi-random number generator only supports 229 dimensions.");
+#endif
+
+  typedef RealType result_type;
+
+  BOOST_STATIC_CONSTANT(std::size_t, dimension_value = Dimension);
+  BOOST_STATIC_CONSTANT(int, qs_base = prime_ge<Dimension>::value);
+
+  // The maximum prime value that fits into vector50_c is 229, which
+  // means that our primes fit into the uint8_t range.
+  // This merits our attention, because binomial values modulo qs_base
+  // will never be bigger than qs_base. We can choose an appropriate
+  // integer type to hold modulo values and shave off memory footprint.
+
+  typedef typename uint_value_t<qs_base>::fast packed_uint_t;
+
+  inline static RealType inv_qs_base()
+  {
+    return static_cast<RealType>(1) / static_cast<RealType>(qs_base);
+  }
+
+  // default copy c-tor is fine
+
+  binomial_coefficients(std::size_t seq) // c-tor to initialize the binomial table
+  {
+    resize(n_elements(seq));
+  }
+
+  void update(std::size_t seq, RealType (&quasi)[Dimension])
+  {
+    const std::size_t hisum = n_elements(seq);
+    if( coeff.size() != size_hint(hisum) )
+      resize(hisum);
+
+    quasi[0] =
+      compute_recip(seq, hisum, ytemp.rbegin(), ytemp.rend());
+
+    // Find components using the Faure method.
+    for( std::size_t k = 1; k < Dimension; ++k )
+    {
+      quasi[k] = RealType();
+      RealType r = inv_qs_base();
+
+      for( std::size_t i = 0; i != hisum; ++i )
+      {
+        RealType ztemp = RealType();
+        for( std::size_t j = i; j != hisum; ++j )
+          ztemp += ytemp[j] * upper_element(i, j, hisum);
+
+        // Sum ( J <= I <= HISUM ) ( old ytemp(i) * binom(i,j) ) mod QS.
+        ytemp[i] = std::fmod(ztemp, static_cast<RealType>(qs_base));
+        quasi[k] += ytemp[i] * r;
+        r *= inv_qs_base();
+      }
+    }
+  }
+
+private:
+  inline static std::size_t n_elements(std::size_t seq)
+  {
+    return integer_log<qs_base>(seq) + 1;
+  }
+
+  inline static std::size_t size_hint(std::size_t n)
+  {
+    return n * (n + 1) / 2;
+  }
+
+  packed_uint_t& upper_element(std::size_t i, std::size_t j, std::size_t dim)
+  {
+    BOOST_ASSERT( i < dim );
+    BOOST_ASSERT( j < dim );
+    BOOST_ASSERT( i <= j );
+    return coeff[(i * (2 * dim - i + 1)) / 2 + j - i];
+  }
+
+  template<typename Iterator>
+  inline static RealType compute_recip(std::size_t seq, std::size_t n,
+      Iterator first, Iterator last)
+  {
+    //BOOST_ASSERT( std::distance(first, last) == n );
+
+    // Here we do
+    //   Sum ( 0 <= J <= HISUM ) YTEMP(J) * QS**J
+    //   Sum ( 0 <= J <= HISUM ) YTEMP(J) / QS**(J+1)
+    // in one go
+    RealType v, r = RealType();
+    std::size_t m, k = integer_pow(qs_base, n - 1);
+    for( ; first != last; ++first, seq = m, k /= qs_base )
+    {
+      m  = seq % k;
+      v  = (seq - m) / k;
+      r += v;
+      r *= inv_qs_base();
+      *first = v; // saves double dereference
+    }
+    return r;
+  }
+
+  void compute_coefficients(const std::size_t n)
+  {
+    // Resize and initialize to zero
+    coeff.resize(size_hint(n));
+    std::fill(coeff.begin(), coeff.end(), 0);
+
+    // The first row and the diagonal is assigned to 1
+    upper_element(0, 0, n) = 1;
+    for( std::size_t i = 1; i < n; ++i )
+    {
+      upper_element(0, i, n) = 1;
+      upper_element(i, i, n) = 1;
+    }
+
+    // Computes binomial coefficients MOD qs_base
+    for( std::size_t i = 1; i < n; ++i )
+    {
+      for( std::size_t j = i + 1; j < n; ++j )
+      {
+        upper_element(i, j, n) = ( upper_element(i, j-1, n) +
+                                   upper_element(i-1, j-1, n) ) % qs_base;
+      }
+    }
+  }
+
+  void resize(std::size_t n)
+  {
+    ytemp.resize(n);
+    compute_coefficients(n);
+  }
+
+private:
+  // here we cache precomputed data; note that binomial coefficients have
+  // to be recomputed iff the integer part of the logarithm of seq changes,
+  // which happens relatively rarely.
+  std::vector<packed_uint_t> coeff; // packed upper (!) triangular matrix
+  std::vector<RealType> ytemp;
+};
+
+}} // namespace detail::fr
+/** @endcond */
+
+//!class template faure implements a quasi-random number generator as described in
+//! \blockquote
+//!Henri Faure,
+//!Discrepance de suites associees a un systeme de numeration (en dimension s),
+//!Acta Arithmetica,
+//!Volume 41, 1982, pages 337-351.
+//! \endblockquote
+//
+//! \blockquote
+//!Bennett Fox,
+//!Algorithm 647:
+//!Implementation and Relative Efficiency of Quasirandom
+//!Sequence Generators,
+//!ACM Transactions on Mathematical Software,
+//!Volume 12, Number 4, December 1986, pages 362-376.
+//! \endblockquote
+//!
+//!\attention\b Important: This implementation supports up to 229 dimensions.
+//!
+//!In the following documentation @c X denotes the concrete class of the template
+//!faure returning objects of type @c RealType, u and v are the values of @c X.
+//!
+//!Some member functions may throw exceptions of type @c std::bad_alloc.
+template<typename RealType, std::size_t Dimension>
+class faure : public detail::qrng_base<
+                       faure<RealType, Dimension>,
+                       detail::fr::binomial_coefficients<RealType, Dimension> >
+{
+
+  typedef faure<RealType, Dimension> self_t;
+  typedef detail::fr::binomial_coefficients<RealType, Dimension> lattice_t;
+  typedef detail::qrng_base<self_t, lattice_t> base_t;
+
+  friend class detail::qrng_base<self_t, lattice_t>;
+
+public:
+  typedef RealType result_type;
+
+  /** @copydoc boost::random::niederreiter_base2::min() */
+  static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () { return static_cast<RealType>(0); }
+
+  /** @copydoc boost::random::niederreiter_base2::max() */
+  static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () { return static_cast<RealType>(1); }
+
+  /** @copydoc boost::random::niederreiter_base2::dimension() */
+  static std::size_t dimension() { return Dimension; }
+
+  //!Effects: Constructs the default Faure quasi-random number generator.
+  //!
+  //!Throws: bad_alloc.
+  faure()
+    : base_t(0) // initialize the binomial table here
+  {}
+
+  //!Effects: Constructs the Faure quasi-random number generator,
+  //!equivalent to X u; u.seed(init).
+  //!
+  //!Throws: bad_alloc.
+  explicit faure(std::size_t init)
+    : base_t(init) // initialize the binomial table here
+  {}
+
+  /** @copydetails boost::random::niederreiter_base2::seed()
+   * Throws: bad_alloc.
+   */
+  void seed()
+  {
+    seed(0);
+  }
+
+  /** @copydetails boost::random::niederreiter_base2::seed(std::size_t)
+   * Throws: bad_alloc.
+   */
+  void seed(std::size_t init)
+  {
+    this->curr_elem = 0;
+    this->seq_count = init;
+    this->lattice.update(this->seq_count, this->quasi_state);
+  }
+
+  //=========================Doxygen needs this!==============================
+
+  //!Requirements: *this is mutable.
+  //!
+  //!Returns: Returns a successive element of an s-dimensional
+  //!(s = X::dimension()) vector at each invocation. When all elements are
+  //!exhausted, X::operator() begins anew with the starting element of a
+  //!subsequent s-dimensional vector.
+  //!
+  //!Throws: bad_alloc.
+
+  // Fixed in Doxygen 1.7.0 -- id 612458: Fixed problem handling @copydoc for function operators.
+  result_type operator()()
+  {
+    return base_t::operator()();
+  }
+
+  /** @copydoc boost::random::niederreiter_base2::discard(std::size_t)
+   * Throws: bad_alloc.
+   */
+  void discard(std::size_t z)
+  {
+    base_t::discard(z);
+  }
+
+  //!Writes a @c faure to a @c std::ostream.
+  BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, faure, f)
+  { return f.write_impl(os); }
+  //!Reads a @c faure from a @c std::istream.
+  BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, faure, f)
+  { return f.read_impl(is); }
+
+  //!Returns true if the two generators will produce identical sequences.
+  BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(faure, x, y)
+  { return x.is_equal(y); }
+  //!Returns true if the two generators will produce different sequences,
+  BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(faure)
+
+private:
+/** @cond hide_private_members */
+  void compute_next_vector()
+  {
+    this->seq_count++;
+    this->lattice.update(this->seq_count, this->quasi_state);
+  }
+/** @endcond */
+};
+
+} // namespace random
+
+//It would have been nice to do something like this, but it seems that doxygen
+//simply won't show those typedefs, so we spell them out by hand.
+
+/*
+#define BOOST_FAURE_TYPEDEF(z, N, text) \
+typedef random::faure<double, N> faure_##N##d; \
+//
+
+BOOST_PP_REPEAT_FROM_TO(1, 229, BOOST_FAURE_TYPEDEF, unused)
+
+#undef BOOST_FAURE_TYPEDEF
+*/
+
+typedef random::faure<double, 1> faure_1d;
+typedef random::faure<double, 2> faure_2d;
+typedef random::faure<double, 3> faure_3d;
+typedef random::faure<double, 4> faure_4d;
+typedef random::faure<double, 5> faure_5d;
+typedef random::faure<double, 6> faure_6d;
+typedef random::faure<double, 7> faure_7d;
+typedef random::faure<double, 8> faure_8d;
+typedef random::faure<double, 9> faure_9d;
+typedef random::faure<double, 10> faure_10d;
+typedef random::faure<double, 11> faure_11d;
+typedef random::faure<double, 12> faure_12d;
+
+} // namespace boost
+
+#endif // BOOST_RANDOM_FAURE_HPP

Added: sandbox/SOC/2010/quasi_random/boost/random/niederreiter_base2.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/boost/random/niederreiter_base2.hpp	2010-07-21 09:27:55 EDT (Wed, 21 Jul 2010)
@@ -0,0 +1,462 @@
+/* boost random/nierderreiter_base2.hpp header file
+ *
+ * Copyright Justinas Vygintas Daugmaudis 2010
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#ifndef BOOST_RANDOM_NIEDERREITER_BASE2_HPP
+#define BOOST_RANDOM_NIEDERREITER_BASE2_HPP
+
+#include <boost/random/detail/gray_coded_qrng_base.hpp>
+#include <boost/random/detail/config.hpp>
+#include <boost/random/detail/operators.hpp>
+
+#include <limits>
+#include <bitset>
+#include <boost/integer/static_log2.hpp>
+
+#include <boost/static_assert.hpp>
+
+#include <boost/cstdint.hpp>
+
+//!\file
+//!Describes the quasi-random number generator class template niederreiter_base2.
+//!
+//!\b Note: it is especially useful in conjunction with class template uniform_real.
+
+namespace boost {
+namespace random {
+
+/** @cond */
+namespace detail {
+namespace nb2 {
+
+/*
+  Primitive polynomials in binary encoding
+  {
+    { 1, 0, 0, 0, 0, 0 },    1
+    { 0, 1, 0, 0, 0, 0 },    x
+    { 1, 1, 0, 0, 0, 0 },    1 + x
+    { 1, 1, 1, 0, 0, 0 },    1 + x + x^2
+    { 1, 1, 0, 1, 0, 0 },    1 + x + x^3
+    { 1, 0, 1, 1, 0, 0 },    1 + x^2 + x^3
+    { 1, 1, 0, 0, 1, 0 },    1 + x + x^4
+    { 1, 0, 0, 1, 1, 0 },    1 + x^3 + x^4
+    { 1, 1, 1, 1, 1, 0 },    1 + x + x^2 + x^3 + x^4
+    { 1, 0, 1, 0, 0, 1 },    1 + x^2 + x^5
+    { 1, 0, 0, 1, 0, 1 },    1 + x^3 + x^5
+    { 1, 1, 1, 1, 0, 1 },    1 + x + x^2 + x^3 + x^5
+    { 1, 1, 1, 0, 1, 1 }     1 + x + x^2 + x^4 + x^5
+  };
+*/
+
+template<std::size_t D>
+struct primitive_polynomial;
+
+#define BOOST_NIEDERREITER_BASE2_PRIM_POLY(D, V) \
+  template<> \
+   struct primitive_polynomial<D> { \
+     static const int value = V; \
+     static const int degree = ::boost::static_log2<V>::value; \
+   } \
+/**/
+
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(0, 1);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(1, 2);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(2, 3);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(3, 7);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(4, 11);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(5, 13);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(6, 19);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(7, 25);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(8, 31);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(9, 37);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(10, 41);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(11, 47);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(12, 55);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(13, 59);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(14, 61);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(15, 67);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(16, 73);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(17, 87);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(18, 91);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(19, 97);
+BOOST_NIEDERREITER_BASE2_PRIM_POLY(20, 103);
+
+#undef BOOST_NIEDERREITER_BASE2_PRIM_POLY
+
+// Return the base 2 logarithm for a given bitset v
+template<typename Bitset>
+inline int ilog2_bitset_trivial(Bitset v)
+{
+  int ret = 0;
+  for( ; v.any(); v >>= 1 )
+    ++ret;
+  return ret - 1;
+}
+
+template<typename Bitset>
+inline int ilog2_bitset(Bitset v)
+{
+  Bitset a = v >> (v.size() / 2);
+  return a.any() ? ilog2_bitset(a) + (v.size() / 2):
+                   ilog2_bitset_trivial(v);
+}
+
+
+// Multiply polynomials over Z_2.
+template<typename Bitset>
+inline Bitset modulo2_multiply(int P, Bitset v)
+{
+  int pos = 0;
+  Bitset pt = 0;
+  for( ; P; P >>= 1, ++pos )
+    if( P & 1 ) pt ^= v << pos;
+  return pt;
+}
+
+
+// Calculate the values of the constants V(J,R) as
+// described in BFN section 3.3.
+//
+// px = appropriate irreducible polynomial for current dimension
+// pb = polynomial defined in section 2.3 of BFN.
+// pb is modified
+template<typename Bitset, typename T, int MaxE>
+inline void calculate_v(int px, Bitset& pb, int& pb_degree, T (&v)[MaxE])
+{
+  const T arbitrary_element = 1;  // arbitray element of Z_2
+
+  // Now choose a value of Kj as defined in section 3.3.
+  // We must have 0 <= Kj < E*J = M.
+  // The limit condition on Kj does not seem very relevant
+  // in this program.
+  int kj = pb_degree;
+
+  // Now multiply B by PX so B becomes PX**J.
+  // In section 2.3, the values of Bi are defined with a minus sign :
+  // don't forget this if you use them later!
+  pb = modulo2_multiply(px, pb);
+  pb_degree = ilog2_bitset(pb);
+
+  // Now choose values of V in accordance with
+  // the conditions in section 3.3.
+  std::fill(v, v + kj, 0);
+
+  // Quoting from BFN: "Our program currently sets each K_q
+  // equal to eq. This has the effect of setting all unrestricted
+  // values of v to 1."
+  // Actually, it sets them to the arbitrary chosen value.
+  // Whatever.
+  for( int r = kj; r < pb_degree; ++r) {
+    v[r] = arbitrary_element;
+  }
+
+  // Calculate the remaining V's using the recursion of section 2.3,
+  // remembering that the B's have the opposite sign.
+  for(int r = 0; r < MaxE - pb_degree; ++r) {
+    T term = 0;
+    Bitset pb_c = pb;
+    for(int k = 0; k < pb_degree; ++k, pb_c >>= 1)
+    {
+      if( pb_c.test(0) )
+        term ^= v[r + k];
+    }
+    v[r+pb_degree] = term;
+  }
+}
+
+
+template<int Degree, int BitCount>
+struct product_degree
+{
+  BOOST_STATIC_CONSTANT(int, value = Degree * ( (BitCount / Degree) + 1 ));
+};
+
+template<std::size_t D>
+struct compute_lattice
+{
+  template<typename T, int BitCount, int MaxE, std::size_t Dimension>
+  void operator()(T (&ci)[BitCount][BitCount], T (&v)[MaxE],
+                  T (&cj)[BitCount][Dimension]) const
+  {
+    enum {
+      iteration  = Dimension - D,
+      poly_index = iteration + 1,
+      px_value   = primitive_polynomial<poly_index>::value,
+      px_degree  = primitive_polynomial<poly_index>::degree,
+      max_degree = product_degree<px_degree, BitCount>::value
+    };
+
+    // For each dimension, we need to calculate powers of an
+    // appropriate irreducible polynomial, see Niederreiter
+    // page 65, just below equation (19).
+    // Copy the appropriate irreducible polynomial into PX,
+    // and its degree into E.  Set polynomial B = PX ** 0 = 1.
+    // M is the degree of B.  Subsequently B will hold higher
+    // powers of PX.
+    int pb_degree = 0;
+    std::bitset<max_degree> pb_value = 1;
+
+    int j = 0;
+    while( j < BitCount )
+    {
+      // If U = 0, we need to set B to the next power of PX
+      // and recalculate V.
+      calculate_v(px_value, pb_value, pb_degree, v);
+
+      // Niederreiter (page 56, after equation (7), defines two
+      // variables Q and U.  We do not need Q explicitly, but we
+      // do need U.
+
+      // Advance Niederreiter's state variables.
+      for( int u = 0; u < px_degree && j < BitCount; ++u, ++j )
+      {
+        // Now C is obtained from V.  Niederreiter
+        // obtains A from V (page 65, near the bottom), and then gets
+        // C from A (page 56, equation (7)).  However this can be done
+        // in one step.  Here CI(J,R) corresponds to
+        // Niederreiter's C(I,J,R).
+        for(int r = 0; r < BitCount; ++r) {
+          ci[r][j] = v[r + u];
+        }
+      }
+    }
+
+    // The array CI now holds the values of C(I,J,R) for this value
+    // of I.  We pack them into array CJ so that CJ(I,R) holds all
+    // the values of C(I,J,R) for J from 1 to NBITS.
+    for(int r = 0; r < BitCount; ++r)
+    {
+      T term = 0;
+      for(int j = 0; j < BitCount; ++j)
+        term = 2*term + ci[r][j];
+      cj[r][iteration] = term;
+    }
+
+    compute_lattice<D - 1> ncj; ncj(ci, v, cj);
+  }
+};
+
+template<>
+struct compute_lattice<0>
+{
+  template<typename T, int BitCount, int MaxE, std::size_t Dimension>
+  void operator()(T (&)[BitCount][BitCount], T (&)[MaxE],
+                  T (&)[BitCount][Dimension]) const
+  {
+    // recursion stop
+  }
+};
+
+} // namespace nb2
+
+template<typename IntType, std::size_t Dimension>
+struct niederreiter_base2_lattice
+{
+#if defined(BOOST_NO_STATIC_ASSERT)
+  BOOST_STATIC_ASSERT( Dimension <= 20 );
+#else
+  static_assert(Dimension <= 20, "The Niederreiter base 2 quasi-random number generator only supports 20 dimensions.");
+#endif
+
+  typedef IntType result_type;
+
+  BOOST_STATIC_CONSTANT(std::size_t, dimension_value = Dimension);
+
+  BOOST_STATIC_CONSTANT(int, bit_count = std::numeric_limits<IntType>::digits);
+  BOOST_STATIC_CONSTANT(int, prim_degree = nb2::primitive_polynomial<Dimension>::degree);
+
+  // Max degree of a polynomial product that will be computed by modulo2_multiply
+  BOOST_STATIC_CONSTANT(int, max_degree = (nb2::product_degree<prim_degree, bit_count>::value));
+
+  // default copy c-tor is fine
+
+  niederreiter_base2_lattice(std::size_t) // c-tor to initialize the lattice
+  {
+    // initial lattice computation
+    IntType ci[bit_count][bit_count];
+    IntType v[prim_degree + max_degree];
+    nb2::compute_lattice<Dimension> compute; compute(ci, v, bits);
+  }
+
+  result_type operator()(int i, int j) const
+  {
+    return bits[i][j];
+  }
+
+  static std::size_t next(std::size_t seq) { return seq++; }
+
+private:
+  IntType bits[bit_count][Dimension];
+};
+
+} // namespace detail
+/** @endcond */
+
+//!class template niederreiter_base2 implements a quasi-random number generator as described in
+//! \blockquote
+//!Bratley, Fox, Niederreiter, ACM Trans. Model. Comp. Sim. 2, 195 (1992).
+//! \endblockquote
+//!
+//!\attention \b Important: This implementation supports up to 20 dimensions.
+//!
+//!In the following documentation @c X denotes the concrete class of the template
+//!niederreiter_base2 returning objects of type @c IntType, u and v are the values of @c X.
+//!
+//!Some member functions may throw exceptions of type std::overflow_error. This
+//!happens when the quasi-random domain is exhausted and the generator cannot produce
+//!any more values. The length of the low discrepancy sequence is given by
+//! \f$L=Dimension \times 2^{digits}\f$, where digits = std::numeric_limits<IntType>::digits.
+template<typename IntType, std::size_t Dimension, IntType c, IntType m>
+class niederreiter_base2 : public detail::gray_coded_qrng_base<
+                                    niederreiter_base2<IntType, Dimension, c, m>,
+                                    detail::niederreiter_base2_lattice<IntType, Dimension> >
+{
+  typedef niederreiter_base2<IntType, Dimension, c, m> self_t;
+  typedef detail::niederreiter_base2_lattice<IntType, Dimension> lattice_t;
+  typedef detail::gray_coded_qrng_base<self_t, lattice_t> base_t;
+
+public:
+  typedef IntType result_type;
+
+  //!Returns: Tight lower bound on the set of values returned by operator().
+  //!
+  //!Throws: nothing.
+  static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () { return c == 0u ? 1u: 0u; }
+
+  //!Returns: Tight upper bound on the set of values returned by operator().
+  //!
+  //!Throws: nothing.
+  static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () { return m - 1u; }
+
+  //!Returns: The dimension of of the quasi-random domain.
+  //!
+  //!Throws: nothing.
+  static std::size_t dimension() { return Dimension; }
+
+  //!Effects: Constructs the default Niederreiter base 2 quasi-random number generator.
+  //!
+  //!Throws: nothing.
+  niederreiter_base2()
+    : base_t(std::size_t()) // initialize lattice here
+  {}
+
+  //!Effects: Constructs the Niederreiter base 2 quasi-random number generator,
+  //!equivalent to X u; u.seed(init).
+  //!
+  //!Throws: overflow_error.
+  explicit niederreiter_base2(std::size_t init)
+    : base_t(init) // initialize lattice here
+  {}
+
+  //!Requirements: *this is mutable.
+  //!
+  //!Effects: Resets the quasi-random number generator state to
+  //!the one given by the default construction. Equivalent to u.seed(0).
+  //!
+  //!\brief Throws: nothing.
+  void seed()
+  {
+    base_t::reset_state();
+  }
+
+  //!Requirements: *this is mutable.
+  //!
+  //!Effects: Effectively sets the quasi-random number generator state to the init-th
+  //!vector in the s-dimensional quasi-random domain, where s == X::dimension().
+  //!\code
+  //!X u, v;
+  //!for(int i = 0; i < N; ++i)
+  //!    for( std::size_t j = 0; j < u.dimension(); ++j )
+  //!        u();
+  //!v.seed(N);
+  //!assert(u() == v());
+  //!\endcode
+  //!
+  //!\brief Throws: overflow_error.
+  void seed(std::size_t init)
+  {
+    base_t::seed(init, "niederreiter_base2::seed");
+  }
+
+  //=========================Doxygen needs this!==============================
+
+  //!Requirements: *this is mutable.
+  //!
+  //!Returns: Returns a successive element of an s-dimensional
+  //!(s = X::dimension()) vector at each invocation. When all elements are
+  //!exhausted, X::operator() begins anew with the starting element of a
+  //!subsequent s-dimensional vector.
+  //!
+  //!Throws: overflow_error.
+  result_type operator()()
+  {
+    return base_t::operator()();
+  }
+
+  //!Requirements: *this is mutable.
+  //!
+  //!Effects: Advances *this state as if z consecutive
+  //!X::operator() invocations were executed.
+  //!\code
+  //!X u = v;
+  //!for(int i = 0; i < N; ++i)
+  //!    u();
+  //!v.discard(N);
+  //!assert(u() == v());
+  //!\endcode
+  //!
+  //!Throws: overflow_error.
+  void discard(std::size_t z)
+  {
+    base_t::discard(z);
+  }
+
+  //!Writes a @c niederreiter_base2 to a @c std::ostream.
+  BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, niederreiter_base2, f)
+  { return f.write_impl(os); }
+  //!Reads a @c niederreiter_base2 from a @c std::istream.
+  BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, niederreiter_base2, f)
+  { return f.read_impl(is); }
+
+  //!Returns true if the two generators will produce identical sequences.
+  BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(niederreiter_base2, x, y)
+  { return x.is_equal(y); }
+  //!Returns true if the two generators will produce different sequences,
+  BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(niederreiter_base2)
+};
+
+} // namespace random
+
+//It would have been nice to do something like this, but it seems that doxygen
+//simply won't show those typedefs, so we spell them out by hand.
+
+/*
+#define BOOST_NIEDERREITER_BASE2_TYPEDEF(z, N, text) \
+typedef random::niederreiter_base2<uint32_t, N, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_##N##d; \
+//
+
+BOOST_PP_REPEAT_FROM_TO(1, 21, BOOST_NIEDERREITER_BASE2_TYPEDEF, unused)
+
+#undef BOOST_NIEDERREITER_BASE2_TYPEDEF
+*/
+
+typedef random::niederreiter_base2<uint32_t, 1, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_1d;
+typedef random::niederreiter_base2<uint32_t, 2, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_2d;
+typedef random::niederreiter_base2<uint32_t, 3, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_3d;
+typedef random::niederreiter_base2<uint32_t, 4, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_4d;
+typedef random::niederreiter_base2<uint32_t, 5, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_5d;
+typedef random::niederreiter_base2<uint32_t, 6, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_6d;
+typedef random::niederreiter_base2<uint32_t, 7, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_7d;
+typedef random::niederreiter_base2<uint32_t, 8, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_8d;
+typedef random::niederreiter_base2<uint32_t, 9, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_9d;
+typedef random::niederreiter_base2<uint32_t, 10, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_10d;
+typedef random::niederreiter_base2<uint32_t, 11, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_11d;
+typedef random::niederreiter_base2<uint32_t, 12, 1, 2*((uint32_t)1 << 31)> niederreiter_base2_12d;
+
+} // namespace boost
+
+#endif // BOOST_RANDOM_NIEDERREITER_BASE2_HPP

Added: sandbox/SOC/2010/quasi_random/boost/random/sobol.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/boost/random/sobol.hpp	2010-07-21 09:27:55 EDT (Wed, 21 Jul 2010)
@@ -0,0 +1,431 @@
+/* boost random/sobol.hpp header file
+ *
+ * Copyright Justinas Vygintas Daugmaudis 2010
+ * Distributed under the Boost Software License, Version 1.0. (See
+ * accompanying file LICENSE_1_0.txt or copy at
+ * http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#ifndef BOOST_RANDOM_SOBOL_HPP
+#define BOOST_RANDOM_SOBOL_HPP
+
+#include <boost/random/detail/gray_coded_qrng_base.hpp>
+#include <boost/random/detail/config.hpp>
+#include <boost/random/detail/operators.hpp>
+
+#include <limits>
+#include <boost/integer/static_log2.hpp>
+
+#include <boost/cstdint.hpp>
+
+#include <boost/static_assert.hpp>
+
+#include <boost/mpl/vector/vector40_c.hpp>
+#include <boost/mpl/at.hpp>
+
+//!\file
+//!Describes the quasi-random number generator class template sobol.
+//!
+//!\b Note: it is especially useful in conjunction with class template uniform_real.
+
+namespace boost {
+namespace random {
+
+/** @cond */
+namespace detail {
+namespace sbl {
+
+// Primitive polynomials in binary encoding
+
+template<std::size_t D>
+struct primitive_polynomial;
+
+#define BOOST_SOBOL_PRIM_POLY(D, V) \
+  template<> \
+   struct primitive_polynomial<D> { \
+     static const int value = V; \
+     static const int degree = ::boost::static_log2<V>::value; \
+   } \
+/**/
+
+BOOST_SOBOL_PRIM_POLY(0, 1);
+BOOST_SOBOL_PRIM_POLY(1, 3);
+BOOST_SOBOL_PRIM_POLY(2, 7);
+BOOST_SOBOL_PRIM_POLY(3, 11);
+BOOST_SOBOL_PRIM_POLY(4, 13);
+BOOST_SOBOL_PRIM_POLY(5, 19);
+BOOST_SOBOL_PRIM_POLY(6, 25);
+BOOST_SOBOL_PRIM_POLY(7, 37);
+BOOST_SOBOL_PRIM_POLY(8, 59);
+BOOST_SOBOL_PRIM_POLY(9, 47);
+BOOST_SOBOL_PRIM_POLY(10, 61);
+BOOST_SOBOL_PRIM_POLY(11, 55);
+BOOST_SOBOL_PRIM_POLY(12, 41);
+BOOST_SOBOL_PRIM_POLY(13, 67);
+BOOST_SOBOL_PRIM_POLY(14, 97);
+BOOST_SOBOL_PRIM_POLY(15, 91);
+BOOST_SOBOL_PRIM_POLY(16, 109);
+BOOST_SOBOL_PRIM_POLY(17, 103);
+BOOST_SOBOL_PRIM_POLY(18, 115);
+BOOST_SOBOL_PRIM_POLY(19, 131);
+BOOST_SOBOL_PRIM_POLY(20, 193);
+BOOST_SOBOL_PRIM_POLY(21, 137);
+BOOST_SOBOL_PRIM_POLY(22, 145);
+BOOST_SOBOL_PRIM_POLY(23, 143);
+BOOST_SOBOL_PRIM_POLY(24, 241);
+BOOST_SOBOL_PRIM_POLY(25, 157);
+BOOST_SOBOL_PRIM_POLY(26, 185);
+BOOST_SOBOL_PRIM_POLY(27, 167);
+BOOST_SOBOL_PRIM_POLY(28, 229);
+BOOST_SOBOL_PRIM_POLY(29, 171);
+BOOST_SOBOL_PRIM_POLY(30, 213);
+BOOST_SOBOL_PRIM_POLY(31, 191);
+BOOST_SOBOL_PRIM_POLY(32, 253);
+BOOST_SOBOL_PRIM_POLY(33, 203);
+BOOST_SOBOL_PRIM_POLY(34, 211);
+BOOST_SOBOL_PRIM_POLY(35, 239);
+BOOST_SOBOL_PRIM_POLY(36, 247);
+BOOST_SOBOL_PRIM_POLY(37, 285);
+BOOST_SOBOL_PRIM_POLY(38, 369);
+BOOST_SOBOL_PRIM_POLY(39, 299);
+
+#undef BOOST_SOBOL_PRIM_POLY
+
+
+// inirandom sets up the random-number generator to produce a Sobol
+// sequence of at most max dims-dimensional quasi-random vectors.
+// Adapted from ACM TOMS algorithm 659, see
+
+// http://doi.acm.org/10.1145/42288.214372
+
+template<std::size_t D>
+struct vinit40;
+
+template<> struct vinit40<1> {
+  typedef mpl::vector40_c<int,
+    0, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+    1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+    1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+    1, 1, 1, 1, 1, 1, 1, 1, 1, 1> type;
+};
+
+template<> struct vinit40<2>
+{
+  typedef mpl::vector40_c<int,
+    0, 0, 1, 3, 1, 3, 1, 3, 3, 1,
+    3, 1, 3, 1, 3, 1, 1, 3, 1, 3,
+    1, 3, 1, 3, 3, 1, 3, 1, 3, 1,
+    3, 1, 1, 3, 1, 3, 1, 3, 1, 3> type;
+};
+
+template<> struct vinit40<3>
+{
+  typedef mpl::vector40_c<int,
+    0, 0, 0, 7, 5, 1, 3, 3, 7, 5,
+    5, 7, 7, 1, 3, 3, 7, 5, 1, 1,
+    5, 3, 3, 1, 7, 5, 1, 3, 3, 7,
+    5, 1, 1, 5, 7, 7, 5, 1, 3, 3> type;
+};
+
+template<> struct vinit40<4>
+{
+  typedef mpl::vector40_c<int,
+    0,  0,  0,  0,  0,  1,  7,  9, 13, 11,
+    1,  3,  7,  9,  5, 13, 13, 11,  3, 15,
+    5,  3, 15,  7,  9, 13,  9,  1, 11,  7,
+    5, 15,  1, 15, 11,  5,  3,  1,  7,  9> type;
+};
+
+template<> struct vinit40<5>
+{
+  typedef mpl::vector40_c<int,
+    0,  0,  0,  0,  0,  0,  0,  9,  3, 27,
+    15, 29, 21, 23, 19, 11, 25,  7, 13, 17,
+     1, 25, 29,  3, 31, 11,  5, 23, 27, 19,
+    21,  5,  1, 17, 13,  7, 15,  9, 31,  9> type;
+};
+
+template<> struct vinit40<6>
+{
+  typedef mpl::vector40_c<int,
+    0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
+    0,  0,  0, 37, 33,  7,  5, 11, 39, 63,
+   27, 17, 15, 23, 29,  3, 21, 13, 31, 25,
+    9, 49, 33, 19, 29, 11, 19, 27, 15, 25> type;
+};
+
+template<> struct vinit40<7>
+{
+  typedef mpl::vector40_c<int,
+    0,   0,  0,  0,  0,  0,    0,  0,  0,   0,
+    0,   0,  0,  0,  0,  0,    0,  0,  0,  13,
+   33, 115, 41, 79, 17,  29, 119, 75, 73, 105,
+    7,  59, 65, 21,  3, 113,  61, 89, 45, 107> type;
+};
+
+template<> struct vinit40<8>
+{
+  typedef mpl::vector40_c<int,
+    0, 0, 0, 0, 0, 0, 0, 0,  0,  0,
+    0, 0, 0, 0, 0, 0, 0, 0,  0,  0,
+    0, 0, 0, 0, 0, 0, 0, 0,  0,  0,
+    0, 0, 0, 0, 0, 0, 0, 7, 23, 39> type;
+};
+
+
+template<std::size_t D, std::size_t Iteration>
+struct leading_elements
+{
+  template<typename T, int BitCount, std::size_t Dimension>
+  static void assign(T (&cj)[BitCount][Dimension])
+  {
+    typedef typename vinit40<D>::type elems_t;
+    typedef typename mpl::at_c<elems_t, Iteration>::type val_t;
+    cj[D - 1][Iteration] = val_t::value;
+    leading_elements<D - 1, Iteration>::assign(cj);
+  }
+};
+
+template<std::size_t Iteration>
+struct leading_elements<0, Iteration>
+{
+  template<typename T, int BitCount, std::size_t Dimension>
+  static void assign(T (&)[BitCount][Dimension])
+  {}
+};
+
+
+template<std::size_t D>
+struct compute_lattice
+{
+  BOOST_STATIC_ASSERT( D != 0 );
+
+  template<typename T, int BitCount, std::size_t Dimension>
+  void operator()(T (&cj)[BitCount][Dimension]) const
+  {
+    enum {
+      iteration  = Dimension - D + 1,
+      px_value   = primitive_polynomial<iteration>::value,
+      degree_i   = primitive_polynomial<iteration>::degree
+    };
+
+    // Leading elements for dimension i come from vinit<>.
+    //for(int k = 0; k < degree_i; ++k )
+    //  cj[k][iteration] = v_init[k][iteration];
+    leading_elements<degree_i, iteration>::assign(cj);
+
+    // Expand the polynomial bit pattern to separate
+    // components of the logical array includ[].
+    int p_i = px_value;
+    bool includ[degree_i];
+    for( int k = degree_i-1; k >= 0; --k, p_i >>= 1 )
+      includ[k] = (p_i & 1);
+
+    // Calculate remaining elements for this dimension,
+    // as explained in Bratley+Fox, section 2.
+    for(int j = degree_i; j < BitCount; ++j)
+    {
+      T p = 2;
+      T w = cj[j - degree_i][iteration];
+      for(int k = 0; k < degree_i; ++k, p <<= 1)
+        if( includ[k] )
+          w ^= (cj[j-k-1][iteration] * p);
+      cj[j][iteration] = w;
+    }
+
+    compute_lattice<D - 1> ncj; ncj(cj);
+  }
+};
+
+template<>
+struct compute_lattice<1>
+{
+  template<typename T, int BitCount, std::size_t Dimension>
+  void operator()(T (&)[BitCount][Dimension]) const
+  {
+    // recursion stop
+  }
+};
+
+} // namespace sbl
+
+template<typename IntType, std::size_t Dimension>
+struct sobol_lattice
+{
+#if defined(BOOST_NO_STATIC_ASSERT)
+  BOOST_STATIC_ASSERT( Dimension <= 40 );
+#else
+  static_assert(Dimension <= 40, "The Sobol quasi-random number generator only supports 40 dimensions.");
+#endif
+
+  typedef IntType result_type;
+  BOOST_STATIC_CONSTANT(std::size_t, dimension_value = Dimension);
+
+  BOOST_STATIC_CONSTANT(int, bit_count = std::numeric_limits<IntType>::digits - 1);
+
+  // default copy c-tor is fine
+
+  sobol_lattice(std::size_t) // c-tor to initialize the lattice
+  {
+    // Initialize direction table in dimension 0.
+    for( int k = 0; k < bit_count; ++k )
+      bits[k][0] = 1;
+
+    // Initialize in remaining dimensions.
+    sbl::compute_lattice<Dimension> compute; compute(bits);
+
+    // Multiply columns of v by appropriate power of 2.
+    IntType p = 2;
+    for(int j = bit_count-1-1; j >= 0; --j, p <<= 1)
+      for(std::size_t k = 0; k < Dimension; ++k )
+        bits[j][k] *= p;
+  }
+
+  result_type operator()(int i, int j) const
+  {
+    return bits[i][j];
+  }
+
+  static std::size_t next(std::size_t seq) { return ++seq; }
+
+private:
+  IntType bits[bit_count][Dimension];
+};
+
+} // namespace detail
+/** @endcond */
+
+//!class template sobol implements a quasi-random number generator as described in
+//! \blockquote
+//![Bratley+Fox, TOMS 14, 88 (1988)]
+//!and [Antonov+Saleev, USSR Comput. Maths. Math. Phys. 19, 252 (1980)]
+//! \endblockquote
+//!
+//!\attention \b Important: This implementation supports up to 40 dimensions.
+//!
+//!In the following documentation @c X denotes the concrete class of the template
+//!sobol returning objects of type @c IntType, u and v are the values of @c X.
+//!
+//!Some member functions may throw exceptions of type @c std::overflow_error. This
+//!happens when the quasi-random domain is exhausted and the generator cannot produce
+//!any more values. The length of the low discrepancy sequence is given by
+//! \f$L=Dimension \times 2^{digits - 1}\f$, where digits = std::numeric_limits<IntType>::digits.
+template<typename IntType, std::size_t Dimension, IntType c, IntType m>
+class sobol : public detail::gray_coded_qrng_base<
+                       sobol<IntType, Dimension, c, m>,
+                       detail::sobol_lattice<IntType, Dimension> >
+{
+  typedef sobol<IntType, Dimension, c, m> self_t;
+  typedef detail::sobol_lattice<IntType, Dimension> lattice_t;
+  typedef detail::gray_coded_qrng_base<self_t, lattice_t> base_t;
+
+public:
+  typedef IntType result_type;
+
+  /** @copydoc boost::random::niederreiter_base2::min() */
+  static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () { return c == 0u ? 1u: 0u; }
+
+  /** @copydoc boost::random::niederreiter_base2::max() */
+  static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () { return m - 1u; }
+
+  /** @copydoc boost::random::niederreiter_base2::dimension() */
+  static std::size_t dimension() { return Dimension; }
+
+  //!Effects: Constructs the default Sobol quasi-random number generator.
+  //!
+  //!Throws: nothing.
+  sobol()
+    : base_t()
+  {}
+
+  //!Effects: Constructs the Sobol quasi-random number generator,
+  //!equivalent to X u; u.seed(init).
+  //!
+  //!Throws: overflow_error.
+  explicit sobol(std::size_t init)
+    : base_t(init)
+  {}
+
+  // default copy c-tor is fine
+
+  /** @copydoc boost::random::niederreiter_base2::seed() */
+  void seed()
+  {
+    base_t::reset_state();
+    base_t::compute_current_vector("sobol::seed");
+  }
+
+  /** @copydoc boost::random::niederreiter_base2::seed(std::size_t) */
+  void seed(std::size_t init)
+  {
+    base_t::seed(init, "sobol::seed");
+  }
+
+  //=========================Doxygen needs this!==============================
+
+  //!Requirements: *this is mutable.
+  //!
+  //!Returns: Returns a successive element of an s-dimensional
+  //!(s = X::dimension()) vector at each invocation. When all elements are
+  //!exhausted, X::operator() begins anew with the starting element of a
+  //!subsequent s-dimensional vector.
+  //!
+  //!Throws: overflow_error.
+
+  // Fixed in Doxygen 1.7.0 -- id 612458: Fixed problem handling @copydoc for function operators.
+  result_type operator()()
+  {
+    return base_t::operator()();
+  }
+
+  /** @copydoc boost::random::niederreiter_base2::discard(std::size_t) */
+  void discard(std::size_t z)
+  {
+    base_t::discard(z);
+  }
+
+  //!Writes a @c sobol to a @c std::ostream.
+  BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, sobol, f)
+  { return f.write_impl(os); }
+  //!Reads a @c niederreiter_base2 from a @c std::istream.
+  BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, sobol, f)
+  { return f.read_impl(is); }
+
+  //!Returns true if the two generators will produce identical sequences.
+  BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(sobol, x, y)
+  { return x.is_equal(y); }
+  //!Returns true if the two generators will produce different sequences,
+  BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(sobol)
+};
+
+} // namespace random
+
+//It would have been nice to do something like this, but it seems that doxygen
+//simply won't show those typedefs, so we spell them out by hand.
+
+/*
+#define BOOST_SOBOL_TYPEDEF(z, N, text) \
+typedef random::sobol<uint32_t, N, 1, (uint32_t)1 << 31> sobol_##N##d; \
+//
+
+BOOST_PP_REPEAT_FROM_TO(1, 21, BOOST_SOBOL_TYPEDEF, unused)
+
+#undef BOOST_SOBOL_TYPEDEF
+*/
+
+typedef random::sobol<uint32_t, 1, 1, (uint32_t)1 << 31> sobol_1d;
+typedef random::sobol<uint32_t, 2, 1, (uint32_t)1 << 31> sobol_2d;
+typedef random::sobol<uint32_t, 3, 1, (uint32_t)1 << 31> sobol_3d;
+typedef random::sobol<uint32_t, 4, 1, (uint32_t)1 << 31> sobol_4d;
+typedef random::sobol<uint32_t, 5, 1, (uint32_t)1 << 31> sobol_5d;
+typedef random::sobol<uint32_t, 6, 1, (uint32_t)1 << 31> sobol_6d;
+typedef random::sobol<uint32_t, 7, 1, (uint32_t)1 << 31> sobol_7d;
+typedef random::sobol<uint32_t, 8, 1, (uint32_t)1 << 31> sobol_8d;
+typedef random::sobol<uint32_t, 9, 1, (uint32_t)1 << 31> sobol_9d;
+typedef random::sobol<uint32_t, 10, 1, (uint32_t)1 << 31> sobol_10d;
+typedef random::sobol<uint32_t, 11, 1, (uint32_t)1 << 31> sobol_11d;
+typedef random::sobol<uint32_t, 12, 1, (uint32_t)1 << 31> sobol_12d;
+
+} // namespace boost
+
+#endif // BOOST_RANDOM_SOBOL_HPP

Added: sandbox/SOC/2010/quasi_random/libs/random/doc/Jamfile.diff
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/libs/random/doc/Jamfile.diff	2010-07-21 09:27:55 EDT (Wed, 21 Jul 2010)
@@ -0,0 +1,34 @@
+--- Jamfile.bak	2010-06-28 20:31:07.000000000 +0300
++++ Jamfile.v2	2010-06-28 19:24:57.000000000 +0300
+@@ -23,6 +23,7 @@
+     discrete_distribution
+     exponential_distribution
+     extreme_value_distribution
++    faure
+     gamma_distribution
+     geometric_distribution
+     inversive_congruential
+@@ -31,11 +32,13 @@
+     linear_feedback_shift
+     lognormal_distribution
+     mersenne_twister
++    niederreiter_base2
+     normal_distribution
+     poisson_distribution
+     random_number_generator
+     ranlux
+     shuffle_order
++    sobol
+     # shuffle_output
+     subtract_with_carry
+     triangle_distribution
+@@ -88,6 +91,9 @@
+         rand48=\"@xmlonly <classname alt=\\\"boost::random::rand48\\\">rand48</classname> @endxmlonly\" \\
+         mt11213b=\"@xmlonly <classname alt=\\\"boost::random::mt11213b\\\">mt11213b</classname> @endxmlonly\" \\
+         mt19937=\"@xmlonly <classname alt=\\\"boost::random::mt19937\\\">mt19937</classname> @endxmlonly\" \\
++        niederreiter_base2=\"@xmlonly <classname alt=\\\"boost::random::niederreiter_base2\\\">niederreiter_base2</classname> @endxmlonly\" \\
++	sobol=\"@xmlonly <classname alt=\\\"boost::random::sobol\\\">sobol</classname> @endxmlonly\" \\
++	faure=\"@xmlonly <classname alt=\\\"boost::random::faure\\\">faure</classname> @endxmlonly\" \\
+         ecuyer1988=\"@xmlonly <classname alt=\\\"boost::ecuyer1988\\\">ecuyer1988</classname> @endxmlonly\" \\
+         lagged_fibonacci607=\"@xmlonly <classname alt=\\\"boost::lagged_fibonacci607\\\">lagged_fibonacci607</classname> @endxmlonly\" \\
+         lagged_fibonacci44497=\"@xmlonly <classname alt=\\\"boost::lagged_fibonacci44497\\\">lagged_fibonacci44497</classname> @endxmlonly\" \\

Added: sandbox/SOC/2010/quasi_random/libs/random/example/qrng_images.cpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/libs/random/example/qrng_images.cpp	2010-07-21 09:27:55 EDT (Wed, 21 Jul 2010)
@@ -0,0 +1,90 @@
+// Copyright Justinas Vygintas Daugmaudis, 2010.
+// Use, modification and distribution is subject to the
+// Boost Software License, Version 1.0. (See accompanying
+// file LICENSE-1.0 or http://www.boost.org/LICENSE-1.0)
+
+
+#include <boost/random/sobol.hpp>
+#include <boost/random/niederreiter_base2.hpp>
+#include <boost/random/faure.hpp>
+
+#include <boost/random/uniform_real.hpp>
+
+#include <boost/gil/image.hpp>
+#include <boost/gil/typedefs.hpp>
+#include <boost/gil/extension/io/png_io.hpp>
+
+using namespace boost::gil;
+
+// Models a Unary Function
+template <typename Qrng, typename Distribution, typename P>   // Models PixelValueConcept
+struct pixelnoise_fn {
+  typedef point2<ptrdiff_t>    point_t;
+
+  typedef pixelnoise_fn        const_t;
+  typedef P                    value_type;
+  typedef value_type           reference;
+  typedef value_type           const_reference;
+  typedef point_t              argument_type;
+  typedef reference            result_type;
+  BOOST_STATIC_CONSTANT(bool, is_mutable=false);
+
+  value_type                    _pnt_color;
+  point_t                       _img_size;
+  std::vector<point_t>          _random_points;
+
+  static const int MAX_ITER=6000;        // max number of iterations
+
+  pixelnoise_fn() {}
+  pixelnoise_fn(const point_t& sz, const value_type& in_color) : _pnt_color(in_color), _img_size(sz) { gen_random(); }
+
+  result_type operator()(const point_t& p) const
+  {
+    // Admitedly silly and slow, but this is straightforward and enough for demonstration purposes
+    return std::find( _random_points.begin(), _random_points.end(), p ) != _random_points.end() ? _pnt_color: value_type(255,255,255);
+  }
+
+private:
+  void gen_random()
+  {
+    Qrng gen;
+    Distribution d;
+    for( int i = 0; i < MAX_ITER; ++i )
+      _random_points.push_back( point_t( _img_size.x * d(gen), _img_size.y * d(gen) ) );
+  }
+};
+
+// To be used in place of uniform_real, where scaling is not necessary
+template<typename T = double>
+struct identity_distribution
+{
+  typedef T result_type;
+  template<typename QRNG>
+  result_type operator()(QRNG& q) const
+  { return q(); }
+};
+
+
+template<typename Qrng, typename Distribution>
+void output_image(const char* img_name)
+{
+  typedef pixelnoise_fn<Qrng, Distribution, rgb8_pixel_t> deref_t;
+  typedef typename deref_t::point_t  point_t;
+  typedef virtual_2d_locator<deref_t,false> locator_t;
+  typedef image_view<locator_t> my_virt_view_t;
+
+  boost::function_requires<PixelLocatorConcept<locator_t> >();
+  gil_function_requires<StepIteratorConcept<typename locator_t::x_iterator> >();
+
+  point_t dims(500,500);
+  my_virt_view_t virt_img(dims, locator_t(point_t(0,0), point_t(1,1), deref_t(dims, rgb8_pixel_t(0,0,0))));
+  png_write_view(img_name, virt_img);
+}
+
+int main()
+{
+  output_image<boost::sobol_4d, boost::uniform_real<> >("out-sobol.png");
+  output_image<boost::niederreiter_base2_4d, boost::uniform_real<> >("out-niederreiter_base2.png");
+  output_image<boost::faure_4d, identity_distribution<> >("out-faure.png");
+  return 0;
+}

Added: sandbox/SOC/2010/quasi_random/libs/random/test/discrepancy.cpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/libs/random/test/discrepancy.cpp	2010-07-21 09:27:55 EDT (Wed, 21 Jul 2010)
@@ -0,0 +1,818 @@
+// Copyright Justinas Vygintas Daugmaudis, 2010.
+// Use, modification and distribution is subject to the
+// Boost Software License, Version 1.0. (See accompanying
+// file LICENSE-1.0 or http://www.boost.org/LICENSE-1.0)
+
+#include "qrng_typedefs.hpp"
+
+#include <boost/random/mersenne_twister.hpp>
+
+#include <boost/random/uniform_real.hpp>
+
+#include <boost/pool/pool.hpp>
+#include <boost/shared_array.hpp>
+
+#define BOOST_TEST_MAIN
+#include <boost/test/unit_test.hpp>
+#include <boost/test/unit_test_log.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+//
+// DESCRIPTION:
+// ~~~~~~~~~~~~
+//
+// The time required to compute the star discrepancy of a sequence of points
+// in a multidimensional unit cube is prohibitive and the best known upper
+// bounds for the star discrepancy of (t,s)-sequences and (t,m,s)-nets are
+// useful only for sample sizes that grow exponentially with the dimension s.
+//
+// Here we implement an algorithm to compute upper bounds for the star discrepancy of an
+// arbitrary set of n points in the s-dimensional unit cube as proposed in
+// Eric Thiemard's paper: "Computing bounds for the star discrepancy", Computing  65 (2000) 2, 169-186.
+// For an integer k > 1, this algorithm computes in O(ns log k+2^sk^s) time
+// and O(k^s) space a bound that is no better than a function depending on s and k.
+//
+// There are 2 tests:
+// 1) test_halton_discrepancy
+// This tests the discrepancy test itself. We give a sequence of points (Halton sequence)
+// for which the discrepancy bounds are known and then check whether the computed values
+// are what we expect. D_n^*(x) between 0.2666670000 and 0.2724280122
+//
+// 2) test_qrng_discrepancy
+// We compare the discrepancy of the Mersenne twister generator against the discrepancy
+// of low discrepancy sequence generators. It must be noted that the Mersenne twister generator
+// is a pseudo-random number generator and the sequence that it produces does not have the
+// property of low discrepancy, thus any quasi-random number generator must have
+// lower discrepancy bounds.
+
+
+// A basic tree node -- POD
+struct tree_node
+{
+  int id;
+  tree_node* node;
+};
+
+struct tree_handle
+{
+  template<typename Pool>
+  explicit tree_handle(Pool& p, int sz)
+    : lower_bound(1)
+  {
+    // allocate 1 bookkeeping node + sz data nodes
+    children = (tree_node*) p.ordered_malloc(1 + sz);
+    if( children == 0 ) throw std::bad_alloc();
+
+    // zero memory
+    tree_node zero = { 0, 0 };
+    std::fill(children, children + 1 + sz, zero);
+
+    children[0].id = sz;
+  }
+
+  int& operator[](int n) { return data()[n].id; }
+  int operator[](int n) const { return data()[n].id; }
+
+  int low() const { return lower_bound - 1; }
+  int high() const { return children_count() - 1; }
+  int children_count() const { return children[0].id; }
+  void low_bound(int sz) { lower_bound = sz; }
+
+  tree_handle subtree(int n) const
+  {
+    return tree_handle(data()[n].node);
+  }
+
+  template<typename Pool, typename Range, typename Comparator>
+  tree_handle get_subtree(Pool& p, Range range, Comparator comp)
+  {
+    if( !has_subtree(range.second) )
+    {
+      int n_elem = range.second - range.first + 1;
+
+      tree_handle tree(p, n_elem);
+      for(int i = 0; i < n_elem; ++i)
+        tree[i] = (*this)[range.first + i]; // copy relevant data
+
+      if( 1 < n_elem )
+        sort_tree(tree, comp);
+
+      attach_subtree(range.second, tree);
+
+      return tree;
+    }
+    else
+    {
+      return subtree(range.second);
+    }
+  }
+
+private:
+   explicit tree_handle(tree_node* ptr)
+    : lower_bound(1)
+    , children(ptr)
+  {}
+
+  tree_node* data() const {
+    return children + 1;  // !!!
+  }
+
+  bool has_subtree(int n) const
+  {
+    return data()[n].node != 0;
+  }
+
+  void attach_subtree(int n, const tree_handle& sub)
+  {
+    data()[n].node = sub.children;
+  }
+
+  // Never throws
+  template<typename Comparator>
+  inline friend void sort_tree(tree_handle& tree, Comparator comp) {
+    std::sort(tree.data(), tree.data() + tree.children_count(), comp);
+  }
+
+private:
+  int lower_bound;
+  tree_node* children;
+};
+
+
+// An s-dimensional table of values that are stored consecutively in the
+// memory area pointed by the given pointer.
+// The table is strictly read-only.
+template<typename T>
+struct value_table
+{
+  value_table(const T* values, std::size_t width)
+    : tbl(values), dim(width) {}
+  std::size_t width() const { return dim; }
+  T operator()(std::size_t i, std::size_t j) const {
+    return tbl[i * dim + j];
+  }
+private:
+  const T* tbl;
+  std::size_t dim;
+};
+
+
+// Sorting predicate
+template<typename T>
+struct branch_comparator_func
+{
+  explicit branch_comparator_func(value_table<T> pt, std::size_t ofs)
+    : table(pt), n_elem(ofs) {}
+  bool operator()(const tree_node& lhs, const tree_node& rhs) const {
+    return table(lhs.id, n_elem) < table(rhs.id, n_elem);
+  }
+private:
+  value_table<T> table;
+  std::size_t n_elem;
+};
+
+template<typename T>
+inline branch_comparator_func<T> branch_compare(value_table<T> pt, std::size_t ofs = 0) {
+  return branch_comparator_func<T>(pt, ofs);
+}
+
+
+// Functors that more or less do forwarding back to the LexTree.
+// Mostly used to customize the behaviour of binary search.
+// (Simpler than member function binding.)
+struct do_nothing
+{
+  template<typename Range>
+  void operator()(Range) const
+  { /* do nothing at all */ }
+};
+
+struct cutoff_tree
+{
+  tree_handle* nd;
+  cutoff_tree(tree_handle& n) : nd(&n) {}
+  template<typename Range>
+  void operator()(Range range)
+  {
+    nd->low_bound(range.first + 1);
+  }
+};
+
+struct count_ranges
+{
+  template<typename LexTree, typename Range, typename BoundVector>
+  int operator()(LexTree&, tree_handle&, Range range, BoundVector) const
+  { return range.second - range.first + 1; }
+};
+
+struct explode_tree
+{
+  template<typename LexTree, typename Range, typename BoundVector>
+  std::size_t operator()(LexTree& lex, tree_handle& nd, Range range, BoundVector rv) const
+  {
+    tree_handle b = subtree(lex, nd, range, rv);
+    lex.insert(b, rv.range());
+    return lex.explore(b, rv);
+  }
+};
+
+struct explore_implicit_tree
+{
+  template<typename LexTree, typename Range, typename BoundVector>
+  std::size_t operator()(LexTree& lex, tree_handle& nd, Range range, BoundVector rv) const
+  {
+    tree_handle b = subtree(lex, nd, range, rv);
+    return lex.explore(b, rv);
+  }
+};
+
+template<typename T>
+class range_tree : private boost::noncopyable
+{
+private:
+  typedef std::vector<T> vector_t;
+
+  struct vector_range_pair
+  {
+    vector_range_pair deepen() const
+    { return vector_range_pair(p_ptr, m_range + 1); }
+    std::size_t range() const { return m_range; }
+    T value() const { return (*p_ptr)[m_range]; }
+    std::size_t next_range() { return ++m_range; }
+    vector_range_pair(const vector_t& p,const std::size_t ofs)
+       : m_range(ofs), p_ptr(&p)
+     {}
+   private:
+    vector_range_pair(const vector_t* p, const std::size_t ofs)
+      : m_range(ofs), p_ptr(p)
+    {}
+    std::size_t m_range;
+    const vector_t* p_ptr;
+  };
+
+public:
+  range_tree(value_table<T> values, const tree_handle& root, boost::pool<>& p)
+    : table(values)
+    , subtotal(0)
+    , lex(values.width())
+    , pool_(p)
+  {
+    insert(root, 0);
+  }
+
+  std::size_t explore(const vector_t& p, const std::size_t ofs)
+  {
+    BOOST_ASSERT( table.width() == p.size() );
+
+    vector_range_pair rv(p, ofs);
+
+    int total = 0;
+    std::size_t sz = lex[ofs].size();
+    if( ofs + 1 == table.width() )
+    {
+      while( sz-- != 0 )
+      {
+        tree_handle tree = lex[ofs][sz];
+        int low  = tree.low();
+        int high = tree.high();
+        if( high < low )
+        {
+          lex[ofs][sz] = lex[ofs].back();
+          lex[ofs].pop_back();
+          subtotal += low;
+        }
+        else
+        {
+          total += low;
+          total += bin_explore(tree, std::make_pair(low, high), rv, count_ranges(),
+              cutoff_tree(lex[ofs][sz]));
+        }
+      }
+      total += subtotal;
+    }
+    else
+    {
+      subtotal = 0;
+      lex[ofs + 1].clear();
+      for(std::size_t i = 0; i != sz; ++i)
+      {
+        tree_handle tree = lex[ofs][i];
+        int start = tree.low();
+        int low   = 0;
+        int high  = tree.high();
+        while( low != start )
+        {
+          int next = ((1 + low) + high) / 2;
+          tree_handle branch = tree.subtree(next);
+          insert(branch, ofs + 1);
+          total += explore(branch, rv.deepen());
+          low = next + 1;
+        }
+        total += bin_explore(tree, std::make_pair(low, high), rv, explode_tree(),
+            cutoff_tree(lex[ofs][i]));
+      }
+    }
+    return total;
+  }
+
+private:
+  friend struct explode_tree;
+  friend struct explore_implicit_tree;
+
+  template<typename Range, typename BoundVector>
+  inline friend tree_handle subtree(range_tree<T>& lex, tree_handle& nd, Range range, BoundVector rv) {
+    // adds a new (or gets the existing) subtree, sorted at the specified level
+    return nd.get_subtree(lex.pool_, range, branch_compare(lex.table, rv.range()));
+  }
+
+  void insert(const tree_handle& nd, std::size_t ofs)
+  {
+    lex[ofs].push_back(nd);
+  }
+
+  template<typename BoundVector>
+  bool in_range(const tree_handle& nd, std::size_t i, BoundVector rv) const
+  {
+    return table(nd[i], rv.range()) < rv.value();
+  }
+
+  template<typename Range, typename BoundVector, typename Explore, typename Right>
+  std::size_t bin_explore(tree_handle& nd, Range range,
+      BoundVector rv, Explore explore, Right on_right)
+  {
+    bool right = true;
+    std::size_t total = 0;
+    while( range.first <= range.second )
+    {
+      int next = ((1 + range.first) + range.second) / 2;
+      if( in_range(nd, next, rv) )
+      {
+        total += explore(*this, nd, std::make_pair(range.first, next), rv.deepen());
+        range.first = next + 1;
+        if( right ) {
+          on_right(range);
+        }
+      }
+      else
+      {
+        right = false;
+        range.second = next - 1;
+      }
+    }
+    return total;
+  }
+
+  template<typename BoundVector, typename Explore>
+  std::size_t bin_explore(tree_handle& nd, BoundVector rv, Explore explore)
+  {
+    return bin_explore(nd, std::make_pair(0, nd.high()), rv, explore, do_nothing());
+  }
+
+  template<typename BoundVector>
+  std::size_t explore(tree_handle& nd, BoundVector rv)
+  {
+    // how many children do we have?
+    const int sz = nd.children_count();
+
+    // quick case
+    if( !in_range(nd, 0, rv) )
+      return 0;
+
+    // a single element case
+    if( sz == 1 )
+    {
+      while( rv.next_range() != table.width() )
+        if( !in_range(nd, 0, rv) )
+          return 0;
+      return 1;
+    }
+
+    // the last column to be explored
+    if( rv.range() + 1 == table.width() )
+    {
+      // another quick case
+      return in_range(nd, sz - 1, rv) ? sz:
+        bin_explore(nd, rv, count_ranges());
+    }
+
+    return bin_explore(nd, rv, explore_implicit_tree());
+  }
+
+private:
+  typedef std::vector<tree_handle> tree_t;
+  typedef std::vector<tree_t> forest_t;
+
+  value_table<T> table; // table of values
+  int subtotal; // cached
+  forest_t lex;
+
+  boost::pool<>& pool_;
+};
+
+
+// D-star state, uses range trees for counting
+template<typename T>
+class star_discrepancy
+{
+private:
+  struct range_tree_pair
+  {
+    range_tree<T> alpha, beta;
+    template<typename Pool>
+    range_tree_pair(value_table<T> values, const tree_handle& root, Pool& p)
+      : alpha(values, root, p)
+      , beta(values, root, p)
+    {}
+  };
+
+public:
+  typedef std::vector<T>    vector_t; // type shortcuts
+  typedef std::pair<T, T>   bounds_t;
+
+  // Initialize point table
+  template<typename QRNG, typename Uniform>
+  star_discrepancy(QRNG& q, Uniform& d, std::size_t dimension, std::size_t n, T eps)
+    : dim(dimension) // for non-QRNGs q.dimension() is malformed
+    , n_elem(n)
+    , epsilon(eps)
+    , points(new T[dim * n_elem])
+  {
+    // All points must be bounded by the unit hypercube.
+    for(std::size_t i = 0; i != dim * n_elem; ++i)
+    {
+      T value = d(q);
+      BOOST_REQUIRE_GE(value, T(0));
+      BOOST_REQUIRE_LE(value, T(1));
+
+      // store the value; it will be used later
+      points[i] = value;
+    }
+  }
+
+  // Initialize point table
+  template<std::size_t Dimension, std::size_t N>
+  star_discrepancy(T (&pt)[N][Dimension], T eps)
+    : dim(Dimension)
+    , n_elem(N)
+    , epsilon(eps)
+    , points(new T[dim * n_elem])
+  {
+    // All points must be bounded by the unit hypercube.
+    for(std::size_t i = 0; i != n_elem; ++i)
+    {
+      for(std::size_t j = 0; j != dim; ++j)
+      {
+        T value = pt[i][j];
+        BOOST_REQUIRE_GE(value, T(0));
+        BOOST_REQUIRE_LE(value, T(1));
+
+        // store the value; it will be used later
+        points[i * dim + j] = value;
+      }
+    }
+  }
+
+  // Computes lower and upper discrepancy bounds
+  bounds_t compute() const
+  {
+    bounds_t bounds = std::make_pair(T(), T());
+
+    boost::pool<> p(sizeof(tree_node));
+    tree_handle root = create_root(p);
+    range_tree_pair tree(get_table(), root, p);
+
+    // Now we carry out sequence decomposition into subintervals
+    vector_t alpha(dim, T(0)), beta(dim, T(1));
+    decompose(alpha, beta, tree, 0, T(1), bounds);
+
+    return bounds;
+  }
+
+private:
+
+  T volume_diff(T volume, std::size_t n) const
+  {
+    return volume - static_cast<T>(n) / static_cast<T>(n_elem);
+  }
+
+  T abs_volume(T volume, std::size_t n) const
+  {
+    return std::abs( volume_diff(volume, n) );
+  }
+
+  static void upper_bound(T volume, bounds_t& bounds)
+  {
+    if( bounds.second < volume )
+    {
+      bounds.second = volume;
+    }
+  }
+
+  T compute_delta(const vector_t& alpha, const vector_t& beta, std::size_t start) const
+  {
+    T pbetamin = T(1);
+    for(std::size_t i = start; i != dim; ++i)
+      pbetamin *= beta[i];
+    T pbeta = pbetamin;
+
+    T palpha = T(1);
+    for(std::size_t i = 0; i != start; ++i) {
+      pbetamin *= beta[i];
+      palpha   *= alpha[i];
+    }
+
+    pbetamin -= epsilon;
+    return std::pow( pbetamin / ( pbeta * palpha ), T(1) / T(dim - start) );
+  }
+
+  template<typename Predicate>
+  T volume_bound(const vector_t& p, Predicate pred, const T init) const
+  {
+    BOOST_ASSERT( p.size() == dim );
+
+    T volume = T(1);
+    for(std::size_t j = 0; j != dim; ++j)
+    {
+      T u = init;
+      for(std::size_t i = 0; i != n_elem; ++i)
+      {
+        T pt = points[i*dim + j];
+        if( pred(u, pt) && !pred(p[j], pt) )
+          u = pt;
+      }
+      volume *= u;
+    }
+    return volume;
+  }
+
+  void lower_bound(T volume, std::size_t npoints, const vector_t& p, bounds_t& b) const
+  {
+    if( b.first < abs_volume(volume, npoints) )
+    {
+      if( volume_diff(volume, npoints) < T() )
+        volume = volume_bound(p, std::less<T>(), T());
+      else
+        volume = volume_bound(p, std::greater<T>(), T(1));
+      // recompute lower bound approximation
+      b.first = abs_volume(volume, npoints);
+    }
+  }
+
+  void evaluate(const vector_t& alpha, const vector_t& beta, range_tree_pair& tree,
+                std::size_t start,
+                bounds_t& b) const
+  {
+    T va = T(1); // volume alpha
+    T vb  = T(1); // volume beta
+    for(std::size_t i = 0; i != dim; ++i) {
+      va *= alpha[i];
+      vb *= beta[i];
+    }
+
+    // Here we evaluate uppper bound
+    std::size_t nalpha = tree.alpha.explore(alpha, start);
+    std::size_t nbeta = tree.beta.explore(beta, start);
+    upper_bound(-volume_diff(va, nbeta), b);
+    upper_bound(volume_diff(vb, nalpha), b);
+
+    // ...and here we evaluate lower bound
+    lower_bound(va, nalpha, alpha, b);
+    lower_bound(vb, nbeta,  beta, b);
+  }
+
+  // In the original Eric Thiemard's paper
+  // "Computing bounds for the star discrepancy", Computing  65 (2000) 2, 169-186.
+  // the decomposition function requires O(k^s) space. Here we implement an optimized
+  // version of the function that requires O(k * (k-1) * (k-2) * .. * (k-s-1)) space.
+  void decompose(vector_t& alpha, vector_t& beta, range_tree_pair& tree,
+                 std::size_t start, // starting element
+                 T threshold, // decomposition threshold
+                 bounds_t& b) const
+  {
+    T delta = compute_delta(alpha, beta, start);
+
+    // Here we preserve only those alpha values that are going to be changed.
+    // This is more efficient memory-wise than doing a full copy of alpha;
+    // the deeper we go the less elements we will have to preserve.
+    const vector_t oalpha(alpha.begin() + start, alpha.end());
+
+    T obeta = beta[start];
+    beta[start] *= delta;
+
+    threshold *= delta;
+    if( epsilon < threshold )
+    {
+      decompose(alpha, beta, tree, start, threshold, b);
+      for(std::size_t i = start + 1; i < dim; ++i)
+      {
+        alpha[i-1] = obeta * delta;
+        beta[i-1]  = obeta;
+        obeta      = beta[i];
+        beta[i]   *= delta;
+        decompose(alpha, beta, tree, i, threshold, b);
+      }
+    }
+    else
+    {
+      std::size_t k = (start == 0) ? 0: start - 1;
+      evaluate(alpha, beta, tree, k++, b);
+      for(std::size_t i = start + 1; i < dim; ++i)
+      {
+        // evaluate member function is not modifying,
+        // so we can save one multiplication here
+        alpha[i-1] = beta[i-1];
+        beta[i-1]  = obeta;
+        obeta      = beta[i];
+        beta[i]   *= delta;
+        evaluate(alpha, beta, tree, k++, b);
+      }
+    }
+
+    alpha[dim-1] = obeta * delta;
+    beta[dim-1] = obeta; // restore the last changed beta value
+    evaluate(alpha, beta, tree, dim - 1, b);
+
+    // Restore alpha values
+    std::copy(oalpha.begin(), oalpha.end(), alpha.begin() + start);
+  }
+
+  template<typename Pool>
+  tree_handle create_root(Pool& p) const
+  {
+    tree_handle tree(p, n_elem);
+
+    for(std::size_t i = 0; i != n_elem; ++i)
+      tree[i] = i;
+
+    sort_tree(tree, branch_compare(get_table()));
+
+    return tree;
+  }
+
+  // The lifetime of the returned table must not exceed the lifetime of d_star object
+  value_table<T> get_table() const
+  {
+    return value_table<T>(points.get(), dim);
+  }
+
+private:
+  const std::size_t dim;
+  const std::size_t n_elem;
+  const T epsilon;
+  boost::shared_array<T> points;
+};
+
+
+
+// This will test the discrepancy test itself. We give a sequence of points (Halton sequence)
+// for which the discrepancy bounds are known and then check whether the computed values
+// are what we expect.
+BOOST_AUTO_TEST_CASE( test_halton_discrepancy )
+{
+  double halton_seq[10][2] =
+  {
+    {         0,           0    },
+    {       0.5,    0.333333    },
+    {      0.25,    0.666667    },
+    {      0.75,    0.111111    },
+    {     0.125,    0.444444    },
+    {     0.625,    0.777778    },
+    {     0.375,    0.222222    },
+    {     0.875,    0.555556    },
+    {    0.0625,    0.888889    },
+    {    0.5625,    0.037037    }
+  };
+
+  star_discrepancy<double> discr(halton_seq, 0.01);
+  std::pair<double, double> bounds = discr.compute();
+
+  BOOST_REQUIRE_CLOSE_FRACTION(bounds.first, 0.2666670000, 0.0001);
+  BOOST_REQUIRE_CLOSE_FRACTION(bounds.second, 0.2724280122, 0.0001);
+}
+
+
+template<typename QRNG, typename DStar, typename Distribution, typename Bounds>
+inline void check_discrepancy(const char* name,
+    Distribution d,
+    std::size_t n_vectors, typename Distribution::result_type eps,
+    Bounds mersenne_twister_discrepancy)
+{
+  QRNG q; // default-construct
+
+  BOOST_TEST_MESSAGE( "Testing " << name << "[" << q.dimension() << "]" );
+
+  DStar ds(q, d, q.dimension(), n_vectors, eps);
+  Bounds qrng_discrepancy = ds.compute();
+  BOOST_CHECK( qrng_discrepancy < mersenne_twister_discrepancy );
+}
+
+template<int D>
+inline void compare_qrng_discrepancy(std::size_t n_vectors, double eps)
+{
+  // Default parameterization of D-dimensional qrng generator templates
+  typedef typename niederreiter_base2_generator<D>::type niederreiter_base2_t;
+  typedef typename sobol_generator<D>::type sobol_t;
+  typedef typename faure_generator<D>::type faure_t;
+
+  // We test with double precision
+  typedef double value_t;
+  typedef star_discrepancy<value_t> measurer_t;
+  typedef typename measurer_t::bounds_t bounds_t;
+
+  boost::uniform_real<value_t> u;
+
+  BOOST_REQUIRE( n_vectors > 0 );
+  BOOST_REQUIRE( eps > 0 && eps < 1 );
+  BOOST_TEST_MESSAGE( "Starting tests in dimension " << D << " with sample size " << n_vectors << ", epsilon=" << eps );
+
+  // Compute the discrepancy of the Mersenne twister pseudo-random number generator
+  bounds_t mt_discr_bounds;
+  {
+    boost::mt19937 mt;
+
+    BOOST_TEST_MESSAGE( "Computing discrepancy bounds for mt19937[" << D << "]" );
+
+    measurer_t discr(mt, u, D, n_vectors, eps);
+    mt_discr_bounds = discr.compute();
+  }
+
+#define UNIT_TEST_CHECK_QRNG_DISCREPANCY(N, U) \
+  check_discrepancy<N, measurer_t>(#N, U, n_vectors, eps, mt_discr_bounds); \
+/**/
+
+  // Compare the discrepancy of quasi-random number generators against the Mersenne twister discrepancy
+  UNIT_TEST_CHECK_QRNG_DISCREPANCY(niederreiter_base2_t, u);
+  UNIT_TEST_CHECK_QRNG_DISCREPANCY(sobol_t, u);
+  UNIT_TEST_CHECK_QRNG_DISCREPANCY(faure_t, identity_distribution<value_t>());
+
+#undef UNIT_TEST_CHECK_QRNG_DISCREPANCY
+}
+
+BOOST_AUTO_TEST_CASE( test_qrng_discrepancy )
+{
+  // Discrepancy computation in high dimensions is quite expensive, so we have
+  // to decrease the sample size and/or decrease subdivision precision (increase epsilon)
+  // accordingly, so that tests finish reasonably fast.
+  // However, sample size and precision cannot be decreased haphazardly. Take too
+  // small a sample and intrinsic random number generator properties will vanish.
+  // Take too rough a precision and those properties will be determined with
+  // a big error margin.
+  // From this follows that there are three possible reasons for discrepancy test
+  // failures:
+  //   1) not enough sample points given;
+  //   2) not enough precision (epsilon is too big);
+  //   3) genuine qrng implementation errors.
+  // All those three points have to be considered carefully, but
+  // first thing to do, if some test fails, is to increase sample point size!
+
+  compare_qrng_discrepancy<2>(  400,    0.01  );
+  compare_qrng_discrepancy<3>(  100,    0.01  );
+  compare_qrng_discrepancy<4>(  100,    0.04  );
+  compare_qrng_discrepancy<5>(  100,    0.06  );
+  compare_qrng_discrepancy<6>(   60,    0.08  );
+  compare_qrng_discrepancy<7>(   60,    0.15  );
+  compare_qrng_discrepancy<8>(   60,    0.4   );
+  compare_qrng_discrepancy<9> (  60,    0.5   );
+  compare_qrng_discrepancy<10>(  80,    0.3   ); // mersenne twister has good discrepancy from dim 10. :-)
+  compare_qrng_discrepancy<11>(  500,   0.4   );
+  compare_qrng_discrepancy<12>(  50,    0.5   );
+
+  // etc.
+}
+
+
+template<int D>
+inline void check_niederreiter_base_2_discrepancy(std::size_t n_vectors, double eps)
+{
+  // Default parameterization of D-dimensional qrng generator templates
+  typedef typename niederreiter_base2_generator<D>::type niederreiter_base2_t;
+
+  BOOST_REQUIRE( n_vectors > 0 );
+  BOOST_REQUIRE( eps > 0 && eps < 1 );
+
+  // D_n^*(x) = O( (log n)^s∕n )
+  double expected_discrepancy = std::pow(log(n_vectors), D) / n_vectors;
+  BOOST_TEST_MESSAGE( "Expected discrepancy O( (log n)^s∕n ) = " << expected_discrepancy );
+
+  niederreiter_base2_t q;
+  boost::uniform_real<double> u;
+  star_discrepancy<double> ds(q, u, q.dimension(), n_vectors, eps);
+  std::pair<double, double> bounds = ds.compute();
+
+  BOOST_CHECK( bounds.first <= expected_discrepancy );
+  BOOST_CHECK( expected_discrepancy <= bounds.second );
+}
+
+BOOST_AUTO_TEST_CASE( test_niederreiter_base_2_discrepancy )
+{
+  // The Niederreiter base 2 discrepancy satisfies
+  // D_n^*(x) = O( (log n)^s∕n ).
+  // Note, that for large dimension s, the theoretical bound (log n)^s∕n may
+  // only be meaningful for extremely large values of n.
+  // The bound in Koksma-Hlawka inequality gives no relevant information
+  // until a very large number of points is used.
+  BOOST_TEST_MESSAGE( "Starting Niederreiter base 2 discrepancy checks" );
+  check_niederreiter_base_2_discrepancy<2>(10000, 0.01);
+  check_niederreiter_base_2_discrepancy<3>(200000, 0.01);
+}
+

Added: sandbox/SOC/2010/quasi_random/libs/random/test/faure_validate.cpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/libs/random/test/faure_validate.cpp	2010-07-21 09:27:55 EDT (Wed, 21 Jul 2010)
@@ -0,0 +1,452 @@
+// Copyright Justinas Vygintas Daugmaudis, 2010.
+// Use, modification and distribution is subject to the
+// Boost Software License, Version 1.0. (See accompanying
+// file LICENSE-1.0 or http://www.boost.org/LICENSE-1.0)
+
+#include <boost/random/faure.hpp>
+
+#include <boost/utility.hpp>
+
+#define BOOST_TEST_MAIN
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+
+//
+// DESCRIPTION:
+// ~~~~~~~~~~~~
+//
+// This file tests the faure quasi-random number generator.
+// These tests compare our results with values produced by the original
+// version of ACM TOMS Algorithm 647, which  is available in the
+// TOMS subdirectory in http://www.netlib.org
+
+
+//Spatial dimension: 2
+//N: 100
+//Base: 2
+//Vectors skipped: 15
+
+double faure_02_100[100][2] =
+{ { 0.937500, 0.062500},
+{ 0.031250, 0.531250},
+{ 0.531250, 0.031250},
+{ 0.281250, 0.281250},
+{ 0.781250, 0.781250},
+{ 0.156250, 0.156250},
+{ 0.656250, 0.656250},
+{ 0.406250, 0.906250},
+{ 0.906250, 0.406250},
+{ 0.093750, 0.468750},
+{ 0.593750, 0.968750},
+{ 0.343750, 0.718750},
+{ 0.843750, 0.218750},
+{ 0.218750, 0.843750},
+{ 0.718750, 0.343750},
+{ 0.468750, 0.093750},
+{ 0.968750, 0.593750},
+{ 0.015625, 0.796875},
+{ 0.515625, 0.296875},
+{ 0.265625, 0.046875},
+{ 0.765625, 0.546875},
+{ 0.140625, 0.421875},
+{ 0.640625, 0.921875},
+{ 0.390625, 0.671875},
+{ 0.890625, 0.171875},
+{ 0.078125, 0.234375},
+{ 0.578125, 0.734375},
+{ 0.328125, 0.984375},
+{ 0.828125, 0.484375},
+{ 0.203125, 0.609375},
+{ 0.703125, 0.109375},
+{ 0.453125, 0.359375},
+{ 0.953125, 0.859375},
+{ 0.046875, 0.265625},
+{ 0.546875, 0.765625},
+{ 0.296875, 0.515625},
+{ 0.796875, 0.015625},
+{ 0.171875, 0.890625},
+{ 0.671875, 0.390625},
+{ 0.421875, 0.140625},
+{ 0.921875, 0.640625},
+{ 0.109375, 0.703125},
+{ 0.609375, 0.203125},
+{ 0.359375, 0.453125},
+{ 0.859375, 0.953125},
+{ 0.234375, 0.078125},
+{ 0.734375, 0.578125},
+{ 0.484375, 0.828125},
+{ 0.984375, 0.328125},
+{ 0.007812, 0.664062},
+{ 0.507812, 0.164062},
+{ 0.257812, 0.414062},
+{ 0.757812, 0.914062},
+{ 0.132812, 0.039062},
+{ 0.632812, 0.539062},
+{ 0.382812, 0.789062},
+{ 0.882812, 0.289062},
+{ 0.070312, 0.351562},
+{ 0.570312, 0.851562},
+{ 0.320312, 0.601562},
+{ 0.820312, 0.101562},
+{ 0.195312, 0.976562},
+{ 0.695312, 0.476562},
+{ 0.445312, 0.226562},
+{ 0.945312, 0.726562},
+{ 0.039062, 0.132812},
+{ 0.539062, 0.632812},
+{ 0.289062, 0.882812},
+{ 0.789062, 0.382812},
+{ 0.164062, 0.507812},
+{ 0.664062, 0.007812},
+{ 0.414062, 0.257812},
+{ 0.914062, 0.757812},
+{ 0.101562, 0.820312},
+{ 0.601562, 0.320312},
+{ 0.351562, 0.070312},
+{ 0.851562, 0.570312},
+{ 0.226562, 0.445312},
+{ 0.726562, 0.945312},
+{ 0.476562, 0.695312},
+{ 0.976562, 0.195312},
+{ 0.023438, 0.398438},
+{ 0.523438, 0.898438},
+{ 0.273438, 0.648438},
+{ 0.773438, 0.148438},
+{ 0.148438, 0.773438},
+{ 0.648438, 0.273438},
+{ 0.398438, 0.023438},
+{ 0.898438, 0.523438},
+{ 0.085938, 0.585938},
+{ 0.585938, 0.085938},
+{ 0.335938, 0.335938},
+{ 0.835938, 0.835938},
+{ 0.210938, 0.210938},
+{ 0.710938, 0.710938},
+{ 0.460938, 0.960938},
+{ 0.960938, 0.460938},
+{ 0.054688, 0.929688},
+{ 0.554688, 0.429688},
+{ 0.304688, 0.179688} };
+
+//Spatial dimension: 7
+//N: 100
+//Base: 7
+//Vectors skipped: 2400
+
+double faure_07_100[100][7] =
+{{ 0.999584, 0.460225, 0.941274, 0.320283, 0.985006, 0.833403, 0.110371},
+{ 0.000059, 0.243708, 0.376569, 0.649372, 0.668531, 0.358244, 0.222883},
+{ 0.142917, 0.386565, 0.519426, 0.792229, 0.811388, 0.501101, 0.365740},
+{ 0.285774, 0.529422, 0.662284, 0.935087, 0.954245, 0.643958, 0.508598},
+{ 0.428631, 0.672279, 0.805141, 0.077944, 0.097102, 0.786815, 0.651455},
+{ 0.571488, 0.815137, 0.947998, 0.220801, 0.239960, 0.929672, 0.794312},
+{ 0.714345, 0.957994, 0.090855, 0.363658, 0.382817, 0.072529, 0.937169},
+{ 0.857202, 0.100851, 0.233712, 0.506515, 0.525674, 0.215386, 0.080026},
+{ 0.020468, 0.406973, 0.682692, 0.098352, 0.260368, 0.092937, 0.100434},
+{ 0.163325, 0.549830, 0.825549, 0.241209, 0.403225, 0.235795, 0.243291},
+{ 0.306182, 0.692688, 0.968406, 0.384066, 0.546082, 0.378652, 0.386149},
+{ 0.449039, 0.835545, 0.111263, 0.526923, 0.688939, 0.521509, 0.529006},
+{ 0.591896, 0.978402, 0.254120, 0.669780, 0.831796, 0.664366, 0.671863},
+{ 0.734753, 0.121259, 0.396977, 0.812638, 0.974653, 0.807223, 0.814720},
+{ 0.877611, 0.264116, 0.539835, 0.955495, 0.117511, 0.950080, 0.957577},
+{ 0.040876, 0.570239, 0.988814, 0.547331, 0.852204, 0.827631, 0.977985},
+{ 0.183733, 0.713096, 0.131671, 0.690189, 0.995062, 0.970488, 0.120843},
+{ 0.326590, 0.855953, 0.274528, 0.833046, 0.137919, 0.113346, 0.263700},
+{ 0.469447, 0.998810, 0.417386, 0.975903, 0.280776, 0.256203, 0.406557},
+{ 0.612304, 0.141667, 0.560243, 0.118760, 0.423633, 0.399060, 0.549414},
+{ 0.755162, 0.284524, 0.703100, 0.261617, 0.566490, 0.541917, 0.692271},
+{ 0.898019, 0.427381, 0.845957, 0.404474, 0.709347, 0.684774, 0.835128},
+{ 0.061284, 0.590647, 0.152079, 0.996311, 0.301184, 0.562325, 0.855536},
+{ 0.204141, 0.733504, 0.294937, 0.139168, 0.444041, 0.705182, 0.998394},
+{ 0.346998, 0.876361, 0.437794, 0.282025, 0.586898, 0.848040, 0.141251},
+{ 0.489855, 0.019218, 0.580651, 0.424882, 0.729755, 0.990897, 0.284108},
+{ 0.632713, 0.162075, 0.723508, 0.567740, 0.872613, 0.133754, 0.426965},
+{ 0.775570, 0.304932, 0.866365, 0.710597, 0.015470, 0.276611, 0.569822},
+{ 0.918427, 0.447790, 0.009222, 0.853454, 0.158327, 0.419468, 0.712679},
+{ 0.081692, 0.753912, 0.458202, 0.302434, 0.893021, 0.154162, 0.590230},
+{ 0.224549, 0.896769, 0.601059, 0.445291, 0.035878, 0.297019, 0.733087},
+{ 0.367406, 0.039626, 0.743916, 0.588148, 0.178735, 0.439876, 0.875945},
+{ 0.510264, 0.182483, 0.886773, 0.731005, 0.321592, 0.582733, 0.018802},
+{ 0.653121, 0.325341, 0.029631, 0.873862, 0.464449, 0.725591, 0.161659},
+{ 0.795978, 0.468198, 0.172488, 0.016719, 0.607306, 0.868448, 0.304516},
+{ 0.938835, 0.611055, 0.315345, 0.159576, 0.750164, 0.011305, 0.447373},
+{ 0.102100, 0.917177, 0.764324, 0.751413, 0.484857, 0.888856, 0.467781},
+{ 0.244957, 0.060035, 0.907182, 0.894270, 0.627715, 0.031713, 0.610638},
+{ 0.387815, 0.202892, 0.050039, 0.037127, 0.770572, 0.174570, 0.753496},
+{ 0.530672, 0.345749, 0.192896, 0.179985, 0.913429, 0.317427, 0.896353},
+{ 0.673529, 0.488606, 0.335753, 0.322842, 0.056286, 0.460284, 0.039210},
+{ 0.816386, 0.631463, 0.478610, 0.465699, 0.199143, 0.603142, 0.182067},
+{ 0.959243, 0.774320, 0.621467, 0.608556, 0.342000, 0.745999, 0.324924},
+{ 0.122508, 0.080443, 0.070447, 0.200393, 0.076694, 0.623550, 0.345332},
+{ 0.265366, 0.223300, 0.213304, 0.343250, 0.219551, 0.766407, 0.488189},
+{ 0.408223, 0.366157, 0.356161, 0.486107, 0.362409, 0.909264, 0.631047},
+{ 0.551080, 0.509014, 0.499018, 0.628964, 0.505266, 0.052121, 0.773904},
+{ 0.693937, 0.651871, 0.641875, 0.771821, 0.648123, 0.194978, 0.916761},
+{ 0.836794, 0.794728, 0.784733, 0.914678, 0.790980, 0.337835, 0.059618},
+{ 0.979651, 0.937586, 0.927590, 0.057536, 0.933837, 0.480693, 0.202475},
+{ 0.002975, 0.409889, 0.889689, 0.917594, 0.977569, 0.993812, 0.307431},
+{ 0.145832, 0.552746, 0.032546, 0.060451, 0.120426, 0.136669, 0.450289},
+{ 0.288689, 0.695603, 0.175403, 0.203308, 0.263283, 0.279526, 0.593146},
+{ 0.431546, 0.838460, 0.318260, 0.346165, 0.406140, 0.422384, 0.736003},
+{ 0.574404, 0.981317, 0.461117, 0.489022, 0.548997, 0.565241, 0.878860},
+{ 0.717261, 0.124174, 0.603975, 0.631880, 0.691855, 0.708098, 0.021717},
+{ 0.860118, 0.267032, 0.746832, 0.774737, 0.834712, 0.850955, 0.164574},
+{ 0.023383, 0.430297, 0.195811, 0.366573, 0.569406, 0.585649, 0.184982},
+{ 0.166240, 0.573154, 0.338668, 0.509431, 0.712263, 0.728506, 0.327840},
+{ 0.309097, 0.716011, 0.481526, 0.652288, 0.855120, 0.871363, 0.470697},
+{ 0.451955, 0.858868, 0.624383, 0.795145, 0.997977, 0.014220, 0.613554},
+{ 0.594812, 0.001725, 0.767240, 0.938002, 0.140834, 0.157077, 0.756411},
+{ 0.737669, 0.144583, 0.910097, 0.080859, 0.283691, 0.299935, 0.899268},
+{ 0.880526, 0.287440, 0.052954, 0.223716, 0.426548, 0.442792, 0.042125},
+{ 0.043791, 0.593562, 0.501934, 0.815553, 0.018385, 0.320343, 0.062533},
+{ 0.186648, 0.736419, 0.644791, 0.958410, 0.161242, 0.463200, 0.205391},
+{ 0.329506, 0.879276, 0.787648, 0.101267, 0.304099, 0.606057, 0.348248},
+{ 0.472363, 0.022134, 0.930505, 0.244124, 0.446957, 0.748914, 0.491105},
+{ 0.615220, 0.164991, 0.073362, 0.386982, 0.589814, 0.891771, 0.633962},
+{ 0.758077, 0.307848, 0.216219, 0.529839, 0.732671, 0.034628, 0.776819},
+{ 0.900934, 0.450705, 0.359077, 0.672696, 0.875528, 0.177486, 0.919676},
+{ 0.064199, 0.756828, 0.808056, 0.264533, 0.610222, 0.055037, 0.940084},
+{ 0.207057, 0.899685, 0.950913, 0.407390, 0.753079, 0.197894, 0.082942},
+{ 0.349914, 0.042542, 0.093770, 0.550247, 0.895936, 0.340751, 0.225799},
+{ 0.492771, 0.185399, 0.236628, 0.693104, 0.038793, 0.483608, 0.368656},
+{ 0.635628, 0.328256, 0.379485, 0.835961, 0.181651, 0.626465, 0.511513},
+{ 0.778485, 0.471113, 0.522342, 0.978818, 0.324508, 0.769322, 0.654370},
+{ 0.921342, 0.613970, 0.665199, 0.121675, 0.467365, 0.912179, 0.797227},
+{ 0.084608, 0.920093, 0.114179, 0.713512, 0.202059, 0.789730, 0.817636},
+{ 0.227465, 0.062950, 0.257036, 0.856369, 0.344916, 0.932588, 0.960493},
+{ 0.370322, 0.205807, 0.399893, 0.999227, 0.487773, 0.075445, 0.103350},
+{ 0.513179, 0.348664, 0.542750, 0.142084, 0.630630, 0.218302, 0.246207},
+{ 0.656036, 0.491521, 0.685607, 0.284941, 0.773487, 0.361159, 0.389064},
+{ 0.798893, 0.634379, 0.828464, 0.427798, 0.916344, 0.504016, 0.531921},
+{ 0.941750, 0.777236, 0.971321, 0.570655, 0.059202, 0.646873, 0.674778},
+{ 0.105016, 0.083358, 0.420301, 0.019635, 0.793895, 0.524424, 0.695187},
+{ 0.247873, 0.226215, 0.563158, 0.162492, 0.936753, 0.667281, 0.838044},
+{ 0.390730, 0.369072, 0.706015, 0.305349, 0.079610, 0.810139, 0.980901},
+{ 0.533587, 0.511930, 0.848872, 0.448206, 0.222467, 0.952996, 0.123758},
+{ 0.676444, 0.654787, 0.991730, 0.591063, 0.365324, 0.095853, 0.266615},
+{ 0.819301, 0.797644, 0.134587, 0.733920, 0.508181, 0.238710, 0.409472},
+{ 0.962159, 0.940501, 0.277444, 0.876778, 0.651038, 0.381567, 0.552329},
+{ 0.125424, 0.246623, 0.583566, 0.468614, 0.385732, 0.259118, 0.429880},
+{ 0.268281, 0.389481, 0.726424, 0.611471, 0.528589, 0.401975, 0.572738},
+{ 0.411138, 0.532338, 0.869281, 0.754329, 0.671446, 0.544833, 0.715595},
+{ 0.553995, 0.675195, 0.012138, 0.897186, 0.814304, 0.687690, 0.858452},
+{ 0.696853, 0.818052, 0.154995, 0.040043, 0.957161, 0.830547, 0.001309},
+{ 0.839710, 0.960909, 0.297852, 0.182900, 0.100018, 0.973404, 0.144166},
+{ 0.982567, 0.103766, 0.440709, 0.325757, 0.242875, 0.116261, 0.287023},
+{ 0.005890, 0.453621, 0.545665, 0.165407, 0.266199, 0.486523, 0.555245} };
+
+
+//Spatial dimension: 16
+//N: 100
+//Base: 17
+//Vectors skipped: 83520
+
+double faure_16_100[100][16] =
+{ { 0.999988, 0.805606, 0.119874, 0.648675, 0.039068, 0.879288, 0.819854, 0.500904, 0.631781, 0.800721, 0.654781, 0.899846, 0.127609, 0.037033, 0.219813, 0.323009},
+{ 0.000001, 0.073934, 0.994601, 0.786222, 0.064715, 0.791814, 0.286264, 0.277967, 0.971074, 0.974378, 0.253685, 0.248238, 0.802326, 0.116438, 0.803439, 0.949630},
+{ 0.058824, 0.132758, 0.053424, 0.845046, 0.123538, 0.850637, 0.345088, 0.336790, 0.029897, 0.033202, 0.312509, 0.307061, 0.861149, 0.175262, 0.862263, 0.008454},
+{ 0.117648, 0.191581, 0.112248, 0.903869, 0.182362, 0.909461, 0.403911, 0.395614, 0.088721, 0.092025, 0.371332, 0.365885, 0.919973, 0.234086, 0.921086, 0.067277},
+{ 0.176471, 0.250405, 0.171071, 0.962693, 0.241186, 0.968284, 0.462735, 0.454437, 0.147544, 0.150849, 0.430156, 0.424708, 0.978796, 0.292909, 0.979910, 0.126101},
+{ 0.235295, 0.309228, 0.229895, 0.021516, 0.300009, 0.027108, 0.521558, 0.513261, 0.206368, 0.209673, 0.488980, 0.483532, 0.037620, 0.351733, 0.038733, 0.184924},
+{ 0.294118, 0.368052, 0.288719, 0.080340, 0.358833, 0.085931, 0.580382, 0.572084, 0.265191, 0.268496, 0.547803, 0.542355, 0.096444, 0.410556, 0.097557, 0.243748},
+{ 0.352942, 0.426875, 0.347542, 0.139163, 0.417656, 0.144755, 0.639205, 0.630908, 0.324015, 0.327320, 0.606627, 0.601179, 0.155267, 0.469380, 0.156381, 0.302571},
+{ 0.411765, 0.485699, 0.406366, 0.197987, 0.476480, 0.203578, 0.698029, 0.689731, 0.382839, 0.386143, 0.665450, 0.660002, 0.214091, 0.528203, 0.215204, 0.361395},
+{ 0.470589, 0.544522, 0.465189, 0.256810, 0.535303, 0.262402, 0.756852, 0.748555, 0.441662, 0.444967, 0.724274, 0.718826, 0.272914, 0.587027, 0.274028, 0.420218},
+{ 0.529412, 0.603346, 0.524013, 0.315634, 0.594127, 0.321225, 0.815676, 0.807378, 0.500486, 0.503790, 0.783097, 0.777649, 0.331738, 0.645850, 0.332851, 0.479042},
+{ 0.588236, 0.662170, 0.582836, 0.374457, 0.652950, 0.380049, 0.874499, 0.866202, 0.559309, 0.562614, 0.841921, 0.836473, 0.390561, 0.704674, 0.391675, 0.537865},
+{ 0.647060, 0.720993, 0.641660, 0.433281, 0.711774, 0.438872, 0.933323, 0.925026, 0.618133, 0.621437, 0.900744, 0.895296, 0.449385, 0.763497, 0.450498, 0.596689},
+{ 0.705883, 0.779817, 0.700483, 0.492104, 0.770597, 0.497696, 0.992146, 0.983849, 0.676956, 0.680261, 0.959568, 0.954120, 0.508208, 0.822321, 0.509322, 0.655512},
+{ 0.764707, 0.838640, 0.759307, 0.550928, 0.829421, 0.556519, 0.050970, 0.042673, 0.735780, 0.739084, 0.018391, 0.012944, 0.567032, 0.881144, 0.568145, 0.714336},
+{ 0.823530, 0.897464, 0.818130, 0.609752, 0.888244, 0.615343, 0.109793, 0.101496, 0.794603, 0.797908, 0.077215, 0.071767, 0.625855, 0.939968, 0.626969, 0.773160},
+{ 0.882354, 0.956287, 0.876954, 0.668575, 0.947068, 0.674166, 0.168617, 0.160320, 0.853427, 0.856731, 0.136038, 0.130591, 0.684679, 0.998791, 0.685792, 0.831983},
+{ 0.941177, 0.015111, 0.935777, 0.727399, 0.005891, 0.732990, 0.227441, 0.219143, 0.912250, 0.915555, 0.194862, 0.189414, 0.743502, 0.057615, 0.744616, 0.890807},
+{ 0.003461, 0.136218, 0.115708, 0.966153, 0.303469, 0.089391, 0.642665, 0.693192, 0.445122, 0.507250, 0.845381, 0.898757, 0.511668, 0.825781, 0.630429, 0.835443},
+{ 0.062284, 0.195041, 0.174532, 0.024976, 0.362293, 0.148215, 0.701489, 0.752015, 0.503946, 0.566074, 0.904204, 0.957580, 0.570492, 0.884605, 0.689253, 0.894267},
+{ 0.121108, 0.253865, 0.233355, 0.083800, 0.421116, 0.207038, 0.760312, 0.810839, 0.562769, 0.624897, 0.963028, 0.016404, 0.629315, 0.943428, 0.748076, 0.953090},
+{ 0.179932, 0.312689, 0.292179, 0.142624, 0.479940, 0.265862, 0.819136, 0.869662, 0.621593, 0.683721, 0.021851, 0.075227, 0.688139, 0.002252, 0.806900, 0.011914},
+{ 0.238755, 0.371512, 0.351002, 0.201447, 0.538763, 0.324686, 0.877960, 0.928486, 0.680416, 0.742544, 0.080675, 0.134051, 0.746963, 0.061075, 0.865723, 0.070737},
+{ 0.297579, 0.430336, 0.409826, 0.260271, 0.597587, 0.383509, 0.936783, 0.987309, 0.739240, 0.801368, 0.139499, 0.192874, 0.805786, 0.119899, 0.924547, 0.129561},
+{ 0.356402, 0.489159, 0.468649, 0.319094, 0.656410, 0.442333, 0.995607, 0.046133, 0.798063, 0.860192, 0.198322, 0.251698, 0.864610, 0.178722, 0.983370, 0.188384},
+{ 0.415226, 0.547983, 0.527473, 0.377918, 0.715234, 0.501156, 0.054430, 0.104956, 0.856887, 0.919015, 0.257146, 0.310521, 0.923433, 0.237546, 0.042194, 0.247208},
+{ 0.474049, 0.606806, 0.586296, 0.436741, 0.774058, 0.559980, 0.113254, 0.163780, 0.915711, 0.977839, 0.315969, 0.369345, 0.982257, 0.296369, 0.101017, 0.306032},
+{ 0.532873, 0.665630, 0.645120, 0.495565, 0.832881, 0.618803, 0.172077, 0.222603, 0.974534, 0.036662, 0.374793, 0.428168, 0.041080, 0.355193, 0.159841, 0.364855},
+{ 0.591696, 0.724453, 0.703943, 0.554388, 0.891705, 0.677627, 0.230901, 0.281427, 0.033358, 0.095486, 0.433616, 0.486992, 0.099904, 0.414016, 0.218664, 0.423679},
+{ 0.650520, 0.783277, 0.762767, 0.613212, 0.950528, 0.736450, 0.289724, 0.340250, 0.092181, 0.154309, 0.492440, 0.545816, 0.158727, 0.472840, 0.277488, 0.482502},
+{ 0.709343, 0.842100, 0.821590, 0.672035, 0.009352, 0.795274, 0.348548, 0.399074, 0.151005, 0.213133, 0.551263, 0.604639, 0.217551, 0.531663, 0.336311, 0.541326},
+{ 0.768167, 0.900924, 0.880414, 0.730859, 0.068175, 0.854097, 0.407371, 0.457898, 0.209828, 0.271956, 0.610087, 0.663463, 0.276374, 0.590487, 0.395135, 0.600149},
+{ 0.826990, 0.959747, 0.939238, 0.789682, 0.126999, 0.912921, 0.466195, 0.516721, 0.268652, 0.330780, 0.668910, 0.722286, 0.335198, 0.649310, 0.453958, 0.658973},
+{ 0.885814, 0.018571, 0.998061, 0.848506, 0.185822, 0.971744, 0.525018, 0.575545, 0.327475, 0.389603, 0.727734, 0.781110, 0.394021, 0.708134, 0.512782, 0.717796},
+{ 0.944637, 0.077394, 0.056885, 0.907329, 0.244646, 0.030568, 0.583842, 0.634368, 0.386299, 0.448427, 0.786557, 0.839933, 0.452845, 0.766958, 0.571605, 0.776620},
+{ 0.006921, 0.198502, 0.177992, 0.146084, 0.542224, 0.386969, 0.999067, 0.108417, 0.919171, 0.040122, 0.437076, 0.549276, 0.221011, 0.593947, 0.457419, 0.721256},
+{ 0.065745, 0.257325, 0.236815, 0.204907, 0.601047, 0.445793, 0.057890, 0.167240, 0.977994, 0.098946, 0.495900, 0.608099, 0.279835, 0.652771, 0.516242, 0.780080},
+{ 0.124568, 0.316149, 0.295639, 0.263731, 0.659871, 0.504616, 0.116714, 0.226064, 0.036818, 0.157769, 0.554723, 0.666923, 0.338658, 0.711594, 0.575066, 0.838903},
+{ 0.183392, 0.374972, 0.354462, 0.322554, 0.718694, 0.563440, 0.175537, 0.284887, 0.095641, 0.216593, 0.613547, 0.725746, 0.397482, 0.770418, 0.633889, 0.897727},
+{ 0.242215, 0.433796, 0.413286, 0.381378, 0.777518, 0.622263, 0.234361, 0.343711, 0.154465, 0.275416, 0.672371, 0.784570, 0.456305, 0.829241, 0.692713, 0.956551},
+{ 0.301039, 0.492619, 0.472110, 0.440201, 0.836341, 0.681087, 0.293184, 0.402534, 0.213288, 0.334240, 0.731194, 0.843393, 0.515129, 0.888065, 0.751536, 0.015374},
+{ 0.359862, 0.551443, 0.530933, 0.499025, 0.895165, 0.739910, 0.352008, 0.461358, 0.272112, 0.393064, 0.790018, 0.902217, 0.573952, 0.946888, 0.810360, 0.074198},
+{ 0.418686, 0.610266, 0.589757, 0.557848, 0.953988, 0.798734, 0.410832, 0.520181, 0.330935, 0.451887, 0.848841, 0.961040, 0.632776, 0.005712, 0.869183, 0.133021},
+{ 0.477509, 0.669090, 0.648580, 0.616672, 0.012812, 0.857557, 0.469655, 0.579005, 0.389759, 0.510711, 0.907665, 0.019864, 0.691599, 0.064535, 0.928007, 0.191845},
+{ 0.536333, 0.727913, 0.707404, 0.675495, 0.071635, 0.916381, 0.528479, 0.637828, 0.448582, 0.569534, 0.966488, 0.078688, 0.750423, 0.123359, 0.986830, 0.250668},
+{ 0.595156, 0.786737, 0.766227, 0.734319, 0.130459, 0.975205, 0.587302, 0.696652, 0.507406, 0.628358, 0.025312, 0.137511, 0.809246, 0.182182, 0.045654, 0.309492},
+{ 0.653980, 0.845561, 0.825051, 0.793143, 0.189282, 0.034028, 0.646126, 0.755475, 0.566230, 0.687181, 0.084135, 0.196335, 0.868070, 0.241006, 0.104477, 0.368315},
+{ 0.712803, 0.904384, 0.883874, 0.851966, 0.248106, 0.092852, 0.704949, 0.814299, 0.625053, 0.746005, 0.142959, 0.255158, 0.926893, 0.299829, 0.163301, 0.427139},
+{ 0.771627, 0.963208, 0.942698, 0.910790, 0.306930, 0.151675, 0.763773, 0.873122, 0.683877, 0.804828, 0.201782, 0.313982, 0.985717, 0.358653, 0.222124, 0.485962},
+{ 0.830451, 0.022031, 0.001521, 0.969613, 0.365753, 0.210499, 0.822596, 0.931946, 0.742700, 0.863652, 0.260606, 0.372805, 0.044540, 0.417477, 0.280948, 0.544786},
+{ 0.889274, 0.080855, 0.060345, 0.028437, 0.424577, 0.269322, 0.881420, 0.990769, 0.801524, 0.922475, 0.319429, 0.431629, 0.103364, 0.476300, 0.339772, 0.603609},
+{ 0.948098, 0.139678, 0.119168, 0.087260, 0.483400, 0.328146, 0.940243, 0.049593, 0.860347, 0.981299, 0.378253, 0.490452, 0.162187, 0.535124, 0.398595, 0.662433},
+{ 0.010381, 0.260785, 0.299099, 0.326015, 0.780978, 0.684547, 0.296645, 0.523641, 0.393219, 0.572994, 0.028772, 0.199795, 0.930354, 0.362113, 0.284408, 0.607070},
+{ 0.069205, 0.319609, 0.357923, 0.384838, 0.839801, 0.743371, 0.355468, 0.582465, 0.452043, 0.631818, 0.087595, 0.258618, 0.989177, 0.420937, 0.343232, 0.665893},
+{ 0.128028, 0.378432, 0.416746, 0.443662, 0.898625, 0.802194, 0.414292, 0.641289, 0.510866, 0.690641, 0.146419, 0.317442, 0.048001, 0.479760, 0.402055, 0.724717},
+{ 0.186852, 0.437256, 0.475570, 0.502485, 0.957449, 0.861018, 0.473115, 0.700112, 0.569690, 0.749465, 0.205242, 0.376265, 0.106824, 0.538584, 0.460879, 0.783540},
+{ 0.245675, 0.496080, 0.534393, 0.561309, 0.016272, 0.919841, 0.531939, 0.758936, 0.628513, 0.808288, 0.264066, 0.435089, 0.165648, 0.597407, 0.519702, 0.842364},
+{ 0.304499, 0.554903, 0.593217, 0.620132, 0.075096, 0.978665, 0.590762, 0.817759, 0.687337, 0.867112, 0.322890, 0.493912, 0.224471, 0.656231, 0.578526, 0.901187},
+{ 0.363323, 0.613727, 0.652040, 0.678956, 0.133919, 0.037488, 0.649586, 0.876583, 0.746160, 0.925935, 0.381713, 0.552736, 0.283295, 0.715054, 0.637349, 0.960011},
+{ 0.422146, 0.672550, 0.710864, 0.737779, 0.192743, 0.096312, 0.708409, 0.935406, 0.804984, 0.984759, 0.440537, 0.611559, 0.342118, 0.773878, 0.696173, 0.018834},
+{ 0.480970, 0.731374, 0.769687, 0.796603, 0.251566, 0.155135, 0.767233, 0.994230, 0.863807, 0.043583, 0.499360, 0.670383, 0.400942, 0.832701, 0.754996, 0.077658},
+{ 0.539793, 0.790197, 0.828511, 0.855426, 0.310390, 0.213959, 0.826056, 0.053053, 0.922631, 0.102406, 0.558184, 0.729207, 0.459765, 0.891525, 0.813820, 0.136481},
+{ 0.598617, 0.849021, 0.887334, 0.914250, 0.369213, 0.272782, 0.884880, 0.111877, 0.981454, 0.161230, 0.617007, 0.788030, 0.518589, 0.950349, 0.872644, 0.195305},
+{ 0.657440, 0.907844, 0.946158, 0.973073, 0.428037, 0.331606, 0.943703, 0.170700, 0.040278, 0.220053, 0.675831, 0.846854, 0.577412, 0.009172, 0.931467, 0.254128},
+{ 0.716264, 0.966668, 0.004981, 0.031897, 0.486860, 0.390429, 0.002527, 0.229524, 0.099102, 0.278877, 0.734654, 0.905677, 0.636236, 0.067996, 0.990291, 0.312952},
+{ 0.775087, 0.025491, 0.063805, 0.090720, 0.545684, 0.449253, 0.061351, 0.288347, 0.157925, 0.337700, 0.793478, 0.964501, 0.695059, 0.126819, 0.049114, 0.371775},
+{ 0.833911, 0.084315, 0.122629, 0.149544, 0.604507, 0.508077, 0.120174, 0.347171, 0.216749, 0.396524, 0.852301, 0.023324, 0.753883, 0.185643, 0.107938, 0.430599},
+{ 0.892734, 0.143138, 0.181452, 0.208367, 0.663331, 0.566900, 0.178998, 0.405994, 0.275572, 0.455347, 0.911125, 0.082148, 0.812706, 0.244466, 0.166761, 0.489423},
+{ 0.951558, 0.201962, 0.240276, 0.267191, 0.722154, 0.625724, 0.237821, 0.464818, 0.334396, 0.514171, 0.969948, 0.140971, 0.871530, 0.303290, 0.225585, 0.548246},
+{ 0.013842, 0.323069, 0.420206, 0.505945, 0.019732, 0.982125, 0.653046, 0.938866, 0.867268, 0.105866, 0.620467, 0.850314, 0.639696, 0.130279, 0.111398, 0.492883},
+{ 0.072665, 0.381893, 0.479030, 0.564769, 0.078556, 0.040948, 0.711870, 0.997690, 0.926091, 0.164690, 0.679291, 0.909137, 0.698520, 0.189103, 0.170221, 0.551706},
+{ 0.131489, 0.440716, 0.537853, 0.623592, 0.137379, 0.099772, 0.770693, 0.056513, 0.984915, 0.223513, 0.738114, 0.967961, 0.757343, 0.247926, 0.229045, 0.610530},
+{ 0.190312, 0.499540, 0.596677, 0.682416, 0.196203, 0.158596, 0.829517, 0.115337, 0.043738, 0.282337, 0.796938, 0.026784, 0.816167, 0.306750, 0.287868, 0.669353},
+{ 0.249136, 0.558363, 0.655501, 0.741239, 0.255026, 0.217419, 0.888340, 0.174160, 0.102562, 0.341160, 0.855762, 0.085608, 0.874990, 0.365573, 0.346692, 0.728177},
+{ 0.307959, 0.617187, 0.714324, 0.800063, 0.313850, 0.276243, 0.947164, 0.232984, 0.161385, 0.399984, 0.914585, 0.144431, 0.933814, 0.424397, 0.405515, 0.787000},
+{ 0.366783, 0.676010, 0.773148, 0.858886, 0.372673, 0.335066, 0.005987, 0.291808, 0.220209, 0.458807, 0.973409, 0.203255, 0.992637, 0.483220, 0.464339, 0.845824},
+{ 0.425606, 0.734834, 0.831971, 0.917710, 0.431497, 0.393890, 0.064811, 0.350631, 0.279032, 0.517631, 0.032232, 0.262079, 0.051461, 0.542044, 0.523163, 0.904647},
+{ 0.484430, 0.793657, 0.890795, 0.976534, 0.490321, 0.452713, 0.123634, 0.409455, 0.337856, 0.576455, 0.091056, 0.320902, 0.110284, 0.600868, 0.581986, 0.963471},
+{ 0.543253, 0.852481, 0.949618, 0.035357, 0.549144, 0.511537, 0.182458, 0.468278, 0.396679, 0.635278, 0.149879, 0.379726, 0.169108, 0.659691, 0.640810, 0.022294},
+{ 0.602077, 0.911304, 0.008442, 0.094181, 0.607968, 0.570360, 0.241281, 0.527102, 0.455503, 0.694102, 0.208703, 0.438549, 0.227931, 0.718515, 0.699633, 0.081118},
+{ 0.660900, 0.970128, 0.067265, 0.153004, 0.666791, 0.629184, 0.300105, 0.585925, 0.514326, 0.752925, 0.267526, 0.497373, 0.286755, 0.777338, 0.758457, 0.139942},
+{ 0.719724, 0.028952, 0.126089, 0.211828, 0.725615, 0.688007, 0.358928, 0.644749, 0.573150, 0.811749, 0.326350, 0.556196, 0.345578, 0.836162, 0.817280, 0.198765},
+{ 0.778547, 0.087775, 0.184912, 0.270651, 0.784438, 0.746831, 0.417752, 0.703572, 0.631974, 0.870572, 0.385173, 0.615020, 0.404402, 0.894985, 0.876104, 0.257589},
+{ 0.837371, 0.146599, 0.243736, 0.329475, 0.843262, 0.805654, 0.476575, 0.762396, 0.690797, 0.929396, 0.443997, 0.673843, 0.463226, 0.953809, 0.934927, 0.316412},
+{ 0.896194, 0.205422, 0.302559, 0.388298, 0.902085, 0.864478, 0.535399, 0.821219, 0.749621, 0.988219, 0.502820, 0.732667, 0.522049, 0.012632, 0.993751, 0.375236},
+{ 0.955018, 0.264246, 0.361383, 0.447122, 0.960909, 0.923301, 0.594223, 0.880043, 0.808444, 0.047043, 0.561644, 0.791490, 0.580873, 0.071456, 0.052574, 0.434059},
+{ 0.017302, 0.385353, 0.541314, 0.685876, 0.258487, 0.279703, 0.009447, 0.295268, 0.341316, 0.638738, 0.212163, 0.500833, 0.349039, 0.898445, 0.938387, 0.378696},
+{ 0.076125, 0.444176, 0.600137, 0.744700, 0.317310, 0.338526, 0.068271, 0.354091, 0.400140, 0.697562, 0.270986, 0.559656, 0.407862, 0.957269, 0.997211, 0.437519},
+{ 0.134949, 0.503000, 0.658961, 0.803523, 0.376134, 0.397350, 0.127094, 0.412915, 0.458963, 0.756385, 0.329810, 0.618480, 0.466686, 0.016092, 0.056035, 0.496343},
+{ 0.193772, 0.561823, 0.717784, 0.862347, 0.434957, 0.456173, 0.185918, 0.471738, 0.517787, 0.815209, 0.388634, 0.677303, 0.525509, 0.074916, 0.114858, 0.555166},
+{ 0.252596, 0.620647, 0.776608, 0.921170, 0.493781, 0.514997, 0.244742, 0.530562, 0.576610, 0.874032, 0.447457, 0.736127, 0.584333, 0.133740, 0.173682, 0.613990},
+{ 0.311419, 0.679471, 0.835431, 0.979994, 0.552604, 0.573820, 0.303565, 0.589385, 0.635434, 0.932856, 0.506281, 0.794950, 0.643156, 0.192563, 0.232505, 0.672814},
+{ 0.370243, 0.738294, 0.894255, 0.038817, 0.611428, 0.632644, 0.362389, 0.648209, 0.694257, 0.991679, 0.565104, 0.853774, 0.701980, 0.251387, 0.291329, 0.731637},
+{ 0.429066, 0.797118, 0.953078, 0.097641, 0.670251, 0.691468, 0.421212, 0.707032, 0.753081, 0.050503, 0.623928, 0.912598, 0.760803, 0.310210, 0.350152, 0.790461},
+{ 0.487890, 0.855941, 0.011902, 0.156464, 0.729075, 0.750291, 0.480036, 0.765856, 0.811904, 0.109327, 0.682751, 0.971421, 0.819627, 0.369034, 0.408976, 0.849284},
+{ 0.546714, 0.914765, 0.070725, 0.215288, 0.787898, 0.809115, 0.538859, 0.824680, 0.870728, 0.168150, 0.741575, 0.030245, 0.878450, 0.427857, 0.467799, 0.908108},
+{ 0.605537, 0.973588, 0.129549, 0.274111, 0.846722, 0.867938, 0.597683, 0.883503, 0.929551, 0.226974, 0.800398, 0.089068, 0.937274, 0.486681, 0.526623, 0.966931},
+{ 0.664361, 0.032412, 0.188372, 0.332935, 0.905545, 0.926762, 0.656506, 0.942327, 0.988375, 0.285797, 0.859222, 0.147892, 0.996097, 0.545504, 0.585446, 0.025755},
+{ 0.723184, 0.091235, 0.247196, 0.391758, 0.964369, 0.985585, 0.715330, 0.001150, 0.047198, 0.344621, 0.918045, 0.206715, 0.054921, 0.604328, 0.644270, 0.084578},
+{ 0.782008, 0.150059, 0.306020, 0.450582, 0.023192, 0.044409, 0.774153, 0.059974, 0.106022, 0.403444, 0.976869, 0.265539, 0.113745, 0.663151, 0.703093, 0.143402} };
+
+// Invokes operator() precisely n times. This is to check that
+// QRNG::discard(n) actually has the same effect.
+template<typename QRNG>
+void trivial_discard(QRNG &q, std::size_t n)
+{
+  for( ; n != 0; --n ) q();
+}
+
+template<typename T, std::size_t Dimension, std::size_t N>
+inline void test_values(T (&dset)[N][Dimension], std::size_t skip)
+{
+  boost::random::faure<T, Dimension> qf;
+
+  qf.seed(skip);
+  for( std::size_t i = 0; i != N; ++i )
+  {
+    for( std::size_t j = 0; j != Dimension; ++j )
+    {
+      T val = qf();
+      // We want to check that quasi-random number generator values differ no
+      // more than 1% of their value.
+      // Such a coarse check is because dset is initialized from fixed point output.
+      BOOST_CHECK_CLOSE_FRACTION(dset[i][j], val, 1);
+    }
+  }
+}
+
+template<typename T, std::size_t Dimension, std::size_t N>
+inline void test_seed(T (&dset)[N][Dimension], std::size_t skip)
+{
+  boost::random::faure<T, Dimension> qf;
+
+  for( std::size_t i = 0; i != N; ++i )
+  {
+    qf.seed(skip + i);
+    for( std::size_t j = 0; j != Dimension; ++j )
+    {
+      T val = qf();
+      // We want to check that quasi-random number generator values differ no
+      // more than 1% of their value.
+      BOOST_CHECK_CLOSE_FRACTION(dset[i][j], val, 1);
+    }
+  }
+}
+
+template<typename T, std::size_t Dimension, std::size_t N>
+inline void test_discard(T (&dset)[N][Dimension], std::size_t skip)
+{
+  boost::random::faure<T, Dimension> qf, qtrivial, qstart;
+
+  const std::size_t element_count = N * Dimension;
+  const T* dset_array = reinterpret_cast<T *>(boost::addressof(dset));
+
+  qstart.seed(skip);
+  for( std::size_t step = 0; step != element_count; ++step )
+  {
+    // Init to the same state
+    qf = qstart;
+    qtrivial = qstart;
+
+    // Discards have to have the same effect
+    qf.discard(step);
+    trivial_discard(qtrivial, step);
+
+    // Therefore, states are equal
+    BOOST_CHECK( qf == qtrivial );
+
+    // Now, let's check whether they really produce the same sequence
+    for( std::size_t k = step; k != element_count; ++k )
+    {
+      T q_val = qf();
+      T t_val = qtrivial();
+      BOOST_CHECK_CLOSE_FRACTION(q_val, t_val, 0.00001);
+      // ~ BOOST_CHECK(q_val == t_val), but those are floating point values,
+      // so strict equality check may fail unnecessarily
+
+      // States remain equal!
+      BOOST_CHECK( qf == qtrivial );
+
+      // We want to check that quasi-random number generator values differ no
+      // more than 1% of their value.
+      BOOST_CHECK_CLOSE_FRACTION(dset_array[k], q_val, 1);
+    }
+  }
+}
+
+BOOST_AUTO_TEST_CASE( validate_faure )
+{
+  test_values(faure_02_100, 15);
+  test_values(faure_07_100, 2400);
+  test_values(faure_16_100, 83520);
+}
+
+BOOST_AUTO_TEST_CASE( validate_faure_seed )
+{
+  test_seed(faure_02_100, 15);
+  test_seed(faure_07_100, 2400);
+  test_seed(faure_16_100, 83520);
+}
+
+BOOST_AUTO_TEST_CASE( validate_faure_discard )
+{
+  test_discard(faure_02_100, 15);
+  test_discard(faure_07_100, 2400);
+  test_discard(faure_16_100, 83520);
+}

Added: sandbox/SOC/2010/quasi_random/libs/random/test/gsl_validate.cpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/libs/random/test/gsl_validate.cpp	2010-07-21 09:27:55 EDT (Wed, 21 Jul 2010)
@@ -0,0 +1,122 @@
+// Copyright Justinas Vygintas Daugmaudis, 2010.
+// Use, modification and distribution is subject to the
+// Boost Software License, Version 1.0. (See accompanying
+// file LICENSE-1.0 or http://www.boost.org/LICENSE-1.0)
+
+#include "qrng_typedefs.hpp"
+
+#include <gsl/gsl_qrng.h>
+
+
+#include <boost/random/uniform_real.hpp>
+
+#include <boost/preprocessor/repeat_from_to.hpp>
+
+#include <boost/shared_ptr.hpp>
+
+#define BOOST_TEST_MAIN
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+//
+// DESCRIPTION:
+// ~~~~~~~~~~~~
+//
+// This file tests the quasi-random number generators.
+// There are two sets of tests:
+// 1) validate_sobol;
+// 2) validate_niederreiter_base2.
+// These tests compare our results with values produced by qrngs in
+// GNU Scientific Library.
+
+// GSL adaptors
+
+struct gsl_niederreiter_base2
+{
+  static const gsl_qrng_type * type() {
+    return ::gsl_qrng_niederreiter_2;
+  }
+};
+
+struct gsl_sobol
+{
+  static const gsl_qrng_type * type() {
+    return ::gsl_qrng_sobol;
+  }
+};
+
+struct gsl_deleter
+{
+  void operator()(gsl_qrng *q) const
+  {
+    if( q != 0 )
+      gsl_qrng_free(q);
+  }
+};
+
+// Test routine
+template<typename GSL, typename QRNG>
+bool test_qrng(const std::size_t D)
+{
+  QRNG test;
+
+  const int validating_sequence_lenght = 10000;
+
+  boost::shared_ptr<gsl_qrng> q(gsl_qrng_alloc(GSL::type(), D), gsl_deleter());
+  // RAII, to ensure that no memory leaks occur
+  BOOST_REQUIRE( q.get() != 0 );
+
+  double val[D];
+  boost::uniform_real<double> rnd;
+
+  for( int i = 0; i < validating_sequence_lenght; ++i )
+  {
+    gsl_qrng_get(q.get(), val);
+    for( std::size_t j = 0; j != D; ++j )
+    {
+      double n = rnd(test);
+
+      // We want to check that quasi-random number generator values differ no
+      // more than 0.0001% of their value.
+      BOOST_REQUIRE_CLOSE_FRACTION(val[j], n, 0.0001);
+    }
+  }
+
+  return true;
+}
+
+
+template<int D>
+void test_sobol()
+{
+  typedef typename sobol_generator<D>::type sobol_t;
+  BOOST_CHECK( (test_qrng<gsl_sobol, sobol_t>(D)) );
+}
+
+template<int D>
+void test_niederreiter_base2()
+{
+  typedef typename niederreiter_base2_generator<D>::type nb2_t;
+  BOOST_CHECK( (test_qrng<gsl_niederreiter_base2, nb2_t>(D)) );
+}
+
+
+BOOST_AUTO_TEST_CASE( validate_sobol )
+{
+#define UNIT_TEST_QRNG_VALIDATE_SOBOL(z, N, text) test_sobol<N>();
+
+  BOOST_PP_REPEAT_FROM_TO(1, 41, UNIT_TEST_QRNG_VALIDATE_SOBOL, unused)
+
+#undef UNIT_TEST_QRNG_VALIDATE_SOBOL
+}
+
+BOOST_AUTO_TEST_CASE( validate_niederreiter_base2 )
+{
+#define UNIT_TEST_QRNG_VALIDATE_NIEDERREITER_BASE2(z, N, text) test_niederreiter_base2<N>();
+
+  // GSL supports only up to 12 dimensions, so we check only that much
+  BOOST_PP_REPEAT_FROM_TO(1, 13, UNIT_TEST_QRNG_VALIDATE_NIEDERREITER_BASE2, unused)
+
+#undef UNIT_TEST_QRNG_VALIDATE_NIEDERREITER_BASE2
+}
+

Added: sandbox/SOC/2010/quasi_random/libs/random/test/qrng_typedefs.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/libs/random/test/qrng_typedefs.hpp	2010-07-21 09:27:55 EDT (Wed, 21 Jul 2010)
@@ -0,0 +1,45 @@
+// Copyright Justinas Vygintas Daugmaudis, 2010.
+// Use, modification and distribution is subject to the
+// Boost Software License, Version 1.0. (See accompanying
+// file LICENSE-1.0 or http://www.boost.org/LICENSE-1.0)
+
+#ifndef HEADER_QRNG_TYPEDEFS_HPP_INCLUDED
+#define HEADER_QRNG_TYPEDEFS_HPP_INCLUDED
+
+#include <boost/random/niederreiter_base2.hpp>
+#include <boost/random/sobol.hpp>
+#include <boost/random/faure.hpp>
+
+
+// Boost.Random adaptors
+
+template<int D>
+struct niederreiter_base2_generator
+{
+  typedef boost::random::niederreiter_base2<boost::uint32_t, D, 1, 2*((boost::uint32_t)1 << 31)> type;
+};
+
+template<int D>
+struct sobol_generator
+{
+  typedef boost::random::sobol<boost::uint32_t, D, 1, (boost::uint32_t)1 << 31> type;
+};
+
+template<int D>
+struct faure_generator
+{
+  typedef boost::random::faure<double, D> type;
+};
+
+
+// To be used in place of uniform_real, where scaling is not necessary
+template<typename T>
+struct identity_distribution
+{
+  typedef T result_type;
+  template<typename QRNG>
+  result_type operator()(QRNG& q) const
+  { return q(); }
+};
+
+#endif // HEADER_QRNG_TYPEDEFS_HPP_INCLUDED