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

486 lines
19 KiB
C++
Raw Blame History

// Copyright (c) 2006-2008 Inria Lorraine (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.
//
// 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: Luis Peñaranda <luis.penaranda@gmx.com>
#ifndef CGAL_RS_REFINE_1_H
#define CGAL_RS_REFINE_1_H
#include <gmp.h>
#include <mpfr.h>
#include <mpfi.h>
#include <CGAL/RS/polynomial_1.h>
#include <CGAL/RS/algebraic_1.h>
#include <CGAL/RS/polynomial_1_utils.h>
#include <CGAL/RS/dyadic.h>
#include <CGAL/RS/sign_1.h>
#include <CGAL/assertions.h>
namespace CGAL{
class Algebraic_1;
// bisects a's isolation interval in point p, returns a
// Comparison_result, which is the comparison between a and p
// precondition: p belongs to a's isolation interval
template <class _Gcd_policy>
Comparison_result bisect_at(const Algebraic_1 &a,mpfr_srcptr p){
typedef _Gcd_policy Gcd;
Sign sl,sp;
int round;
sp=RSSign::signat(sfpart_1<Gcd>()(a.pol()),p);
if(sp==ZERO){
// set a=[p,p]
mpfr_set_prec(&(a.mpfi()->left),mpfr_get_prec(p));
round=mpfr_set(&(a.mpfi()->left),p,GMP_RNDN);
CGAL_assertion(!round);
mpfr_set_prec(&(a.mpfi()->right),mpfr_get_prec(p));
round=mpfr_set(&(a.mpfi()->right),p,GMP_RNDN);
CGAL_assertion(!round);
a.set_lefteval(ZERO);
return EQUAL;
}
sl=a.lefteval();
if(sp==sl){
// set a=[p,a_right]
mpfr_set_prec(&(a.mpfi()->left),mpfr_get_prec(p));
round=mpfr_set(&(a.mpfi()->left),p,GMP_RNDN);
CGAL_assertion(!round);
return LARGER;
}else{
// set a=[a_left,p]
mpfr_set_prec(&(a.mpfi()->right),mpfr_get_prec(p));
round=mpfr_set(&(a.mpfi()->right),p,GMP_RNDN);
CGAL_assertion(!round);
return SMALLER;
}
}
// refines b until having either only one point in common with a or
// all points inside a; the result will be zero when all points of
// b result to be in a, negative when a<b and positive when a>b
// precondition: a and b overlap
template <class _Gcd_policy>
int bisect_at_endpoints(const Algebraic_1 &a,Algebraic_1 &b){
typedef _Gcd_policy Gcd;
CGAL_precondition(a.overlaps(b));
if(mpfr_cmp(&(a.mpfi()->left),&(b.mpfi()->left))>0){
Comparison_result refinement=
bisect_at<Gcd>(b,&(a.mpfi()->left));
if(refinement==EQUAL)
return 0;
if(refinement==SMALLER)
return 1;
}
if(mpfr_cmp(&(a.mpfi()->right),&(b.mpfi()->right))<0 &&
bisect_at<Gcd>(b,&(a.mpfi()->right))==LARGER)
return -1;
return 0;
}
// bisect n times an interval; this function returns the number of
// refinements made
template <class _Gcd_policy>
int bisect_n(const Algebraic_1 &a,unsigned long n=1){
typedef _Gcd_policy Gcd;
Sign sl,sc;
mp_prec_t pl,pc;
mpfr_t center;
unsigned long i;
int round;
sl=a.lefteval();
if(sl==ZERO)
return 0;
pl=mpfr_get_prec(&(a.mpfi()->left));
pc=mpfr_get_prec(&(a.mpfi()->right));
pc=(pl>pc?pl:pc)+(mp_prec_t)n;
mpfr_init2(center,pc);
round=mpfr_prec_round((mpfr_ptr)&(a.mpfi()->left),pc,GMP_RNDN);
CGAL_assertion(!round);
round=mpfr_prec_round((mpfr_ptr)&(a.mpfi()->right),pc,GMP_RNDN);
CGAL_assertion(!round);
for(i=0;i<n;++i){
round=mpfr_add(
center,
&(a.mpfi()->left),
&(a.mpfi()->right),
GMP_RNDN);
CGAL_assertion(!round);
round=mpfr_div_2ui(center,center,1,GMP_RNDN);
CGAL_assertion(!round);
sc=RSSign::signat(sfpart_1<Gcd>()(a.pol()),center);
if(sc==ZERO){ // we have a root
round=mpfr_set((mpfr_ptr)&(a.mpfi()->left),
center,
GMP_RNDN);
CGAL_assertion(!round);
mpfr_swap((mpfr_ptr)&(a.mpfi()->right),center);
a.set_lefteval(ZERO);
break;
}
if(sc==sl)
mpfr_swap((mpfr_ptr)&(a.mpfi()->left),center);
else
mpfr_swap((mpfr_ptr)&(a.mpfi()->right),center);
}
mpfr_clear(center);
return i;
}
// The same bisection function, but with dyadic numbers. The difference
// is that dyadics will use the exact amount of bits needed, without
// allocating a big amount of memory. Note that the dyadic numbers are
// implemented as mpfrs, this implies that no conversion is needed.
template <class _Gcd_policy>
int bisect_n_dyadic(const Algebraic_1 &a,unsigned long n=1){
typedef _Gcd_policy Gcd;
Sign sl,sc;
CGALRS_dyadic_t center;
unsigned long i;
sl=a.lefteval();
if(sl==ZERO)
return 0;
CGALRS_dyadic_init(center);
for(i=0;i<n;++i){
CGALRS_dyadic_midpoint(center,a.left(),a.right());
sc=RSSign::signat(sfpart_1<Gcd>()(a.pol()),center);
if(sc==ZERO){ // we have a root
CGALRS_dyadic_set((CGALRS_dyadic_ptr)a.left(),center);
CGALRS_dyadic_swap((CGALRS_dyadic_ptr)a.right(),center);
a.set_lefteval(ZERO);
break;
}
if(sc==sl)
CGALRS_dyadic_swap((CGALRS_dyadic_ptr)a.left(),center);
else
CGALRS_dyadic_swap((CGALRS_dyadic_ptr)a.right(),center);
}
CGALRS_dyadic_clear(center);
return i;
}
// refine an interval, by bisecting it, until having a size smaller than 2^(-s)
template <class _Gcd_policy>
int bisect(const Algebraic_1 &a,mp_exp_t s){
typedef _Gcd_policy Gcd;
long ed;
mpfr_t d; // interval size
mpfr_init(d);
mpfr_sub(d,a.right(),a.left(),GMP_RNDU);
mpfr_log2(d,d,GMP_RNDU);
ed=1+mpfr_get_si(d,GMP_RNDU);
// we found ed such that the interval size is between 2^(ed-1) and 2^ed
mpfr_clear(d);
// following tests, dyadic bisection is a bit faster than mpfr one
return bisect_n_dyadic<Gcd>(a,s-ed);
}
// the four following functions are used to apply quadratic interval
// refinement (Abbott, 2006)
// calculate the value of kappa
//--------------------------------------------------
// unsigned long calc_kappa(const RS_polynomial_1 &p,
// CGALRS_dyadic_ptr f_lo,CGALRS_dyadic_ptr f_hi,
// unsigned n){
// unsigned long ret;
// CGALRS_dyadic_t temp;
// mpfr_t numerator,denominator;
// mpfr_inits2(MPFR_PREC_MIN,numerator,denominator,NULL);
// CGALRS_dyadic_init(temp);
//
// CGALRS_dyadic_sub(temp,f_lo,f_hi);
// CGALRS_dyadic_get_exactfr(denominator,temp);
// CGALRS_dyadic_mul_2exp(temp,f_lo,n);
// CGALRS_dyadic_get_exactfr(numerator,temp);
//
// mpfr_div(numerator,numerator,denominator,GMP_RNDN);
// ret=1+mpfr_get_ui(numerator,GMP_RNDN);
//
// CGALRS_dyadic_clear(temp);
// mpfr_clears(numerator,denominator,NULL);
// return ret;
// }
//--------------------------------------------------
unsigned long calc_kappa(const RS_polynomial_1 &p,
CGALRS_dyadic_ptr f_lo,CGALRS_dyadic_ptr f_hi,
unsigned n){
unsigned long ret;
//--------------------------------------------------
// CGALRS_dyadic_t numerator,denominator;
// CGALRS_dyadic_init(numerator);
// CGALRS_dyadic_init(denominator);
//--------------------------------------------------
CGALRS_dyadic_sub(f_hi,f_lo,f_hi);
CGALRS_dyadic_mul_2exp(f_lo,f_lo,n);
mpfr_div(f_lo,f_lo,f_hi,GMP_RNDN);
ret=mpfr_get_ui(f_lo,GMP_RNDN);
//--------------------------------------------------
// CGALRS_dyadic_clear(numerator);
// CGALRS_dyadic_clear(denominator);
//--------------------------------------------------
return ret;
}
// calculate kappa using doubles; faster but less accurate
//--------------------------------------------------
// inline unsigned calc_kappa(const RS_polynomial_1 &p,
// double f_lo,double f_hi,
// unsigned n){
// return(1+(int)nearbyint(f_lo*pow(2,n)/(f_lo-f_hi)));
// }
//--------------------------------------------------
// returns 0 for success, -1 for failure and 1 if it finds the root;
// the "N" of the paper is represented here as 2^n
int refine_interval_by_factor(CGALRS_dyadic_ptr x_lo,CGALRS_dyadic_ptr x_hi,
CGALRS_dyadic_ptr f_lo,CGALRS_dyadic_ptr f_hi,
const RS_polynomial_1 &p,unsigned n){
Sign sl=RSSign::signat(p,x_lo);
unsigned long kappa;
CGALRS_dyadic_t xkappa,xside,w;
Sign s1,s2;
kappa=calc_kappa(p,f_lo,f_hi,n);
//kappa=calc_kappa(p,
// CGALRS_dyadic_get_d(f_lo),
// CGALRS_dyadic_get_d(f_hi),
// n);
if(n==2){ // N==4: bisect twice
unsigned b[2],inter;
// we will use first xkappa as bisection point
CGALRS_dyadic_init(xkappa);
for(int i=0;i<2;++i){
CGALRS_dyadic_midpoint(xkappa,x_lo,x_hi);
if((s1=RSSign::signat(p,xkappa))==ZERO){ //exact
CGALRS_dyadic_set(x_lo,xkappa);
CGALRS_dyadic_swap(x_hi,xkappa);
CGALRS_dyadic_clear(xkappa);
return 1;
}
if(sl==s1){ // we take the right half
b[i]=2-i;
CGALRS_dyadic_swap(x_lo,xkappa);
}else{ // we take the left half
b[i]=0;
CGALRS_dyadic_swap(x_hi,xkappa);
}
}
CGALRS_dyadic_clear(xkappa);
inter=b[0]+b[1];
switch(inter){
case 0:p.inexact_eval_mpfr(f_hi,x_hi);break;
case 3:p.inexact_eval_mpfr(f_lo,x_lo);break;
default:p.inexact_eval_mpfr(f_hi,x_hi);
p.inexact_eval_mpfr(f_lo,x_lo);
break;
}
if(inter+1==kappa)
return 0; // success
else
return -2; // failure (but we refined anyway)
}else{ // N!=4
// calculate w, the width of every interval
CGALRS_dyadic_init(w);
CGALRS_dyadic_sub(w,x_hi,x_lo);
CGALRS_dyadic_div_2exp(w,w,n);
// calculate x[kappa]
CGALRS_dyadic_init_set(xkappa,x_lo);
CGALRS_dyadic_addmul_ui(xkappa,w,kappa);
if((s1=RSSign::signat(p,xkappa))==ZERO){
CGALRS_dyadic_set(x_lo,xkappa);
CGALRS_dyadic_swap(x_hi,xkappa);
CGALRS_dyadic_clear(xkappa);
CGALRS_dyadic_clear(w);
return 1; // exact root
}
if(s1==sl){
// calculate x[kappa+1]
CGALRS_dyadic_init(xside);
CGALRS_dyadic_add(xside,xkappa,w);
if((s2=RSSign::signat(p,xside))==ZERO){
CGALRS_dyadic_set(x_lo,xside);
CGALRS_dyadic_swap(x_hi,xside);
CGALRS_dyadic_clear(xkappa);
CGALRS_dyadic_clear(w);
CGALRS_dyadic_clear(xside);
return 1; // exact root
}
if(s2==sl){
CGALRS_dyadic_clear(xkappa);
CGALRS_dyadic_clear(w);
CGALRS_dyadic_clear(xside);
return -1; // failure
}else{ // s2!=sl && s2!=0
CGALRS_dyadic_set(x_lo,xkappa);
CGALRS_dyadic_swap(x_hi,xside);
p.inexact_eval_mpfr(f_lo,x_lo);
p.inexact_eval_mpfr(f_hi,x_hi);
CGALRS_dyadic_clear(xkappa);
CGALRS_dyadic_clear(w);
CGALRS_dyadic_clear(xside);
return 0; // success
}
}else{ // s1!=sl && s1!=ZERO
// we calculate x[kappa-1];
CGALRS_dyadic_init(xside);
CGALRS_dyadic_sub(xside,xkappa,w);
if((s2=RSSign::signat(p,xside))==ZERO){
CGALRS_dyadic_set(x_lo,xside);
CGALRS_dyadic_swap(x_hi,xside);
CGALRS_dyadic_clear(xkappa);
CGALRS_dyadic_clear(w);
CGALRS_dyadic_clear(xside);
return 1; // exact root
}
if(s2==sl){
CGALRS_dyadic_swap(x_lo,xside);
CGALRS_dyadic_swap(x_hi,xkappa);
p.inexact_eval_mpfr(f_lo,x_lo);
p.inexact_eval_mpfr(f_hi,x_hi);
CGALRS_dyadic_clear(xkappa);
CGALRS_dyadic_clear(w);
CGALRS_dyadic_clear(xside);
return 0; // success
}else{ // s2!=sl && s2!=ZERO
CGALRS_dyadic_clear(xkappa);
CGALRS_dyadic_clear(w);
CGALRS_dyadic_clear(xside);
return -1; // failure
}
}
}
CGAL_assertion_msg(false,"never reached");
return 2;
}
// applies qir a given number of times
template <class _Gcd_policy>
int qir_n(const Algebraic_1 &a,unsigned t=1){
typedef _Gcd_policy Gcd;
unsigned count=0;
unsigned n=2; // N=2^n
CGALRS_dyadic_t x_lo,x_hi,f_lo,f_hi; // endpoints and their evaluations
if(a.is_point())
return 0;
CGALRS_dyadic_init_set_fr(x_lo,a.left());
CGALRS_dyadic_init_set_fr(x_hi,a.right());
CGALRS_dyadic_init(f_lo);
CGALRS_dyadic_init(f_hi);
sfpart_1<Gcd>()(a.pol()).eval_dyadic(f_lo,x_lo);
sfpart_1<Gcd>()(a.pol()).eval_dyadic(f_hi,x_hi);
if(!CGALRS_dyadic_sgn(f_lo)){
CGALRS_dyadic_get_exactfr(&(a.mpfi()->right),x_lo);
return 0;
}
if(!CGALRS_dyadic_sgn(f_hi)){
CGALRS_dyadic_get_exactfr(&(a.mpfi()->left),x_hi);
return 0;
}
CGAL_assertion(CGALRS_dyadic_sgn(f_lo)!=CGALRS_dyadic_sgn(f_hi));
while(t>count){
switch(refine_interval_by_factor(x_lo,x_hi,f_lo,f_hi,
sfpart_1<Gcd>()(a.pol()),n))
{
case 0:n*=2;++count;break; // N=N^2
case -1:if(n>2)n/=2;break; // N=sqrt(N)
case -2:++count;break;
default:++count;goto endloop; // we found the root
}
}
endloop:
CGALRS_dyadic_get_exactfr(&(a.mpfi()->left),x_lo);
CGALRS_dyadic_get_exactfr(&(a.mpfi()->right),x_hi);
CGALRS_dyadic_clear(x_lo);
CGALRS_dyadic_clear(x_hi);
CGALRS_dyadic_clear(f_lo);
CGALRS_dyadic_clear(f_hi);
a.set_lefteval(RSSign::signat(sfpart_1<Gcd>()(a.pol()),a.left()));
return count;
}
// applies qir until having an interval of size less than 2^(-t)
template <class _Gcd_policy>
int qir(const Algebraic_1 &a,mp_exp_t t){
typedef _Gcd_policy Gcd;
int count=0;
long ed;
unsigned n=2; // N=2^n
CGALRS_dyadic_t d,f_lo,f_hi;
if(a.is_point())
return 0;
CGALRS_dyadic_init(f_lo);
CGALRS_dyadic_init(f_hi);
sfpart_1<Gcd>()
(a.pol()).eval_dyadic(f_lo,(CGALRS_dyadic_ptr)(a.left()));
sfpart_1<Gcd>()
(a.pol()).eval_dyadic(f_hi,(CGALRS_dyadic_ptr)(a.right()));
// we make sure we have a nice interval
if(!CGALRS_dyadic_sgn(f_lo)){
CGALRS_dyadic_set((CGALRS_dyadic_ptr)a.right(),
(CGALRS_dyadic_srcptr)a.left());
return 0;
}
if(!CGALRS_dyadic_sgn(f_hi)){
CGALRS_dyadic_set((CGALRS_dyadic_ptr)a.left(),
(CGALRS_dyadic_srcptr)a.right());
return 0;
}
CGAL_assertion(CGALRS_dyadic_sgn(f_lo)!=CGALRS_dyadic_sgn(f_hi));
// calculate ed, such that 2^(ed-1)<diam<2^ed
CGALRS_dyadic_init(d);
CGALRS_dyadic_sub(d,
(CGALRS_dyadic_ptr)a.right(),
(CGALRS_dyadic_ptr)a.left());
mpfr_log2((mpfr_ptr)d,(mpfr_ptr)d,GMP_RNDU);
ed=mpfr_get_si((mpfr_ptr)d,GMP_RNDU);
while(ed<t){
++count;
switch(refine_interval_by_factor(
(CGALRS_dyadic_ptr)a.left(),
(CGALRS_dyadic_ptr)a.right(),
f_lo,
f_hi,
sfpart_1<Gcd>()(a.pol()),
n)){
case 0:t-=n;n*=2;break; // N=N^2
case -1:if(n>2)n/=2;break; // N=sqrt(N)
case -2:t-=2;break;
default:t=0; // end the loop, we found the root
}
}
CGAL_assertion(CGALRS_dyadic_sgn(f_lo)!=CGALRS_dyadic_sgn(f_hi));
a.set_lefteval(
CGALRS_dyadic_sgn(f_lo)==0?
ZERO:
(CGALRS_dyadic_sgn(f_lo)<0?NEGATIVE:POSITIVE));
CGALRS_dyadic_clear(d);
CGALRS_dyadic_clear(f_lo);
CGALRS_dyadic_clear(f_hi);
return count;
}
} // namespace CGAL
#endif // CGAL_RS_REFINE_1_H
Reference in New Issue View Git Blame Copy Permalink