mirror of https://github.com/CGAL/cgal
1267 lines
43 KiB
C++
1267 lines
43 KiB
C++
// Copyright (c) 2005 Stanford University (USA).
|
|
// All rights reserved.
|
|
//
|
|
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
|
|
// modify it under the terms of the GNU Lesser General Public License as
|
|
// published by the Free Software Foundation; version 2.1 of the License.
|
|
// See the file LICENSE.LGPL distributed with CGAL.
|
|
//
|
|
// Licensees holding a valid commercial license may use this file in
|
|
// accordance with the commercial license agreement provided with the software.
|
|
//
|
|
// 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) : Daniel Russel <drussel@alumni.princeton.edu>
|
|
|
|
#ifndef CGAL_STURM_LAZY_SOLVER_H
|
|
#define CGAL_STURM_LAZY_SOLVER_H
|
|
|
|
#include <CGAL/Polynomial/basic.h>
|
|
//#include <CGAL/Polynomial/utilities.h>
|
|
//#include <CGAL/Polynomial/make_simple.h>
|
|
|
|
// the Standard sequence
|
|
#include <CGAL/Polynomial/internal/Rational/Standard_sequence.h>
|
|
|
|
// representations for roots and intervals
|
|
#include <CGAL/Polynomial/internal/Sturm_root_rep.h>
|
|
#include <CGAL/Polynomial/internal/Sturm_isolating_interval.h>
|
|
|
|
// the root counters
|
|
//#include <CGAL/Polynomial/internal/Rational/Sturm_root_counter.h>
|
|
//#include <CGAL/Polynomial/internal/Rational/Descartes_root_counter.h>
|
|
//#include <CGAL/Polynomial/internal/Rational/Bezier_root_counter.h>
|
|
|
|
#include <list>
|
|
|
|
#include <CGAL/Polynomial/Kernel.h>
|
|
|
|
CGAL_POLYNOMIAL_BEGIN_NAMESPACE
|
|
|
|
//================================================================
|
|
//================================================================
|
|
// ************** the lazy Sturm solver **************************
|
|
//================================================================
|
|
//================================================================
|
|
|
|
#if 0
|
|
namespace Polynomiali
|
|
{
|
|
|
|
template<class Tag>
|
|
struct Which_root_counter_chooser
|
|
{
|
|
template<class P, class M>
|
|
struct Which
|
|
{
|
|
typedef
|
|
CGAL::Polynomial::internal::Sturm_root_counter<P,M>
|
|
Root_count;
|
|
};
|
|
};
|
|
|
|
template<>
|
|
struct Which_root_counter_chooser<Descartes_tag>
|
|
{
|
|
template<class P, class M>
|
|
struct Which
|
|
{
|
|
typedef
|
|
CGAL::Polynomial::internal::Descartes_root_counter<P>
|
|
Root_count;
|
|
};
|
|
};
|
|
|
|
template<>
|
|
struct Which_root_counter_chooser<Bezier_tag>
|
|
{
|
|
template<class P, class M>
|
|
struct Which
|
|
{
|
|
typedef
|
|
CGAL::Polynomial::internal::Bezier_root_counter<P>
|
|
Root_count;
|
|
};
|
|
};
|
|
|
|
template<class STag, class P, class M>
|
|
struct Which_root_counter
|
|
{
|
|
typedef Which_root_counter_chooser<STag> Which;
|
|
typedef typename Which::template Which<P,M>::Root_count
|
|
Root_count;
|
|
};
|
|
|
|
} // namespace Polynomiali
|
|
#endif
|
|
|
|
struct Sturm_tag {};
|
|
|
|
template<class T, class S = Sturm_tag, class M = Field_tag,
|
|
bool Drop_even_roots = false>
|
|
class Sturm_root_stack
|
|
{
|
|
public:
|
|
typedef T Traits;
|
|
typedef typename Traits::Function Polynomial;
|
|
typedef typename Traits::Isolating_interval Interval;
|
|
typedef typename Traits::NT NT;
|
|
typedef typename Traits::Sign_at Sign_at;
|
|
|
|
typedef M Method_tag;
|
|
typedef S Subdivision_tag;
|
|
|
|
enum { DROP_EVEN_ROOTS = Drop_even_roots };
|
|
|
|
typedef typename Traits::Standard_sequence Standard_sequence;
|
|
|
|
protected:
|
|
// typedef CGAL::Sign Sign;
|
|
// typedef internal::Sign_at<Polynomial> Sign_at;
|
|
typedef std::pair<unsigned int, unsigned int> uint_pair;
|
|
|
|
// typedef internal::Sturm_isolating_interval<NT> Interval;
|
|
typedef std::list<Interval> Interval_container;
|
|
typedef std::pair<Interval,Interval> Interval_pair;
|
|
typedef typename Traits::Root_count Root_count;
|
|
// typedef typename
|
|
// Polynomiali::Which_root_counter<S,P,M>::Root_count
|
|
// Root_count;
|
|
|
|
typedef Sturm_root_stack<T,S,M,DROP_EVEN_ROOTS> Self;
|
|
|
|
public:
|
|
typedef internal::Sturm_root_rep<Self,Interval> Root;
|
|
|
|
protected:
|
|
//====================
|
|
// PROTECTED METHODS
|
|
//====================
|
|
|
|
#if 0
|
|
// ACCESS METHODS TO THE INTERVAL LIST AND SIGN VARIATIONS LIST
|
|
void pop_front( Descartes_tag ) const
|
|
{
|
|
i_list.pop_front();
|
|
nroots.pop_front();
|
|
}
|
|
#endif
|
|
void pop_front( Sturm_tag ) const
|
|
{
|
|
i_list.pop_front();
|
|
s_variations.pop_front();
|
|
}
|
|
|
|
unsigned int num_roots_in_first_interval( Sturm_tag ) const
|
|
{
|
|
uint_pair ip = s_variations.front();
|
|
return ip.first - ip.second;
|
|
}
|
|
#if 0
|
|
unsigned int num_roots_in_first_interval( Descartes_tag ) const
|
|
{
|
|
return nroots.front();
|
|
}
|
|
#endif
|
|
|
|
void push_front(const Interval& ivl, unsigned int sva,
|
|
unsigned int svb, Sturm_tag) const
|
|
{
|
|
i_list.push_front(ivl);
|
|
s_variations.push_front( uint_pair(sva, svb) );
|
|
}
|
|
#if 0
|
|
void push_front(const Interval& ivl,
|
|
unsigned int nr, Descartes_tag) const
|
|
{
|
|
i_list.push_front(ivl);
|
|
nroots.push_front(nr);
|
|
}
|
|
#endif
|
|
void push_front(const Interval& ivl, Sturm_tag) const
|
|
{
|
|
i_list.push_front(ivl);
|
|
s_variations.push_front( uint_pair(1, 0) );
|
|
}
|
|
#if 0
|
|
void push_front(const Interval& ivl, Descartes_tag) const
|
|
{
|
|
i_list.push_front(ivl);
|
|
nroots.push_front(1);
|
|
}
|
|
#endif
|
|
|
|
//------------------------------------------------------------------
|
|
|
|
//------------------------------------------------------------------
|
|
// the method for subdivision; it guarantees that the first interval
|
|
// contains only one root
|
|
template<class SA, class RC>
|
|
void subdivide(Sturm_tag tag, const SA& sign_at_p,
|
|
const RC& rc) const
|
|
{
|
|
if ( i_list.size() == 0 ) { return; }
|
|
|
|
if ( num_roots_in_first_interval(tag) == 1 ) { return; }
|
|
|
|
while ( num_roots_in_first_interval(tag) > 1 ) {
|
|
Interval ivl = i_list.front();
|
|
|
|
#if 1
|
|
Interval_pair ivl_pair = ivl.split();
|
|
Interval first_half = ivl_pair.first;
|
|
Interval second_half = ivl_pair.second;
|
|
|
|
unsigned int sv_a = s_variations.front().first;
|
|
unsigned int sv_b = s_variations.front().second;
|
|
|
|
Sign s_mid =
|
|
first_half.apply_to_endpoint(sign_at_p, Interval::UPPER);
|
|
|
|
if ( s_mid == CGAL::ZERO ) {
|
|
Interval mid_ivl = ivl;
|
|
|
|
Interval mid = first_half.endpoint_interval(Interval::UPPER);
|
|
|
|
Sign s_midp, s_midm;
|
|
unsigned int sv_midm(0), sv_midp(0);
|
|
|
|
do {
|
|
mid_ivl = mid_ivl.middle_half();
|
|
s_midm = mid_ivl.apply_to_endpoint(sign_at_p, Interval::LOWER);
|
|
s_midp = mid_ivl.apply_to_endpoint(sign_at_p, Interval::UPPER);
|
|
|
|
if ( s_midm != CGAL::ZERO && s_midp != CGAL::ZERO ) {
|
|
sv_midm = mid_ivl.apply_to_endpoint(rc, Interval::LOWER);
|
|
sv_midp = mid_ivl.apply_to_endpoint(rc, Interval::UPPER);
|
|
}
|
|
|
|
} while ( s_midm == CGAL::ZERO || s_midp == CGAL::ZERO ||
|
|
sv_midm != sv_midp + 1 );
|
|
|
|
pop_front( tag );
|
|
if ( sv_midp > sv_b ) {
|
|
Interval ivl_ = ivl.split_at( mid_ivl.upper_bound() ).second;
|
|
|
|
push_front(ivl_, sv_midp, sv_b, tag);
|
|
}
|
|
|
|
push_front(mid, tag);
|
|
|
|
if ( sv_a > sv_midm ) {
|
|
Interval ivl_ = ivl.split_at( mid_ivl.lower_bound() ).first;
|
|
|
|
push_front(ivl_, sv_a, sv_midm, tag);
|
|
}
|
|
}
|
|
else {
|
|
unsigned int sv_mid =
|
|
first_half.apply_to_endpoint(rc, Interval::UPPER);
|
|
|
|
// num roots in (a, mid)
|
|
unsigned int n1 = sv_a - sv_mid;
|
|
// num roots in (mid, b)
|
|
unsigned int n2 = sv_mid - sv_b;
|
|
|
|
CGAL_assertion( n1 + n2 == num_roots_in_first_interval(tag) );
|
|
|
|
pop_front( tag );
|
|
if ( n2 > 0 ) {
|
|
push_front(second_half, sv_mid, sv_b, tag);
|
|
}
|
|
if ( n1 > 0 ) {
|
|
push_front(first_half, sv_a, sv_mid, tag);
|
|
}
|
|
} // end-if
|
|
#else
|
|
// Interval_pair ivl_pair = ivl.split();
|
|
// Interval first_half = ivl_pair.first;
|
|
// Interval second_half = ivl_pair.second;
|
|
|
|
typename Interval::NT a = ivl.lower_bound();
|
|
typename Interval::NT b = ivl.upper_bound();
|
|
|
|
typename Interval::NT mid = ivl.middle();
|
|
|
|
typename Interval::NT a1 = a, b1 = b;
|
|
|
|
unsigned int sv_a = s_variations.front().first;
|
|
unsigned int sv_b = s_variations.front().second;
|
|
|
|
Sign s_mid = apply(sign_at_p, mid);
|
|
|
|
if ( s_mid == CGAL::ZERO ) {
|
|
Sign s_midp, s_midm;
|
|
unsigned int sv_midm(0), sv_midp(0);
|
|
|
|
do {
|
|
a1 = midpoint(a1, mid);
|
|
b1 = midpoint(b1, mid);
|
|
s_midm = apply(sign_at_p, a1);
|
|
s_midp = apply(sign_at_p, b1);
|
|
|
|
if ( s_midm != CGAL::ZERO && s_midp != CGAL::ZERO ) {
|
|
sv_midm = apply(rc, a1);
|
|
sv_midp = apply(rc, b1);
|
|
}
|
|
|
|
} while ( s_midm == CGAL::ZERO || s_midp == CGAL::ZERO ||
|
|
sv_midm != sv_midp + 1 );
|
|
|
|
pop_front( tag );
|
|
if ( sv_midp > sv_b ) {
|
|
Interval ivl_(b1, b);
|
|
push_front(ivl_, sv_midp, sv_b, tag);
|
|
}
|
|
|
|
Interval imid(mid);
|
|
push_front(imid, tag);
|
|
|
|
if ( sv_a > sv_midm ) {
|
|
Interval ivl_(a, a1);
|
|
push_front(ivl_, sv_a, sv_midm, tag);
|
|
}
|
|
}
|
|
else {
|
|
unsigned int sv_mid = apply(rc, mid);
|
|
|
|
// num roots in (a, mid)
|
|
unsigned int n1 = sv_a - sv_mid;
|
|
// num roots in (mid, b)
|
|
unsigned int n2 = sv_mid - sv_b;
|
|
|
|
CGAL_assertion( n1 + n2 == num_roots_in_first_interval(tag) );
|
|
|
|
pop_front( tag );
|
|
if ( n2 > 0 ) {
|
|
Interval ivl_(mid, b);
|
|
push_front(ivl_, sv_mid, sv_b, tag);
|
|
}
|
|
if ( n1 > 0 ) {
|
|
Interval ivl_(a, mid);
|
|
push_front(ivl_, sv_a, sv_mid, tag);
|
|
}
|
|
} // end-if
|
|
#endif
|
|
|
|
} // end-while
|
|
} // end subdivide method
|
|
#if 0
|
|
//------------------------------------------------------------------
|
|
// the method for subdivision; it guarantees that the first interval
|
|
// contains only one root
|
|
template<class SA, class RC>
|
|
void subdivide(Descartes_tag tag, const SA& sign_at_p,
|
|
const RC& rc) const
|
|
{
|
|
if ( i_list.size() == 0 ) { return; }
|
|
|
|
if ( num_roots_in_first_interval(tag) == 1 ) { return; }
|
|
|
|
// MK: WRITE THIS METHOD AS IN THE CASE OF Sturm_tag
|
|
|
|
while ( num_roots_in_first_interval(tag) > 1 ) {
|
|
// std::cout << "." << std::flush;
|
|
Interval ivl = i_list.front();
|
|
Interval_pair ivl_pair = ivl.split();
|
|
Interval first_half = ivl_pair.first;
|
|
Interval second_half = ivl_pair.second;
|
|
|
|
// Sign s_mid = sign_at_p(mid);
|
|
Sign s_mid =
|
|
ivl_pair.first.apply_to_bound(sign_at_p, Interval::UPPER);
|
|
|
|
if ( s_mid == CGAL::ZERO ) {
|
|
Interval mid_ivl = ivl;
|
|
|
|
Sign s_midp, s_midm;
|
|
unsigned int nr(0);
|
|
|
|
do {
|
|
mid_ivl = mid_ivl.middle_half();
|
|
s_midm = mid_ivl.apply_to_bound(sign_at_p, Interval::LOWER);
|
|
s_midp = mid_ivl.apply_to_bound(sign_at_p, Interval::UPPER);
|
|
|
|
// s_midm = sign_at_p(midm);
|
|
// s_midp = sign_at_p(midp);
|
|
|
|
if ( s_midm != CGAL::ZERO && s_midp != CGAL::ZERO ) {
|
|
nr = mid_ivl.apply_to_interval(rc);
|
|
}
|
|
} while ( s_midm == CGAL::ZERO || s_midp == CGAL::ZERO ||
|
|
nr != 1 );
|
|
|
|
// Interval_NT a = ivl.lower_bound();
|
|
// Interval_NT b = ivl.upper_bound();
|
|
|
|
Interval ivl_upper =
|
|
ivl.split_at(mid_ivl, Interval::UPPER).second;
|
|
// Interval::create_from_endpoints(mid_ivl, Interval::UPPER,
|
|
// ivl, Interval::UPPER);
|
|
|
|
Interval ivl_lower =
|
|
ivl.split_at(mid_ivl, Interval::LOWER).first;
|
|
// Interval::create_from_endpoints(ivl, Interval::LOWER,
|
|
// mid_ivl, Interval::LOWER);
|
|
|
|
unsigned int np = ivl_upper.apply_to_interval(rc);
|
|
unsigned int nm = ivl_lower.apply_to_interval(rc);
|
|
|
|
pop_front( tag );
|
|
|
|
if ( np > 0 ) {
|
|
push_front(ivl_upper, np, tag);
|
|
}
|
|
|
|
Interval mid = first_half.collapse_to_bound(Interval::UPPER);
|
|
push_front(mid, tag);
|
|
|
|
if ( nm > 0 ) {
|
|
push_front(ivl_lower, nm, tag);
|
|
}
|
|
|
|
}
|
|
else {
|
|
// Interval_NT a = ivl.lower_bound();
|
|
// Interval_NT b = ivl.upper_bound();
|
|
|
|
int n1 = first_half.apply_to_interval(rc);
|
|
int n2 = second_half.apply_to_interval(rc);
|
|
|
|
pop_front( tag );
|
|
if ( n2 > 0 ) {
|
|
push_front(second_half, n2, tag);
|
|
}
|
|
if ( n1 > 0 ) {
|
|
push_front(first_half, n1, tag);
|
|
}
|
|
|
|
if ( n1 + n2 == 0 ) { break; }
|
|
} // end-if
|
|
} // end-while
|
|
} // end subdivide method
|
|
#endif
|
|
|
|
//------------------------------------------------------------------
|
|
|
|
// the method for initializing the interval list
|
|
template<class RC, class SA>
|
|
void initialize_interval_list(const Interval initial_ivl,
|
|
Sturm_tag tag, const SA& sign_at_p,
|
|
const RC& rc) {
|
|
#if 1
|
|
Sign s_a = initial_ivl.apply_to_endpoint(sign_at_p, Interval::LOWER);
|
|
Sign s_b = initial_ivl.apply_to_endpoint(sign_at_p, Interval::UPPER);
|
|
|
|
if ( s_a != CGAL::ZERO && s_b != CGAL::ZERO ) {
|
|
|
|
unsigned int sv_a = initial_ivl.apply_to_endpoint(rc, Interval::LOWER);
|
|
unsigned int sv_b = initial_ivl.apply_to_endpoint(rc, Interval::UPPER);
|
|
|
|
if ( sv_a > sv_b ) {
|
|
push_front(initial_ivl, sv_a, sv_b, tag);
|
|
}
|
|
return;
|
|
}
|
|
|
|
if ( s_a == CGAL::ZERO && s_b != CGAL::ZERO ) {
|
|
Interval ia = initial_ivl.interval_around_endpoint(Interval::LOWER);
|
|
Interval iam(ia.lower_bound(), initial_ivl.lower_bound());
|
|
Interval iap = initial_ivl.split_at(ia.upper_bound()).first;
|
|
|
|
Sign s_ap, s_am;
|
|
unsigned int sv_ap(0), sv_am(0);
|
|
|
|
do {
|
|
iap = iap.split().first;
|
|
iam = iam.split().second;
|
|
|
|
s_ap = iap.apply_to_endpoint(sign_at_p, Interval::UPPER);
|
|
s_am = iam.apply_to_endpoint(sign_at_p, Interval::LOWER);
|
|
|
|
if ( s_ap != CGAL::ZERO ) {
|
|
sv_ap = iap.apply_to_endpoint(rc, Interval::UPPER);
|
|
}
|
|
if ( s_am != CGAL::ZERO ) {
|
|
sv_am = iam.apply_to_endpoint(rc, Interval::LOWER);
|
|
}
|
|
|
|
} while ( s_ap == CGAL::ZERO || s_am == CGAL::ZERO ||
|
|
sv_ap + 1 != sv_am );
|
|
|
|
unsigned int sv_b =
|
|
initial_ivl.apply_to_endpoint(rc, Interval::UPPER);
|
|
|
|
Interval ivl_upper =
|
|
initial_ivl.split_at(iap.upper_bound()).second;
|
|
|
|
push_front(ivl_upper, sv_ap, sv_b, tag);
|
|
|
|
Interval ivl_a = initial_ivl.endpoint_interval(Interval::LOWER);
|
|
push_front(ivl_a, tag);
|
|
return;
|
|
}
|
|
|
|
if ( s_a != CGAL::ZERO && s_b == CGAL::ZERO ) {
|
|
Interval ib = initial_ivl.interval_around_endpoint(Interval::UPPER);
|
|
|
|
Interval ibm =
|
|
initial_ivl.split_at(ib.lower_bound()).second;
|
|
|
|
Interval ibp(initial_ivl.upper_bound(), ib.upper_bound());
|
|
|
|
Sign s_bp, s_bm;
|
|
unsigned int sv_bp(0), sv_bm(0);
|
|
do {
|
|
ibp = ibp.split().first;
|
|
ibm = ibm.split().second;
|
|
s_bp = ibp.apply_to_endpoint(sign_at_p, Interval::UPPER);
|
|
s_bm = ibm.apply_to_endpoint(sign_at_p, Interval::LOWER);
|
|
|
|
if ( s_bp != CGAL::ZERO ) {
|
|
sv_bp = ibp.apply_to_endpoint(rc, Interval::UPPER);
|
|
}
|
|
if ( s_bm != CGAL::ZERO ) {
|
|
sv_bm = ibm.apply_to_endpoint(rc, Interval::LOWER);
|
|
}
|
|
} while ( s_bp == CGAL::ZERO || s_bm == CGAL::ZERO ||
|
|
sv_bp != sv_bm - 1 );
|
|
|
|
unsigned int sv_a =
|
|
initial_ivl.apply_to_endpoint(rc, Interval::LOWER);
|
|
|
|
Interval ivl_b = initial_ivl.endpoint_interval(Interval::UPPER);
|
|
push_front(ivl_b, tag);
|
|
|
|
Interval ivl_lower =
|
|
initial_ivl.split_at( ibm.lower_bound() ).first;
|
|
|
|
push_front(ivl_lower, sv_a, sv_bm, tag);
|
|
return;
|
|
}
|
|
|
|
if ( s_a == CGAL::ZERO && s_b == CGAL::ZERO ) {
|
|
// both a and b are roots
|
|
Interval ia =
|
|
initial_ivl.interval_around_endpoint(Interval::LOWER);
|
|
Interval ib =
|
|
initial_ivl.interval_around_endpoint(Interval::UPPER);
|
|
|
|
Interval iam(ia.lower_bound(), initial_ivl.lower_bound());
|
|
Interval iap = initial_ivl.split_at(ia.upper_bound()).first;
|
|
Interval ibm = initial_ivl.split_at(ib.lower_bound()).second;
|
|
Interval ibp(initial_ivl.upper_bound(), ib.upper_bound());
|
|
|
|
Sign s_ap, s_am, s_bp, s_bm;
|
|
unsigned int sv_am(0), sv_ap(0), sv_bm(0), sv_bp(0);
|
|
unsigned int nbig(0), nsmall(0);
|
|
|
|
do {
|
|
iap = iap.split().first;
|
|
iam = iam.split().second;
|
|
s_ap = iap.apply_to_endpoint(sign_at_p, Interval::UPPER);
|
|
s_am = iam.apply_to_endpoint(sign_at_p, Interval::LOWER);
|
|
|
|
ibp = ibp.split().first;
|
|
ibm = ibm.split().second;
|
|
s_bp = ibp.apply_to_endpoint(sign_at_p, Interval::UPPER);
|
|
s_bm = ibm.apply_to_endpoint(sign_at_p, Interval::LOWER);
|
|
|
|
if ( s_ap != CGAL::ZERO && s_bm != CGAL::ZERO ) {
|
|
sv_ap = iap.apply_to_endpoint(rc, Interval::UPPER);
|
|
sv_bm = ibm.apply_to_endpoint(rc, Interval::LOWER);
|
|
|
|
nsmall = sv_ap - sv_bm; // num roots in (ap, bm)
|
|
}
|
|
if ( s_am != CGAL::ZERO && s_bp != CGAL::ZERO) {
|
|
sv_am = iam.apply_to_endpoint(rc, Interval::LOWER);
|
|
sv_bp = ibp.apply_to_endpoint(rc, Interval::UPPER);
|
|
|
|
nbig = sv_am - sv_bp; // num roots in (am, bp)
|
|
}
|
|
} while ( s_ap == CGAL::ZERO || s_am == CGAL::ZERO ||
|
|
s_bp == CGAL::ZERO || s_bm == CGAL::ZERO ||
|
|
nbig != nsmall + 2 );
|
|
|
|
Interval ivl_b = initial_ivl.endpoint_interval(Interval::UPPER);
|
|
push_front(ivl_b, tag);
|
|
|
|
Interval ivl_interior(iap.upper_bound(), ibm.lower_bound());
|
|
|
|
push_front(ivl_interior, sv_ap, sv_bm, tag);
|
|
|
|
Interval ivl_a = initial_ivl.endpoint_interval(Interval::LOWER);
|
|
push_front(ivl_a, tag);
|
|
return;
|
|
}
|
|
#else
|
|
typename Interval::NT a = initial_ivl.lower_bound();
|
|
typename Interval::NT b = initial_ivl.upper_bound();
|
|
|
|
Sign s_a = apply(sign_at_p, a);
|
|
Sign s_b = apply(sign_at_p, b);
|
|
|
|
if ( s_a != CGAL::ZERO && s_b != CGAL::ZERO ) {
|
|
|
|
unsigned int sv_a = apply(rc, a);
|
|
unsigned int sv_b = apply(rc, b);
|
|
|
|
if ( sv_a > sv_b ) {
|
|
push_front(initial_ivl, sv_a, sv_b, tag);
|
|
}
|
|
return;
|
|
}
|
|
|
|
if ( s_a == CGAL::ZERO && s_b != CGAL::ZERO ) {
|
|
Interval ia =
|
|
initial_ivl.interval_around_endpoint(Interval::LOWER);
|
|
typename Interval::NT am = ia.lower_bound();
|
|
typename Interval::NT ap = ia.upper_bound();
|
|
|
|
Sign s_ap, s_am;
|
|
unsigned int sv_ap(0), sv_am(0);
|
|
|
|
do {
|
|
ap = midpoint(ap, a);
|
|
am = midpoint(am, a);
|
|
|
|
s_ap = apply(sign_at_p, ap);
|
|
s_am = apply(sign_at_p, am);
|
|
|
|
if ( s_ap != CGAL::ZERO ) {
|
|
sv_ap = apply(rc, ap);
|
|
}
|
|
if ( s_am != CGAL::ZERO ) {
|
|
sv_am = apply(rc, am);
|
|
}
|
|
|
|
} while ( s_ap == CGAL::ZERO || s_am == CGAL::ZERO ||
|
|
sv_ap + 1 != sv_am );
|
|
|
|
unsigned int sv_b = apply(rc, b);
|
|
|
|
Interval ivl_upper(ap, b);
|
|
push_front(ivl_upper, sv_ap, sv_b, tag);
|
|
|
|
Interval ivl_a(a);
|
|
push_front(ivl_a, tag);
|
|
return;
|
|
}
|
|
|
|
if ( s_a != CGAL::ZERO && s_b == CGAL::ZERO ) {
|
|
Interval ib =
|
|
initial_ivl.interval_around_endpoint(Interval::UPPER);
|
|
typename Interval::NT bm = ib.lower_bound();
|
|
typename Interval::NT bp = ib.upper_bound();
|
|
|
|
Sign s_bp, s_bm;
|
|
unsigned int sv_bp(0), sv_bm(0);
|
|
do {
|
|
bp = midpoint(bp, b);
|
|
bm = midpoint(bm, b);
|
|
s_bp = apply(sign_at_p, bp);
|
|
s_bm = apply(sign_at_p, bm);
|
|
|
|
if ( s_bp != CGAL::ZERO ) {
|
|
sv_bp = apply(rc, bp);
|
|
}
|
|
if ( s_bm != CGAL::ZERO ) {
|
|
sv_bm = apply(rc, bm);
|
|
}
|
|
} while ( s_bp == CGAL::ZERO || s_bm == CGAL::ZERO ||
|
|
sv_bp != sv_bm - 1 );
|
|
|
|
unsigned int sv_a = apply(rc, a);
|
|
|
|
Interval ivl_b(b);
|
|
push_front(ivl_b, tag);
|
|
|
|
Interval ivl_lower(a, bm);
|
|
push_front(ivl_lower, sv_a, sv_bm, tag);
|
|
return;
|
|
}
|
|
|
|
if ( s_a == CGAL::ZERO && s_b == CGAL::ZERO ) {
|
|
// both a and b are roots
|
|
Interval ia =
|
|
initial_ivl.interval_around_endpoint(Interval::LOWER);
|
|
Interval ib =
|
|
initial_ivl.interval_around_endpoint(Interval::UPPER);
|
|
|
|
typename Interval::NT am = ia.lower_bound();
|
|
typename Interval::NT ap = ia.upper_bound();
|
|
typename Interval::NT bm = ib.lower_bound();
|
|
typename Interval::NT bp = ib.upper_bound();
|
|
|
|
Sign s_ap, s_am, s_bp, s_bm;
|
|
unsigned int sv_am(0), sv_ap(0), sv_bm(0), sv_bp(0);
|
|
unsigned int nbig(0), nsmall(0);
|
|
|
|
do {
|
|
ap = midpoint(ap, a);
|
|
am = midpoint(am, a);
|
|
bp = midpoint(bp, b);
|
|
bm = midpoint(bm, b);
|
|
|
|
s_ap = apply(sign_at_p, ap);
|
|
s_am = apply(sign_at_p, am);
|
|
s_bp = apply(sign_at_p, bp);
|
|
s_bm = apply(sign_at_p, bm);
|
|
|
|
if ( s_ap != CGAL::ZERO && s_bm != CGAL::ZERO ) {
|
|
sv_ap = apply(rc, ap);
|
|
sv_bm = apply(rc, bm);
|
|
|
|
nsmall = sv_ap - sv_bm; // num roots in (ap, bm)
|
|
}
|
|
if ( s_am != CGAL::ZERO && s_bp != CGAL::ZERO) {
|
|
sv_am = apply(rc, am);
|
|
sv_bp = apply(rc, bp);
|
|
|
|
nbig = sv_am - sv_bp; // num roots in (am, bp)
|
|
}
|
|
} while ( s_ap == CGAL::ZERO || s_am == CGAL::ZERO ||
|
|
s_bp == CGAL::ZERO || s_bm == CGAL::ZERO ||
|
|
nbig != nsmall + 2 );
|
|
|
|
Interval ivl_b(b);
|
|
push_front(ivl_b, tag);
|
|
|
|
Interval ivl_interior(ap, bm);
|
|
push_front(ivl_interior, sv_ap, sv_bm, tag);
|
|
|
|
Interval ivl_a(a);
|
|
push_front(ivl_a, tag);
|
|
return;
|
|
}
|
|
#endif
|
|
//bool this_line_should_not_have_been_reached = false;
|
|
CGAL_postcondition_msg(0,"this_line_should_not_have_been_reached" );
|
|
//if (0) this_line_should_not_have_been_reached=false;
|
|
}
|
|
|
|
#if 0
|
|
//------------------------------------------------------------------
|
|
template<class RC, class SA>
|
|
void initialize_interval_list(const Interval& initial_ivl,
|
|
Bezier_tag, const SA& sign_at_p,
|
|
const RC& rc) {
|
|
initialize_interval_list(initial_ivl, Descartes_tag(),
|
|
sign_at_p, rc);
|
|
}
|
|
|
|
//------------------------------------------------------------------
|
|
|
|
// the method for initializing the interval list
|
|
template<class RC, class SA>
|
|
void initialize_interval_list(const Interval& initial_ivl,
|
|
Descartes_tag tag,
|
|
const SA& sign_at_p,
|
|
const RC& rc) {
|
|
// MK: WRITE THIS METHOD AS IN THE CASE OF Sturm_tag
|
|
|
|
Sign s_a = initial_ivl.apply_to_bound(sign_at_p, Interval::LOWER);
|
|
Sign s_b = initial_ivl.apply_to_bound(sign_at_p, Interval::UPPER);
|
|
|
|
if ( s_a != CGAL::ZERO && s_b != CGAL::ZERO ) {
|
|
int nr = initial_ivl.apply_to_interval(rc);
|
|
|
|
if ( nr > 0 ) {
|
|
push_front(initial_ivl, nr, tag);
|
|
}
|
|
return;
|
|
}
|
|
|
|
if ( s_a == CGAL::ZERO && s_b != CGAL::ZERO ) {
|
|
Interval ia = initial_ivl.interval_around_bound(Interval::LOWER);
|
|
Interval iam(ia.lower_bound(), initial_ivl.lower_bound());
|
|
Interval iap = initial_ivl.split_at(ia, Interval::UPPER).first;
|
|
// Interval_NT ap = a_plus(a, b);
|
|
// Interval_NT am = a - NT(1);
|
|
|
|
Sign s_ap, s_am;
|
|
int np(0), nm(0);
|
|
|
|
do {
|
|
iap = iap.split().first;
|
|
iam = iam.split().second;
|
|
|
|
s_ap = iap.apply_to_bound(sign_at_p, Interval::UPPER);
|
|
s_am = iam.apply_to_bound(sign_at_p, Interval::LOWER);
|
|
|
|
// ap = Interval::midpoint(a, ap);
|
|
// am = Interval::midpoint(a, am);
|
|
// s_ap = sign_at_p(ap);
|
|
// s_am = sign_at_p(am);
|
|
|
|
if ( s_ap != CGAL::ZERO ) {
|
|
np = iap.apply_to_interval(rc);
|
|
// np = root_counter_(ap, b);
|
|
}
|
|
if ( s_am != CGAL::ZERO ) {
|
|
nm = iam.apply_to_interval(rc);
|
|
// nm = root_counter_(am, b);
|
|
}
|
|
|
|
} while ( s_ap == ZERO || s_am == CGAL::ZERO || np + 1 != nm );
|
|
|
|
Interval ivl_upper =
|
|
initial_ivl.split_at(iap, Interval::UPPER).second;
|
|
push_front(ivl_upper, np, tag);
|
|
|
|
Interval ivl_a = initial_ivl.collapse_to_bound(Interval::LOWER);
|
|
push_front(ivl_a, tag);
|
|
|
|
// push_front(ap, b, np, tag);
|
|
// push_front(a, tag);
|
|
return;
|
|
}
|
|
|
|
if ( s_a != CGAL::ZERO && s_b == CGAL::ZERO ) {
|
|
Interval ib = initial_ivl.interval_around_bound(Interval::UPPER);
|
|
|
|
Interval ibm =
|
|
initial_ivl.split_at(ib, Interval::LOWER).second;
|
|
Interval ibp(initial_ivl.upper_bound(), ib.upper_bound());
|
|
// Interval_NT bp = b + NT(1);
|
|
// Interval_NT bm = b_minus(a, b);
|
|
|
|
Sign s_bp, s_bm;
|
|
int np(0), nm(0);
|
|
do {
|
|
ibp = ibp.split().first;
|
|
ibm = ibm.split().second;
|
|
// bp = Interval::midpoint(b, bp);
|
|
// bm = Interval::midpoint(b, bm);
|
|
// s_bp = sign_at_p(bp);
|
|
// s_bm = sign_at_p(bm);
|
|
s_bp = ibp.apply_to_bound(sign_at_p, Interval::UPPER);
|
|
s_bm = ibm.apply_to_bound(sign_at_p, Interval::LOWER);
|
|
|
|
if ( s_bp != CGAL::ZERO ) {
|
|
np = ibp.apply_to_interval(rc);
|
|
// np = root_counter_(a, bp, check);
|
|
}
|
|
if ( s_bm != CGAL::ZERO ) {
|
|
nm = ibm.apply_to_interval(rc);
|
|
// nm = root_counter_(a, bm, check);
|
|
}
|
|
|
|
} while ( s_bp == CGAL::ZERO || s_bm == CGAL::ZERO ||
|
|
np != nm + 1 );
|
|
|
|
Interval ivl_b = initial_ivl.collapse_to_bound(Interval::UPPER);
|
|
push_front(ivl_b, tag);
|
|
|
|
Interval ivl_lower =
|
|
initial_ivl.split_at(ibm, Interval::LOWER).first;
|
|
// Interval::create_from_endpoints(initial_ivl, Interval::LOWER,
|
|
// ibm, Interval::LOWER);
|
|
|
|
push_front(ivl_lower, nm, tag);
|
|
// push_front(b, tag);
|
|
// push_front(a, bm, nm, tag);
|
|
return;
|
|
}
|
|
|
|
if ( s_a == CGAL::ZERO && s_b == CGAL::ZERO ) {
|
|
// both a and b are roots
|
|
// Interval_NT ap = a_plus(a, b);
|
|
// Interval_NT am = a - NT(1);
|
|
// Interval_NT bp = b + NT(1);
|
|
// Interval_NT bm = b_minus(a, b);
|
|
Interval ia =
|
|
initial_ivl.interval_around_bound(Interval::LOWER);
|
|
Interval ib =
|
|
initial_ivl.interval_around_bound(Interval::UPPER);
|
|
|
|
// Interval_NT ap = Interval::midpoint(a, b);
|
|
// Interval_NT am = a - NT(1);
|
|
// Interval_NT bp = b + NT(1);
|
|
// Interval_NT bm = Interval::midpoint(a, b);
|
|
Interval iam(ia.lower_bound(), initial_ivl.lower_bound());
|
|
Interval iap = initial_ivl.split_at(ia, Interval::UPPER).first;
|
|
Interval ibm = initial_ivl.split_at(ib, Interval::LOWER).second;
|
|
Interval ibp(initial_ivl.upper_bound(), ib.lower_bound());
|
|
Interval::create_from_endpoints(initial_ivl, Interval::UPPER,
|
|
ib, Interval::UPPER);
|
|
|
|
Sign s_ap, s_am, s_bp, s_bm;
|
|
int nbig(0), nsmall(0);
|
|
|
|
do {
|
|
iap = iap.split().first;
|
|
iam = iam.split().second;
|
|
s_ap = iap.apply_to_bound(sign_at_p, Interval::UPPER);
|
|
s_am = iam.apply_to_bound(sign_at_p, Interval::LOWER);
|
|
|
|
ibp = ibp.split().first;
|
|
ibm = ibm.split().second;
|
|
s_bp = ibp.apply_to_bound(sign_at_p, Interval::UPPER);
|
|
s_bm = ibm.apply_to_bound(sign_at_p, Interval::LOWER);
|
|
|
|
// ap = Interval::midpoint(a, ap);
|
|
// am = Interval::midpoint(a, am);
|
|
// s_ap = sign_at_p(ap);
|
|
// s_am = sign_at_p(am);
|
|
|
|
// bp = Interval::midpoint(b, bp);
|
|
// bm = Interval::midpoint(b, bm);
|
|
// s_bp = sign_at_p(bp);
|
|
// s_bm = sign_at_p(am);
|
|
|
|
if ( s_ap != CGAL::ZERO && s_bm != CGAL::ZERO ) {
|
|
Interval ismall(iap.upper_bound(), ibm.lower_bound());
|
|
nsmall = ismall.apply_to_interval(rc);
|
|
// nsmall = root_counter_(ap, bm, check);
|
|
}
|
|
if ( s_am != CGAL::ZERO && s_bp != CGAL::ZERO) {
|
|
Interval ibig(iam.lower_bound(), ibp.upper_bound());
|
|
nbig = ibig.apply_to_interval(rc);
|
|
// nbig = root_counter_(am, bp, check);
|
|
}
|
|
} while ( s_ap == CGAL::ZERO || s_am == CGAL::ZERO ||
|
|
s_bp == CGAL::ZERO || s_bm == CGAL::ZERO ||
|
|
nbig != nsmall + 2 );
|
|
|
|
Interval ivl_b = initial_ivl.collapse_to_bound(Interval::UPPER);
|
|
push_front(ivl_b, tag);
|
|
|
|
Interval ivl_interior(iam.upper_bound(), ibm.lower_bound());
|
|
// push_front(ap, bm, sv_ap, sv_bm, tag);
|
|
push_front(ivl_interior, nsmall, tag);
|
|
|
|
Interval ivl_a = initial_ivl.collapse_to_bound(Interval::LOWER);
|
|
push_front(ivl_a, tag);
|
|
|
|
// push_front(b, tag);
|
|
// push_front(ap, bm, nsmall, tag);
|
|
// push_front(a, tag);
|
|
return;
|
|
}
|
|
|
|
bool this_line_should_not_have_been_reached = false;
|
|
CGAL_assertion( this_line_should_not_have_been_reached );
|
|
}
|
|
#endif
|
|
//------------------------------------------------------------------
|
|
|
|
void initialize_root_counter(Sturm_tag) {
|
|
root_counter_ = Root_count(sseq, traits_);
|
|
}
|
|
#if 0
|
|
void initialize_root_counter(Bezier_tag) {
|
|
initialize_root_counter(Descartes_tag());
|
|
}
|
|
|
|
void initialize_root_counter(Descartes_tag) {
|
|
root_counter_ = Root_count(p);
|
|
}
|
|
#endif
|
|
|
|
Interval compute_initial_interval(Method_tag tag) {
|
|
typedef typename Traits::Root_bound RBE;
|
|
|
|
// MK: I MAY NEED TO CHANGE THIS SO THAT THE ROOT BOUND EVALUATOR
|
|
// TAKES A FUNCTION HANDLE; I MAY ALSO WANT TO USE THE
|
|
// INTERVAL VERSION OF THE FUNCTION TO DO THE COMPUTATIONS
|
|
RBE rbe = traits_.root_bound_object(false, tag);
|
|
|
|
typename RBE::result_type root_bound = rbe(p);
|
|
|
|
Root B_lower = -root_bound;
|
|
Root B_upper = root_bound;
|
|
|
|
Root lower = start_;
|
|
Root upper = end_;
|
|
|
|
if ( start_ < B_lower ) {
|
|
lower = B_lower;
|
|
}
|
|
|
|
if ( B_upper < end_ ) {
|
|
upper = B_upper;
|
|
}
|
|
|
|
typedef typename Interval::NT INT;
|
|
|
|
std::pair<double,double> ilower =
|
|
CGAL_POLYNOMIAL_TO_INTERVAL(lower.lower_bound());
|
|
std::pair<double,double> iupper =
|
|
CGAL_POLYNOMIAL_TO_INTERVAL(upper.upper_bound());
|
|
|
|
double low = std::min(ilower.first, iupper.first);
|
|
double high = std::max(ilower.second, iupper.second);
|
|
|
|
if ( CGAL::is_finite(low) && CGAL::is_finite(high) ) {
|
|
CGAL_Polynomial_assertion( low <= high );
|
|
return Interval(INT(low), INT(high));
|
|
}
|
|
|
|
INT low_nt = std::min(lower.lower_bound(), upper.lower_bound());
|
|
INT high_nt = std::max(lower.upper_bound(), upper.upper_bound());
|
|
|
|
CGAL_Polynomial_assertion( low_nt <= high_nt );
|
|
|
|
return Interval(low_nt, high_nt);
|
|
}
|
|
|
|
// initialization method: unused check_if_simple
|
|
void initialize(bool, Method_tag tag) {
|
|
if ( p.is_zero() ) { return; }
|
|
|
|
if ( end_ == -infinity() || start_ == infinity() ) {
|
|
return;
|
|
}
|
|
|
|
Subdivision_tag subdivision_tag;
|
|
|
|
sseq = Standard_sequence(p);
|
|
|
|
#if 0
|
|
normalize(p_simple, tag, success);
|
|
#endif
|
|
|
|
Sign_at sign_at_p(p);
|
|
|
|
initialize_root_counter( subdivision_tag );
|
|
|
|
Interval initial_ivl = compute_initial_interval(tag);
|
|
|
|
initialize_interval_list( initial_ivl, subdivision_tag, sign_at_p,
|
|
root_counter_ );
|
|
}
|
|
|
|
public:
|
|
//===================================
|
|
// POSITIVE INFINITY
|
|
//===================================
|
|
static Root infinity() { return Root::infinity(); }
|
|
|
|
Sturm_root_stack() {}
|
|
|
|
public:
|
|
//==============
|
|
// CONSTRUCTORS
|
|
//==============
|
|
Sturm_root_stack(const Polynomial& p,
|
|
const Root& start = -infinity(),
|
|
const Root& end = infinity(),
|
|
const Traits& tr = Traits(),
|
|
bool check_if_simple = true)
|
|
: start_(start), end_(end), traits_(tr), root_counter_(),
|
|
p(p), root_idx(0) {
|
|
initialize( check_if_simple, Method_tag() );
|
|
pop();
|
|
}
|
|
|
|
//===============
|
|
// DESTRUCTOR
|
|
//===============
|
|
virtual ~Sturm_root_stack() {}
|
|
|
|
//=============
|
|
// METHODS
|
|
//=============
|
|
protected:
|
|
#if 0
|
|
const Root& top(Sturm_tag tag) const
|
|
{
|
|
return current_;
|
|
#if 0
|
|
if ( i_list.size() == 0 ) {
|
|
return infinity();
|
|
}
|
|
else {
|
|
Interval ivl = i_list.front();
|
|
return Root(ivl, p, sseq, root_idx+1);
|
|
}
|
|
#endif
|
|
}
|
|
#endif
|
|
|
|
void pop(Sturm_tag tag) const
|
|
{
|
|
Sign_at sign_at_p(p);
|
|
|
|
do {
|
|
|
|
bool is_even_parity(false);
|
|
do {
|
|
if ( i_list.size() == 0 ) {
|
|
current_ = infinity();
|
|
return;
|
|
}
|
|
else {
|
|
subdivide( tag, sign_at_p, root_counter_ );
|
|
|
|
Interval ivl = i_list.front();
|
|
pop_front( tag );
|
|
// now I pass the simple polynomial and the gcd; I do not
|
|
// need to compute the simple polynomial in the solver?
|
|
// maybe for the descartes solver...
|
|
// cur = Root(ivl, p_simple, p_orig, ++root_idx);
|
|
current_ = Root(ivl, p, sseq, ++root_idx, traits_);
|
|
|
|
} // end-if
|
|
|
|
is_even_parity = current_.is_even_multiplicity();
|
|
|
|
if ( DROP_EVEN_ROOTS && is_even_parity ) {
|
|
root_idx--;
|
|
}
|
|
|
|
} while ( DROP_EVEN_ROOTS && current_ != infinity() &&
|
|
is_even_parity );
|
|
|
|
if ( current_ <= start_ ) {
|
|
root_idx--;
|
|
}
|
|
|
|
} while ( current_ != infinity() && current_ <= start_ );
|
|
|
|
#if 0
|
|
if ( current_ == infinity() || current_ >= end_ ) {
|
|
return infinity();
|
|
}
|
|
|
|
return cur;
|
|
#endif
|
|
}
|
|
|
|
#if 0
|
|
Root next_root(Bezier_tag) const
|
|
{
|
|
return next_root(Descartes_tag());
|
|
}
|
|
|
|
Root next_root(Descartes_tag tag) const
|
|
{
|
|
Root cur;
|
|
Sign_at sign_at_p(p);
|
|
|
|
do {
|
|
bool is_even_parity(false);
|
|
do {
|
|
if ( i_list.size() == 0 ) {
|
|
cur = infinity();
|
|
}
|
|
else {
|
|
subdivide( tag, sign_at_p, root_counter_ );
|
|
|
|
if ( i_list.size() == 0 ) {
|
|
cur = infinity();
|
|
}
|
|
else {
|
|
Interval ivl = i_list.front();
|
|
pop_front( tag );
|
|
// cur = Root(ivl, p_simple, p_orig, ++root_idx);
|
|
cur = Root(ivl, p, sseq[sseq.size() - 1], ++root_idx);
|
|
}
|
|
} // end-if
|
|
|
|
is_even_parity = cur.is_even_parity();
|
|
|
|
if ( DROP_EVEN_ROOTS && is_even_parity ) {
|
|
root_idx--;
|
|
}
|
|
|
|
} while ( DROP_EVEN_ROOTS && cur != infinity() &&
|
|
is_even_parity );
|
|
|
|
if ( cur <= start_ ) {
|
|
root_idx--;
|
|
}
|
|
|
|
} while ( cur != infinity() && cur <= start_ );
|
|
|
|
if ( cur == infinity() || cur >= end_ ) {
|
|
return infinity();
|
|
}
|
|
|
|
return cur;
|
|
}
|
|
#endif
|
|
public:
|
|
const Root& top() const
|
|
{
|
|
return current_;
|
|
}
|
|
|
|
bool empty() const
|
|
{
|
|
// return ( i_list.size() == 0 || current() >= end_ );
|
|
return ( top() >= end_ );
|
|
}
|
|
|
|
void pop() const
|
|
{
|
|
return pop( Subdivision_tag() );
|
|
}
|
|
|
|
const Root& start() const { return start_; }
|
|
const Root& end() const { return end_; }
|
|
|
|
const Standard_sequence& standard_sequence() const { return sseq; }
|
|
Standard_sequence& standard_sequence() { return sseq; }
|
|
|
|
const Root_count& root_count() const { return root_counter_; }
|
|
Root_count& root_count() { return root_counter_; }
|
|
|
|
std::list<uint_pair>& sv_list() { return s_variations; }
|
|
std::list<unsigned int>& nr_list() { return nroots; }
|
|
|
|
Interval_container& ivl_list() { return i_list; }
|
|
const Interval_container& ivl_list() const { return i_list; }
|
|
|
|
int& root_index() { return root_idx; }
|
|
|
|
template<class Stream>
|
|
Stream& write(Stream& os) const
|
|
{
|
|
for (unsigned int i = 0; i < sseq.size(); i++) {
|
|
os << sseq[i] << std::endl;
|
|
}
|
|
|
|
typename Interval_container::iterator ivl_it;
|
|
typename std::list<uint_pair>::iterator nr_it;
|
|
for (ivl_it = i_list.begin(), nr_it = s_variations.begin();
|
|
ivl_it != i_list.end(); ++ivl_it, ++nr_it) {
|
|
os << "{";
|
|
ivl_it->write(os);
|
|
int nr = nr_it->first - nr_it->second;
|
|
os << ", " << nr << "} ";
|
|
}
|
|
os << std::endl;
|
|
return os;
|
|
}
|
|
|
|
protected:
|
|
Root start_, end_;
|
|
Traits traits_;
|
|
Root_count root_counter_;
|
|
Polynomial p;
|
|
Standard_sequence sseq;
|
|
mutable Root current_;
|
|
mutable std::list<uint_pair> s_variations;
|
|
mutable std::list<unsigned int> nroots;
|
|
mutable Interval_container i_list;
|
|
mutable int root_idx;
|
|
};
|
|
|
|
template<class Stream, class T, class S, class M, bool D>
|
|
Stream& operator<<(Stream& os,
|
|
const Sturm_root_stack<T,S,M,D>& solver)
|
|
{
|
|
return solver.write(os);
|
|
}
|
|
|
|
|
|
CGAL_POLYNOMIAL_END_NAMESPACE
|
|
#endif // CGAL_STURM_LAZY_SOLVER_H
|