From 73e4618db3ad365635e729b8dca1ebb66b78bb03 Mon Sep 17 00:00:00 2001 From: Sebastian Limbach Date: Thu, 15 Mar 2007 10:07:42 +0000 Subject: [PATCH] Imported 'Descartes.h' from EXACUS and adapted to CGAL, first version. --- .gitattributes | 1 + .../CGAL/Algebraic_kernel_d/Descartes.h | 520 ++++++++++++++++++ 2 files changed, 521 insertions(+) create mode 100644 Algebraic_kernel_d/include/CGAL/Algebraic_kernel_d/Descartes.h diff --git a/.gitattributes b/.gitattributes index b070cedd045..03eb47e4ba9 100644 --- a/.gitattributes +++ b/.gitattributes @@ -26,6 +26,7 @@ Algebraic_kernel_GBRS/test/Gbrs_polynomial/Parser_polys.C -text Algebraic_kernel_GBRS/test/Gbrs_polynomial/parsers.h -text Algebraic_kernel_d/doc_tex/Algebraic_kernel_d_ref/figures/cpvl.eps -text svneol=unset#application/postscript Algebraic_kernel_d/doc_tex/Algebraic_kernel_d_ref/figures/cpvl.fig -text svneol=unset#application/octet-stream +Algebraic_kernel_d/include/CGAL/Algebraic_kernel_d/Descartes.h -text Alpha_shapes_2/demo/Alpha_shapes_2/alpha_shapes_2.vcproj eol=crlf Alpha_shapes_2/demo/Alpha_shapes_2/data/m30f.jpg -text svneol=unset#image/jpeg Alpha_shapes_2/demo/Alpha_shapes_2/help/index.html svneol=native#text/html diff --git a/Algebraic_kernel_d/include/CGAL/Algebraic_kernel_d/Descartes.h b/Algebraic_kernel_d/include/CGAL/Algebraic_kernel_d/Descartes.h new file mode 100644 index 00000000000..6d1b3f9463d --- /dev/null +++ b/Algebraic_kernel_d/include/CGAL/Algebraic_kernel_d/Descartes.h @@ -0,0 +1,520 @@ +// 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) : +// +// ============================================================================ + +// TODO: The comments are all original EXACUS comments and aren't adapted. So +// they may be wrong now. + +/*! \file NiX/Descartes.h + \brief Defines class NiX::Descartes. + + Isolate real roots of polynomials. + + This file provides a class to isolate real roots of polynomials, + using the algorithm based on the method of Descartes. + + The polynomial has to be a univariat polynomial over any number + type which is contained in the real numbers. + +*/ + +#ifndef CGAL_ALGEBRAIC_KERNEL_D_DESCARTES_H +#define CGAL_ALGEBRAIC_KERNEL_D_DESCARTES_H + +#include +#include + +#include + +#ifdef CGAL_USE_LEDA +#include +#include +#endif +#ifdef CGAL_USE_CORE +#include +#include +#endif + +CGAL_BEGIN_NAMESPACE + +namespace CGALi { + +// TODO: Copied from EXACUS, not part of any concept!? + +#ifdef CGAL_USE_LEDA +// Constructs 2^e from an integer e. Needed in Descartes +inline void construct_binary(const ::leda::integer& e, ::leda::integer& x) { + typedef ::leda::integer Integer; + x = Integer(1) << e.to_long(); +} + +// Constructs m*2^e from two integers m,e. Needed in Descartes +inline void construct_binary(const ::leda::integer& m, const ::leda::integer& e, + ::leda::rational& x) { + + // TODO: Default implementation should use the construct_binary function above + typedef ::leda::integer Integer; + typedef ::leda::rational Rational; + + Integer den(1); + Integer num(m); + if(e>0) { + num <<= e.to_long(); + } + else { + den <<= (-e).to_long(); + } + x = Rational(num, den); +} + +#endif // CGAL_USE_LEDA + +#ifdef CGAL_USE_CORE +// Constructs 2^e from an integer e. Needed in Descartes +inline void construct_binary(const ::CORE::BigInt& e, ::CORE::BigInt& x) { + typedef ::CORE::BigInt Integer; + x = Integer(1) << ::CORE::ulongValue(e); +} + +// Constructs m*2^e from two integers m,e. Needed in Descardes +inline void construct_binary(const ::CORE::BigInt& m, const ::CORE::BigInt& e, + ::CORE::BigRat& x) { + typedef ::CORE::BigInt Integer; + typedef ::CORE::BigRat Rational; + + Integer den(1); + Integer num(m); + if(e>0) { + num <<= ::CORE::ulongValue(e); + } + else { + den <<= ::CORE::ulongValue(-e); + } + x = Rational(num, den); +} + +#endif // CGAL_USE_CORE + + +/*! \ingroup NiX_Algebraic_real + * \brief A model of concept RealRootIsolator. + */ +template +class Descartes { + typedef CGAL::Fraction_traits FT_poly; + typedef Fraction_traits FT_rat; +public: + //! First template parameter + typedef Polynomial_ Polynomial; + //! Second template parameter + typedef Rational_ Rational; + //! Boundary type of the isolating intervals + typedef Rational_ Boundary; + // Integer or Numerator/Denominator type of boundary. + typedef typename CGAL::Fraction_traits::Numerator_type Integer; +private: + typedef typename Polynomial::NT Coeff; + typedef Integer IT; + + Polynomial poly_; + int number_of_real_roots_; + IT* numerator; + IT* denominator_exponent; + bool* is_exact; + IT LEFT,SCALE,DENOM; + bool is_strong_; + int k; + bool interval_given; + +public: + /*! \brief Constructor from univariate square free polynomial. + + The RealRootIsolator provides isolating intervals for the real + roots of the polynomial. + \pre the polynomial is square free + */ + Descartes(const Polynomial& P = Polynomial(Coeff(0)), + bool is_strong = false, + int kk = 2) + : poly_(P) , + is_strong_(is_strong), + k(kk), + interval_given(false) { + + numerator = new IT[P.degree()]; + denominator_exponent = new IT[P.degree()]; + is_exact = new bool[P.degree()]; + number_of_real_roots_ = 0; + if(P.degree() == 0) + { + if(P.is_zero()) number_of_real_roots_ = -1; + return; + } + + intern_decompose(poly_,typename FT_poly::Is_fraction()); + } + + // constructor for coefficient types \c Coeff with given interval + // (experimental) + Descartes(const Polynomial& P, + const Rational& left, + const Rational& right, + bool is_strong = false, + int kk = 2) + : poly_(P) , + is_strong_(is_strong), + k(kk), + interval_given(true) { + + numerator = new IT[P.degree()]; + denominator_exponent = new IT[P.degree()]; + is_exact = new bool[P.degree()]; + number_of_real_roots_ = 0; + if(P.degree() == 0) + { + if(P.is_zero()) number_of_real_roots_ = -1; + return; + } + typename FT_rat::Decompose decompose; + typedef typename FT_rat::Numerator Numerator; + typedef typename FT_rat::Denominator Denominator; + Numerator numleft, numright; + Denominator denleft, denright; + + decompose(left,numleft,denleft); + decompose(right,numright,denright); + + LEFT = numleft * denright; + SCALE = numright * denleft - LEFT; + DENOM = denleft * denright; + poly_.scale_down(denleft*denright); + + intern_decompose(poly_,typename FT_poly::Is_decomposable()); + } + + //! copy constructor + Descartes(const Descartes& D) + : poly_(D.poly_), + number_of_real_roots_(D.number_of_real_roots_), + LEFT(D.LEFT), + SCALE(D.SCALE), + DENOM(D.DENOM), + is_strong_(D.is_strong_), + k(D.k), + interval_given(D.interval_given) { + + numerator = new IT[poly_.degree()]; + denominator_exponent = new IT[poly_.degree()]; + is_exact = new bool[poly_.degree()]; + for(int i=0; i= 0 && i < number_of_real_roots_); + construct_binary(denominator_exponent[i], denominator_); + numerator_= SCALE * numerator[i] + LEFT * denominator_; + denominator_ = denominator_ * DENOM; + } + + + void right_boundary(int i,IT& numerator_, IT& denominator_) const { + CGAL_assertion(i >= 0 && i < number_of_real_roots_); + if(is_exact[i]){ + return left_boundary(i,numerator_,denominator_); + } + else{ + construct_binary(denominator_exponent[i],denominator_); + numerator_= SCALE * (numerator[i]+1) + LEFT * denominator_; + denominator_ = denominator_ * DENOM; + } + } +public: + + /*! \brief returns \f${l_i}\f$ the left boundary of the isolating interval + for root \f$root_{i}\f$. + + In case is_exact_root(i) is true, \f$l_i = root_{i}\f$,\n + otherwise: \f$l_i < root_{i}\f$. + + If \f$i-1>=0\f$, then \f$l_i > root_{i-1}\f$. \n + If \f$i-1>=0\f$, then \f$l_i >= r_{i-1}\f$, + the right boundary of \f$root_{i-1}\f$\n + + \pre 0 <= i < number_of_real_roots() + */ + Rational left_boundary(int i) const { + IT numerator_, denominator_; + left_boundary(i,numerator_,denominator_); + return Rational(numerator_) / Rational(denominator_); + } + + /*! \brief returns \f${r_i}\f$ the right boundary of the isolating interval + for root \f$root_{i}\f$. + + In case is_exact_root(i) is true, \f$r_i = root_{i}\f$,\n + otherwise: \f$r_i > root_{i}\f$. + + If \f$i+1< n \f$, then \f$r_i < root_{i+1}\f$, + where \f$n\f$ is number of real roots.\n + If \f$i+1< n \f$, then \f$r_i <= l_{i+1}\f$, + the left boundary of \f$root_{i+1}\f$\n + + + \pre 0 <= i < number_of_real_roots() + */ + Rational right_boundary(int i) const { + IT numerator_, denominator_; + right_boundary(i,numerator_,denominator_); + return Rational(numerator_) / Rational(denominator_); + } + +private: + void intern_decompose( Polynomial P_, ::CGAL::Tag_true){ + typename FT_poly::Decompose decompose; + typedef typename FT_poly::Numerator_type Numerator_poly; + typedef typename Numerator_poly::NT Coeff; + typename FT_poly::Numerator_type NumP; + typename FT_poly::Denominator_type dummy; + + decompose(P_,NumP,dummy); + init_with(NumP); + } + + void intern_decompose( Polynomial P, ::CGAL::Tag_false){ + init_with(P); + } + + + template + void init_with(const Polynomial__& P){ + typedef typename Polynomial__::NT Coeff; + if(!interval_given) + { + LEFT = -weak_upper_root_bound(P); + SCALE = - LEFT * IT(2); + DENOM = IT(1); + } + Polynomial__ R = ::CGAL::translate(P,Coeff(LEFT)); + Polynomial__ Q = ::CGAL::scale_up(R,Coeff(SCALE)); + zero_one_descartes(Q,0,0); + } + + + //! returns the polynomial $(1 + x)^n P(1/(1 + x))$. + template + CGAL::Polynomial + variation_transformation(const CGAL::Polynomial& P) { + CGAL::Polynomial R = reversal(P); + return translate_by_one(R); + } + + //! Returns an upper bound on the absolute value of all roots of $P$. + /*! The upper bound is a power of two. Only works for univariate + * polynomials. + */ + template + IT weak_upper_root_bound(const CGAL::Polynomial& P) { + + typename Real_embeddable_traits::Abs abs; + const int n = P.degree(); + IT r(1); // return value + Coeff__ x(1); // needed to "evaluate" the polynomial + Coeff__ val; + for (;;) { + val = -abs(P[n]); + for (int i = n-1; i >= 0; i--) { + val = val*x + abs(P[i]); + } + if (val < Coeff__(0)) return r; + r *= IT(2); + x = Coeff__(r); + } + } + + //! tests if the polynomial has no root in the interval. + template + bool not_zero_in_interval(const CGAL::Polynomial& P) + { + if(P.degree() == 0) return true; + if(INTERN_POLYNOMIAL::sign_variations(variation_transformation(P)) != 0) + return false; + return (P[0] != Coeff__(0) && P.evaluate(Coeff__(1)) != Coeff__(0)); + } + //! Descartes algoritm to determine isolating intervals for the roots + //! lying in the interval (0,1). + // The parameters $(i,D)$ describe the interval $(i/2^D, (i+1)/2^D)$. + // Here $0\leq i < 2^D$. + template + void zero_one_descartes(const CGAL::Polynomial& P, + IT i, IT D) { + // Determine the number of sign variations of the transformed + // polynomial $(1+x)^nP(1/(1+x))$. This gives the number of + // roots of $P$ in $(0,1)$. + + CGAL::Polynomial R = variation_transformation(P); + int descarte = sign_variations(R); + + // no root + if ( descarte == 0 ) return; + + // exactly one root + // Note the termination criterion $P(0)\neq 0$ and $P(1)\neq 0$. + // This ensures that the given interval is an isolating interval. + if ( descarte == 1 + && P[0] != Coeff__(0) + && P.evaluate(Coeff__(1)) != Coeff__(0) ) { + if(is_strong_) { + strong_zero_one_descartes(P,i,D); + return; + } + else { + numerator[number_of_real_roots_] = i; + denominator_exponent[number_of_real_roots_] = D; + is_exact[number_of_real_roots_] = false; + number_of_real_roots_++; + return; + } + } + + // more than one root + // Refine the interval. + i = 2*i; D = D+1; + + // Transform the polynomial such that the first half of the interval + // is mapped to the unit interval. + CGAL::Polynomial Q = scale_down(P,Coeff__(2)); + + // Consider the first half of the interval. + zero_one_descartes(Q,i,D); + + // Test if the polynomial is zero at the midpoint of the interval + CGAL::Polynomial S = translate_by_one(Q); + if ( S[0] == Coeff__(0) ) { + numerator[number_of_real_roots_] = i + 1; + denominator_exponent[number_of_real_roots_] = D; + is_exact[number_of_real_roots_] = true; + number_of_real_roots_++; + } + + // Consider the second half of the interval. + zero_one_descartes(S,i+1,D); + } + + + //! Strong Descartes algoritm to determine isolating intervals for the + //! roots lying in the interval (0,1), where the first + //! derivative have no sign change. \pre $P$ has only one root in the + //! interval given by $(i,D)$. + // The parameters $(i,D)$ describe the interval $(i/2^D, (i+1)/2^D)$. + // Here $0\leq i < D$. + template + void strong_zero_one_descartes(const CGAL::Polynomial& P, + IT i, IT D) { + + // Test if the polynomial P' has no roots in the + // interval. For further use in Newton, the interval should be not + // too large. + + // test if isolating interval is smaller than epsilon + // [l,r] -> r-l < epsilon + // l = (r-l) * i/2^D + l + // r = (r-l) * (i+1)/2^D + l + // r-l = (r-l) * 1/2^D + // r-l < epsilon = 2^(-k) + // <=> (r-l) * 1/2^D < 2^(-k) + // <=> 2^D > (r-l) / 2^(-k) + // <=> 2^D > (r-l) * 2^k + + CGAL::Polynomial PP = diff(P); + if(not_zero_in_interval(PP)) { // P' + IT tmp; + construct_binary(D-k, tmp); // tmp = 2^{D-k} + if(tmp * DENOM > SCALE ) { + numerator[number_of_real_roots_] = i; + denominator_exponent[number_of_real_roots_] = D; + is_exact[number_of_real_roots_] = false; + number_of_real_roots_++; + return; + } + } + + // either $P'$ fails the test, + // or the interval is too large + // Refine the interval. + i = 2*i; D = D+1; + + // Transform the polynomial such that the first half of the interval + // is mapped to the unit interval. + CGAL::Polynomial Q = scale_down(P,Coeff__(2)); + + // Test if the polynomial is zero at the midpoint of the interval + CGAL::Polynomial S = translate_by_one(Q); + if ( S[0] == Coeff__(0) ) { + numerator[number_of_real_roots_] = i + 1; + denominator_exponent[number_of_real_roots_] = D; + is_exact[number_of_real_roots_] = true; + number_of_real_roots_++; + return; + } + + // Consider the first half of the interval. + if(sign_variations(variation_transformation(Q)) == 1) { + strong_zero_one_descartes(Q,i,D); + return; + } + + // Consider the second half of the interval. + strong_zero_one_descartes(S,i+1,D); + return; + } +}; + +} // namespace CGALi + +CGAL_END_NAMESPACE + +#endif // CGAL_ALGEBRAIC_KERNEL_D_DESCARTES_H