mirror of https://github.com/CGAL/cgal
276 lines
8.1 KiB
C++
276 lines
8.1 KiB
C++
#include <iostream>
|
|
#include <cassert>
|
|
#include <cstdlib>
|
|
#include <CGAL/Polynomial/Kernel.h>
|
|
#include <CGAL/Polynomial/Numeric_root_stack.h>
|
|
#include <CGAL/Polynomial/Polynomial.h>
|
|
#include <CGAL/Polynomial/Root_stack_default_traits.h>
|
|
//#include <CGAL/Polynomial/Upper_bound_root_stack_Descartes_traits.h>
|
|
//#include <CGAL/Polynomial/Upper_bound_root_stack.h>
|
|
#include <CGAL/Polynomial/internal/numeric_solvers.h>
|
|
|
|
#ifdef CGAL_POLYNOMIAL_USE_GSL
|
|
#include <CGAL/Polynomial/internal/GSL_numeric_solver.h>
|
|
#endif
|
|
|
|
#include "Check_solver.h"
|
|
bool verbose=true;
|
|
typedef CGAL_POLYNOMIAL_NS::Polynomial<double> Pd;
|
|
typedef CGAL_POLYNOMIAL_NS::Polynomial<CGAL::POLYNOMIAL::Default_field_nt> Pe;
|
|
typedef CGAL_POLYNOMIAL_NS::Root_stack_default_traits<Pd> Dt;
|
|
|
|
struct Interval_root_stack {
|
|
typedef CGAL::POLYNOMIAL::Interval_polynomial Pi;
|
|
typedef CGAL::POLYNOMIAL::Default_field_nt ENT;
|
|
Interval_root_stack(){}
|
|
typedef CGAL::POLYNOMIAL::Interval_nt Root;
|
|
typedef Root Interval;
|
|
typedef Dt Traits;
|
|
Interval_root_stack(const Pd &p, Root lb, Root ub, Traits): lb_(lb),
|
|
ub_(ub) {
|
|
CGAL_POLYNOMIAL_NS::Polynomial_converter<Pd, Pi, CGAL_POLYNOMIAL_NS::To_interval<double> > pc;
|
|
CGAL_POLYNOMIAL_NS::Polynomial_converter<Pd, Pe, CGAL::NT_converter<double, ENT> > pce;
|
|
p_=pc(p);
|
|
pe_= pce(p);
|
|
ok_=false;
|
|
stack_.push_back(Root(std::numeric_limits<double>::infinity()));
|
|
stack_.push_back(Root(lb.inf(), ub.sup()));
|
|
refine();
|
|
}
|
|
|
|
const Root &top() const {
|
|
CGAL_assertion(ok_);
|
|
return stack_.back();
|
|
}
|
|
|
|
void pop() {
|
|
if (stack_.size()>1) {
|
|
CGAL_assertion(!stack_.empty());
|
|
stack_.pop_back();
|
|
ok_=false;
|
|
if (stack_.size() >1) refine();
|
|
}
|
|
}
|
|
|
|
bool empty() const {
|
|
return stack_.size()==1;
|
|
}
|
|
|
|
static bool is_unknown_sign(Interval rv) {
|
|
return rv.inf() <=0 && rv.sup() >=0;
|
|
}
|
|
|
|
bool is_odd(double lb, double ub) const {
|
|
// put in filtered evaluation
|
|
Interval lbi= p_(Interval(lb));
|
|
Interval ubi= p_(Interval(ub));
|
|
bool ret=false;
|
|
if (is_unknown_sign(lbi)
|
|
|| is_unknown_sign(ubi)){
|
|
ret= CGAL::sign(pe_(ENT(lb))) != CGAL::sign(pe_(ENT(ub)))
|
|
|| CGAL::sign(pe_(ENT(lb))) == CGAL::ZERO;
|
|
} else {
|
|
ret= CGAL::sign(lbi) != CGAL::sign(ubi)
|
|
|| CGAL::sign(lbi) == CGAL::ZERO;
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
bool no_root(Interval r) const {
|
|
Interval rv= p_(r);
|
|
return !is_unknown_sign(rv);
|
|
}
|
|
|
|
static bool is_small(Interval r) {
|
|
return r.sup()- r.inf() < .001;
|
|
}
|
|
|
|
double mp(double lb, double ub) const {
|
|
if (lb == -std::numeric_limits<double>::infinity()
|
|
&& ub == std::numeric_limits<double>::infinity()){
|
|
return 0;
|
|
} else if (lb == -std::numeric_limits<double>::infinity()) {
|
|
return ub-1000000;
|
|
} else if (ub == std::numeric_limits<double>::infinity()){
|
|
return lb + 1000000;
|
|
} else {
|
|
return (lb+ub)/2.0;
|
|
}
|
|
}
|
|
|
|
|
|
|
|
void refine() {
|
|
CGAL_POLYNOMIAL_NS::Interval_arithmetic_guard gd;
|
|
CGAL_assertion(stack_.size()>1);
|
|
while (stack_.size() >1 ){
|
|
if (! is_small(stack_.back())) {
|
|
//std::cout << "Splitting " << stack_.back() << std::endl;
|
|
Interval b= stack_.back();
|
|
stack_.pop_back();
|
|
Interval fh= Interval(b.inf(), mp(b.inf(),b.sup()));
|
|
Interval sh(fh.sup(), b.sup());
|
|
stack_.push_back(sh);
|
|
stack_.push_back(fh);
|
|
}
|
|
bool changed=false;
|
|
do {
|
|
changed=false;
|
|
while(stack_.size() >1&&no_root(stack_.back()) ){
|
|
CGAL_assertion(!stack_.empty());
|
|
//std::cout << "Discarding " << stack_.back() << std::endl;
|
|
stack_.pop_back();
|
|
changed =true;
|
|
}
|
|
|
|
while (stack_.size() >1 && is_small(stack_.back())
|
|
&& !is_odd(stack_.back().inf(), stack_.back().sup())) {
|
|
CGAL_assertion(!stack_.empty());
|
|
//std::cout << "Discarding " << stack_.back() << std::endl;
|
|
stack_.pop_back();
|
|
changed =true;
|
|
}
|
|
} while (changed);
|
|
if (is_small(stack_.back())) break;
|
|
}
|
|
|
|
if (!empty()) {
|
|
Interval b= stack_.back();
|
|
CGAL_assertion(is_small(b));
|
|
CGAL_assertion(is_odd(b.inf(),b.sup()));
|
|
}
|
|
ok_=true;
|
|
}
|
|
|
|
bool ok_;
|
|
std::vector<Root> stack_;
|
|
Root lb_, ub_;
|
|
CGAL::POLYNOMIAL::Interval_polynomial p_;
|
|
Pe pe_;
|
|
};
|
|
|
|
|
|
namespace std {
|
|
template <>
|
|
struct numeric_limits<CGAL_POLYNOMIAL_NS::Interval_nt>: public numeric_limits<double> {
|
|
typedef numeric_limits<double> P;
|
|
typedef CGAL_POLYNOMIAL_NS::Interval_nt T;
|
|
static const bool is_specialized = true;
|
|
static T min BOOST_PREVENT_MACRO_SUBSTITUTION () throw() {return T((P::min)());}
|
|
static T max BOOST_PREVENT_MACRO_SUBSTITUTION () throw() {return T((P::max)());}
|
|
static T infinity() throw() {return P::infinity();}
|
|
};
|
|
}
|
|
|
|
|
|
int main(int argc, char* argv[])
|
|
{
|
|
//assert(std::numeric_limits<double>::has_infinity());
|
|
if ( argc > 1 ) {
|
|
int is_verbose = std::atoi(argv[1]);
|
|
if ( is_verbose == 0 ) {
|
|
verbose = false;
|
|
} else verbose = true;
|
|
}
|
|
|
|
|
|
#ifdef CGAL_USE_TNT
|
|
{
|
|
if (verbose) std::cout <<"JAMA______________________________________\n";
|
|
else std::cout << "JAMA &\t";
|
|
typedef CGAL_POLYNOMIAL_NS::Numeric_root_stack<Dt,
|
|
CGAL_POLYNOMIAL_NS::internal::JAMA_numeric_solver> NRE;
|
|
typedef CGAL_POLYNOMIAL_NS::Kernel<Pd, NRE> K;
|
|
K k;
|
|
Check_solver<K > cg(k,verbose);
|
|
cg.all();
|
|
std::cout << std::endl;
|
|
}
|
|
#endif
|
|
if (0) {
|
|
if (verbose) std::cout <<"Inter______________________________________\n";
|
|
else std::cout << "Interval &\t";
|
|
typedef CGAL_POLYNOMIAL_NS::Kernel<Pd, Interval_root_stack> K;
|
|
K k;
|
|
Check_solver<K > cg(k,verbose);
|
|
cg.all();
|
|
std::cout << std::endl;
|
|
}
|
|
#ifdef CGAL_POLYNOMIAL_USE_GSL
|
|
{
|
|
if (verbose) std::cout <<"GSL________________________________________\n";
|
|
else std::cout << "GSL &\t";
|
|
typedef CGAL_POLYNOMIAL_NS::Numeric_root_stack<Dt, CGAL_POLYNOMIAL_NS::internal::GSL_numeric_solver> NRE;
|
|
typedef CGAL_POLYNOMIAL_NS::Kernel<Pd, NRE> K;
|
|
K k;
|
|
Check_solver<K > cg(k,verbose);
|
|
cg.all();
|
|
std::cout << std::endl;
|
|
}
|
|
#endif
|
|
{
|
|
if (verbose) std::cout <<"Turk______________________________________\n";
|
|
else std::cout << "Turk &\t";
|
|
typedef CGAL_POLYNOMIAL_NS::Numeric_root_stack<Dt,
|
|
CGAL_POLYNOMIAL_NS::internal::Turkowski_numeric_solver> NRE;
|
|
typedef CGAL_POLYNOMIAL_NS::Kernel<Pd, NRE> K;
|
|
K k;
|
|
Check_solver<K > cg(k,verbose);
|
|
cg.all();
|
|
std::cout << std::endl;
|
|
}
|
|
|
|
#ifdef CGAL_POLYNOMIAL_USE_GSL
|
|
{
|
|
if (verbose) std::cout <<"CleanGSL__________________________________\n";
|
|
else std::cout << "CleanGSL &\t";
|
|
typedef CGAL_POLYNOMIAL_NS::Numeric_root_stack<Dt,
|
|
CGAL_POLYNOMIAL_NS::internal::GSL_cleaned_numeric_solver> NRE;
|
|
typedef CGAL_POLYNOMIAL_NS::Kernel<Pd, NRE> K;
|
|
K k;
|
|
Check_solver<K > cg(k,verbose);
|
|
cg.cleaned();
|
|
std::cout << std::endl;
|
|
}
|
|
#endif
|
|
|
|
#ifdef CGAL_USE_TNT
|
|
{
|
|
if (verbose) std::cout <<"CleanJAMA________________________________\n";
|
|
else std::cout << "CleanJAMA &\t";
|
|
typedef CGAL_POLYNOMIAL_NS::Numeric_root_stack<Dt,
|
|
CGAL_POLYNOMIAL_NS::internal::JAMA_cleaned_numeric_solver> NRE;
|
|
typedef CGAL_POLYNOMIAL_NS::Kernel<Pd, NRE> K;
|
|
K k;
|
|
Check_solver<K > cg(k,verbose);
|
|
cg.cleaned();
|
|
std::cout << std::endl;
|
|
}
|
|
#endif
|
|
{
|
|
if (verbose) std::cout <<"CleanTurk________________________________\n";
|
|
else std::cout << "CleanTurk &\t";
|
|
typedef CGAL_POLYNOMIAL_NS::Numeric_root_stack<Dt,
|
|
CGAL_POLYNOMIAL_NS::internal::Turkowski_cleaned_numeric_solver> NRE;
|
|
typedef CGAL_POLYNOMIAL_NS::Kernel<Pd, NRE> K;
|
|
K k;
|
|
Check_solver<K > cg(k,verbose);
|
|
cg.cleaned();
|
|
std::cout << std::endl;
|
|
}
|
|
|
|
/*{
|
|
if (verbose) std::cout <<"Descartes__________________________________\n";
|
|
else std::cout << "Descartes &\t";
|
|
typedef CGAL_POLYNOMIAL_NS::Upper_bound_enumerator_Descartes_traits<Pd> Dt;
|
|
typedef CGAL_POLYNOMIAL_NS::Upper_bound_root_enumerator<Dt> NRE;
|
|
typedef CGAL_POLYNOMIAL_NS::Kernel<Pd, NRE> K;
|
|
K k;
|
|
Check_solver<K > cg(k,verbose);
|
|
cg.all();
|
|
std::cout << std::endl;
|
|
}*/
|
|
|
|
return 0;
|
|
}
|