Worked on issue 131

Update boost
Spent a mutex for UUIDgen

1 files changed, 384 insertions, 384 deletions

diff --git a/contrib/src/boost/random/detail/polynomial.hpp b/contrib/src/boost/random/detail/polynomial.hpp
index a8c4b26..cc4ecd9 100644
--- a/contrib/src/boost/random/detail/polynomial.hpp
+++ b/contrib/src/boost/random/detail/polynomial.hpp
@@ -1,384 +1,384 @@
-/* boost random/detail/polynomial.hpp header file
- *
- * Copyright Steven Watanabe 2014
- * 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)
- *
- * See http://www.boost.org for most recent version including documentation.
- *
- * $Id$
- */
-
-#ifndef BOOST_RANDOM_DETAIL_POLYNOMIAL_HPP
-#define BOOST_RANDOM_DETAIL_POLYNOMIAL_HPP
-
-#include <cstddef>
-#include <limits>
-#include <vector>
-#include <algorithm>
-#include <boost/assert.hpp>
-#include <boost/cstdint.hpp>
-
-namespace boost {
-namespace random {
-namespace detail {
-
-class polynomial_ops {
-public:
-    typedef unsigned long digit_t;
-
-    static void add(std::size_t size, const digit_t * lhs,
-                       const digit_t * rhs, digit_t * output)
-    {
-        for(std::size_t i = 0; i < size; ++i) {
-            output[i] = lhs[i] ^ rhs[i];
-        }
-    }
-
-    static void add_shifted_inplace(std::size_t size, const digit_t * lhs,
-                                    digit_t * output, std::size_t shift)
-    {
-        if(shift == 0) {
-            add(size, lhs, output, output);
-            return;
-        }
-        std::size_t bits = std::numeric_limits<digit_t>::digits;
-        digit_t prev = 0;
-        for(std::size_t i = 0; i < size; ++i) {
-            digit_t tmp = lhs[i];
-            output[i] ^= (tmp << shift) | (prev >> (bits-shift));
-            prev = tmp;
-        }
-        output[size] ^= (prev >> (bits-shift));
-    }
-
-    static void multiply_simple(std::size_t size, const digit_t * lhs,
-                                   const digit_t * rhs, digit_t * output)
-    {
-        std::size_t bits = std::numeric_limits<digit_t>::digits;
-        for(std::size_t i = 0; i < 2*size; ++i) {
-            output[i] = 0;
-        }
-        for(std::size_t i = 0; i < size; ++i) {
-            for(std::size_t j = 0; j < bits; ++j) {
-                if((lhs[i] & (digit_t(1) << j)) != 0) {
-                    add_shifted_inplace(size, rhs, output + i, j);
-                }
-            }
-        }
-    }
-
-    // memory requirements: (size - cutoff) * 4 + next_smaller
-    static void multiply_karatsuba(std::size_t size,
-                               const digit_t * lhs, const digit_t * rhs,
-                               digit_t * output)
-    {
-        if(size < 64) {
-            multiply_simple(size, lhs, rhs, output);
-            return;
-        }
-        // split in half
-        std::size_t cutoff = size/2;
-        multiply_karatsuba(cutoff, lhs, rhs, output);
-        multiply_karatsuba(size - cutoff, lhs + cutoff, rhs + cutoff,
-                              output + cutoff*2);
-        std::vector<digit_t> local1(size - cutoff);
-        std::vector<digit_t> local2(size - cutoff);
-        // combine the digits for the inner multiply
-        add(cutoff, lhs, lhs + cutoff, &local1[0]);
-        if(size & 1) local1[cutoff] = lhs[size - 1];
-        add(cutoff, rhs + cutoff, rhs, &local2[0]);
-        if(size & 1) local2[cutoff] = rhs[size - 1];
-        std::vector<digit_t> local3((size - cutoff) * 2);
-        multiply_karatsuba(size - cutoff, &local1[0], &local2[0], &local3[0]);
-        add(cutoff * 2, output, &local3[0], &local3[0]);
-        add((size - cutoff) * 2, output + cutoff*2, &local3[0], &local3[0]);
-        // Finally, add the inner result
-        add((size - cutoff) * 2, output + cutoff, &local3[0], output + cutoff);
-    }
-    
-    static void multiply_add_karatsuba(std::size_t size,
-                                       const digit_t * lhs, const digit_t * rhs,
-                                       digit_t * output)
-    {
-        std::vector<digit_t> buf(size * 2);
-        multiply_karatsuba(size, lhs, rhs, &buf[0]);
-        add(size * 2, &buf[0], output, output);
-    }
-
-    static void multiply(const digit_t * lhs, std::size_t lhs_size,
-                         const digit_t * rhs, std::size_t rhs_size,
-                         digit_t * output)
-    {
-        std::fill_n(output, lhs_size + rhs_size, digit_t(0));
-        multiply_add(lhs, lhs_size, rhs, rhs_size, output);
-    }
-
-    static void multiply_add(const digit_t * lhs, std::size_t lhs_size,
-                             const digit_t * rhs, std::size_t rhs_size,
-                             digit_t * output)
-    {
-        // split into pieces that can be passed to
-        // karatsuba multiply.
-        while(lhs_size != 0) {
-            if(lhs_size < rhs_size) {
-                std::swap(lhs, rhs);
-                std::swap(lhs_size, rhs_size);
-            }
-            
-            multiply_add_karatsuba(rhs_size, lhs, rhs, output);
-            
-            lhs += rhs_size;
-            lhs_size -= rhs_size;
-            output += rhs_size;
-        }
-    }
-
-    static void copy_bits(const digit_t * x, std::size_t low, std::size_t high,
-                   digit_t * out)
-    {
-        const std::size_t bits = std::numeric_limits<digit_t>::digits;
-        std::size_t offset = low/bits;
-        x += offset;
-        low -= offset*bits;
-        high -= offset*bits;
-        std::size_t n = (high-low)/bits;
-        if(low == 0) {
-            for(std::size_t i = 0; i < n; ++i) {
-                out[i] = x[i];
-            }
-        } else {
-            for(std::size_t i = 0; i < n; ++i) {
-                out[i] = (x[i] >> low) | (x[i+1] << (bits-low));
-            }
-        }
-        if((high-low)%bits) {
-            digit_t low_mask = (digit_t(1) << ((high-low)%bits)) - 1;
-            digit_t result = (x[n] >> low);
-            if(low != 0 && (n+1)*bits < high) {
-                result |= (x[n+1] << (bits-low));
-            }
-            out[n] = (result & low_mask);
-        }
-    }
-
-    static void shift_left(digit_t * val, std::size_t size, std::size_t shift)
-    {
-        const std::size_t bits = std::numeric_limits<digit_t>::digits;
-        BOOST_ASSERT(shift > 0);
-        BOOST_ASSERT(shift < bits);
-        digit_t prev = 0;
-        for(std::size_t i = 0; i < size; ++i) {
-            digit_t tmp = val[i];
-            val[i] = (prev >> (bits - shift)) | (val[i] << shift);
-            prev = tmp;
-        }
-    }
-
-    static digit_t sqr(digit_t val) {
-        const std::size_t bits = std::numeric_limits<digit_t>::digits;
-        digit_t mask = (digit_t(1) << bits/2) - 1;
-        for(std::size_t i = bits; i > 1; i /= 2) {
-            val = ((val & ~mask) << i/2) | (val & mask);
-            mask = mask & (mask >> i/4);
-            mask = mask | (mask << i/2);
-        }
-        return val;
-    }
-
-    static void sqr(digit_t * val, std::size_t size)
-    {
-        const std::size_t bits = std::numeric_limits<digit_t>::digits;
-        digit_t mask = (digit_t(1) << bits/2) - 1;
-        for(std::size_t i = 0; i < size; ++i) {
-            digit_t x = val[size - i - 1];
-            val[(size - i - 1) * 2] = sqr(x & mask);
-            val[(size - i - 1) * 2 + 1] = sqr(x >> bits/2);
-        }
-    }
-
-    // optimized for the case when the modulus has few bits set.
-    struct sparse_mod {
-        sparse_mod(const digit_t * divisor, std::size_t divisor_bits)
-        {
-            const std::size_t bits = std::numeric_limits<digit_t>::digits;
-            _remainder_bits = divisor_bits - 1;
-            for(std::size_t i = 0; i < divisor_bits; ++i) {
-                if(divisor[i/bits] & (digit_t(1) << i%bits)) {
-                    _bit_indices.push_back(i);
-                }
-            }
-            BOOST_ASSERT(_bit_indices.back() == divisor_bits - 1);
-            _bit_indices.pop_back();
-            if(_bit_indices.empty()) {
-                _block_bits = divisor_bits;
-                _lower_bits = 0;
-            } else {
-                _block_bits = divisor_bits - _bit_indices.back() - 1;
-                _lower_bits = _bit_indices.back() + 1;
-            }
-            
-            _partial_quotient.resize((_block_bits + bits - 1)/bits);
-        }
-        void operator()(digit_t * dividend, std::size_t dividend_bits)
-        {
-            const std::size_t bits = std::numeric_limits<digit_t>::digits;
-            while(dividend_bits > _remainder_bits) {
-                std::size_t block_start = (std::max)(dividend_bits - _block_bits, _remainder_bits);
-                std::size_t block_size = (dividend_bits - block_start + bits - 1) / bits;
-                copy_bits(dividend, block_start, dividend_bits, &_partial_quotient[0]);
-                for(std::size_t i = 0; i < _bit_indices.size(); ++i) {
-                    std::size_t pos = _bit_indices[i] + block_start - _remainder_bits;
-                    add_shifted_inplace(block_size, &_partial_quotient[0], dividend + pos/bits, pos%bits);
-                }
-                add_shifted_inplace(block_size, &_partial_quotient[0], dividend + block_start/bits, block_start%bits);
-                dividend_bits = block_start;
-            }
-        }
-        std::vector<digit_t> _partial_quotient;
-        std::size_t _remainder_bits;
-        std::size_t _block_bits;
-        std::size_t _lower_bits;
-        std::vector<std::size_t> _bit_indices;
-    };
-
-    // base should have the same number of bits as mod
-    // base, and mod should both be able to hold a power
-    // of 2 >= mod_bits.  out needs to be twice as large.
-    static void mod_pow_x(boost::uintmax_t exponent, const digit_t * mod, std::size_t mod_bits, digit_t * out)
-    {
-        const std::size_t bits = std::numeric_limits<digit_t>::digits;
-        const std::size_t n = (mod_bits + bits - 1) / bits;
-        const std::size_t highbit = mod_bits - 1;
-        if(exponent == 0) {
-            out[0] = 1;
-            std::fill_n(out + 1, n - 1, digit_t(0));
-            return;
-        }
-        boost::uintmax_t i = std::numeric_limits<boost::uintmax_t>::digits - 1;
-        while(((boost::uintmax_t(1) << i) & exponent) == 0) {
-            --i;
-        }
-        out[0] = 2;
-        std::fill_n(out + 1, n - 1, digit_t(0));
-        sparse_mod m(mod, mod_bits);
-        while(i--) {
-            sqr(out, n);
-            m(out, 2 * mod_bits - 1);
-            if((boost::uintmax_t(1) << i) & exponent) {
-                shift_left(out, n, 1);
-                if(out[highbit / bits] & (digit_t(1) << highbit%bits))
-                    add(n, out, mod, out);
-            }
-        }
-    }
-};
-
-class polynomial
-{
-    typedef polynomial_ops::digit_t digit_t;
-public:
-    polynomial() : _size(0) {}
-    class reference {
-    public:
-        reference(digit_t &value, int idx)
-            : _value(value), _idx(idx) {}
-        operator bool() const { return (_value & (digit_t(1) << _idx)) != 0; }
-        reference& operator=(bool b)
-        {
-            if(b) {
-                _value |= (digit_t(1) << _idx);
-            } else {
-                _value &= ~(digit_t(1) << _idx);
-            }
-            return *this;
-        }
-        reference &operator^=(bool b)
-        {
-            _value ^= (digit_t(b) << _idx);
-            return *this;
-        }
-
-        reference &operator=(const reference &other)
-        {
-            return *this = static_cast<bool>(other);
-        }
-    private:
-        digit_t &_value;
-        int _idx;
-    };
-    reference operator[](std::size_t i)
-    {
-        static const std::size_t bits = std::numeric_limits<digit_t>::digits;
-        ensure_bit(i);
-        return reference(_storage[i/bits], i%bits);
-    }
-    bool operator[](std::size_t i) const
-    {
-        static const std::size_t bits = std::numeric_limits<digit_t>::digits;
-        if(i < size())
-            return (_storage[i/bits] & (digit_t(1) << (i%bits))) != 0;
-        else
-            return false;
-    }
-    std::size_t size() const
-    {
-        return _size;
-    }
-    void resize(std::size_t n)
-    {
-        static const std::size_t bits = std::numeric_limits<digit_t>::digits;
-        _storage.resize((n + bits - 1)/bits);
-        // clear the high order bits in case we're shrinking.
-        if(n%bits) {
-            _storage.back() &= ((digit_t(1) << (n%bits)) - 1);
-        }
-        _size = n;
-    }
-    friend polynomial operator*(const polynomial &lhs, const polynomial &rhs);
-    friend polynomial mod_pow_x(boost::uintmax_t exponent, polynomial mod);
-private:
-    std::vector<polynomial_ops::digit_t> _storage;
-    std::size_t _size;
-    void ensure_bit(std::size_t i)
-    {
-        if(i >= size()) {
-            resize(i + 1);
-        }
-    }
-    void normalize()
-    {
-        while(size() && (*this)[size() - 1] == 0)
-            resize(size() - 1);
-    }
-};
-
-inline polynomial operator*(const polynomial &lhs, const polynomial &rhs)
-{
-    polynomial result;
-    result._storage.resize(lhs._storage.size() + rhs._storage.size());
-    polynomial_ops::multiply(&lhs._storage[0], lhs._storage.size(),
-                             &rhs._storage[0], rhs._storage.size(),
-                             &result._storage[0]);
-    result._size = lhs._size + rhs._size;
-    return result;
-}
-
-inline polynomial mod_pow_x(boost::uintmax_t exponent, polynomial mod)
-{
-    polynomial result;
-    mod.normalize();
-    std::size_t mod_size = mod.size();
-    result._storage.resize(mod._storage.size() * 2);
-    result._size = mod.size() * 2;
-    polynomial_ops::mod_pow_x(exponent, &mod._storage[0], mod_size, &result._storage[0]);
-    result.resize(mod.size() - 1);
-    return result;
-}
-
-}
-}
-}
-
-#endif // BOOST_RANDOM_DETAIL_POLYNOMIAL_HPP
+/* boost random/detail/polynomial.hpp header file

+ *

+ * Copyright Steven Watanabe 2014

+ * 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)

+ *

+ * See http://www.boost.org for most recent version including documentation.

+ *

+ * $Id$

+ */

