New default microstepper implementation

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 cc4ecd9..a8c4b26 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