// Copyright (c) 2006-2009 Max-Planck-Institute Saarbruecken (Germany). // All rights reserved. // // This file is part of CGAL (www.cgal.org) // // $URL: https://github.com/CGAL/cgal/blob/v5.2/Algebraic_kernel_d/include/CGAL/Algebraic_kernel_d/Real_roots.h $ // $Id: Real_roots.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial // // // Author(s) : Michael Hemmer // // ============================================================================ // TODO: The comments are all original EXACUS comments and aren't adapted. So // they may be wrong now. /*! \file NiX/Real_roots.h \brief This file defines the class NiX::Real_roots. */ #ifndef CGAL_ALGEBRAIC_KERNEL_D_REAL_ROOTS_H #define CGAL_ALGEBRAIC_KERNEL_D_REAL_ROOTS_H #include #include #include #include #include namespace CGAL { namespace internal { /*! \ingroup NiX_Real_roots * \brief This class provides operators for a comfortable construction of * AlgebraicReal from Polynomials using a specific RealRootIsolator. * * A valid template argument for AlgebraicReal is NiX::Algebraic_real. */ template < class AlgebraicReal , class RealRootIsolator > class Real_roots{ public: //! The Real_roots type it self. typedef Real_roots Self; //! First template argument. typedef AlgebraicReal Algebraic_real; //! Second template argument. typedef RealRootIsolator Real_root_isolator; //! The Polnomial type used by Algebraic_real and the Real_root_isolator. typedef typename Algebraic_real::Polynomial_1 Polynomial; private: typedef typename AlgebraicReal::Coefficient Coefficient; typedef typename AlgebraicReal::Rational Rational; private: template < class PolynomialIterator, class IntIterator, class AlgebraicRealOutputIterator, class IntOutputIterator> int gen_agebraic_reals_with_mults( PolynomialIterator fac, PolynomialIterator fac_end, IntIterator mul, IntIterator CGAL_precondition_code(mul_end), AlgebraicRealOutputIterator oi_root, IntOutputIterator oi_mult){ Self real_roots; typedef std::pair PAIR; // find zeroes of each factor and sort them in ascending order std::priority_queue< PAIR,std::vector,std::greater > pqueue; std::vector tmp; while(fac != fac_end){ CGAL_assertion(mul != mul_end); tmp.clear(); real_roots(*fac, std::back_inserter(tmp)); for (int j = 0; j < static_cast(tmp.size()); j++) { pqueue.push(PAIR(tmp[j], *mul)); } fac++; mul++; } // output factors and multiplicities int n = 0; while (!pqueue.empty()) { *oi_root++ = pqueue.top().first; *oi_mult++ = pqueue.top().second; n++; pqueue.pop(); } return n; } private: template < class PolynomialConstIterator, class IntConstIterator, class PolynomialOutputIterator> void write_factors_by_multiplicity(PolynomialConstIterator fac, PolynomialConstIterator fac_end, IntConstIterator mul, IntConstIterator CGAL_assertion_code(mul_end), PolynomialOutputIterator oi_poly){ // output table such that table[m] contains square-free factor of // multiplicity m (or else constant poly 1) int m = 0; while (fac != fac_end) { CGAL_assertion(mul != mul_end); while (m < *mul) { *oi_poly++ = Polynomial(Coefficient(1)); m++; } *oi_poly++ = *fac; m++; ++fac; ++mul; } } public: /*! \brief computes all roots of the square free polynomial P in * ascending order and returns the number of real roots. */ template int operator()(const Polynomial& poly , AlgebraicRealOutputIterator it){ CGAL_precondition_msg( typename CGAL::Polynomial_traits_d< Polynomial >::Is_square_free()(poly), "P not square free."); return (*this)(Real_root_isolator(poly),it); } public: /*! \brief computes all roots of the polynomial P in ascending order * and their multiplicity and returns the number of real roots. * * This operator returns the number \e n of distinct real roots of \c poly. * Each root is represented as an object of an instance AlgebraicReal. * The operator writes these \e n real zeroes in ascending order to \c * oi_root. * It writes the multiplicities of the zeroes in the same order to * \c oi_mult . */ template int operator()(const Polynomial& poly, AlgebraicRealOutputIterator oi_root, IntOutputIterator oi_mult){ CGAL_precondition(CGAL::degree(poly) >= 0); // fast exit if (CGAL::degree(poly) == 0) return (poly.is_zero())?-1:0; std::list sqffac; std::list facmul; filtered_square_free_factorize_utcf(poly, std::back_inserter(sqffac), std::back_inserter(facmul)); int number_of_real_roots= gen_agebraic_reals_with_mults(sqffac.begin(),sqffac.end(), facmul.begin(),facmul.end(), oi_root, oi_mult); return number_of_real_roots; } public: /*! \brief computes all roots defined by the Real_root_isolator object in * ascending order */ template int operator()(const Real_root_isolator& isolator, AlgebraicRealOutputIterator it){ Polynomial poly = isolator.polynomial(); //cout << "P: "< FT; typename FT::Numerator_type num; typename FT::Denominator_type denom; typename FT::Decompose decomp; decomp(root,num,denom); Polynomial linear_factor(Coefficient(-num), Coefficient(denom)); poly=CGAL::integral_division(poly,linear_factor); } } //cout << "P_without_exact: "< conjugated_roots; std::back_insert_iterator > con_it = std::back_inserter(conjugated_roots); // construct AlgebraicReal for(int j = 0 ; j < isolator.number_of_real_roots(); j++){ if(isolator.is_exact_root(j)){ // exact roots (Rational Rational root=isolator.left_bound(j); CGAL::simplify(root); *it++=AlgebraicReal(root); }else{ // other roots Rational left = isolator.left_bound(j); Rational right= isolator.right_bound(j); CGAL::simplify(left); CGAL::simplify(right); AlgebraicReal tmp(poly,left,right); *it++=tmp; *con_it++=tmp; } } AlgebraicReal::conjugate(conjugated_roots.begin(), conjugated_roots.end()); return isolator.number_of_real_roots(); } /*! \brief factor \c p by multiplicities, return both factors and their roots * * This operator is for those users who have an * Polynomial \c poly which is not necessarily square-free, and who * want to get both its square-free factorization and its real roots with * their respective multiplicities. * * This operator returns the number \e n of distinct real roots of \c poly . * Each root is represented as an object of an instance AlgebraicReal. * The operator writes these \e n real zeroes in ascending order to \c * oi_root. * It writes the multiplicities of the roots in the same order to * \c oi_mult . * Finally, it writes the square-free factors of \c p * to \c oi_poly such that the factor #k written is the factor * with exponent \e k in the square-free factorization of \c p , or * the constant polynomial 1 if a factor of multiplicity \e k does * not occur. Yes, this means that the first factor written, which is * factor #0, will always be 1. * * The data types involved are determined the \c AlgebraicReal * template argument. In particular, \c p must be of type * \c NiX::Polynomial, as are its factors; * the roots are of type \c AlgebraicReal ; and the multiplicities are * of type \c int . * */ template int operator()( const Polynomial& poly, AlgebraicRealOutputIterator oi_root, IntOutputIterator oi_mult, PolynomialOutputIterator oi_poly) { CGAL_precondition(CGAL::degree(poly) >= 0); // fast exit if (CGAL::degree(poly) == 0) return (poly.is_zero())?-1:0; std::list sqffac; std::list facmul; filtered_square_free_factorize_utcf(poly, std::back_inserter(sqffac), std::back_inserter(facmul)); write_factors_by_multiplicity(sqffac.begin(),sqffac.end(), facmul.begin(),facmul.end(), oi_poly); int numer_of_real_roots = gen_agebraic_reals_with_mults(sqffac.begin(),sqffac.end(), facmul.begin(),facmul.end(), oi_root, oi_mult); return numer_of_real_roots; } }; } // namespace internal } //namespace CGAL #endif // CGAL_ALGEBRAIC_KERNEL_D_REAL_ROOTS_ROOTS_H