+

+#ifndef BOOST_RANDOM_DETAIL_POLYNOMIAL_HPP

+#define BOOST_RANDOM_DETAIL_POLYNOMIAL_HPP

+

+#include <cstddef>

+#include <limits>

+#include <vector>

+#include <algorithm>

+#include <boost/assert.hpp>

+#include <boost/cstdint.hpp>

+

+namespace boost {

+namespace random {

+namespace detail {

+

+class polynomial_ops {

+public:

+    typedef unsigned long digit_t;

+

+    static void add(std::size_t size, const digit_t * lhs,

+                       const digit_t * rhs, digit_t * output)

+    {

+        for(std::size_t i = 0; i < size; ++i) {

+            output[i] = lhs[i] ^ rhs[i];

+        }

+    }

+

+    static void add_shifted_inplace(std::size_t size, const digit_t * lhs,

+                                    digit_t * output, std::size_t shift)

+    {

+        if(shift == 0) {

+            add(size, lhs, output, output);

+            return;

+        }

+        std::size_t bits = std::numeric_limits<digit_t>::digits;

+        digit_t prev = 0;

+        for(std::size_t i = 0; i < size; ++i) {

+            digit_t tmp = lhs[i];

+            output[i] ^= (tmp << shift) | (prev >> (bits-shift));

+            prev = tmp;

+        }

+        output[size] ^= (prev >> (bits-shift));

+    }

+

+    static void multiply_simple(std::size_t size, const digit_t * lhs,

+                                   const digit_t * rhs, digit_t * output)

+    {

+        std::size_t bits = std::numeric_limits<digit_t>::digits;

+        for(std::size_t i = 0; i < 2*size; ++i) {

+            output[i] = 0;

+        }

+        for(std::size_t i = 0; i < size; ++i) {

+            for(std::size_t j = 0; j < bits; ++j) {

+                if((lhs[i] & (digit_t(1) << j)) != 0) {

+                    add_shifted_inplace(size, rhs, output + i, j);

+                }

+            }

+        }

+    }

+

+    // memory requirements: (size - cutoff) * 4 + next_smaller

+    static void multiply_karatsuba(std::size_t size,

+                               const digit_t * lhs, const digit_t * rhs,

+                               digit_t * output)

+    {

+        if(size < 64) {

+            multiply_simple(size, lhs, rhs, output);

+            return;

+        }

+        // split in half

+        std::size_t cutoff = size/2;

+        multiply_karatsuba(cutoff, lhs, rhs, output);

+        multiply_karatsuba(size - cutoff, lhs + cutoff, rhs + cutoff,

+                              output + cutoff*2);

+        std::vector<digit_t> local1(size - cutoff);

+        std::vector<digit_t> local2(size - cutoff);

+        // combine the digits for the inner multiply

+        add(cutoff, lhs, lhs + cutoff, &local1[0]);

+        if(size & 1) local1[cutoff] = lhs[size - 1];

+        add(cutoff, rhs + cutoff, rhs, &local2[0]);

+        if(size & 1) local2[cutoff] = rhs[size - 1];

+        std::vector<digit_t> local3((size - cutoff) * 2);

+        multiply_karatsuba(size - cutoff, &local1[0], &local2[0], &local3[0]);

+        add(cutoff * 2, output, &local3[0], &local3[0]);

+        add((size - cutoff) * 2, output + cutoff*2, &local3[0], &local3[0]);

+        // Finally, add the inner result

+        add((size - cutoff) * 2, output + cutoff, &local3[0], output + cutoff);

+    }

+    

+    static void multiply_add_karatsuba(std::size_t size,

+                                       const digit_t * lhs, const digit_t * rhs,

+                                       digit_t * output)

+    {

+        std::vector<digit_t> buf(size * 2);

+        multiply_karatsuba(size, lhs, rhs, &buf[0]);

+        add(size * 2, &buf[0], output, output);

+    }

+

+    static void multiply(const digit_t * lhs, std::size_t lhs_size,

+                         const digit_t * rhs, std::size_t rhs_size,

+                         digit_t * output)

+    {

+        std::fill_n(output, lhs_size + rhs_size, digit_t(0));

+        multiply_add(lhs, lhs_size, rhs, rhs_size, output);

+    }

+

+    static void multiply_add(const digit_t * lhs, std::size_t lhs_size,

+                             const digit_t * rhs, std::size_t rhs_size,

+                             digit_t * output)

+    {

+        // split into pieces that can be passed to

+        // karatsuba multiply.

+        while(lhs_size != 0) {

+            if(lhs_size < rhs_size) {

+                std::swap(lhs, rhs);

+                std::swap(lhs_size, rhs_size);

+            }

+            

+            multiply_add_karatsuba(rhs_size, lhs, rhs, output);

+            

+            lhs += rhs_size;

+            lhs_size -= rhs_size;

+            output += rhs_size;

+        }

+    }

+

+    static void copy_bits(const digit_t * x, std::size_t low, std::size_t high,

+                   digit_t * out)

+    {

+        const std::size_t bits = std::numeric_limits<digit_t>::digits;

+        std::size_t offset = low/bits;

+        x += offset;

+        low -= offset*bits;

+        high -= offset*bits;

+        std::size_t n = (high-low)/bits;

+        if(low == 0) {

+            for(std::size_t i = 0; i < n; ++i) {

+                out[i] = x[i];

+            }

+        } else {

+            for(std::size_t i = 0; i < n; ++i) {

+                out[i] = (x[i] >> low) | (x[i+1] << (bits-low));

+            }

+        }

+        if((high-low)%bits) {

+            digit_t low_mask = (digit_t(1) << ((high-low)%bits)) - 1;

+            digit_t result = (x[n] >> low);

+            if(low != 0 && (n+1)*bits < high) {

+                result |= (x[n+1] << (bits-low));

+            }

+            out[n] = (result & low_mask);

+        }

+    }

+

+    static void shift_left(digit_t * val, std::size_t size, std::size_t shift)

+    {

+        const std::size_t bits = std::numeric_limits<digit_t>::digits;

+        BOOST_ASSERT(shift > 0);

+        BOOST_ASSERT(shift < bits);

+        digit_t prev = 0;

+        for(std::size_t i = 0; i < size; ++i) {

+            digit_t tmp = val[i];

+            val[i] = (prev >> (bits - shift)) | (val[i] << shift);

+            prev = tmp;

+        }

+    }

+

+    static digit_t sqr(digit_t val) {

+        const std::size_t bits = std::numeric_limits<digit_t>::digits;

+        digit_t mask = (digit_t(1) << bits/2) - 1;

+        for(std::size_t i = bits; i > 1; i /= 2) {

+            val = ((val & ~mask) << i/2) | (val & mask);

+            mask = mask & (mask >> i/4);

+            mask = mask | (mask << i/2);

+        }

+        return val;

+    }

+

+    static void sqr(digit_t * val, std::size_t size)

+    {

+        const std::size_t bits = std::numeric_limits<digit_t>::digits;

+        digit_t mask = (digit_t(1) << bits/2) - 1;

+        for(std::size_t i = 0; i < size; ++i) {

+            digit_t x = val[size - i - 1];

+            val[(size - i - 1) * 2] = sqr(x & mask);

+            val[(size - i - 1) * 2 + 1] = sqr(x >> bits/2);

+        }

+    }

+

+    // optimized for the case when the modulus has few bits set.

+    struct sparse_mod {

+        sparse_mod(const digit_t * divisor, std::size_t divisor_bits)

+        {

+            const std::size_t bits = std::numeric_limits<digit_t>::digits;

+            _remainder_bits = divisor_bits - 1;

+            for(std::size_t i = 0; i < divisor_bits; ++i) {

+                if(divisor[i/bits] & (digit_t(1) << i%bits)) {

+                    _bit_indices.push_back(i);

+                }

+            }

+            BOOST_ASSERT(_bit_indices.back() == divisor_bits - 1);

+            _bit_indices.pop_back();

+            if(_bit_indices.empty()) {

+                _block_bits = divisor_bits;

+                _lower_bits = 0;

+            } else {

+                _block_bits = divisor_bits - _bit_indices.back() - 1;

+                _lower_bits = _bit_indices.back() + 1;

+            }

+            

+            _partial_quotient.resize((_block_bits + bits - 1)/bits);

+        }

+        void operator()(digit_t * dividend, std::size_t dividend_bits)

+        {

+            const std::size_t bits = std::numeric_limits<digit_t>::digits;

+            while(dividend_bits > _remainder_bits) {

+                std::size_t block_start = (std::max)(dividend_bits - _block_bits, _remainder_bits);

+                std::size_t block_size = (dividend_bits - block_start + bits - 1) / bits;

+                copy_bits(dividend, block_start, dividend_bits, &_partial_quotient[0]);

+                for(std::size_t i = 0; i < _bit_indices.size(); ++i) {

+                    std::size_t pos = _bit_indices[i] + block_start - _remainder_bits;

+                    add_shifted_inplace(block_size, &_partial_quotient[0], dividend + pos/bits, pos%bits);

+                }

+                add_shifted_inplace(block_size, &_partial_quotient[0], dividend + block_start/bits, block_start%bits);

+                dividend_bits = block_start;

+            }

+        }

+        std::vector<digit_t> _partial_quotient;

+        std::size_t _remainder_bits;

+        std::size_t _block_bits;

+        std::size_t _lower_bits;

+        std::vector<std::size_t> _bit_indices;

+    };

+

+    // base should have the same number of bits as mod

+    // base, and mod should both be able to hold a power

+    // of 2 >= mod_bits.  out needs to be twice as large.

+    static void mod_pow_x(boost::uintmax_t exponent, const digit_t * mod, std::size_t mod_bits, digit_t * out)

+    {

+        const std::size_t bits = std::numeric_limits<digit_t>::digits;

+        const std::size_t n = (mod_bits + bits - 1) / bits;

+        const std::size_t highbit = mod_bits - 1;

+        if(exponent == 0) {

+            out[0] = 1;

+            std::fill_n(out + 1, n - 1, digit_t(0));

+            return;

+        }

+        boost::uintmax_t i = std::numeric_limits<boost::uintmax_t>::digits - 1;

+        while(((boost::uintmax_t(1) << i) & exponent) == 0) {

+            --i;

+        }

+        out[0] = 2;

+        std::fill_n(out + 1, n - 1, digit_t(0));

+        sparse_mod m(mod, mod_bits);

+        while(i--) {

+            sqr(out, n);

+            m(out, 2 * mod_bits - 1);

+            if((boost::uintmax_t(1) << i) & exponent) {

+                shift_left(out, n, 1);

+                if(out[highbit / bits] & (digit_t(1) << highbit%bits))

+                    add(n, out, mod, out);

+            }

+        }

+    }

+};

+

+class polynomial

+{

+    typedef polynomial_ops::digit_t digit_t;

+public:

+    polynomial() : _size(0) {}

+    class reference {

+    public:

+        reference(digit_t &value, int idx)

+            : _value(value), _idx(idx) {}

+        operator bool() const { return (_value & (digit_t(1) << _idx)) != 0; }

+        reference& operator=(bool b)

+        {

+            if(b) {

+                _value |= (digit_t(1) << _idx);

+            } else {

+                _value &= ~(digit_t(1) << _idx);

+            }

+            return *this;

+        }

+        reference &operator^=(bool b)

+        {

+            _value ^= (digit_t(b) << _idx);

+            return *this;

+        }

+

+        reference &operator=(const reference &other)

+        {

+            return *this = static_cast<bool>(other);

+        }

+    private:

+        digit_t &_value;

+        int _idx;

+    };

+    reference operator[](std::size_t i)

+    {

+        static const std::size_t bits = std::numeric_limits<digit_t>::digits;

+        ensure_bit(i);

+        return reference(_storage[i/bits], i%bits);

+    }

+    bool operator[](std::size_t i) const

+    {

+        static const std::size_t bits = std::numeric_limits<digit_t>::digits;

+        if(i < size())

+            return (_storage[i/bits] & (digit_t(1) << (i%bits))) != 0;

+        else

+            return false;

+    }

+    std::size_t size() const

+    {

+        return _size;

+    }

+    void resize(std::size_t n)

+    {

+        static const std::size_t bits = std::numeric_limits<digit_t>::digits;

+        _storage.resize((n + bits - 1)/bits);

+        // clear the high order bits in case we're shrinking.

+        if(n%bits) {

+            _storage.back() &= ((digit_t(1) << (n%bits)) - 1);

+        }

+        _size = n;

+    }

+    friend polynomial operator*(const polynomial &lhs, const polynomial &rhs);

+    friend polynomial mod_pow_x(boost::uintmax_t exponent, polynomial mod);

+private:

+    std::vector<polynomial_ops::digit_t> _storage;

+    std::size_t _size;

+    void ensure_bit(std::size_t i)

+    {

+        if(i >= size()) {

+            resize(i + 1);

+        }

+    }

+    void normalize()

+    {

+        while(size() && (*this)[size() - 1] == 0)

+            resize(size() - 1);

+    }

+};

+

+inline polynomial operator*(const polynomial &lhs, const polynomial &rhs)

+{

+    polynomial result;

+    result._storage.resize(lhs._storage.size() + rhs._storage.size());

+    polynomial_ops::multiply(&lhs._storage[0], lhs._storage.size(),

+                             &rhs._storage[0], rhs._storage.size(),

+                             &result._storage[0]);

+    result._size = lhs._size + rhs._size;

+    return result;

+}

+

+inline polynomial mod_pow_x(boost::uintmax_t exponent, polynomial mod)

+{

+    polynomial result;

+    mod.normalize();

+    std::size_t mod_size = mod.size();

+    result._storage.resize(mod._storage.size() * 2);

+    result._size = mod.size() * 2;

+    polynomial_ops::mod_pow_x(exponent, &mod._storage[0], mod_size, &result._storage[0]);

+    result.resize(mod.size() - 1);

+    return result;

+}

+

+}

+}

+}

+

+#endif // BOOST_RANDOM_DETAIL_POLYNOMIAL_HPP