songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

added rational kernel 

Given that the integer kernel works like a charm, I put the rational
one.
This commit is contained in:
Luis Peñaranda 2013-11-23 02:44:59 -02:00
parent ed04ff7c18
commit 112385eef0
9 changed files with 245 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
Algebraic_kernel_d/include/CGAL/Algebraic_kernel_rs_gmpq_d_1.h
View File

@ -1,9 +1,9 @@
// Copyright (c) 2006-2010 Inria Lorraine (France). All rights reserved.
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). 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; either version 3 of the License,
// or (at your option) any later version.
// 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.
@ -19,17 +19,31 @@
#ifndef CGAL_ALGEBRAIC_KERNEL_RS_GMPQ_1
#define CGAL_ALGEBRAIC_KERNEL_RS_GMPQ_1
#include <CGAL/config.h>
#ifdef CGAL_USE_RS
#include <CGAL/Gmpq.h>
#include <CGAL/RS/Algebraic_kernel_rs_1.h>
#include <CGAL/Gmpfr.h>
#include <CGAL/Polynomial.h>
#include "RS/functors_1.h"
#include "RS/rs2_isolator_1.h"
#include "RS/bisection_refiner_1.h"
#include "RS/ak_1.h"
namespace CGAL{
typedef Algebraic_kernel_rs_1<Gmpq> Algebraic_kernel_rs_gmpq_d_1;
// Choice of the isolator: RS default.
typedef CGAL::RS2::RS2_isolator_1<CGAL::Polynomial<CGAL::Gmpq>,CGAL::Gmpfr>
Isolator;
// Choice of the refiner: bisection.
typedef CGAL::Bisection_refiner_1<CGAL::Polynomial<CGAL::Gmpq>,CGAL::Gmpfr>
Refiner;
typedef CGAL::RS_AK1::Algebraic_kernel_1<
CGAL::Polynomial<CGAL::Gmpq>,
CGAL::Gmpfr,
Isolator,
Refiner>
Algebraic_kernel_rs_gmpq_d_1;
}
#endif
#endif // CGAL_ALGEBRAIC_KERNEL_RS_GMPQ_1

3
Algebraic_kernel_d/include/CGAL/Algebraic_kernel_rs_gmpz_d_1.h
View File

@ -1,5 +1,4 @@
// Copyright (c) 2011 National and Kapodistrian University of Athens (Greece).
// All rights reserved.
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). 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

40
Algebraic_kernel_d/include/CGAL/RS/algebraic_1.h
View File

