// TODO: Add licence // // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // // $URL$ // $Id$ // // // Author(s) : Arno Eigenwillig // // ============================================================================ // TODO: The comments are all original EXACUS comments and aren't adapted. So // they may be wrong now. /*! \file NiX/prs_resultant.h * \brief Resultant computation via polynomial remainder sequences (PRS) * */ #include #include #ifndef CGAL_PRS_RESULTANT_H #define CGAL_PRS_RESULTANT_H CGAL_BEGIN_NAMESPACE template inline NT prs_resultant_ufd(Polynomial A, Polynomial B) { // implemented using the subresultant algorithm for resultant computation // see [Cohen, 1993], algorithm 3.3.7 if (A.is_zero() || B.is_zero()) return NT(0); int signflip; if (A.degree() < B.degree()) { Polynomial T = A; A = B; B = T; signflip = (A.degree() & B.degree() & 1); } else { signflip = 0; } NT a = A.content(), b = B.content(); NT g(1), h(1), t = ipower(a, B.degree()) * ipower(b, A.degree()); Polynomial Q, R; NT d; int delta; A /= a; B /= b; do { signflip ^= (A.degree() & B.degree() & 1); Polynomial::pseudo_division(A, B, Q, R, d); delta = A.degree() - B.degree(); CGAL_expensive_assertion( typename CGAL::Algebraic_structure_traits::Is_exact == CGAL_Tag_false || d == ipower(B.lcoeff(), delta + 1) ); A = B; B = R / (g * ipower(h, delta)); g = A.lcoeff(); // h = h^(1-delta) * g^delta INTERN_POLYNOMIAL::hgdelta_update(h, g, delta); } while (B.degree() > 0); // h = h^(1-deg(A)) * lcoeff(B)^deg(A) delta = A.degree(); g = B.lcoeff(); INTERN_POLYNOMIAL::hgdelta_update(h, g, delta); h = signflip ? -(t*h) : t*h; typename Algebraic_structure_traits::Simplify simplify; simplify(h); return h; } template inline NT prs_resultant_field(Polynomial A, Polynomial B) { // implemented using the Euclidean algorithm for resultant computation // compare [Cox et al, 1997], p.157 if (A.is_zero() || B.is_zero()) return NT(0); int signflip; if (A.degree() < B.degree()) { Polynomial T = A; A = B; B = T; signflip = (A.degree() & B.degree() & 1); } else { signflip = 0; } NT res(1); Polynomial Q, R; while (B.degree() > 0) { signflip ^= (A.degree() & B.degree() & 1); Polynomial::euclidean_division(A, B, Q, R); res *= ipower(B.lcoeff(), A.degree() - R.degree()); A = B; B = R; } res = ipower(B.lcoeff(), A.degree()) * (signflip ? -res : res); typename Algebraic_structure_traits::Simplify simplify; simplify(res); return res; } // definition follows below template inline NT prs_resultant_decompose(Polynomial A, Polynomial B); namespace INTERN_PRS_RESULTANT { template inline NT prs_resultant_(Polynomial A, Polynomial B, ::CGAL::Tag_false) { return prs_resultant_field(A, B); } template inline NT prs_resultant_(Polynomial A, Polynomial B, ::CGAL::Tag_true) { return prs_resultant_decompose(A, B); } template inline NT prs_resultant_(Polynomial A, Polynomial B, Field_tag) { typedef typename Fraction_traits::Is_fraction Is_decomposable; return prs_resultant_(A, B, Is_decomposable()); } template inline NT prs_resultant_(Polynomial A, Polynomial B, Unique_factorization_domain_tag) { return prs_resultant_ufd(A, B); } } // namespace Intern template inline NT prs_resultant_decompose(Polynomial A, Polynomial B){ typedef Polynomial POLY; typedef typename Fraction_traits::Numerator_type INTPOLY; typedef typename Fraction_traits::Denominator_type DENOM; typedef typename INTPOLY::NT RES; DENOM a, b; A.simplify_coefficients(); B.simplify_coefficients(); INTPOLY A0 = integralize_polynomial(A, a); INTPOLY B0 = integralize_polynomial(B, b); DENOM c = ipower(a, B.degree()) * ipower(b, A.degree()); typedef typename Algebraic_structure_traits::Algebraic_category Algebraic_category; RES res0 = INTERN_PRS_RESULTANT::prs_resultant_(A0, B0, Algebraic_category()); typename Fraction_traits::Compose comp_frac; NT res = comp_frac(res0, c); typename Algebraic_structure_traits::Simplify simplify; simplify(res); return res; } /*! \ingroup NiX_Polynomial * \relates NiX::Polynomial * \brief compute the resultant of polynomials \c A and \c B * * The resultant of two polynomials is computed from their * polynomial remainder sequence (PRS), in the Euclidean or * subresultant version. This depends on the coefficient type: * If \c NT is a \c UFDomain , the subresultant PRS is formed. * If \c NT is a \c Field that is not decomposable (see * \c NiX::Fraction_traits ), then a Euclidean PRS is formed. * If \c NT is a \c Field that is decomposable, then the * \c Numerator must be a \c UFDomain, and the subresultant * PRS is formed for the decomposed polynomials. * * Using \c NiX::hybrid_bezout_subresultant() may be faster in some cases * and works for non-UFDomains, too. * Using \c NiX::resultant() from \c NiX/resultant.h * chooses automatically among these alternative methods of resultant * computation for you. * * For the benefit of those who want to do their own template * metaprogramming to choose the method of resultant computation, * the three variants of resultant computation from a PRS * can be called directly as \c prs_resultant_field() , * \c prs_resultant_ufd() and \c prs_resultant_decompose() . * Do not use them directly unless you know what you are doing! * */ template inline NT prs_resultant(Polynomial A, Polynomial B) { typedef typename Algebraic_structure_traits::Algebraic_category Algebraic_category; return INTERN_PRS_RESULTANT::prs_resultant_(A, B, Algebraic_category()); } CGAL_END_NAMESPACE #endif // CGAL_PRS_RESULTANT_H // EOF