1 /** @file heur_gcd_uvar.h
3 * Heuristic GCD code. */
6 * GiNaC Copyright (C) 1999-2014 Johannes Gutenberg University Mainz, Germany
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 #ifndef GINAC_UPOLY_HEUR_GCD_H
24 #define GINAC_UPOLY_HEUR_GCD_H
27 #include "ring_traits.h"
28 #include "normalize.h"
29 #include "remainder.h"
30 #include "eval_uvar.h"
31 #include "interpolate_padic_uvar.h"
37 /// Compute GCD of primitive univariate polynomials.
38 template<typename T> static bool
39 heur_gcd_z_pp(T& g, const T& a, const T& b, unsigned max_tries = 66)
41 typedef typename T::value_type ring_t;
42 const ring_t n73794 = get_ring_elt(b[0], 73794);
43 const ring_t n27011 = get_ring_elt(b[0], 27011);
44 const std::size_t maxdeg = std::max(degree(a), degree(b));
49 // find the evaluation point
50 ring_t xi = (std::min(max_coeff(a), max_coeff(b)) + 1) << 1;
53 const ring_t av = eval(a, xi);
54 const ring_t bv = eval(b, xi);
55 const ring_t gamma = gcd(av, bv);
56 interpolate(gg, gamma, xi, maxdeg);
57 normalize_in_ring(gg);
58 remainder_in_ring(r, a, gg);
63 // next evaluation point
64 xi = truncate1(xi*isqrt(isqrt(xi))*n73794, n27011);
65 } while (--max_tries != 0);
69 template<typename T> static bool
70 heur_gcd_z_priv(T& g, const T& a, const T& b, const unsigned max_tries = 66)
72 typedef typename T::value_type ring_t;
75 normalize_in_ring(a_, &acont);
76 normalize_in_ring(b_, &bcont);
77 const ring_t gc = gcd(acont, bcont);
78 bool found = heur_gcd_z_pp(g, a_, b_, max_tries);
88 #endif // ndef GINAC_UPOLY_HEUR_GCD_H