diff options
Diffstat (limited to 'contrib/src/boost/random/detail/polynomial.hpp')
| -rw-r--r-- | contrib/src/boost/random/detail/polynomial.hpp | 768 |
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 |