mirror of https://github.com/CGAL/cgal
Imported 'Descartes.h' from EXACUS and adapted to CGAL, first version.
This commit is contained in:
Sebastian Limbach
parent
cfdb45f40b
commit
73e4618db3
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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 <CGAL/basic.h>
|
||||
#include <CGAL/Polynomial.h>
|
||||
|
||||
#include <CGAL/Polynomial/univariate_polynomial_utils.h>
|
||||
|
||||
#ifdef CGAL_USE_LEDA
|
||||
#include <CGAL/leda_integer.h>
|
||||
#include <CGAL/leda_rational.h>
|
||||
#endif
|
||||
#ifdef CGAL_USE_CORE
|
||||
#include <CGAL/CORE_BigInt.h>
|
||||
#include <CGAL/CORE_BigRat.h>
|
||||
#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 Polynomial_, class Rational_>
|
||||
class Descartes {
|
||||
typedef CGAL::Fraction_traits<Polynomial_> FT_poly;
|
||||
typedef Fraction_traits<Rational_> 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<Rational>::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<number_of_real_roots(); i++)
|
||||
{
|
||||
numerator[i] = D.numerator[i];
|
||||
denominator_exponent[i] = D.denominator_exponent[i];
|
||||
is_exact[i] = D.is_exact[i];
|
||||
}
|
||||
}
|
||||
|
||||
// destructor
|
||||
~Descartes() {
|
||||
delete[] numerator;
|
||||
delete[] denominator_exponent;
|
||||
delete[] is_exact;
|
||||
}
|
||||
|
||||
public: // functions
|
||||
|
||||
/*! \brief returns the defining polynomial*/
|
||||
Polynomial polynomial() const { return poly_; }
|
||||
|
||||
//! returns the number of real roots
|
||||
int number_of_real_roots() const { return number_of_real_roots_; }
|
||||
|
||||
/*! \brief returns true if the isolating interval is degenerated to a
|
||||
single point.
|
||||
|
||||
If is_exact_root(i) is true,
|
||||
then left_boundary(int i) equals \f$root_i\f$. \n
|
||||
If is_exact_root(i) is true,
|
||||
then right_boundary(int i) equals \f$root_i\f$. \n
|
||||
*/
|
||||
bool is_exact_root(int i) const { return is_exact[i]; }
|
||||
|
||||
|
||||
|
||||
public:
|
||||
|
||||
|
||||
void left_boundary(int i, IT& numerator_, IT& denominator_) const {
|
||||
CGAL_assertion(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<class Polynomial__>
|
||||
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 <class Coeff__>
|
||||
CGAL::Polynomial<Coeff__>
|
||||
variation_transformation(const CGAL::Polynomial<Coeff__>& P) {
|
||||
CGAL::Polynomial<Coeff__> 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 <class Coeff__>
|
||||
IT weak_upper_root_bound(const CGAL::Polynomial<Coeff__>& P) {
|
||||
|
||||
typename Real_embeddable_traits<Coeff__>::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 <class Coeff__>
|
||||
bool not_zero_in_interval(const CGAL::Polynomial<Coeff__>& 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 <class Coeff__>
|
||||
void zero_one_descartes(const CGAL::Polynomial<Coeff__>& 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<Coeff__> 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<Coeff__> 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<Coeff__> 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 <class Coeff__>
|
||||
void strong_zero_one_descartes(const CGAL::Polynomial<Coeff__>& 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<Coeff__> 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<Coeff__> Q = scale_down(P,Coeff__(2));
|
||||
|
||||
// Test if the polynomial is zero at the midpoint of the interval
|
||||
CGAL::Polynomial<Coeff__> 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
|
||||
Loading…
Reference in New Issue