math/tools/polynomial_gcd.hpp

0001 //  (C) Copyright Jeremy William Murphy 2016.
0002 //  (C) Copyright Matt Borland 2021.
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP
0008 #define BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP
0009
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #endif
0013
0014 #include <algorithm>
0015 #include <type_traits>
0016 #include <boost/math/tools/is_standalone.hpp>
0017 #include <boost/math/tools/polynomial.hpp>
0018
0019 #ifndef BOOST_MATH_STANDALONE
0020 #include <boost/integer/common_factor_rt.hpp>
0021
0022 #else
0023 #include <numeric>
0024 #include <utility>
0025 #include <iterator>
0026 #include <boost/math/tools/assert.hpp>
0027 #include <boost/math/tools/config.hpp>
0028
0029 namespace boost { namespace integer {
0030
0031 namespace gcd_detail {
0032
0033 template <typename EuclideanDomain>
0034 inline EuclideanDomain Euclid_gcd(EuclideanDomain a, EuclideanDomain b) noexcept(std::is_arithmetic<EuclideanDomain>::value)
0035 {
0036     using std::swap;
0037     while (b != EuclideanDomain(0))
0038     {
0039         a %= b;
0040         swap(a, b);
0041     }
0042     return a;
0043 }
0044
0045 enum method_type
0046 {
0047     method_euclid = 0,
0048     method_binary = 1,
0049     method_mixed = 2
0050 };
0051
0052 } // gcd_detail
0053
0054 template <typename Iter, typename T = typename std::iterator_traits<Iter>::value_type>
0055 std::pair<T, Iter> gcd_range(Iter first, Iter last) noexcept(std::is_arithmetic<T>::value)
0056 {
0057     BOOST_MATH_ASSERT(first != last);
0058
0059     T d = *first;
0060     ++first;
0061     while (d != T(1) && first != last)
0062     {
0063         #ifdef BOOST_MATH_HAS_CXX17_NUMERIC
0064         d = std::gcd(d, *first);
0065         #else
0066         d = gcd_detail::Euclid_gcd(d, *first);
0067         #endif
0068         ++first;
0069     }
0070     return std::make_pair(d, first);
0071 }
0072
0073 }} // namespace boost::integer
0074 #endif
0075
0076 namespace boost{
0077
0078    namespace integer {
0079
0080       namespace gcd_detail {
0081
0082          template <class T>
0083          struct gcd_traits;
0084
0085          template <class T>
0086          struct gcd_traits<boost::math::tools::polynomial<T> >
0087          {
0088             inline static const boost::math::tools::polynomial<T>& abs(const boost::math::tools::polynomial<T>& val) { return val; }
0089
0090             static const method_type method = method_euclid;
0091          };
0092
0093       }
0094 }
0095
0096 namespace math{ namespace tools{
0097
0098 /* From Knuth, 4.6.1:
0099 *
0100 * We may write any nonzero polynomial u(x) from R[x] where R is a UFD as
0101 *
0102 *      u(x) = cont(u) . pp(u(x))
0103 *
0104 * where cont(u), the content of u, is an element of S, and pp(u(x)), the primitive
0105 * part of u(x), is a primitive polynomial over S.
0106 * When u(x) = 0, it is convenient to define cont(u) = pp(u(x)) = O.
0107 */
0108
0109 template <class T>
0110 T content(polynomial<T> const &x)
0111 {
0112     return x ? boost::integer::gcd_range(x.data().begin(), x.data().end()).first : T(0);
0113 }
0114
0115 // Knuth, 4.6.1
0116 template <class T>
0117 polynomial<T> primitive_part(polynomial<T> const &x, T const &cont)
0118 {
0119     return x ? x / cont : polynomial<T>();
0120 }
0121
0122
0123 template <class T>
0124 polynomial<T> primitive_part(polynomial<T> const &x)
0125 {
0126     return primitive_part(x, content(x));
0127 }
0128
0129
0130 // Trivial but useful convenience function referred to simply as l() in Knuth.
0131 template <class T>
0132 T leading_coefficient(polynomial<T> const &x)
0133 {
0134     return x ? x.data().back() : T(0);
0135 }
0136
0137
0138 namespace detail
0139 {
0140     /* Reduce u and v to their primitive parts and return the gcd of their
0141     * contents. Used in a couple of gcd algorithms.
0142     */
0143     template <class T>
0144     T reduce_to_primitive(polynomial<T> &u, polynomial<T> &v)
0145     {
0146         T const u_cont = content(u), v_cont = content(v);
0147         u /= u_cont;
0148         v /= v_cont;
0149
0150         #ifdef BOOST_MATH_HAS_CXX17_NUMERIC
0151         return std::gcd(u_cont, v_cont);
0152         #else
0153         return boost::integer::gcd_detail::Euclid_gcd(u_cont, v_cont);
0154         #endif
0155     }
0156 }
0157
0158
0159 /**
0160 * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
0161 * Algorithm 4.6.1C: Greatest common divisor over a unique factorization domain.
0162 *
0163 * The subresultant algorithm by George E. Collins [JACM 14 (1967), 128-142],
0164 * later improved by W. S. Brown and J. F. Traub [JACM 18 (1971), 505-514].
0165 *
0166 * Although step C3 keeps the coefficients to a "reasonable" size, they are
0167 * still potentially several binary orders of magnitude larger than the inputs.
0168 * Thus, this algorithm should only be used where T is a multi-precision type.
0169 *
0170 * @tparam   T   Polynomial coefficient type.
0171 * @param    u   First polynomial.
0172 * @param    v   Second polynomial.
0173 * @return       Greatest common divisor of polynomials u and v.
0174 */
0175 template <class T>
0176 typename std::enable_if< std::numeric_limits<T>::is_integer, polynomial<T> >::type
0177 subresultant_gcd(polynomial<T> u, polynomial<T> v)
0178 {
0179     using std::swap;
0180     BOOST_MATH_ASSERT(u || v);
0181
0182     if (!u)
0183         return v;
0184     if (!v)
0185         return u;
0186
0187     typedef typename polynomial<T>::size_type N;
0188
0189     if (u.degree() < v.degree())
0190         swap(u, v);
0191
0192     T const d = detail::reduce_to_primitive(u, v);
0193     T g = 1, h = 1;
0194     polynomial<T> r;
0195     while (true)
0196     {
0197         BOOST_MATH_ASSERT(u.degree() >= v.degree());
0198         // Pseudo-division.
0199         r = u % v;
0200         if (!r)
0201             return d * primitive_part(v); // Attach the content.
0202         if (r.degree() == 0)
0203             return d * polynomial<T>(T(1)); // The content is the result.
0204         N const delta = u.degree() - v.degree();
0205         // Adjust remainder.
0206         u = v;
0207         v = r / (g * detail::integer_power(h, delta));
0208         g = leading_coefficient(u);
0209         T const tmp = detail::integer_power(g, delta);
0210         if (delta <= N(1))
0211             h = tmp * detail::integer_power(h, N(1) - delta);
0212         else
0213             h = tmp / detail::integer_power(h, delta - N(1));
0214     }
0215 }
0216
0217
0218 /**
0219  * @brief GCD for polynomials with unbounded multi-precision integral coefficients.
0220  *
0221  * The multi-precision constraint is enforced via numeric_limits.
0222  *
0223  * Note that intermediate terms in the evaluation can grow arbitrarily large, hence the need for
0224  * unbounded integers, otherwise numeric overflow would break the algorithm.
0225  *
0226  * @tparam  T   A multi-precision integral type.
0227  */
0228 template <typename T>
0229 typename std::enable_if<std::numeric_limits<T>::is_integer && !std::numeric_limits<T>::is_bounded, polynomial<T> >::type
0230 gcd(polynomial<T> const &u, polynomial<T> const &v)
0231 {
0232     return subresultant_gcd(u, v);
0233 }
0234 // GCD over bounded integers is not currently allowed:
0235 template <typename T>
0236 typename std::enable_if<std::numeric_limits<T>::is_integer && std::numeric_limits<T>::is_bounded, polynomial<T> >::type
0237 gcd(polynomial<T> const &u, polynomial<T> const &v)
0238 {
0239    static_assert(sizeof(v) == 0, "GCD on polynomials of bounded integers is disallowed due to the excessive growth in the size of intermediate terms.");
0240    return subresultant_gcd(u, v);
0241 }
0242 // GCD over polynomials of floats can go via the Euclid algorithm:
0243 template <typename T>
0244 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
0245 gcd(polynomial<T> const &u, polynomial<T> const &v)
0246 {
0247     return boost::integer::gcd_detail::Euclid_gcd(u, v);
0248 }
0249
0250 }
0251 //
0252 // Using declaration so we overload the default implementation in this namespace:
0253 //
0254 using boost::math::tools::gcd;
0255
0256 }
0257
0258 namespace integer
0259 {
0260    //
0261    // Using declaration so we overload the default implementation in this namespace:
0262    //
0263    using boost::math::tools::gcd;
0264 }
0265
0266 } // namespace boost::math::tools
0267
0268 #endif