269 lines · plain
1// (C) Copyright Jeremy William Murphy 2016.2// (C) Copyright Matt Borland 2021.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP8#define BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP9 10#ifdef _MSC_VER11#pragma once12#endif13 14#include <algorithm>15#include <type_traits>16#include <boost/math/tools/is_standalone.hpp>17#include <boost/math/tools/polynomial.hpp>18 19#ifndef BOOST_MATH_STANDALONE20#include <boost/integer/common_factor_rt.hpp>21 22#else23#include <numeric>24#include <utility>25#include <iterator>26#include <boost/math/tools/assert.hpp>27#include <boost/math/tools/config.hpp>28 29namespace boost { namespace integer {30 31namespace gcd_detail {32 33template <typename EuclideanDomain>34inline EuclideanDomain Euclid_gcd(EuclideanDomain a, EuclideanDomain b) noexcept(std::is_arithmetic<EuclideanDomain>::value)35{36 using std::swap;37 while (b != EuclideanDomain(0))38 {39 a %= b;40 swap(a, b);41 }42 return a;43}44 45enum method_type46{47 method_euclid = 0,48 method_binary = 1,49 method_mixed = 250};51 52} // gcd_detail53 54template <typename Iter, typename T = typename std::iterator_traits<Iter>::value_type>55std::pair<T, Iter> gcd_range(Iter first, Iter last) noexcept(std::is_arithmetic<T>::value)56{57 BOOST_MATH_ASSERT(first != last);58 59 T d = *first;60 ++first;61 while (d != T(1) && first != last)62 {63 #ifdef BOOST_MATH_HAS_CXX17_NUMERIC64 d = std::gcd(d, *first);65 #else66 d = gcd_detail::Euclid_gcd(d, *first);67 #endif68 ++first;69 }70 return std::make_pair(d, first);71}72 73}} // namespace boost::integer74#endif75 76namespace boost{77 78 namespace integer {79 80 namespace gcd_detail {81 82 template <class T>83 struct gcd_traits;84 85 template <class T>86 struct gcd_traits<boost::math::tools::polynomial<T> >87 {88 inline static const boost::math::tools::polynomial<T>& abs(const boost::math::tools::polynomial<T>& val) { return val; }89 90 static const method_type method = method_euclid;91 };92 93 }94}95 96namespace math{ namespace tools{97 98/* From Knuth, 4.6.1:99*100* We may write any nonzero polynomial u(x) from R[x] where R is a UFD as101*102* u(x) = cont(u) . pp(u(x))103*104* where cont(u), the content of u, is an element of S, and pp(u(x)), the primitive105* part of u(x), is a primitive polynomial over S.106* When u(x) = 0, it is convenient to define cont(u) = pp(u(x)) = O.107*/108 109template <class T>110T content(polynomial<T> const &x)111{112 return x ? boost::integer::gcd_range(x.data().begin(), x.data().end()).first : T(0);113}114 115// Knuth, 4.6.1116template <class T>117polynomial<T> primitive_part(polynomial<T> const &x, T const &cont)118{119 return x ? x / cont : polynomial<T>();120}121 122 123template <class T>124polynomial<T> primitive_part(polynomial<T> const &x)125{126 return primitive_part(x, content(x));127}128 129 130// Trivial but useful convenience function referred to simply as l() in Knuth.131template <class T>132T leading_coefficient(polynomial<T> const &x)133{134 return x ? x.data().back() : T(0);135}136 137 138namespace detail139{140 /* Reduce u and v to their primitive parts and return the gcd of their141 * contents. Used in a couple of gcd algorithms.142 */143 template <class T>144 T reduce_to_primitive(polynomial<T> &u, polynomial<T> &v)145 {146 T const u_cont = content(u), v_cont = content(v);147 u /= u_cont;148 v /= v_cont;149 150 #ifdef BOOST_MATH_HAS_CXX17_NUMERIC151 return std::gcd(u_cont, v_cont);152 #else153 return boost::integer::gcd_detail::Euclid_gcd(u_cont, v_cont);154 #endif155 }156}157 158 159/**160* Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998161* Algorithm 4.6.1C: Greatest common divisor over a unique factorization domain.162*163* The subresultant algorithm by George E. Collins [JACM 14 (1967), 128-142],164* later improved by W. S. Brown and J. F. Traub [JACM 18 (1971), 505-514].165*166* Although step C3 keeps the coefficients to a "reasonable" size, they are167* still potentially several binary orders of magnitude larger than the inputs.168* Thus, this algorithm should only be used where T is a multi-precision type.169*170* @tparam T Polynomial coefficient type.171* @param u First polynomial.172* @param v Second polynomial.173* @return Greatest common divisor of polynomials u and v.174*/175template <class T>176typename std::enable_if< std::numeric_limits<T>::is_integer, polynomial<T> >::type177subresultant_gcd(polynomial<T> u, polynomial<T> v)178{179 using std::swap;180 BOOST_MATH_ASSERT(u || v);181 182 if (!u)183 return v;184 if (!v)185 return u;186 187 typedef typename polynomial<T>::size_type N;188 189 if (u.degree() < v.degree())190 swap(u, v);191 192 T const d = detail::reduce_to_primitive(u, v);193 T g = 1, h = 1;194 polynomial<T> r;195 while (true)196 {197 BOOST_MATH_ASSERT(u.degree() >= v.degree());198 // Pseudo-division.199 r = u % v;200 if (!r)201 return d * primitive_part(v); // Attach the content.202 if (r.degree() == 0)203 return d * polynomial<T>(T(1)); // The content is the result.204 N const delta = u.degree() - v.degree();205 // Adjust remainder.206 u = v;207 v = r / (g * detail::integer_power(h, delta));208 g = leading_coefficient(u);209 T const tmp = detail::integer_power(g, delta);210 if (delta <= N(1))211 h = tmp * detail::integer_power(h, N(1) - delta);212 else213 h = tmp / detail::integer_power(h, delta - N(1));214 }215}216 217 218/**219 * @brief GCD for polynomials with unbounded multi-precision integral coefficients.220 *221 * The multi-precision constraint is enforced via numeric_limits.222 *223 * Note that intermediate terms in the evaluation can grow arbitrarily large, hence the need for224 * unbounded integers, otherwise numeric overflow would break the algorithm.225 *226 * @tparam T A multi-precision integral type.227 */228template <typename T>229typename std::enable_if<std::numeric_limits<T>::is_integer && !std::numeric_limits<T>::is_bounded, polynomial<T> >::type230gcd(polynomial<T> const &u, polynomial<T> const &v)231{232 return subresultant_gcd(u, v);233}234// GCD over bounded integers is not currently allowed:235template <typename T>236typename std::enable_if<std::numeric_limits<T>::is_integer && std::numeric_limits<T>::is_bounded, polynomial<T> >::type237gcd(polynomial<T> const &u, polynomial<T> const &v)238{239 static_assert(sizeof(v) == 0, "GCD on polynomials of bounded integers is disallowed due to the excessive growth in the size of intermediate terms.");240 return subresultant_gcd(u, v);241}242// GCD over polynomials of floats can go via the Euclid algorithm:243template <typename T>244typename 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> >::type245gcd(polynomial<T> const &u, polynomial<T> const &v)246{247 return boost::integer::gcd_detail::Euclid_gcd(u, v);248}249 250}251//252// Using declaration so we overload the default implementation in this namespace:253//254using boost::math::tools::gcd;255 256}257 258namespace integer259{260 //261 // Using declaration so we overload the default implementation in this namespace:262 //263 using boost::math::tools::gcd;264}265 266} // namespace boost::math::tools267 268#endif269