diff options
Diffstat (limited to 'boost/math/tools/polynomial_gcd.hpp')
-rw-r--r-- | boost/math/tools/polynomial_gcd.hpp | 209 |
1 files changed, 209 insertions, 0 deletions
diff --git a/boost/math/tools/polynomial_gcd.hpp b/boost/math/tools/polynomial_gcd.hpp new file mode 100644 index 0000000000..fdbafda6ca --- /dev/null +++ b/boost/math/tools/polynomial_gcd.hpp @@ -0,0 +1,209 @@ +// (C) Copyright Jeremy William Murphy 2016. + +// Use, modification and distribution are subject to 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) + +#ifndef BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP +#define BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include <boost/math/tools/polynomial.hpp> +#include <boost/math/common_factor_rt.hpp> +#include <boost/type_traits/is_pod.hpp> + + +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{ + +/* From Knuth, 4.6.1: +* +* We may write any nonzero polynomial u(x) from R[x] where R is a UFD as +* +* u(x) = cont(u) ยท pp(u(x)) +* +* where cont(u), the content of u, is an element of S, and pp(u(x)), the primitive +* part of u(x), is a primitive polynomial over S. +* When u(x) = 0, it is convenient to define cont(u) = pp(u(x)) = O. +*/ + +template <class T> +T content(polynomial<T> const &x) +{ + return x ? gcd_range(x.data().begin(), x.data().end()).first : T(0); +} + +// Knuth, 4.6.1 +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)); +} + + +// Trivial but useful convenience function referred to simply as l() in Knuth. +template <class T> +T leading_coefficient(polynomial<T> const &x) +{ + return x ? x.data().back() : T(0); +} + + +namespace detail +{ + /* Reduce u and v to their primitive parts and return the gcd of their + * contents. Used in a couple of gcd algorithms. + */ + template <class T> + T reduce_to_primitive(polynomial<T> &u, polynomial<T> &v) + { + using boost::math::gcd; + T const u_cont = content(u), v_cont = content(v); + u /= u_cont; + v /= v_cont; + return gcd(u_cont, v_cont); + } +} + + +/** +* Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998 +* Algorithm 4.6.1C: Greatest common divisor over a unique factorization domain. +* +* The subresultant algorithm by George E. Collins [JACM 14 (1967), 128-142], +* later improved by W. S. Brown and J. F. Traub [JACM 18 (1971), 505-514]. +* +* Although step C3 keeps the coefficients to a "reasonable" size, they are +* still potentially several binary orders of magnitude larger than the inputs. +* Thus, this algorithm should only be used where T is a multi-precision type. +* +* @tparam T Polynomial coefficient type. +* @param u First polynomial. +* @param v Second polynomial. +* @return Greatest common divisor of polynomials u and v. +*/ +template <class T> +typename enable_if_c< std::numeric_limits<T>::is_integer, polynomial<T> >::type +subresultant_gcd(polynomial<T> u, polynomial<T> v) +{ + using std::swap; + BOOST_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_ASSERT(u.degree() >= v.degree()); + // Pseudo-division. + r = u % v; + if (!r) + return d * primitive_part(v); // Attach the content. + if (r.degree() == 0) + return d * polynomial<T>(T(1)); // The content is the result. + N const delta = u.degree() - v.degree(); + // Adjust remainder. + 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)); + } +} + + +/** + * @brief GCD for polynomials with unbounded multi-precision integral coefficients. + * + * The multi-precision constraint is enforced via numeric_limits. + * + * Note that intermediate terms in the evaluation can grow arbitrarily large, hence the need for + * unbounded integers, otherwise numeric loverflow would break the algorithm. + * + * @tparam T A multi-precision integral type. + */ +template <typename T> +typename enable_if_c<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); +} +// GCD over bounded integers is not currently allowed: +template <typename T> +typename enable_if_c<std::numeric_limits<T>::is_integer && std::numeric_limits<T>::is_bounded, polynomial<T> >::type +gcd(polynomial<T> const &u, polynomial<T> const &v) +{ + BOOST_STATIC_ASSERT_MSG(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); +} +// GCD over polynomials of floats can go via the Euclid algorithm: +template <typename T> +typename enable_if_c<!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 declaration so we overload the default implementation in this namespace: +// +using boost::math::tools::gcd; + +} + +namespace integer +{ + // + // Using declaration so we overload the default implementation in this namespace: + // + using boost::math::tools::gcd; +} + +} // namespace boost::math::tools + +#endif |