mirror of https://github.com/CGAL/cgal
Horner's method to evaluate poly_1's and easier refine_and_compare algorithm
This commit is contained in:
Luis Peñaranda
parent
ebf91a88f9
commit
8ed87d7a05
|
|
@ -70,7 +70,7 @@ class Rational_polynomial_1 {
|
|||
void set_root(const Algebraic_1&);*/
|
||||
//CGAL::Algebraic_1 eval_alg(const CGAL::Algebraic_1&)const;
|
||||
CGAL::Gmpz eval (const CGAL::Gmpz &) const;
|
||||
void eval_mpfr(mpfr_ptr,mpfr_srcptr,mp_prec_t)const;
|
||||
void eval_mpfr(mpfr_ptr,mpfr_srcptr)const;
|
||||
void eval_mpfi(mpfi_ptr,mpfi_srcptr)const;
|
||||
Rational_polynomial_1 derive () const;
|
||||
std::ostream& show (std::ostream &) const;
|
||||
|
|
|
|||
|
|
@ -153,20 +153,49 @@ CGAL::Gmpz Rational_polynomial_1::eval(const CGAL::Gmpz &x)const{
|
|||
return ret;
|
||||
};
|
||||
|
||||
void Rational_polynomial_1::eval_mpfr
|
||||
(mpfr_ptr result,mpfr_srcptr x,mp_prec_t prec)const{
|
||||
mpfr_t x_pow,temp;
|
||||
mp_prec_t prec_r=mpfr_get_prec(result);
|
||||
mpfr_inits2(prec<prec_r?prec_r:prec,x_pow,temp,NULL);
|
||||
mpfr_set_ui(x_pow,1,GMP_RNDN);
|
||||
mpfr_set_ui(result,0,GMP_RNDN);
|
||||
for(int i=0;i<=degree;++i){ // mpfr_fma?
|
||||
mpfr_mul_z(temp,x_pow,coef[i],GMP_RNDN);
|
||||
mpfr_add(result,temp,result,GMP_RNDN);
|
||||
mpfr_mul(x_pow,x_pow,x,GMP_RNDN);
|
||||
// This implementation uses the Horner's method. The precision paramenter
|
||||
// is not used anymore because calculations are exact (this function was
|
||||
// based on the signat in solve_1 but results are after calculations
|
||||
// converted to mpfr).
|
||||
void Rational_polynomial_1::eval_mpfr(mpfr_ptr result,mpfr_srcptr xcoord)const{
|
||||
// we have numbers H0=h0*2^e0 and X=x*2^ex
|
||||
mpz_t h0,x,temp;
|
||||
mp_exp_t e0,ex;
|
||||
// we convert the evaluation point (mpfr) to fit our needs
|
||||
mpz_init(x);
|
||||
ex=mpfr_get_z_exp(x,xcoord);
|
||||
// data from the polynomial
|
||||
mpz_t *c=get_coefs();
|
||||
int d=get_degree();
|
||||
// we start the algorithm by setting H0=c[d]*2^0
|
||||
mpz_init_set(h0,c[d]);
|
||||
e0=0;
|
||||
// the iteration is H0=H0*X+c[d-i]
|
||||
mpz_init(temp);
|
||||
for(int i=1;i<d+1;++i){
|
||||
mpz_mul(h0,h0,x);
|
||||
e0+=ex; // at this point, we did H0*=X
|
||||
// now, we have to add H0+c[d-i]
|
||||
if(e0>0){ // shift h0 and add
|
||||
mpz_mul_2exp(h0,h0,e0);
|
||||
mpz_add(h0,h0,c[d-i]);
|
||||
e0=0;
|
||||
}else{
|
||||
if(e0<0){ // shift c[d-i] and add
|
||||
mpz_mul_2exp(temp,c[d-i],-e0);
|
||||
mpz_add(h0,h0,temp);
|
||||
e0=0;
|
||||
}else // we are lucky, e0=0
|
||||
mpz_add(h0,h0,c[d-i]);
|
||||
}
|
||||
// at this point, H0 is the evaluation of the polynomial
|
||||
}
|
||||
mpfr_clears(x_pow,temp,NULL);
|
||||
return;
|
||||
mpfr_set_prec(result,mpz_sizeinbase(h0,2));
|
||||
mpfr_set_z(result,h0,GMP_RNDN);
|
||||
mpfr_mul_2si(result,result,e0,GMP_RNDN);
|
||||
mpz_clear(h0);
|
||||
mpz_clear(x);
|
||||
mpz_clear(temp);
|
||||
};
|
||||
|
||||
// I think RS should do this
|
||||
|
|
|
|||
|
|
@ -299,140 +299,15 @@ int get_root (mpfi_ptr x, int n) {
|
|||
return 0; // false: we couldn't copy the root
|
||||
}
|
||||
|
||||
int refine_1_rs(Algebraic_1 &a){
|
||||
rs_reset_all ();
|
||||
create_rs_upoly (a.pol().get_coefs (), a.pol().get_degree (),
|
||||
rs_get_default_up ());
|
||||
int newprec = a.rsprec() * CGAL_RS_PREC_FACTOR;
|
||||
a.set_rsprec (newprec);
|
||||
set_rs_precisol (newprec);
|
||||
set_rs_verbose (CGAL_RS_VERB);
|
||||
rs_run_algo ("UISOLE");
|
||||
return get_root (a.mpfi(), a.nr());
|
||||
}
|
||||
|
||||
// TODO: rewrite this awful function
|
||||
// TODO: test cases where RS is called
|
||||
Comparison_result refine_and_compare_1(Algebraic_1 &r1,Algebraic_1 &r2){
|
||||
mpfr_t left1,right1,center1,evall1,evalc1,evalr1,
|
||||
left2,right2,center2,evall2,evalc2,evalr2;
|
||||
mp_prec_t prec1,prec2,prec;
|
||||
if((prec1=r1.rsprec())<(prec2=r2.rsprec()))
|
||||
prec=(prec2*=2);
|
||||
else
|
||||
prec=(prec1*=2);
|
||||
mp_prec_t local_prec=prec;
|
||||
mpfr_inits2(local_prec,left1,right1,center1,evall1,evalc1,evalr1,
|
||||
left2,right2,center2,evall2,evalc2,evalr2,NULL);
|
||||
bool flag1,flag2;
|
||||
int res1,res2;
|
||||
int prec_changes=1;
|
||||
do{
|
||||
// refine r1
|
||||
flag1=true;
|
||||
r1.get_endpoints(left1,right1);
|
||||
mpfi_get_fr(center1,r1.mpfi());
|
||||
r1.pol().eval_mpfr(evall1,left1,prec);
|
||||
r1.pol().eval_mpfr(evalc1,center1,prec);
|
||||
r1.pol().eval_mpfr(evalr1,right1,prec);
|
||||
int sl1=mpfr_sgn(evall1);
|
||||
int sr1=mpfr_sgn(evalr1);
|
||||
// if the following assertion fails, it means that the
|
||||
// precision of the mpfr's was not correctly calculated
|
||||
//std::cout<<"sl1="<<sl1<<", sr1="<<sr1<<std::endl;
|
||||
CGAL_assertion(sl1&&sr1&&(sl1!=sr1));
|
||||
int sc1=mpfr_sgn(evalc1);
|
||||
while(r1.rsprec()<(int)prec){
|
||||
if(!sc1){
|
||||
refine_1_rs(r1);
|
||||
r1.get_endpoints(left1,right1);
|
||||
flag1=false;
|
||||
break;
|
||||
}else{
|
||||
if(sl1==sc1){
|
||||
mpfr_set(left1,center1,GMP_RNDN);
|
||||
mpfr_set(evall1,evalc1,GMP_RNDN);
|
||||
sl1=sc1;
|
||||
}else{
|
||||
mpfr_set(right1,center1,GMP_RNDN);
|
||||
mpfr_set(evalr1,evalc1,GMP_RNDN);
|
||||
sr1=sc1;
|
||||
}
|
||||
mpfr_add(center1,left1,right1,GMP_RNDN);
|
||||
mpfr_div_ui(center1,center1,2,GMP_RNDN);
|
||||
r1.pol().eval_mpfr(evalc1,center1,prec);
|
||||
sc1=mpfr_sgn(evalc1);
|
||||
r1.set_rsprec(r1.rsprec()+1);
|
||||
}
|
||||
}
|
||||
// refine r2
|
||||
flag2=true;
|
||||
r2.get_endpoints(left2,right2);
|
||||
mpfi_get_fr(center2,r2.mpfi());
|
||||
r2.pol().eval_mpfr(evall2,left2,prec);
|
||||
r2.pol().eval_mpfr(evalc2,center2,prec);
|
||||
r2.pol().eval_mpfr(evalr2,right2,prec);
|
||||
int sl2=mpfr_sgn(evall2);
|
||||
int sr2=mpfr_sgn(evalr2);
|
||||
CGAL_assertion(sl2&&sr2&&(sl2!=sr2));
|
||||
int sc2=mpfr_sgn(evalc2);
|
||||
while(r2.rsprec()<(int)prec){
|
||||
if(!sc2){
|
||||
refine_1_rs(r2);
|
||||
r2.get_endpoints(left2,right2);
|
||||
flag2=false;
|
||||
break;
|
||||
}else{
|
||||
if(sl2==sc2){
|
||||
mpfr_set(left2,center2,GMP_RNDN);
|
||||
mpfr_set(evall2,evalc2,GMP_RNDN);
|
||||
sl2=sc2;
|
||||
}else{
|
||||
mpfr_set(right2,center2,GMP_RNDN);
|
||||
mpfr_set(evalr2,evalc2,GMP_RNDN);
|
||||
sr2=sc2;
|
||||
}
|
||||
mpfr_add(center2,left2,right2,GMP_RNDN);
|
||||
mpfr_div_ui(center2,center2,2,GMP_RNDN);
|
||||
r2.pol().eval_mpfr(evalc2,center2,prec);
|
||||
sc2=mpfr_sgn(evalc2);
|
||||
r2.set_rsprec(r2.rsprec()+1);
|
||||
}
|
||||
}
|
||||
res1=mpfr_less_p(right1,left2);
|
||||
res2=mpfr_less_p(right2,left1);
|
||||
if((prec*=2)>local_prec){
|
||||
++prec_changes;
|
||||
local_prec*=(prec_changes*prec_changes);
|
||||
// change al the precisions
|
||||
mpfr_set_prec(left1,local_prec);
|
||||
mpfr_set_prec(center1,local_prec);
|
||||
mpfr_set_prec(right1,local_prec);
|
||||
mpfr_set_prec(evall1,local_prec);
|
||||
mpfr_set_prec(evalc1,local_prec);
|
||||
mpfr_set_prec(evalr1,local_prec);
|
||||
mpfr_set_prec(left2,local_prec);
|
||||
mpfr_set_prec(center2,local_prec);
|
||||
mpfr_set_prec(right2,local_prec);
|
||||
mpfr_set_prec(evall2,local_prec);
|
||||
mpfr_set_prec(evalc2,local_prec);
|
||||
mpfr_set_prec(evalr2,local_prec);
|
||||
}
|
||||
}while(!res1&&!res2);
|
||||
if(flag1)
|
||||
mpfi_interv_fr(r1.mpfi(),left1,right1);
|
||||
if(flag2)
|
||||
mpfi_interv_fr(r2.mpfi(),left2,right2);
|
||||
mpfr_clears(left1,right1,center1,evall1,evalc1,evalr1,
|
||||
left2,right2,center2,evall2,evalc2,evalr2,NULL);
|
||||
if(res1)
|
||||
return SMALLER;
|
||||
if(res2)
|
||||
return LARGER;
|
||||
CGAL_assertion_msg(false,"this point should never be reached");
|
||||
return EQUAL;
|
||||
refine(r1,5);
|
||||
refine(r2,5);
|
||||
}while(r1.overlaps(r2));
|
||||
return(r1<r2?SMALLER:LARGER);
|
||||
}
|
||||
|
||||
// TODO: rewrite using gcd
|
||||
Comparison_result compare_1(Algebraic_1 &r1,Algebraic_1 &r2){
|
||||
// there is a strange RS result where (r1<r2)&&(r2<r1)&&(r1!=r2)
|
||||
try{
|
||||
|
|
|
|||
Loading…
Reference in New Issue