@ -91,7 +91,7 @@ boost::totally_ordered<Algebraic_1<Polynomial_,
}
Algebraic_1(const Algebraic &a):
pol(a.pol),left(a.left),right(a.right){}
// This assumes that Gmpq is constructible from bound type and
// XXX: This assumes that Gmpq is constructible from bound type and
// that polynomial coefficient type is constructible from mpz_t.
Algebraic_1(const Bound &b):left(b),right(b){
typedef typename Ptraits::Shift shift;
@ -101,22 +101,9 @@ boost::totally_ordered<Algebraic_1<Polynomial_,
Coefficient(mpq_numref(q.mpq()));
CGAL_assertion(left==right&&left==b);
}
// TODO: make this constructor generic, the coefficient type is
// assumed to be constructible from mpz_t (rewrite in terms of
// Gmpq/Gmpz)
Algebraic_1(double d){
typedef typename Ptraits::Shift shift;
Gmpq q(d);
pol=Coefficient(mpq_denref(q.mpq()))*
shift()(Polynomial(1),1,0)-
Coefficient(mpq_numref(q.mpq()));
left=Bound(d/*,std::round_toward_neg_infinity*/);
right=Bound(d/*,std::round_toward_infinity*/);
CGAL_assertion((left==right&&left==d)||(left<d&&right>d));
}
// TODO: Constructors from types such as int, unsigned and long. This
// implementation assumes that the bound type is Gmpfr and that T can
// be exactly converted to Gmpq.
// XXX: This implementation assumes that the bound type is Gmpfr
// and that T can be exactly converted to Gmpq. This constructor
// can handle types such as int, unsigned and long.
template <class T>
Algebraic_1(const T &t){
typedef typename Ptraits::Shift shift;
@ -128,6 +115,23 @@ boost::totally_ordered<Algebraic_1<Polynomial_,
right=Bound(t,std::round_toward_infinity);
CGAL_assertion(left<=t&&right>=t);
}
// XXX: This constructor assumes the bound is a Gmpfr.
Algebraic_1(const CGAL::Gmpq &q){
typedef typename Ptraits::Shift shift;
pol=Coefficient(mpq_denref(q.mpq()))*
shift()(Polynomial(1),1,0)-
Coefficient(mpq_numref(q.mpq()));
left=Bound();
right=Bound();
mpfr_t b;
mpfr_init(b);
mpfr_set_q(b,q.mpq(),GMP_RNDD);
mpfr_swap(b,left.fr());
mpfr_set_q(b,q.mpq(),GMP_RNDU);
mpfr_swap(b,right.fr());
mpfr_clear(b);
CGAL_assertion(left<=q&&right>=q);
}
~Algebraic_1(){}
Algebraic_1& operator=(const Algebraic &a){
@ -316,11 +320,11 @@ std::ostream& operator<<(std::ostream &o,
a.get_right()<<']');
}
// XXX: This function works, but it will be nice to rewrite it cleanly.
template <class P,class B,class R,class C,class T>
inline
std::istream& operator>>(std::istream &i,
RS_AK1::Algebraic_1<P,B,R,C,T> &a){
// TODO: cleanly write this function
std::istream::int_type c;
P pol;
B lb,rb;

84
Algebraic_kernel_d/include/CGAL/RS/bisection_refiner_1.h
View File

@ -33,7 +33,7 @@ struct Bisection_refiner_1{
void operator()(const Polynomial_&,Bound_&,Bound_&,int);
}; // class Bisection_refiner_1
// TODO: Write in a generic way, if possible.
// TODO: Write in a generic way, if possible (see next function).
template <class Polynomial_,class Bound_>
void
Bisection_refiner_1<Polynomial_,Bound_>::
@ -42,7 +42,9 @@ operator()(const Polynomial_&,Bound_&,Bound_&,int){
return;
}
// TODO: Optimize this function.
// This works with any type of polynomial, but only for Gmpfr bounds.
// TODO: Beyond writing generically, optimize this function. This would
// remove the need for the next function, which essentially the same.
template<>
void
Bisection_refiner_1<Polynomial<Gmpz>,Gmpfr>::
@ -121,6 +123,84 @@ operator()(const Polynomial<Gmpz> &pol,Gmpfr &left,Gmpfr &right,int prec){
return;
}
template<>
void
Bisection_refiner_1<Polynomial<Gmpq>,Gmpfr>::
operator()(const Polynomial<Gmpq> &pol,Gmpfr &left,Gmpfr &right,int prec){
typedef Polynomial<Gmpq> Polynomial;
typedef Polynomial_traits_d<Polynomial> Ptraits;
typedef Ptraits::Make_square_free Sfpart;
typedef CGAL::RS_AK1::Signat_1<Polynomial,Gmpfr>
Signat;
CGAL_precondition(left<=right);
// TODO: add precondition to check whether the interval is a point
// or the evaluations on its endpoints have different signs
//std::cout<<"refining ["<<left<<","<<right<<"]"<<std::endl;
#ifndef CGAL_GMPFR_NO_REFCOUNT
// Make sure the endpoints do not share references. If some of them
// does, copy it.
if(!left.is_unique()){
Gmpfr new_left(0,left.get_precision());
mpfr_set(new_left.fr(),left.fr(),GMP_RNDN);
left=new_left;
CGAL_assertion_code(new_left=Gmpfr();)
CGAL_assertion(left.is_unique());
}
if(!right.is_unique()){
Gmpfr new_right(0,right.get_precision());
mpfr_set(new_right.fr(),right.fr(),GMP_RNDN);
right=new_right;
CGAL_assertion_code(new_right=Gmpfr();)
CGAL_assertion(right.is_unique());
}
#endif // CGAL_GMPFR_NO_REFCOUNT
Polynomial sfpp=Sfpart()(pol);
Signat signof(sfpp);
CGAL::Sign sl,sc;
mp_prec_t pl,pc;
mpfr_t center;
sl=signof(left);
if(sl==ZERO)
return;
pl=left.get_precision();
pc=right.get_precision();
pc=(pl>pc?pl:pc)+(mp_prec_t)prec;
mpfr_init2(center,pc);
CGAL_assertion_code(int round=)
mpfr_prec_round(left.fr(),pc,GMP_RNDN);
CGAL_assertion(!round);
CGAL_assertion_code(round=)
mpfr_prec_round(right.fr(),pc,GMP_RNDN);
CGAL_assertion(!round);
for(int i=0;i<prec;++i){
CGAL_assertion_code(round=)
mpfr_add(center,left.fr(),right.fr(),GMP_RNDN);
CGAL_assertion(!round);
CGAL_assertion_code(round=)
mpfr_div_2ui(center,center,1,GMP_RNDN);
CGAL_assertion(!round);
sc=signof(Gmpfr(center));
if(sc==ZERO){ // we have a root
CGAL_assertion_code(round=)
mpfr_set(left.fr(),center,GMP_RNDN);
CGAL_assertion(!round);
mpfr_swap(right.fr(),center);
break;
}
if(sc==sl)
mpfr_swap(left.fr(),center);
else
mpfr_swap(right.fr(),center);
}
mpfr_clear(center);
CGAL_postcondition(left<=right);
//std::cout<<"ref root is ["<<left<<","<<right<<"]"<<std::endl;
return;
}
} // namespace CGAL
#endif // CGAL_RS_BISECTION_REFINER_1_H

