boost/boost/random/detail/polynomial.hpp - nest-cam/4320010/boost - Git at Google

 /* 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