ceph/src/boost/boost/math/tools/polynomial_gcd.hpp

   1 //  (C) Copyright Jeremy William Murphy 2016.
   2
   3 //  Use, modification and distribution are subject to the
   4 //  Boost Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP
   8 #define BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP
   9
  10 #ifdef _MSC_VER
  11 #pragma once
  12 #endif
  13
  14 #include <boost/math/tools/polynomial.hpp>
  15 #include <boost/math/common_factor_rt.hpp>
  16 #include <boost/type_traits/is_pod.hpp>
  17
  18
  19 namespace boost{
  20
  21    namespace integer {
  22
  23       namespace gcd_detail {
  24
  25          template <class T>
  26          struct gcd_traits;
  27
  28          template <class T>
  29          struct gcd_traits<boost::math::tools::polynomial<T> >
  30          {
  31             inline static const boost::math::tools::polynomial<T>& abs(const boost::math::tools::polynomial<T>& val) { return val; }
  32
  33             static const method_type method = method_euclid;
  34          };
  35
  36       }
  37 }
  38
  39
  40
  41 namespace math{ namespace tools{
  42
  43 /* From Knuth, 4.6.1:
  44 *
  45 * We may write any nonzero polynomial u(x) from R[x] where R is a UFD as
  46 *
  47 *      u(x) = cont(u) . pp(u(x))
  48 *
  49 * where cont(u), the content of u, is an element of S, and pp(u(x)), the primitive
  50 * part of u(x), is a primitive polynomial over S.
  51 * When u(x) = 0, it is convenient to define cont(u) = pp(u(x)) = O.
  52 */
  53
  54 template <class T>
  55 T content(polynomial<T> const &x)
  56 {
  57     return x ? gcd_range(x.data().begin(), x.data().end()).first : T(0);
  58 }
  59
  60 // Knuth, 4.6.1
  61 template <class T>
  62 polynomial<T> primitive_part(polynomial<T> const &x, T const &cont)
  63 {
  64     return x ? x / cont : polynomial<T>();
  65 }
  66
  67
  68 template <class T>
  69 polynomial<T> primitive_part(polynomial<T> const &x)
  70 {
  71     return primitive_part(x, content(x));
  72 }
  73
  74
  75 // Trivial but useful convenience function referred to simply as l() in Knuth.
  76 template <class T>
  77 T leading_coefficient(polynomial<T> const &x)
  78 {
  79     return x ? x.data().back() : T(0);
  80 }
  81
  82
  83 namespace detail
  84 {
  85     /* Reduce u and v to their primitive parts and return the gcd of their
  86     * contents. Used in a couple of gcd algorithms.
  87     */
  88     template <class T>
  89     T reduce_to_primitive(polynomial<T> &u, polynomial<T> &v)
  90     {
  91         using boost::math::gcd;
  92         T const u_cont = content(u), v_cont = content(v);
  93         u /= u_cont;
  94         v /= v_cont;
  95         return gcd(u_cont, v_cont);
  96     }
  97 }
  98
  99
 100 /**
 101 * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
 102 * Algorithm 4.6.1C: Greatest common divisor over a unique factorization domain.
 103 *
 104 * The subresultant algorithm by George E. Collins [JACM 14 (1967), 128-142],
 105 * later improved by W. S. Brown and J. F. Traub [JACM 18 (1971), 505-514].
 106 *
 107 * Although step C3 keeps the coefficients to a "reasonable" size, they are
 108 * still potentially several binary orders of magnitude larger than the inputs.
 109 * Thus, this algorithm should only be used where T is a multi-precision type.
 110 *
 111 * @tparam   T   Polynomial coefficient type.
 112 * @param    u   First polynomial.
 113 * @param    v   Second polynomial.
 114 * @return       Greatest common divisor of polynomials u and v.
 115 */
 116 template <class T>
 117 typename enable_if_c< std::numeric_limits<T>::is_integer, polynomial<T> >::type
 118 subresultant_gcd(polynomial<T> u, polynomial<T> v)
 119 {
 120     using std::swap;
 121     BOOST_ASSERT(u || v);
 122
 123     if (!u)
 124         return v;
 125     if (!v)
 126         return u;
 127
 128     typedef typename polynomial<T>::size_type N;
 129
 130     if (u.degree() < v.degree())
 131         swap(u, v);
 132
 133     T const d = detail::reduce_to_primitive(u, v);
 134     T g = 1, h = 1;
 135     polynomial<T> r;
 136     while (true)
 137     {
 138         BOOST_ASSERT(u.degree() >= v.degree());
 139         // Pseudo-division.
 140         r = u % v;
 141         if (!r)
 142             return d * primitive_part(v); // Attach the content.
 143         if (r.degree() == 0)
 144             return d * polynomial<T>(T(1)); // The content is the result.
 145         N const delta = u.degree() - v.degree();
 146         // Adjust remainder.
 147         u = v;
 148         v = r / (g * detail::integer_power(h, delta));
 149         g = leading_coefficient(u);
 150         T const tmp = detail::integer_power(g, delta);
 151         if (delta <= N(1))
 152             h = tmp * detail::integer_power(h, N(1) - delta);
 153         else
 154             h = tmp / detail::integer_power(h, delta - N(1));
 155     }
 156 }
 157
 158
 159 /**
 160  * @brief GCD for polynomials with unbounded multi-precision integral coefficients.
 161  *
 162  * The multi-precision constraint is enforced via numeric_limits.
 163  *
 164  * Note that intermediate terms in the evaluation can grow arbitrarily large, hence the need for
 165  * unbounded integers, otherwise numeric loverflow would break the algorithm.
 166  *
 167  * @tparam  T   A multi-precision integral type.
 168  */
 169 template <typename T>
 170 typename enable_if_c<std::numeric_limits<T>::is_integer && !std::numeric_limits<T>::is_bounded, polynomial<T> >::type
 171 gcd(polynomial<T> const &u, polynomial<T> const &v)
 172 {
 173     return subresultant_gcd(u, v);
 174 }
 175 // GCD over bounded integers is not currently allowed:
 176 template <typename T>
 177 typename enable_if_c<std::numeric_limits<T>::is_integer && std::numeric_limits<T>::is_bounded, polynomial<T> >::type
 178 gcd(polynomial<T> const &u, polynomial<T> const &v)
 179 {
 180    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.");
 181    return subresultant_gcd(u, v);
 182 }
 183 // GCD over polynomials of floats can go via the Euclid algorithm:
 184 template <typename T>
 185 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
 186 gcd(polynomial<T> const &u, polynomial<T> const &v)
 187 {
 188    return boost::integer::gcd_detail::Euclid_gcd(u, v);
 189 }
 190
 191 }
 192 //
 193 // Using declaration so we overload the default implementation in this namespace:
 194 //
 195 using boost::math::tools::gcd;
 196
 197 }
 198
 199 namespace integer
 200 {
 201    //
 202    // Using declaration so we overload the default implementation in this namespace:
 203    //
 204    using boost::math::tools::gcd;
 205 }
 206
 207 } // namespace boost::math::tools
 208
 209 #endif