34
Algebraic_kernel_d/include/CGAL/RS/rs2_isolator_1.h
View File

@ -51,6 +51,7 @@ RS2_isolator_1(const Polynomial_ &p){
CGAL_error_msg("not implemented for these polynomial/bound types");
}
// Isolator of Gmpz polynomials and Gmpfr bounds (the best case for RS).
template <>
RS2_isolator_1<CGAL::Polynomial<CGAL::Gmpz>,Gmpfr>::
RS2_isolator_1(const CGAL::Polynomial<CGAL::Gmpz> &p):_polynomial(p){
@ -69,6 +70,39 @@ RS2_isolator_1(const CGAL::Polynomial<CGAL::Gmpz> &p):_polynomial(p){
free(intervals_mpfi);
}
// Isolator of Gmpq polynomials and Gmpfr bounds. The polynomial must be
// converted to a Gmpz one with the same roots.
template <>
RS2_isolator_1<CGAL::Polynomial<CGAL::Gmpq>,Gmpfr>::
RS2_isolator_1(const CGAL::Polynomial<CGAL::Gmpq> &p):_polynomial(p){
unsigned degree=p.degree();
// Compute first the LCM of the divisors of the coefficients.
mpz_t lcm;
mpz_init(lcm);
mpz_lcm(lcm,mpq_denref(p[0].mpq()),mpq_denref(p[degree].mpq()));
for(unsigned i=1;i<degree;++i)
mpz_lcm(lcm,lcm,mpq_denref(p[i].mpq()));
// Create an array of mpz_t, to store the new integer polynomial.
mpz_t *coeffs=(mpz_t*)malloc((degree+1)*sizeof(mpz_t));
for(unsigned i=0;i<=degree;++i){
mpz_init(coeffs[i]);
mpz_div(coeffs[i],lcm,mpq_denref(p[i].mpq()));
mpz_mul(coeffs[i],coeffs[i],mpq_numref(p[i].mpq()));
}
// Call RS to solve the computed integer polynomial.
RS2::RS2_calls::init_solver();
RS2::RS2_calls::create_rs_upoly(coeffs,degree,rs_get_default_up());
set_rs_precisol(0);
set_rs_verbose(0);
rs_run_algo((char*)"UISOLE");
RS2::RS2_calls::insert_roots(std::back_inserter(_real_roots));
// Done. Clear the mpz_t's and the free the array.
mpz_clear(lcm);
for(unsigned i=0;i<=degree;++i)
mpz_clear(coeffs[i]);
free(coeffs);
}
template <class Polynomial_,class Bound_>
Polynomial_
RS2_isolator_1<Polynomial_,Bound_>::polynomial()const{

22
Algebraic_kernel_d/include/CGAL/RS/signat_1.h
View File

@ -77,6 +77,28 @@ Signat_1<Polynomial<Gmpz>,Gmpfr>::operator()(const Gmpfr &x)const{
return h.sign();
}
// This is the same code as above.
template <>
inline CGAL::Sign
Signat_1<Polynomial<Gmpq>,Gmpfr>::operator()(const Gmpfr &x)const{
typedef Signat_1<Polynomial,Gmpq> Exact_sign;
int d=Degree()(pol);
if(d==0)
return pol[0].sign();
Gmpfi h(pol[d],x.get_precision()+2*d);
Uncertain<CGAL::Sign> indet=Uncertain<CGAL::Sign>::indeterminate();
if(h.sign().is_same(indet))
return Exact_sign(pol)(x);
for(int i=1;i<=d;++i){
h*=x;
h+=pol[d-i];
if(h.sign().is_same(indet))
return Exact_sign(pol)(x);
}
CGAL_assertion(!h.sign().is_same(indet));
return h.sign();
}
} // namespace RS_AK1
} // namespace CGAL

