123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268 |
- #ifndef BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP
- #define BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP
- #ifdef _MSC_VER
- #pragma once
- #endif
- #include <algorithm>
- #include <type_traits>
- #include <boost/math/tools/is_standalone.hpp>
- #include <boost/math/tools/polynomial.hpp>
- #ifndef BOOST_MATH_STANDALONE
- #include <boost/integer/common_factor_rt.hpp>
- #else
- #include <numeric>
- #include <utility>
- #include <iterator>
- #include <boost/math/tools/assert.hpp>
- #include <boost/math/tools/config.hpp>
- namespace boost { namespace integer {
- namespace gcd_detail {
- template <typename EuclideanDomain>
- inline EuclideanDomain Euclid_gcd(EuclideanDomain a, EuclideanDomain b) noexcept(std::is_arithmetic<EuclideanDomain>::value)
- {
- using std::swap;
- while (b != EuclideanDomain(0))
- {
- a %= b;
- swap(a, b);
- }
- return a;
- }
- enum method_type
- {
- method_euclid = 0,
- method_binary = 1,
- method_mixed = 2
- };
- }
- template <typename Iter, typename T = typename std::iterator_traits<Iter>::value_type>
- std::pair<T, Iter> gcd_range(Iter first, Iter last) noexcept(std::is_arithmetic<T>::value)
- {
- BOOST_MATH_ASSERT(first != last);
- T d = *first;
- ++first;
- while (d != T(1) && first != last)
- {
- #ifdef BOOST_MATH_HAS_CXX17_NUMERIC
- d = std::gcd(d, *first);
- #else
- d = gcd_detail::Euclid_gcd(d, *first);
- #endif
- ++first;
- }
- return std::make_pair(d, first);
- }
- }}
- #endif
- namespace boost{
- namespace integer {
- namespace gcd_detail {
- template <class T>
- struct gcd_traits;
- template <class T>
- struct gcd_traits<boost::math::tools::polynomial<T> >
- {
- inline static const boost::math::tools::polynomial<T>& abs(const boost::math::tools::polynomial<T>& val) { return val; }
- static const method_type method = method_euclid;
- };
- }
- }
- namespace math{ namespace tools{
- template <class T>
- T content(polynomial<T> const &x)
- {
- return x ? boost::integer::gcd_range(x.data().begin(), x.data().end()).first : T(0);
- }
- template <class T>
- polynomial<T> primitive_part(polynomial<T> const &x, T const &cont)
- {
- return x ? x / cont : polynomial<T>();
- }
- template <class T>
- polynomial<T> primitive_part(polynomial<T> const &x)
- {
- return primitive_part(x, content(x));
- }
- template <class T>
- T leading_coefficient(polynomial<T> const &x)
- {
- return x ? x.data().back() : T(0);
- }
- namespace detail
- {
-
- template <class T>
- T reduce_to_primitive(polynomial<T> &u, polynomial<T> &v)
- {
- T const u_cont = content(u), v_cont = content(v);
- u /= u_cont;
- v /= v_cont;
- #ifdef BOOST_MATH_HAS_CXX17_NUMERIC
- return std::gcd(u_cont, v_cont);
- #else
- return boost::integer::gcd_detail::Euclid_gcd(u_cont, v_cont);
- #endif
- }
- }
- template <class T>
- typename std::enable_if< std::numeric_limits<T>::is_integer, polynomial<T> >::type
- subresultant_gcd(polynomial<T> u, polynomial<T> v)
- {
- using std::swap;
- BOOST_MATH_ASSERT(u || v);
- if (!u)
- return v;
- if (!v)
- return u;
- typedef typename polynomial<T>::size_type N;
- if (u.degree() < v.degree())
- swap(u, v);
- T const d = detail::reduce_to_primitive(u, v);
- T g = 1, h = 1;
- polynomial<T> r;
- while (true)
- {
- BOOST_MATH_ASSERT(u.degree() >= v.degree());
-
- r = u % v;
- if (!r)
- return d * primitive_part(v);
- if (r.degree() == 0)
- return d * polynomial<T>(T(1));
- N const delta = u.degree() - v.degree();
-
- u = v;
- v = r / (g * detail::integer_power(h, delta));
- g = leading_coefficient(u);
- T const tmp = detail::integer_power(g, delta);
- if (delta <= N(1))
- h = tmp * detail::integer_power(h, N(1) - delta);
- else
- h = tmp / detail::integer_power(h, delta - N(1));
- }
- }
- template <typename T>
- typename std::enable_if<std::numeric_limits<T>::is_integer && !std::numeric_limits<T>::is_bounded, polynomial<T> >::type
- gcd(polynomial<T> const &u, polynomial<T> const &v)
- {
- return subresultant_gcd(u, v);
- }
- template <typename T>
- typename std::enable_if<std::numeric_limits<T>::is_integer && std::numeric_limits<T>::is_bounded, polynomial<T> >::type
- gcd(polynomial<T> const &u, polynomial<T> const &v)
- {
- static_assert(sizeof(v) == 0, "GCD on polynomials of bounded integers is disallowed due to the excessive growth in the size of intermediate terms.");
- return subresultant_gcd(u, v);
- }
- template <typename T>
- typename std::enable_if<!std::numeric_limits<T>::is_integer && (std::numeric_limits<T>::min_exponent != std::numeric_limits<T>::max_exponent) && !std::numeric_limits<T>::is_exact, polynomial<T> >::type
- gcd(polynomial<T> const &u, polynomial<T> const &v)
- {
- return boost::integer::gcd_detail::Euclid_gcd(u, v);
- }
- }
- using boost::math::tools::gcd;
- }
- namespace integer
- {
-
-
-
- using boost::math::tools::gcd;
- }
- }
- #endif
|