102
Algebraic_kernel_d/test/Algebraic_kernel_d/Algebraic_kernel_rs_gmpq_d_1.cpp
View File

@ -1,9 +1,9 @@
// Copyright (c) 2009,2010 Inria Lorraine (France). All rights reserved.
// Copyright (c) 2006-2013 INRIA Nancy-Grand Est (France). 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; either version 3 of the License,
// or (at your option) any later version.
// 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.
@ -16,112 +16,124 @@
//
// Author: Luis Peñaranda <luis.penaranda@gmx.com>
#include <CGAL/basic.h>
#include <CGAL/config.h>
#if defined(CGAL_USE_GMP) && defined(CGAL_USE_MPFI) && defined(CGAL_USE_RS)
#include <CGAL/Algebraic_kernel_rs_gmpq_d_1.h>
#include "include/CGAL/_test_algebraic_kernel_1.h"
int main(){
#include <CGAL/Algebraic_kernel_rs_gmpq_d_1.h>
typedef CGAL::Algebraic_kernel_rs_gmpq_d_1 AK;
typedef AK::Polynomial_1 Polynomial_1;
typedef AK::Coefficient Coefficient;
typedef AK::Bound Bound;
typedef AK::Algebraic_real_1 Algebraic_real_1;
typedef AK::Multiplicity_type Multiplicity_type;
template <class AK_>
int test_ak(){
typedef AK_ AK;
typedef typename AK::Polynomial_1 Polynomial_1;
//typedef typename AK::Coefficient Coefficient;
typedef typename AK::Bound Bound;
typedef typename AK::Algebraic_real_1 Algebraic_real_1;
typedef typename AK::Multiplicity_type Multiplicity_type;
AK ak;
CGAL::test_algebraic_kernel_1<AK>(ak);
AK ak; // an algebraic kernel object
CGAL::test_algebraic_kernel_1<AK>(ak); // we run standard tests first
AK::Solve_1 solve_1 = ak.solve_1_object();
Polynomial_1 x = CGAL::shift(AK::Polynomial_1(1),1);
typename AK::Solve_1 solve_1 = ak.solve_1_object();
Polynomial_1 x = CGAL::shift(Polynomial_1(1),1);
int returnvalue=0;
// variant using a bool indicating a square free polynomial
// multiplicities are not computed
std::vector<Algebraic_real_1> roots;
solve_1(CGAL::Gmpq(1,2)*x*x-1,true, std::back_inserter(roots));
solve_1(x*x-2,true, std::back_inserter(roots));
if(roots.size()!=2){
returnvalue-=1;
std::cerr<<"error 1: the number of roots of"<<
" (1/2)x^2-1 must be 2"<<std::endl;
std::cerr<<"error 1: the number of roots of x^2-2 must be 2"<<
std::endl;
}
if(-1.42>=roots[0] || -1.41<=roots[0] ||
1.41>=roots[1] || 1.42<=roots[1]){
returnvalue-=2;
std::cerr<<"error 2: the roots of (1/2)x^2-1 are wrong"<<
std::endl;
std::cerr<<"error 2: the roots of x^2-2 are wrong"<<std::endl;
}
roots.clear();
// variant for roots in a given range of a square free polynomial
solve_1((CGAL::Gmpq(1,4)*x*x-CGAL::Gmpq(1,2))*(x*x-3),
true,
Bound(0),
Bound(10),
solve_1((x*x-2)*(x*x-3),true, Bound(0),Bound(10),
std::back_inserter(roots));
if(roots.size()!=2){
returnvalue-=4;
std::cerr<<"error 3: the number of roots of"<<
" ((1/4)x^2-(1/2))*(x^2-3) between 0 and 10 must be 2"
<<std::endl;
std::cerr<<"error 3: the number of roots of (x^2-2)*(x^2-3)"<<
" between 0 and 10 must be 2"<<std::endl;
}
if(1.41>=roots[0] || 1.42<=roots[0] ||
1.73>=roots[1] || 1.74<=roots[1]){
returnvalue-=8;
std::cerr<<"error 4: the roots of ((1/4)x^2-(1/2))*(x^2-3)"<<
std::cerr<<"error 4: the roots of (x^2-2)*(x^2-3)"<<
" between 0 and 10 are wrong"<<std::endl;
}
roots.clear();
// variant computing all roots with multiplicities
std::vector<std::pair<Algebraic_real_1,Multiplicity_type> > mroots;
solve_1(CGAL::Gmpq(1,4)*x*x-CGAL::Gmpq(1,2),std::back_inserter(mroots));
solve_1((x*x-2), std::back_inserter(mroots));
if(mroots.size()!=2){
returnvalue-=16;
std::cerr<<"error 5: the number of roots of"<<
" (1/4)x^2-(1/2) must be 2"<<std::endl;
std::cerr<<"error 5: the number of roots of x^2-2 must be 2"<<
std::endl;
}
if(-1.42>=mroots[0].first || -1.41<=mroots[0].first ||
1.41>=mroots[1].first || 1.42<=mroots[1].first){
returnvalue-=32;
std::cerr<<"error 6: the roots of (1/4)x^2-(1/2) are wrong"<<
std::endl;
std::cerr<<"error 6: the roots of x^2-2 are wrong"<<std::endl;
}
if(mroots[0].second!=1 && mroots[1].second!=1){
returnvalue-=64;
std::cerr<<"error 7: the multiplicities of the"<<
" roots of (1/4)x^2-(1/2) are wrong"<<std::endl;
" roots of x^2-2 are wrong"<<std::endl;
}
mroots.clear();
// variant computing roots with multiplicities for a range
solve_1((x*x-2)*(CGAL::Gmpq(1,2)*x*x-CGAL::Gmpq(3,2)),
Bound(0),
Bound(10),
std::back_inserter(mroots));
solve_1((x*x-2)*(x*x-3),Bound(0),Bound(10),std::back_inserter(mroots));
if(mroots.size()!=2){
returnvalue-=128;
std::cerr<<"error 8: the number of roots of"<<
" (x^2-2)*((1/2)x^2-(3/2)) between 0 and 10 must be 2"
<<std::endl;
std::cerr<<"error 8: the number of roots of (x^2-2)*(x^2-3)"<<
" between 0 and 10 must be 2"<<std::endl;
}
if(1.41>=mroots[0].first || 1.42<=mroots[0].first ||
1.73>=mroots[1].first || 1.74<=mroots[1].first){
returnvalue-=256;
std::cerr<<"error 9: the roots of (x^2-2)*((1/2)x^2-(3/2))"<<
" are wrong"<<std::endl;
std::cerr<<"error 9: the roots of (x^2-2)*(x^2-3) are wrong"<<
std::endl;
}
if(mroots[0].second!=1 && mroots[1].second!=1){
returnvalue-=512;
std::cerr<<"error 10: the multiplicities of the roots of"<<
" (x^2-2)*((1/2)x^2-(3/2)) are wrong"<<std::endl;
" (x^2-2)*(x^2-3) are wrong"<<std::endl;
}
typename AK::Number_of_solutions_1 nos_1 = ak.number_of_solutions_1_object();
if(nos_1(x*x*x-2)!=1){
returnvalue-=1024;
std::cerr<<"error 11: x^3-2 must have only one root"<<std::endl;
}
return returnvalue;
}
int main(){
typedef CGAL::Algebraic_kernel_rs_gmpq_d_1 AK_rational;
typedef CGAL::Polynomial<CGAL::Gmpq> Polynomial;
typedef CGAL::Gmpfr Bound;
// test and return the result
long result=0;
std::cerr<<"*** testing rational RS AK_1:";
result+=test_ak<AK_rational>();
std::cerr<<"*** result of the tests (should be 0): "<<result<<std::endl;
return result;
}
#else
int main(){
return 0;

6
Algebraic_kernel_d/test/Algebraic_kernel_d/Algebraic_kernel_rs_gmpz_d_1.cpp
View File

@ -166,12 +166,12 @@ int main(){
// test all and return the result
long result=0;
std::cerr<<"testing default RS AK_1:";
std::cerr<<"*** testing default RS AK_1:";
result+=test_ak<AK_default>();
std::cerr<<"testing RS2 AK_1:";
std::cerr<<"*** testing RS2 AK_1:";
result+=(2048*test_ak<AK_RS2>());
#ifdef CGAL_USE_RS3
std::cerr<<"testing RS2/RS3 k_AK_1:";
std::cerr<<"*** testing RS2/RS3 k_AK_1:";
result+=(4096*test_ak<AK_RS2_RS3>());
#endif // CGAL_USE_RS3
std::cerr<<"*** result of the tests (should be 0): "<<result<<std::endl;

2
Algebraic_kernel_d/test/Algebraic_kernel_d/CMakeLists.txt
View File

@ -65,7 +65,7 @@ create_single_source_cgal_program( "Curve_pair_analysis_2.cpp" )
create_single_source_cgal_program( "Descartes.cpp" )
create_single_source_cgal_program( "Real_embeddable_traits_extension.cpp" )
if(RS_FOUND)
#create_single_source_cgal_program( "Algebraic_kernel_rs_gmpq_d_1.cpp" )
create_single_source_cgal_program( "Algebraic_kernel_rs_gmpq_d_1.cpp" )
create_single_source_cgal_program( "Algebraic_kernel_rs_gmpz_d_1.cpp" )
else()
message(STATUS "NOTICE: Some tests require the RS library, and will not be compiled.")