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

removed old files

This commit is contained in:
Luis Peñaranda 2013-11-19 16:08:29 -02:00
parent cbdca2c35d
commit 8a58557f40
38 changed files with 0 additions and 6252 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

101
Algebraic_kernel_d/include/CGAL/RS/Algebraic_kernel_rs_1.h
View File

@ -1,101 +0,0 @@
// Copyright (c) 2006-2010 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_ALGEBRAIC_KERNEL_RS_1
#define CGAL_RS_ALGEBRAIC_KERNEL_RS_1
#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Polynomial_type_generator.h>
#include <CGAL/RS/functors_1.h>
template <class CoefficientType,
class GcdFunction=
#ifdef CGAL_RS_USE_UGCD
CGAL::Modgcd_1
#else
CGAL::Rsgcd_1
#endif
>
struct Algebraic_kernel_rs_1{
typedef CoefficientType Coefficient;
typedef GcdFunction Gcd;
typedef CGAL::Polynomial<Coefficient> Polynomial_1;
typedef RSFunctors::Algebraic Algebraic_real_1;
typedef RSFunctors::Bound Bound;
typedef RSFunctors::Multiplicity Multiplicity_type;
// constructor: we must initialize RS just a time, so this is a
// good time to do it
Algebraic_kernel_rs_1(){CGAL::init_solver();};
~Algebraic_kernel_rs_1(){CGAL::reset_solver();};
typedef RSFunctors::Construct_alg_1<Polynomial_1,Coefficient,Gcd>
Construct_algebraic_real_1;
typedef RSFunctors::Compute_polynomial_1<Polynomial_1>
Compute_polynomial_1;
typedef RSFunctors::Isolate_1<Polynomial_1,Gcd> Isolate_1;
typedef RSFunctors::Is_square_free_1<Polynomial_1,Gcd>
Is_square_free_1;
typedef RSFunctors::Make_square_free_1<Polynomial_1,Gcd>
Make_square_free_1;
typedef RSFunctors::Square_free_factorize_1<Polynomial_1,Gcd>
Square_free_factorize_1;
typedef RSFunctors::Is_coprime_1<Polynomial_1,Gcd>
Is_coprime_1;
typedef RSFunctors::Make_coprime_1<Polynomial_1,Gcd>
Make_coprime_1;
typedef RSFunctors::Solve_1<Polynomial_1,Gcd> Solve_1;
typedef RSFunctors::Number_of_solutions_1<Polynomial_1>
Number_of_solutions_1;
typedef RSFunctors::Sign_at_1<Polynomial_1,Gcd> Sign_at_1;
typedef RSFunctors::Is_zero_at_1<Polynomial_1,Gcd>
Is_zero_at_1;
typedef RSFunctors::Compare_1<Gcd> Compare_1;
typedef RSFunctors::Bound_between_1<Gcd> Bound_between_1;
typedef RSFunctors::Approximate_absolute_1 Approximate_absolute_1;
typedef RSFunctors::Approximate_relative_1 Approximate_relative_1;
#define CGALRS_CREATE_FUNCTION_OBJECT(T,N) \
T N##_object()const{return T();}
CGALRS_CREATE_FUNCTION_OBJECT(Construct_algebraic_real_1,
construct_algebraic_real_1)
CGALRS_CREATE_FUNCTION_OBJECT(Compute_polynomial_1,compute_polynomial_1)
CGALRS_CREATE_FUNCTION_OBJECT(Isolate_1,isolate_1)
CGALRS_CREATE_FUNCTION_OBJECT(Is_square_free_1,is_square_free_1)
CGALRS_CREATE_FUNCTION_OBJECT(Make_square_free_1,make_square_free_1)
CGALRS_CREATE_FUNCTION_OBJECT(Square_free_factorize_1,
square_free_factorize_1)
CGALRS_CREATE_FUNCTION_OBJECT(Is_coprime_1,is_coprime_1)
CGALRS_CREATE_FUNCTION_OBJECT(Make_coprime_1,make_coprime_1)
CGALRS_CREATE_FUNCTION_OBJECT(Solve_1,solve_1)
CGALRS_CREATE_FUNCTION_OBJECT(Number_of_solutions_1,
number_of_solutions_1)
CGALRS_CREATE_FUNCTION_OBJECT(Sign_at_1,sign_at_1)
CGALRS_CREATE_FUNCTION_OBJECT(Is_zero_at_1,is_zero_at_1)
CGALRS_CREATE_FUNCTION_OBJECT(Compare_1,compare_1)
CGALRS_CREATE_FUNCTION_OBJECT(Bound_between_1,bound_between_1)
CGALRS_CREATE_FUNCTION_OBJECT(Approximate_absolute_1,
approximate_absolute_1)
CGALRS_CREATE_FUNCTION_OBJECT(Approximate_relative_1,
approximate_relative_1)
#undef CGALRS_CREATE_FUNCTION_OBJECT
}; // Algebraic_kernel_d_1_RS
#endif // CGAL_RS_ALGEBRAIC_KERNEL_RS_1

60
Algebraic_kernel_d/include/CGAL/RS/algebraic_1_comparisons.h
View File

@ -1,60 +0,0 @@
// 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_ALGEBRAIC_1_COMPARISONS_H
#define CGAL_RS_ALGEBRAIC_1_COMPARISONS_H
#include <boost/config.hpp>
#include <CGAL/RS/compare_1.h>
#include <CGAL/RS/polynomial_1_utils.h>
namespace CGAL{
#if CGAL_USE_RS3
inline
bool operator<(const Algebraic_1 &n1,const Algebraic_1 &n2){
typedef CGAL::Rsgcd_1 Gcd;
return(CGAL::RS_COMPARE::compare_1<Gcd>(n1,n2)==SMALLER);
}
#else
inline
bool operator<(const Algebraic_1 &n1,const Algebraic_1 &n2){
typedef CGAL::Cgalgcd_1 Gcd;
return(CGAL::RS_COMPARE::compare_1<Gcd>(n1,n2)==SMALLER);
}
#endif
inline
bool operator==(const Algebraic_1 &n1,const Algebraic_1 &n2){
typedef CGAL::Rsgcd_1 Gcd;
return(CGAL::RS_COMPARE::compare_1<Gcd>(n1,n2)==EQUAL);
}
inline
Algebraic_1 min BOOST_PREVENT_MACRO_SUBSTITUTION (const Algebraic_1 &a,const Algebraic_1 &b){
return (a<b?a:b);
}
inline
Algebraic_1 max BOOST_PREVENT_MACRO_SUBSTITUTION (const Algebraic_1 &a,const Algebraic_1 &b){
return (a>b?a:b);
}
} // namespace CGAL
#endif // CGAL_RS_ALGEBRAIC_1_COMPARISONS_H

227
Algebraic_kernel_d/include/CGAL/RS/algebraic_1_constructors.h
View File

@ -1,227 +0,0 @@
// 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_ALGEBRAIC_1_CONSTRUCTORS_H
#define CGAL_RS_ALGEBRAIC_1_CONSTRUCTORS_H
#include <CGAL/RS/polynomial_1.h>
#include <CGAL/RS/polynomial_1_utils.h>
#include <CGAL/RS/sign_1.h>
namespace CGAL{
inline
Algebraic_1::Algebraic_1(const Algebraic_1 &i,mpfr_prec_t pl,mpfr_prec_t pr){
if(pl>pr)
mpfi_init2(mpfi(),pl);
else
mpfi_init2(mpfi(),pr);
set_mpfi(i.mpfi());
set_pol(i.pol());
set_nr(i.nr());
set_mult(i.mult());
set_lefteval(i.lefteval());
}
inline
Algebraic_1::Algebraic_1(){
mpfi_init(mpfi());
}
inline
Algebraic_1::Algebraic_1(int i){
mpq_t temp;
mpfi_init(mpfi());
mpfi_set_si(mpfi(),(long int)i);
mpq_init(temp);
mpq_set_si(temp,(long int)i,1);
RS_polynomial_1 *p=new RS_polynomial_1(temp);
mpq_clear(temp);
set_pol(*p);
set_nr(0);
set_lefteval(CGAL::NEGATIVE);
}
inline
Algebraic_1::Algebraic_1(unsigned i){
mpq_t temp;
mpfi_init(mpfi());
mpfi_set_ui(mpfi(),i);
mpq_init(temp);
mpq_set_ui(temp,(unsigned)i,1);
RS_polynomial_1 *p=new RS_polynomial_1(temp);
mpq_clear(temp);
set_pol(*p);
set_nr(0);
set_lefteval(CGAL::NEGATIVE);
}
inline
Algebraic_1::Algebraic_1(long i){
mpq_t temp;
mpfi_init(mpfi());
mpfi_set_si(mpfi(),i);
mpq_init(temp);
mpq_set_si(temp,i,1);
RS_polynomial_1 *p=new RS_polynomial_1(temp);
mpq_clear(temp);
set_pol(*p);
set_nr(0);
set_lefteval(CGAL::NEGATIVE);
}
inline
Algebraic_1::Algebraic_1(unsigned long i){
mpq_t temp;
mpfi_init(mpfi());
mpfi_set_ui(mpfi(),i);
mpq_init(temp);
mpq_set_ui(temp,i,1);
RS_polynomial_1 *p=new RS_polynomial_1(temp);
mpq_clear(temp);
set_pol(*p);
set_nr(0);
set_lefteval(CGAL::NEGATIVE);
}
inline
Algebraic_1::Algebraic_1(double d){
mpq_t temp;
mpfi_init(mpfi());
mpfi_set_d(mpfi(),d);
mpq_init(temp);
mpq_set_d(temp,d);
RS_polynomial_1 *p=new RS_polynomial_1(temp);
mpq_clear(temp);
set_pol(*p);
set_nr(0);
set_lefteval(CGAL::NEGATIVE);
}
inline
Algebraic_1::Algebraic_1(mpz_srcptr z){
mpq_t temp;
mpfi_init(mpfi());
mpfi_set_z(mpfi(),z);
mpq_init(temp);
mpq_set_z(temp,z);
RS_polynomial_1 *p=new RS_polynomial_1(temp);
mpq_clear(temp);
set_pol(*p);
set_nr(0);
set_lefteval(CGAL::NEGATIVE);
}
inline
Algebraic_1::Algebraic_1(mpq_srcptr q){
mpfi_init(mpfi());
mpfi_set_q(mpfi(),q);
RS_polynomial_1 *p=new RS_polynomial_1(q);
set_pol(*p);
set_nr(0);
set_lefteval(CGAL::NEGATIVE);
}
inline
Algebraic_1::Algebraic_1(mpfr_srcptr src){
Gmpfr r(src);
mpfi_init2(mpfi(),r.get_precision());
mpfi_set_fr(mpfi(),r.fr());
CGAL_assertion(mpfr_equal_p(r.fr(),&mpfi()->left)!=0);
CGAL_assertion(mpfr_equal_p(r.fr(),&mpfi()->right)!=0);
RS_polynomial_1 *rsp=new RS_polynomial_1(Gmpq(r).mpq());
set_pol(*rsp);
set_nr(0);
set_lefteval(CGAL::NEGATIVE);
}
inline
Algebraic_1::Algebraic_1(mpfi_srcptr i){
mpfi_init(mpfi());
set_mpfi(i);
}
inline
Algebraic_1::Algebraic_1(const Gmpz &z){
mpq_t temp;
mpfi_init(mpfi());
mpfi_set_z(mpfi(),z.mpz());
mpq_init(temp);
mpq_set_z(temp,z.mpz());
RS_polynomial_1 *p=new RS_polynomial_1(temp);
mpq_clear(temp);
set_pol(*p);
set_nr(0);
set_lefteval(CGAL::NEGATIVE);
}
inline
Algebraic_1::Algebraic_1(const Gmpq &q){
mpfi_init(mpfi());
mpfi_set_q(mpfi(),q.mpq());
RS_polynomial_1 *p=new RS_polynomial_1(q.mpq());
set_pol(*p);
set_nr(0);
set_lefteval(CGAL::NEGATIVE);
}
inline
Algebraic_1::Algebraic_1(const Gmpfr &r){
mpfi_init2(mpfi(),r.get_precision());
mpfi_set_fr(mpfi(),r.fr());
CGAL_assertion(mpfr_equal_p(r.fr(),&mpfi()->left)!=0);
CGAL_assertion(mpfr_equal_p(r.fr(),&mpfi()->right)!=0);
RS_polynomial_1 *rsp=new RS_polynomial_1(Gmpq(r).mpq());
set_pol(*rsp);
set_nr(0);
set_lefteval(CGAL::NEGATIVE);
}
// interesting constructor
inline
Algebraic_1::Algebraic_1(
mpfi_srcptr i,
const RS_polynomial_1 &p,
int n,
int m){
mpfi_init(mpfi());
set_mpfi_ptr(i);
set_pol(p);
set_nr(n);
set_mult(m);
// we don't evaluate in the sf part of p, since p is sf
// TODO: add assertion
set_lefteval(RSSign::signat(p,&i->left));
}
// another interesting constructor, where we have already calculated the
// left evaluation
inline
Algebraic_1::Algebraic_1(mpfi_srcptr i,const RS_polynomial_1 &p,int n,int m,
CGAL::Sign s){
mpfi_init(mpfi());
set_mpfi_ptr(i);
set_pol(p);
set_nr(n);
set_mult(m);
set_lefteval(s);
}
} // namespace CGAL
#endif // CGAL_RS_ALGEBRAIC_1_CONSTRUCTORS_H

213
Algebraic_kernel_d/include/CGAL/RS/algebraic_1_member.h
View File

@ -1,213 +0,0 @@
// 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_ALGEBRAIC_1_MEMBER_H
#define CGAL_RS_ALGEBRAIC_1_MEMBER_H
#include <CGAL/RS/refine_1_rs.h>
#include <CGAL/assertions.h>
namespace CGAL{
inline
mpfi_srcptr Algebraic_1::mpfi()const{
return Ptr()->interval;
}
inline
mpfi_ptr Algebraic_1::mpfi(){
return ptr()->interval;
}
inline
Gmpfi Algebraic_1::interval()const{
return Gmpfi(Ptr()->interval);
}
inline
Gmpfr Algebraic_1::inf()const{
return Gmpfr(&(mpfi()->left));
}
inline
Gmpfr Algebraic_1::sup()const{
return Gmpfr(&(mpfi()->right));
}
inline
const RS_polynomial_1& Algebraic_1::pol()const{
return *(Ptr()->polynomial);
}
inline
int Algebraic_1::nr()const{
return ptr()->numroot;
}
inline
int Algebraic_1::mult()const{
return ptr()->multiplicity;
}
inline
void Algebraic_1::set_mpfi_ptr(mpfi_srcptr x){
// *mpfi()=*x;
// mpfi_set(mpfi(),x);
set_mpfi(x);
}
inline
void Algebraic_1::clear_pol(){
ptr()->polynomial=NULL;
}
inline
void Algebraic_1::set_pol(const RS_polynomial_1 &p){
ptr()->polynomial=const_cast<RS_polynomial_1*>(&p);
}
inline
void Algebraic_1::set_nr(int n){
ptr()->numroot=n;
}
inline
void Algebraic_1::set_mult(int m){
ptr()->multiplicity=m;
}
inline
void Algebraic_1::set_prec(mp_prec_t p){
mpfi_round_prec(mpfi(),p);
}
inline
void Algebraic_1::set_lefteval(Sign s)const{
Ptr()->lefteval=s;
}
inline
mp_prec_t Algebraic_1::get_prec()const{
return mpfi_get_prec(mpfi());
}
inline
mpfr_srcptr Algebraic_1::left()const{
return &(mpfi()->left);
}
inline
mpfr_srcptr Algebraic_1::right()const{
return &(mpfi()->right);
}
inline
Sign Algebraic_1::lefteval()const{
return ptr()->lefteval;
}
inline
bool Algebraic_1::is_consistent()const{
return(&pol()==NULL?false:true);
}
inline
bool Algebraic_1::is_point()const{
return(mpfr_equal_p(&(mpfi()->left),&(mpfi()->right))!=0);
}
inline
bool Algebraic_1::contains(int n)const{
return(mpfr_cmp_si(&(mpfi()->left),n)<=0 &&
mpfr_cmp_si(&(mpfi()->right),n)>=0);
}
inline
bool Algebraic_1::contains(mpfr_srcptr n)const{
return(mpfr_lessequal_p(&(mpfi()->left),n) &&
mpfr_greaterequal_p(&(mpfi()->right),n));
}
inline
bool Algebraic_1::contains(const Gmpz &n)const{
return(mpfr_cmp_z(&(mpfi()->left),n.mpz())<=0 &&
mpfr_cmp_z(&(mpfi()->right),n.mpz())>=0);
}
inline
std::pair<double,double> Algebraic_1::to_interval()const{
return std::make_pair(
mpfr_get_d(left(),GMP_RNDD),
mpfr_get_d(right(),GMP_RNDU));
}
inline
void Algebraic_1::set_mpfi(mpfi_srcptr x){
mp_prec_t xp;
xp=mpfi_get_prec(x);
if(xp>mpfr_get_prec(left()) || xp>mpfr_get_prec(right()))
mpfi_set_prec(mpfi(),xp);
mpfi_set(mpfi(),x);
}
inline
bool Algebraic_1::overlaps(const Algebraic_1&a)const{
if(mpfr_lessequal_p(left(),a.left()))
return (mpfr_lessequal_p(a.left(),right())!=0);
else
return (mpfr_lessequal_p(left(),a.right())!=0);
}
inline
bool Algebraic_1::is_valid()const{
return (mpfi_nan_p(mpfi())==0);
}
inline
bool Algebraic_1::is_finite()const{
return (mpfi_bounded_p(mpfi())!=0);
}
//template <class GcdPolicy>
inline
double Algebraic_1::to_double()const{
//typedef GcdPolicy Gcd;
//while(mpfr_get_d(left(),GMP_RNDU)!=mpfr_get_d(right(),GMP_RNDU))
// bisect_n<Gcd>(*this,33);
RS3::refine_1(*this,100);
CGAL_assertion(mpfr_get_d(left(),GMP_RNDD)==
mpfr_get_d(right(),GMP_RNDD));
CGAL_assertion(mpfr_get_d(left(),GMP_RNDU)==
mpfr_get_d(right(),GMP_RNDU));
CGAL_assertion(mpfr_get_d(left(),GMP_RNDN)==
mpfr_get_d(right(),GMP_RNDN));
return mpfr_get_d(right(),GMP_RNDU);
}
inline
Algebraic_1 Algebraic_1::sqrt()const{
mpfi_t s;
mpfi_init(s);
mpfi_sqrt(s,mpfi());
Algebraic_1 ret(s);
return ret;
}
} // namespace CGAL
#endif // CGAL_RS_ALGEBRAIC_1_MEMBER_H

44
Algebraic_kernel_d/include/CGAL/RS/algebraic_1_operators.h
View File

@ -1,44 +0,0 @@
// 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_ALGEBRAIC_1_OPERATORS_H
#define CGAL_RS_ALGEBRAIC_1_OPERATORS_H
namespace CGAL{
inline
Algebraic_1 Algebraic_1::operator+()const{
return *this;
}
inline
Algebraic_1 Algebraic_1::operator-()const{
mpfi_t inv;
mpfi_init2(inv,mpfi_get_prec(mpfi()));
mpfi_neg(inv,mpfi());
Algebraic_1 *inverse=new Algebraic_1(inv,
pol().minusx(),
nr(),
mult(),
-lefteval());
return *inverse;
}
} // namespace CGAL
#endif // CGAL_RS_ALGEBRAIC_1_OPERATORS_H

114
Algebraic_kernel_d/include/CGAL/RS/algebraic_1_other.h
View File

@ -1,114 +0,0 @@
// 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_ALGEBRAIC_1_OTHER_H
#define CGAL_RS_ALGEBRAIC_1_OTHER_H
namespace CGAL{
inline
bool is_valid(const Algebraic_1 &n){
return n.is_valid();
}
inline
bool is_finite(const Algebraic_1 &n){
return n.is_finite();
}
//template <class GcdPolicy>
inline
double to_double(Algebraic_1 &n){
//typedef GcdPolicy Gcd;
//return n.to_double<Gcd>();
return n.to_double();
}
inline
std::pair<double,double>to_interval(const Algebraic_1 &n){
return n.to_interval();
}
inline
Algebraic_1 sqrt(const Algebraic_1 &ntval){
return ntval.sqrt();
}
inline
std::ostream& operator<<(std::ostream &os,const Algebraic_1 &a){
// (interval,polynomial,number_of_root,multiplicity,left_sign)
return os<<'('<<a.interval()<<','<<a.pol()<<','<<
a.nr()<<','<<a.mult()<<','<<a.lefteval()<<')';
}
inline
std::istream& operator>>(std::istream &is,Algebraic_1 &a){
Gmpfi i;
RS_polynomial_1 pol;
int nr,mult,eval;
std::istream::int_type c;
std::ios::fmtflags old_flags=is.flags();
is.unsetf(std::ios::skipws);
gmpz_eat_white_space(is);
c=is.get();
if(c!='(')
goto is_fail_ret;
gmpz_eat_white_space(is);
// read the Gmpfi
is>>i;
c=is.get();
if(c!=',')
goto is_fail_ret;
is>>pol;
gmpz_eat_white_space(is);
c=is.get();
if(c!=',')
goto is_fail_ret;
is>>nr;
gmpz_eat_white_space(is);
c=is.get();
if(c!=',')
goto is_fail_ret;
is>>mult;
gmpz_eat_white_space(is);
c=is.get();
if(c!=',')
goto is_fail_ret;
is>>eval;
gmpz_eat_white_space(is);
c=is.get();
if(c!=')')
goto is_fail_ret;
a=Algebraic_1(i.mpfi(), // interval
*(new RS_polynomial_1(pol)), // polynomial
nr, // number of root
mult, // multiplicity
//(mpfi_ptr)NULL, // previous root
//(mpfi_ptr)NULL, // next root
(CGAL::Sign)eval);// evaluation on the left bound
goto is_ret;
is_fail_ret:
is.setstate(std::ios_base::failbit);
is_ret:
is.flags(old_flags);
return is;
}
} // namespace CGAL
#endif // CGAL_RS_ALGEBRAIC_1_OTHER_H

88
Algebraic_kernel_d/include/CGAL/RS/algebraic_1_real_embeddable.h
View File

@ -1,88 +0,0 @@
// 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(s) : Michael Hemmer <hemmer@mpi-inf.mpg.de>
//
// ============================================================================
namespace CGAL {
template<>
class Real_embeddable_traits< Algebraic_1 >
: public INTERN_RET::Real_embeddable_traits_base< Algebraic_1 , CGAL::Tag_true > {
public:
typedef INTERN_RET::Real_embeddable_traits_base< Algebraic_1 , CGAL::Tag_true > Base;
typedef CGAL::Tag_true Is_real_embeddable;
typedef bool Boolean;
typedef CGAL::Sign Sign;
typedef CGAL::Comparison_result Comparison_result;
typedef Algebraic_1 Type;
typedef Base::Compare Compare; // todo: get a more efficient impl
class Sgn
: public std::unary_function< Type, CGAL::Sign > {
public:
CGAL::Sign operator()( const Type& a ) const {
return Compare()(a, Type(0));
}
};
class To_double
: public std::unary_function< Type, double > {
public:
double operator()(const Type& a) const {
return a.to_double();
}
};
class To_interval
: public std::unary_function< Type, std::pair<double, double> > {
public:
std::pair<double, double> operator()(const Type& a) const {
return a.to_interval();
}
};
class Is_zero
: public std::unary_function< Type, Boolean> {
public:
bool operator()(const Type& a) const {
return Sgn()(a) == CGAL::ZERO;
}
};
class Is_finite
:public std::unary_function< Type, Boolean> {
public:
bool operator()(const Type& ) const {
return true;
}
};
class Abs
:public std::unary_function< Type, Type> {
public:
Type operator()(const Type& a) const {
return Sgn()(a)==CGAL::NEGATIVE?-a:a;
}
};
};
} // namespace CGAL

79
Algebraic_kernel_d/include/CGAL/RS/basic.h
View File

@ -1,79 +0,0 @@
// Copyright (c) 2006-2010 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_BASIC_H
#define CGAL_RS_BASIC_H
#include <CGAL/basic.h>
#include <CGAL/enum.h>
#ifndef CGAL_RS_VERB
#ifdef CGAL_RS_DEBUG
#define CGAL_RS_VERB 2
#else
#define CGAL_RS_VERB 0
#endif
#endif
// the default precision of RS to calculate a root (precision is 2^n)
#ifndef CGAL_RS_DEF_PREC
#define CGAL_RS_DEF_PREC 0
#endif
// the minimum, used when calculating a sign
#ifndef CGAL_RS_MIN_PREC
#define CGAL_RS_MIN_PREC 0
#endif
#if ( defined CGAL_HAS_THREADS && !defined CGAL_RS_NO_TLS )
# ifdef _MSC_VER
# ifdef _WINDLL
# error "Can't build CGAL_RS as thread safe."
# define CGALRS_THREAD_ATTR
# else
# define CGALRS_THREAD_ATTR __declspec(thread)
# endif
# else
# define CGALRS_THREAD_ATTR __thread
# endif
#else
# define CGALRS_THREAD_ATTR
#endif
namespace RS{
enum rs_sign{
RS_NEGATIVE= CGAL::NEGATIVE,
RS_ZERO= CGAL::ZERO,
RS_POSITIVE= CGAL::POSITIVE,
RS_UNKNOWN
};
// this function must only be called when s is not RS_UNKNOWN
inline CGAL::Sign convert_rs_sign(rs_sign s){
CGAL_precondition(s!=RS_UNKNOWN);
switch(s){
case RS_NEGATIVE:return CGAL::NEGATIVE;break;
case RS_POSITIVE:return CGAL::POSITIVE;break;
default:return CGAL::ZERO;
}
}
} // namespace RS
#endif // CGAL_RS_BASIC_H

90
Algebraic_kernel_d/include/CGAL/RS/compare_1.h
View File

@ -1,90 +0,0 @@
// 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_COMPARE_1_H
#define CGAL_RS_COMPARE_1_H
#include <mpfr.h>
#include <CGAL/RS/polynomial_1.h>
#include <CGAL/RS/algebraic_1.h>
#include <CGAL/RS/polynomial_1_utils.h>
#include <CGAL/RS/sign_1.h>
#include <CGAL/RS/refine_1_rs.h>
// Default sign function on this file.
#define CGALRS_SIGNAT(P,M) RSSign::signat(P,M)
namespace CGAL{
namespace RS_COMPARE{
// compare two algebraic numbers, knowing they are not equal
inline
Comparison_result
compare_1_unequal(const Algebraic_1 &r1,const Algebraic_1 &r2){
mp_prec_t prec1=r1.get_prec();
mp_prec_t prec2=r2.get_prec();
unsigned refsteps=prec1>prec2?prec1:prec2;
RS3::refine_1(r1,refsteps);
RS3::refine_1(r2,refsteps);
while(r1.overlaps(r2)){
refsteps*=2;
RS3::refine_1(r1,refsteps);
RS3::refine_1(r2,refsteps);
}
return(mpfr_less_p(r1.right(),r2.left())?SMALLER:LARGER);
}
template <class GcdPolicy>
Comparison_result
compare_1(const Algebraic_1 &r1,const Algebraic_1 &r2){
typedef GcdPolicy Gcd;
//if(r1.pol()==r2.pol())
// return(r1.nr()!=r2.nr()?(r1.nr()<r2.nr()?SMALLER:LARGER):EQUAL);
if(mpfr_lessequal_p(r1.left(),r2.left())){
if(mpfr_less_p(r1.right(),r2.left()))
return SMALLER;
}else{
if(mpfr_less_p(r2.right(),r1.left()))
return LARGER;
}
RS_polynomial_1 gcd=Gcd()(r1.pol(),r2.pol());
if(!gcd.get_degree())
return RS_COMPARE::compare_1_unequal(r1,r2);
Sign sleft,sright;
if(mpfr_greater_p(r1.left(),r2.left()))
sleft=CGALRS_SIGNAT(gcd,r1.left());
else
sleft=CGALRS_SIGNAT(gcd,r2.left());
if(sleft==ZERO)
return EQUAL;
if(mpfr_less_p(r1.right(),r2.right()))
sright=CGALRS_SIGNAT(gcd,r1.right());
else
sright=CGALRS_SIGNAT(gcd,r2.right());
if(sleft!=sright)
return EQUAL;
else
return RS_COMPARE::compare_1_unequal(r1,r2);
}
} // namespace RS_COMPARE
} // namespace CGAL
#undef CGALRS_SIGNAT
#endif // CGAL_RS_COMPARE_1_H

438
Algebraic_kernel_d/include/CGAL/RS/dyadic.h
View File

@ -1,438 +0,0 @@
// Copyright (c) 2007-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_DYADIC_H
#define CGAL_RS_DYADIC_H
#include <stdio.h>
#include <math.h>
#include <gmp.h>
#include <mpfr.h>
#include <CGAL/assertions.h>
// for c++, compile with -lgmpxx
#ifdef __cplusplus
#include <iostream>
#endif
#define CGALRS_dyadic_struct __mpfr_struct
#define CGALRS_dyadic_t mpfr_t
#define CGALRS_dyadic_ptr mpfr_ptr
#define CGALRS_dyadic_srcptr mpfr_srcptr
// some auxiliary defines
#define CGALRS_dyadic_set_prec(D,P) \
( mpfr_set_prec( (D), (P)>MPFR_PREC_MIN?(P):MPFR_PREC_MIN) )
#define CGALRS_dyadic_prec_round(D,P) \
( mpfr_prec_round( (D), (P)>MPFR_PREC_MIN?(P):MPFR_PREC_MIN, GMP_RNDN) )
#define CGALRS_dyadic_set_exp(D,E) \
( CGAL_assertion( (E) <= mpfr_get_emax() && \
(E) >= mpfr_get_emin() ) ,\
mpfr_set_exp(D,E) )
// init functions
#define CGALRS_dyadic_init(D) mpfr_init2(D,MPFR_PREC_MIN)
#define CGALRS_dyadic_init2(D,P) mpfr_init2(D,P)
#define CGALRS_dyadic_clear(D) mpfr_clear(D)
inline void CGALRS_dyadic_set(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op){
if(rop!=op){
CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op));
mpfr_set(rop,op,GMP_RNDN);
}
CGAL_assertion(mpfr_equal_p(rop,op)!=0);
}
inline void CGALRS_dyadic_set_z(CGALRS_dyadic_ptr rop,mpz_srcptr z){
size_t prec;
prec=mpz_sizeinbase(z,2)-(mpz_tstbit(z,0)?0:mpz_scan1(z,0));
CGALRS_dyadic_set_prec(rop,prec);
mpfr_set_z(rop,z,GMP_RNDN);
CGAL_assertion(!mpfr_cmp_z(rop,z));
}
inline void CGALRS_dyadic_set_si(CGALRS_dyadic_ptr rop,long s){
CGALRS_dyadic_set_prec(rop,sizeof(long));
mpfr_set_si(rop,s,GMP_RNDN);
CGAL_assertion(!mpfr_cmp_si(rop,s));
}
inline void CGALRS_dyadic_set_ui(CGALRS_dyadic_ptr rop,unsigned long u){
CGALRS_dyadic_set_prec(rop,sizeof(unsigned long));
mpfr_set_ui(rop,u,GMP_RNDN);
CGAL_assertion(!mpfr_cmp_ui(rop,u));
}
inline void CGALRS_dyadic_set_fr(CGALRS_dyadic_ptr rop,mpfr_srcptr op){
if(rop!=op){
CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op));
mpfr_set(rop,op,GMP_RNDN);
CGAL_assertion(mpfr_equal_p(rop,op)!=0);
}
}
#define CGALRS_dyadic_init_set(R,D) \
( CGALRS_dyadic_init(R), CGALRS_dyadic_set((R), (D)) )
#define CGALRS_dyadic_init_set_z(R,Z) \
( CGALRS_dyadic_init(R), CGALRS_dyadic_set_z((R), (Z)) )
#define CGALRS_dyadic_init_set_si(R,I) \
( CGALRS_dyadic_init(R), CGALRS_dyadic_set_si((R), (I)) )
#define CGALRS_dyadic_init_set_ui(R,I) \
( CGALRS_dyadic_init(R), CGALRS_dyadic_set_ui((R), (I)) )
#define CGALRS_dyadic_init_set_fr(R,F) \
( CGALRS_dyadic_init(R), CGALRS_dyadic_set_fr((R), (F)) )
#define CGALRS_dyadic_get_fr(M,D) mpfr_set(M,D,GMP_RNDN)
#define CGALRS_dyadic_get_d(D,RM) mpfr_get_d(D,RM)
inline void CGALRS_dyadic_get_exactfr(mpfr_ptr rop,CGALRS_dyadic_srcptr op){
if(rop!=op){
CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op));
mpfr_set(rop,op,GMP_RNDN);
CGAL_assertion(mpfr_equal_p(rop,op)!=0);
}
}
#define CGALRS_dyadic_canonicalize(D) ()
// comparison functions
#define CGALRS_dyadic_sgn(D) mpfr_sgn(D)
#define CGALRS_dyadic_zero(D) mpfr_zero_p(D)
#define CGALRS_dyadic_cmp(D,E) mpfr_cmp(D,E)
// arithmetic functions
#define CGALRS_dyadic_add(R,D,E) CGALRS_dyadic_ll_add(R,D,E,0)
#define CGALRS_dyadic_sub(R,D,E) CGALRS_dyadic_ll_sub(R,D,E,0)
#define CGALRS_dyadic_mul(R,D,E) CGALRS_dyadic_ll_mul(R,D,E,0)
inline void CGALRS_dyadic_neg(CGALRS_dyadic_ptr rop,CGALRS_dyadic_srcptr op){
if(rop!=op)
CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op));
mpfr_neg(rop,op,GMP_RNDN);
CGAL_assertion(
rop==op||
(!mpfr_cmpabs(rop,op)&&
((CGALRS_dyadic_zero(op)&&CGALRS_dyadic_zero(rop))||
(CGALRS_dyadic_sgn(op)!=CGALRS_dyadic_sgn(rop)))));
}
// low-level addition:
// add op1 and op2 and reserve b bits for future lowlevel operations
inline void CGALRS_dyadic_ll_add(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
CGALRS_dyadic_srcptr op2,
mp_prec_t b){
mp_exp_t l,r,temp1,temp2;
mp_prec_t rop_prec;
if(mpfr_zero_p(op1)){
if(rop!=op2)
CGALRS_dyadic_set(rop,op2);
return;
}
if(mpfr_zero_p(op2)){
if(rop!=op1)
CGALRS_dyadic_set(rop,op1);
return;
}
l=mpfr_get_exp(op1)>mpfr_get_exp(op2)?
mpfr_get_exp(op1):
mpfr_get_exp(op2);
temp1=mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1);
temp2=mpfr_get_exp(op2)-(mp_exp_t)mpfr_get_prec(op2);
r=temp1>temp2?temp2:temp1;
CGAL_assertion(l>r);
rop_prec=b+1+(mp_prec_t)(l-r);
CGAL_assertion(rop_prec>=mpfr_get_prec(op1)&&
rop_prec>=mpfr_get_prec(op2));
if(rop==op1||rop==op2)
CGALRS_dyadic_prec_round(rop,rop_prec);
else
CGALRS_dyadic_set_prec(rop,rop_prec);
CGAL_assertion_code(int round=)
mpfr_add(rop,op1,op2,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_add_z(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
mpz_srcptr z){
mp_exp_t l,r;
mp_prec_t rop_prec;
if(mpfr_zero_p(op1)){
CGALRS_dyadic_set_z(rop,z);
return;
}
if(!mpz_sgn(z)){
if(rop!=op1)
CGALRS_dyadic_set(rop,op1);
return;
}
l=mpfr_get_exp(op1)>(mp_exp_t)mpz_sizeinbase(z,2)?
mpfr_get_exp(op1):
mpz_sizeinbase(z,2);
r=mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1)<0?
mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1):
0;
CGAL_assertion(l>r);
rop_prec=1+(mp_prec_t)(l-r);
CGAL_assertion(rop_prec>=mpfr_get_prec(op1)&&
rop_prec>=(mp_prec_t)mpz_sizeinbase(z,2));
if(rop==op1)
CGALRS_dyadic_prec_round(rop,rop_prec);
else
CGALRS_dyadic_set_prec(rop,rop_prec);
CGAL_assertion_code(int round=)
mpfr_add_z(rop,op1,z,GMP_RNDN);
CGAL_assertion(!round);
}
// low-level subtraction:
// subtract op2 to op1 and reserve b bits for future lowlevel operations
inline void CGALRS_dyadic_ll_sub(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
CGALRS_dyadic_srcptr op2,
mp_prec_t b){
mp_exp_t l,r,temp1,temp2;
mp_prec_t rop_prec;
if(mpfr_zero_p(op1)){
CGALRS_dyadic_neg(rop,op2);
return;
}
if(mpfr_zero_p(op2)){
if(rop!=op1)
CGALRS_dyadic_set(rop,op1);
return;
}
l=mpfr_get_exp(op1)>mpfr_get_exp(op2)?
mpfr_get_exp(op1):
mpfr_get_exp(op2);
temp1=mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1);
temp2=mpfr_get_exp(op2)-(mp_exp_t)mpfr_get_prec(op2);
r=temp1>temp2?temp2:temp1;
CGAL_assertion(l>r);
rop_prec=b+1+(mp_prec_t)(l-r);
CGAL_assertion(rop_prec>=mpfr_get_prec(op1)&&
rop_prec>=mpfr_get_prec(op2));
if(rop==op1||rop==op2)
CGALRS_dyadic_prec_round(rop,rop_prec);
else
CGALRS_dyadic_set_prec(rop,rop_prec);
CGAL_assertion_code(int round=)
mpfr_sub(rop,op1,op2,GMP_RNDN);
CGAL_assertion(!round);
}
// low-level multiplication:
// multiply op1 and op2 and reserve b bits for future lowlevel operations
inline void CGALRS_dyadic_ll_mul(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
CGALRS_dyadic_srcptr op2,
mp_prec_t b){
if(rop==op1||rop==op2)
CGALRS_dyadic_prec_round(rop,b+mpfr_get_prec(op1)+mpfr_get_prec(op2));
else
CGALRS_dyadic_set_prec(rop,b+mpfr_get_prec(op1)+mpfr_get_prec(op2));
CGAL_assertion_code(int round=)
mpfr_mul(rop,op1,op2,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_mul_z(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
mpz_srcptr z){
if(rop==op1)
CGALRS_dyadic_prec_round(
rop,
mpfr_get_prec(op1)+mpz_sizeinbase(z,2));
else
CGALRS_dyadic_set_prec(
rop,
mpfr_get_prec(op1)+mpz_sizeinbase(z,2));
CGAL_assertion_code(int round=)
mpfr_mul_z(rop,op1,z,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_mul_si(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
long s){
if(rop==op1)
CGALRS_dyadic_prec_round(rop,mpfr_get_prec(op1)+sizeof(long));
else
CGALRS_dyadic_set_prec(rop,mpfr_get_prec(op1)+sizeof(long));
CGAL_assertion_code(int round=)
mpfr_mul_si(rop,op1,s,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_mul_ui(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
unsigned long u){
if(rop==op1)
CGALRS_dyadic_prec_round(
rop,
mpfr_get_prec(op1)+sizeof(unsigned long));
else
CGALRS_dyadic_set_prec(
rop,
mpfr_get_prec(op1)+sizeof(unsigned long));
CGAL_assertion_code(int round=)
mpfr_mul_ui(rop,op1,u,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_pow_ui(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
unsigned long u){
if(!u){
CGAL_assertion_msg(!mpfr_zero_p(op1),"0^0");
CGALRS_dyadic_set_ui(rop,1);
return;
}
if(u==1){
if(rop!=op1)
CGALRS_dyadic_set(rop,op1);
return;
}
if(mpfr_zero_p(op1)){
CGAL_assertion_msg(u!=0,"0^0");
CGALRS_dyadic_set_ui(rop,0);
return;
}
if(!mpfr_cmp_ui(op1,1)){
if(rop!=op1)
CGALRS_dyadic_set(rop,op1);
return;
}
if(rop==op1)
CGALRS_dyadic_prec_round(rop,u*mpfr_get_prec(op1));
else
CGALRS_dyadic_set_prec(rop,u*mpfr_get_prec(op1));
CGAL_assertion_code(int round=)
mpfr_pow_ui(rop,op1,u,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_addmul(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
CGALRS_dyadic_srcptr op2){
CGALRS_dyadic_t temp;
CGALRS_dyadic_init(temp);
CGALRS_dyadic_mul(temp,op1,op2);
CGALRS_dyadic_add(rop,rop,temp);
CGALRS_dyadic_clear(temp);
}
inline void CGALRS_dyadic_addmul_si(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
long op2){
CGALRS_dyadic_t temp;
CGALRS_dyadic_init(temp);
CGALRS_dyadic_mul_si(temp,op1,op2);
CGALRS_dyadic_add(rop,rop,temp);
CGALRS_dyadic_clear(temp);
}
inline void CGALRS_dyadic_addmul_ui(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
unsigned long u){
//CGALRS_dyadic_t temp;
//CGALRS_dyadic_init(temp);
//CGALRS_dyadic_mul_ui(temp,op1,u);
//CGALRS_dyadic_add(rop,rop,temp);
//CGALRS_dyadic_clear(temp);
CGALRS_dyadic_t temp;
mp_exp_t l,r,temp1,temp2;
mp_prec_t rop_prec;
if(u==0||mpfr_zero_p(op1))
return;
if(u==1){
CGALRS_dyadic_add(rop,rop,op1);
return;
}
// TODO: if(op1==1)
// calculate temp=op1*u
mpfr_init2(temp,mpfr_get_prec(op1)+sizeof(unsigned int));
CGAL_assertion_code(int round1=)
mpfr_mul_ui(temp,op1,u,GMP_RNDN);
CGAL_assertion(!round1);
// calculate the precision needed for rop
l=mpfr_get_exp(op1)>0?mpfr_get_exp(op1):0;
temp1=mpfr_get_exp(op1)-(mp_exp_t)mpfr_get_prec(op1);
temp2=sizeof(unsigned long);
r=temp1>temp2?temp2:temp1;
CGAL_assertion(l>r);
rop_prec=sizeof(unsigned long)+1+(mp_prec_t)(l-r);
CGAL_assertion(rop_prec>=mpfr_get_prec(op1)&&
rop_prec>=mpfr_get_prec(rop));
// set precision and add
CGALRS_dyadic_prec_round(rop,rop_prec);
CGAL_assertion_code(int round2=)
mpfr_add(rop,rop,temp,GMP_RNDN);
CGAL_assertion(!round2);
}
inline void CGALRS_dyadic_mul_2exp(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
unsigned long ui){
// mpfr_mul_2ui does change the mantissa!
if(rop==op1)
CGALRS_dyadic_prec_round(
rop,
sizeof(unsigned long)+mpfr_get_prec(op1));
else
CGALRS_dyadic_set_prec(
rop,
sizeof(unsigned long)+mpfr_get_prec(op1));
CGAL_assertion_code(int round=)
mpfr_mul_2ui(rop,op1,ui,GMP_RNDN);
CGAL_assertion(!round);
}
inline void CGALRS_dyadic_div_2exp(CGALRS_dyadic_ptr rop,
CGALRS_dyadic_srcptr op1,
unsigned long ui){
// mpfr_div_2ui does not change the mantissa... am I sure?
CGAL_assertion_code(int round=)
mpfr_div_2ui(rop,op1,ui,GMP_RNDN);
CGAL_assertion(!round);
}
// miscelaneous functions
#define CGALRS_dyadic_midpoint(R,D,E) \
( CGALRS_dyadic_ll_add(R,D,E,1) , mpfr_div_2ui(R,R,1,GMP_RNDN) )
#define CGALRS_dyadic_swap(D,E) mpfr_swap(D,E)
// I/O functions
#define CGALRS_dyadic_out_str(F,D) mpfr_out_str(F,10,0,D,GMP_RNDN)
#ifdef __cplusplus
inline std::ostream& operator<<(std::ostream &s,CGALRS_dyadic_srcptr op){
mp_exp_t exponent;
mpz_t mantissa;
mpz_init(mantissa);
exponent=mpfr_get_z_exp(mantissa,op);
s<<"["<<mantissa<<","<<exponent<<"]";
return s;
}
#endif
#endif // CGAL_RS_DYADIC_H

157
Algebraic_kernel_d/include/CGAL/RS/isolator_1.h
View File

@ -1,157 +0,0 @@
// Copyright (c) 2006-2009, 2011 Max-Planck-Institute Saarbruecken (Germany).
// 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(s) : Eric Berberich <eric@mpi-inf.mpg.de>
// Michael Hemmer <hemmer@mpi-inf.mpg.de>
// Luis Peñaranda <luis.penaranda@gmx.com>
//
// ============================================================================
/// \file RS/isolator.h
/// Defines class CGAL::RS_real_root_isolator
///
/// Isolate real roots of polynomials with Fabrice Roullier's Rs.
///
/// The polynomial has to be a univariate polynomial over any number type
/// which is contained in the real numbers. The polynomial does not need to
/// be square-free.
///
#ifndef CGAL_ALGEBRAIC_KERNEL_D_RS_ISOLATOR_H
#define CGAL_ALGEBRAIC_KERNEL_D_RS_ISOLATOR_H
#include <CGAL/config.h>
#include <CGAL/Gmpfr.h>
#include <CGAL/Gmpfi.h>
#include <CGAL/RS/isole_1.h>
namespace CGAL {
namespace internal {
/// A model of concept RealRootIsolator.
///
/// Polynomial_ must be Polynomial<Gmpz>, and Bound_ must be Gmpfr.
template <class Polynomial_, class Bound_>
class RS_real_root_isolator {
public:
//! First template parameter
typedef Polynomial_ Polynomial;
//! Second template parameter
typedef Bound_ Bound;
private:
//! Coefficient type of polynomial
typedef typename Polynomial::NT Coefficient;
typedef Gmpfi Interval;
public:
/// Constructor from univariate square free polynomial.
/// The RealRootIsolator provides isolating intervals for the real
/// roots of the polynomial
RS_real_root_isolator(const Polynomial& p = Polynomial(Coefficient(0))) :
input_polynomial(p)
//, interval_given(false)
{
real_roots=RS::isolator<Polynomial>()(p);
}
// @LUIS: add constructor from interval (maybe some time in the future)
public: // functions
//! Returns the defining polynomial.
Polynomial polynomial() const {
return input_polynomial;
}
//! Returns the number of real roots.
int number_of_real_roots() const {
return real_roots.size();
}
/// Returns true if the isolating interval is degenerated to a
/// single point.
///
/// If is_exact_root(i) is true,
/// then left_bound(int i) equals \f$root_i\f$. \n
/// If is_exact_root(i) is true,
/// then right_bound(int i) equals \f$root_i\f$. \n
bool is_exact_root(int i) const {
return(real_roots[i].inf()==real_roots[i].sup());
}
public:
/// Returns \f${l_i}\f$ the left bound of the isolating interval
/// for root \f$root_{i}\f$.
///
/// In case is_exact_root(i) is true, \f$l_i = root_{i}\f$,\n
/// otherwise: \f$l_i < root_{i}\f$.
///
/// If \f$i-1>=0\f$, then \f$l_i > root_{i-1}\f$. \n
/// If \f$i-1>=0\f$, then \f$l_i >= r_{i-1}\f$,
/// the right bound of \f$root_{i-1}\f$\n
///
/// \pre 0 <= i < number_of_real_roots()
Bound left_bound(int i) const {
CGAL_assertion(i >= 0);
CGAL_assertion(i < this->number_of_real_roots());
return real_roots[i].inf();
}
/// Returns \f${r_i}\f$ the right bound of the isolating interval
/// for root \f$root_{i}\f$.
///
/// In case is_exact_root(i) is true, \f$r_i = root_{i}\f$,\n
/// otherwise: \f$r_i > root_{i}\f$.
///
/// If \f$i+1< n \f$, then \f$r_i < root_{i+1}\f$,
/// where \f$n\f$ is number of real roots.\n
/// If \f$i+1< n \f$, then \f$r_i <= l_{i+1}\f$,
/// the left bound of \f$root_{i+1}\f$\n
///
/// \pre 0 <= i < number_of_real_roots()
Bound right_bound(int i) const {
CGAL_assertion(i >= 0);
CGAL_assertion(i < this->number_of_real_roots());
return real_roots[i].sup();
}
private:
//! the input polynomial
Polynomial input_polynomial;
//! the solutions
std::vector<Interval> real_roots;
//! restricted interval?
// TODO bool interval_given;
};
} // namespace internal
} //namespace CGAL
#endif // CGAL_ALGEBRAIC_KERNEL_D_RS_ISOLATOR_H

79
Algebraic_kernel_d/include/CGAL/RS/isole_1.h
View File

@ -1,79 +0,0 @@
// Copyright (c) 2011 National and Kapodistrian University of Athens (Greece).
// 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_ISOLE_1_H
#define CGAL_RS_ISOLE_1_H
#include <CGAL/RS/basic.h>
#include <CGAL/RS/rs_calls_1.h>
#include <CGAL/Gmpfi.h>
#include <vector>
namespace CGAL{
namespace RS{
// CGAL::RS::isolator<P> solves a polynomial of type P and returns a vector
// of Gmpfi's containing the solutions. Currently, the only implemented
// isolator solves polynomials of type CGAL::Polynomial<CGAL::Gmpz>.
template <class PolynomialType>
struct isolator{
typedef PolynomialType PolynomialT;
inline std::vector<Gmpfi>
operator()(const PolynomialT&,unsigned int prec=CGAL_RS_DEF_PREC);
};
template <class PolynomialType>
inline std::vector<Gmpfi>
isolator<PolynomialType>::operator()(const PolynomialType &p,unsigned int prec){
CGAL_error_msg(
"isolator not implemented for this type of polynomials");
return std::vector<Gmpfi>();
}
template <>
inline std::vector<Gmpfi>
isolator<Polynomial<Gmpz> >::operator()(const Polynomial<Gmpz> &p,
unsigned int prec){
int numsols;
unsigned int degree=p.degree();
mpz_t *coeffs=(mpz_t*)malloc((degree+1)*sizeof(mpz_t));
mpfi_ptr *intervals_mpfi=(mpfi_ptr*)malloc(degree*sizeof(mpfi_ptr));
std::vector<Gmpfi> intervals;
for(unsigned int i=0;i<=degree;++i)
coeffs[i][0]=*(p[i].mpz());
init_solver();
create_rs_upoly(coeffs,degree,rs_get_default_up());
free(coeffs);
set_rs_precisol(prec);
set_rs_verbose(CGAL_RS_VERB);
rs_run_algo(CGALRS_CSTR("UISOLE"));
numsols=affiche_sols_eqs(intervals_mpfi);
for(int j=0;j<numsols;++j)
intervals.push_back(Gmpfi(intervals_mpfi[j]));
free(intervals_mpfi);
return intervals;
}
} // namespace RS
} // namespace CGAL
#endif // CGAL_RS_ISOLE_1_H

115
Algebraic_kernel_d/include/CGAL/RS/memory.h
View File

@ -1,115 +0,0 @@
// Copyright (c) 2009 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_MEMORY_H
#define CGAL_RS_MEMORY_H
#include <CGAL/RS/basic.h>
#include <gmp.h>
#include <cstdlib>
#include <rs_exports.h>
#ifdef CGAL_USE_OLD_RS3
extern "C"{
#include <rs_gc.h>
}
#endif
namespace CGAL{
extern "C"{
extern void* RS_gmpalloc(size_t);
extern void* RS_gmprealloc(void*,size_t,size_t);
extern void RS_gmpfree(void*,size_t);
}
#ifdef CGAL_USE_OLD_RS3
inline void * rs3_rs_alloc(size_t s) {
return(rs_alloc(s,&(rs_default_gc_session.default_heap[0])));
}
inline void * rs3_rs_realloc(void * p,size_t s) {
assert(1==0);
return(NULL);
}
inline void rs3_rs_free(void * p) {
assert(1==0);
}
inline void * rs3_rs_gmp_realloc(void * old_p,size_t old_size,size_t new_size){
return(rs_realloc(old_p,((RS_UI)old_size),((RS_UI)new_size),&(rs_default_gc_session.default_heap[0])));
}
inline void rs3_rs_gmp_free(void * p,size_t s) {}
#endif // CGAL_USE_OLD_RS3
//--------------------------------------------------
// extern void * (*cgalrs_allocate_func) (size_t);
// extern void * (*cgalrs_reallocate_func) (void *, size_t, size_t);
// extern void (*cgalrs_free_func) (void *, size_t);
//--------------------------------------------------
inline void* cgalrs_default_allocate(size_t s){
return malloc(s);
}
inline void* cgalrs_default_reallocate(void *a,size_t o,size_t n){
return realloc(a,n);
}
inline void cgalrs_default_free(void *a,size_t s){
return free(a);
}
CGALRS_THREAD_ATTR void * (*cgalrs_allocate_func) (size_t) =
cgalrs_default_allocate;
CGALRS_THREAD_ATTR void * (*cgalrs_reallocate_func) (void *, size_t, size_t) =
cgalrs_default_reallocate;
CGALRS_THREAD_ATTR void (*cgalrs_free_func) (void *, size_t) =
cgalrs_default_free;
inline void cgalrs_dummy_free(void *p,size_t s){}
inline void cgalrs_set_memory_functions(
void *(*alloc_func) (size_t),
void *(*realloc_func) (void *, size_t, size_t),
void (*free_func) (void *, size_t)){
cgalrs_allocate_func=
alloc_func?alloc_func:cgalrs_default_allocate;
cgalrs_reallocate_func=
realloc_func?realloc_func:cgalrs_default_reallocate;
cgalrs_free_func=
free_func?free_func:cgalrs_default_free;
}
inline void cgalrs_get_memory_functions(
void *(**alloc_func) (size_t),
void *(**realloc_func) (void *, size_t, size_t),
void (**free_func) (void *, size_t)){
if(alloc_func!=NULL)
*alloc_func=cgalrs_allocate_func;
if(realloc_func!=NULL)
*realloc_func=cgalrs_reallocate_func;
if(free_func!=NULL)
*free_func=cgalrs_free_func;
}
} // namespace CGAL
#endif // CGAL_RS_MEMORY_H

144
Algebraic_kernel_d/include/CGAL/RS/polynomial_1.h
View File

@ -1,144 +0,0 @@
// Copyright (c) 2006-2010 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_POLYNOMIAL_1_H
#define CGAL_RS_POLYNOMIAL_1_H
#include <iostream>
#include <vector>
#include <CGAL/RS/basic.h>
#include <CGAL/Gmpz.h>
#include <CGAL/Gmpfr.h>
#include <CGAL/RS/dyadic.h>
#include <gmp.h>
#include <mpfr.h>
#include <mpfi.h>
#include <boost/operators.hpp>
#include <boost/shared_ptr.hpp>
namespace CGAL{
class Algebraic_1;
class RS_polynomial_1;
typedef boost::shared_ptr<RS_polynomial_1> polyptr;
typedef std::pair<RS_polynomial_1,int> polypow;
typedef std::vector<polypow> sqfrvec;
typedef boost::shared_ptr<sqfrvec> sqfrptr;
class RS_polynomial_1:
boost::addable1<RS_polynomial_1,
boost::subtractable1<RS_polynomial_1
> >
{
private:
int capacity;
mutable int polynomial_degree;
mpz_t* coefficient_array;
mutable bool is_square_free;
mutable polyptr square_free_part;
mutable sqfrptr square_free_factorization;
void create_storage(int);
void free_storage();
// fetch_gmp_functions gathers the memory functions used by
// gmp at the object creation and stores them in
// alloc_function, realloc_function and free_function
void *(*alloc_function)(size_t);
void *(*realloc_function)(void*,size_t,size_t);
void (*free_function)(void*,size_t);
void fetch_gmp_functions();
public:
// copy constructor and copy assignement operator
RS_polynomial_1(const RS_polynomial_1&);
RS_polynomial_1& operator=(const RS_polynomial_1&);
// other constructors and destructor
RS_polynomial_1();
RS_polynomial_1(unsigned int);
RS_polynomial_1(int);
RS_polynomial_1(std::string&);
RS_polynomial_1(mpq_srcptr);
RS_polynomial_1(mpz_t**,int);
~RS_polynomial_1();
// member functions
void set_degree(int);
void force_degree(int);
int resize(int);
void set_coef(int,mpz_srcptr);
void set_coef(int,const CGAL::Gmpz&);
void set_coef_ui(int,unsigned long);
void set_coef_si(int,long);
int get_degree()const;
int get_degree_static()const;
bool has_sfpart()const;
const RS_polynomial_1& sfpart()const;
void set_sfpart(RS_polynomial_1*)const;
void set_sfpart(const polyptr&)const;
void set_sf()const;
bool has_sqfr()const;
sqfrvec& sqfr()const;
void set_sqfr(sqfrvec*)const;
void set_sqfr(const sqfrptr&)const;
mpz_ptr leading_coefficient()const;
int first_non_zero()const;
mpz_t* get_coefs()const;
mpz_ptr coef(int)const;
RS_polynomial_1& derive()const;
RS_polynomial_1& minusx()const;
void get_lower_bound(mpfr_ptr)const;
void get_upper_bound(mpfr_ptr)const;
RS_polynomial_1 times_monomial(mpz_srcptr,int)const;
// member evaluation and sign functions
void eval_dyadic(CGALRS_dyadic_ptr,CGALRS_dyadic_srcptr)const;
void eval_mpfr(mpfr_ptr,mpfr_srcptr)const;
void inexact_eval_mpfr(mpfr_ptr,mpfr_srcptr)const;
void eval_mpfi(mpfi_ptr,mpfi_srcptr)const;
Sign sign_dyadic(CGALRS_dyadic_srcptr)const;
Sign sign_mpfr(mpfr_srcptr)const;
RS::rs_sign sign_mpfi(mpfi_srcptr)const;
double operator()(double)const;
CGAL::Gmpz operator()(int)const;
RS_polynomial_1 operator-()const;
RS_polynomial_1& operator+=(const RS_polynomial_1&);
RS_polynomial_1& operator-=(const RS_polynomial_1&);
RS_polynomial_1 operator*(const RS_polynomial_1&)const;
RS_polynomial_1& operator*=(const RS_polynomial_1&);
RS_polynomial_1& operator*=(mpz_srcptr);
RS_polynomial_1& operator*=(const CGAL::Gmpz &);
// division is always assumed to be exact in this class
RS_polynomial_1& operator/=(mpz_srcptr);
RS_polynomial_1& operator/=(const CGAL::Gmpz&);
bool operator==(const RS_polynomial_1&)const;
// template members, including constructor
template<class InputIt>RS_polynomial_1(InputIt,InputIt);
template<class T>T operator()(const T&)const;
template<class T>Sign sign_at(const T&)const;
template<class T>RS_polynomial_1 operator*(const T&)const;
template<class T>RS_polynomial_1& operator*=(const T&);
template<class T>RS_polynomial_1& operator/=(const T&);
template<class T>RS_polynomial_1 operator/(const T&)const;
};
} // namespace CGAL
#include <CGAL/RS/polynomial_1_constructors.h>
#include <CGAL/RS/polynomial_1_eval.h>
#include <CGAL/RS/polynomial_1_member.h>
#include <CGAL/RS/polynomial_1_operators.h>
#include <CGAL/RS/polynomial_1_impl.h>
#include <CGAL/RS/polynomial_1_io.h>
#endif // CGAL_GRBS_POLYNOMIAL_1_H

169
Algebraic_kernel_d/include/CGAL/RS/polynomial_1_constructors.h
View File

@ -1,169 +0,0 @@
// Copyright (c) 2006-2009 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_POLYNOMIAL_1_CONSTRUCTORS_H
#define CGAL_RS_POLYNOMIAL_1_CONSTRUCTORS_H
#include <CGAL/assertions.h>
#include <CGAL/RS/polynomial_1_parser.h>
namespace CGAL{
// auxiliar function to parse a c++ string and return a vector
// containing the coefficients as Gmpz's
inline
std::vector<Gmpz> parse_string(std::string tstr){
Polynomial_parser_1 parser;
parser.parse(tstr);
CGAL_assertion(parser.is_correct());
std::vector<Gmpz> Coeff;
parser.result(std::back_inserter(Coeff),CGAL::Convert_to_Gmpz());
return Coeff;
}
// copy constructor
inline
RS_polynomial_1::RS_polynomial_1(const RS_polynomial_1 &p){
polynomial_degree=p.polynomial_degree;
fetch_gmp_functions();
create_storage(polynomial_degree+1);
// we have to copy the contents, not just the pointer
for(int i=0;i<polynomial_degree+1;++i)
mpz_init_set(coef(i),p.coef(i));
is_square_free=p.is_square_free;
square_free_part=p.square_free_part;
square_free_factorization=p.square_free_factorization;
}
// copy assignement operator
inline
RS_polynomial_1& RS_polynomial_1::operator=(const RS_polynomial_1 &p){
square_free_part=p.square_free_part;
square_free_factorization=p.square_free_factorization;
if(capacity<p.get_degree()+1){
// destroy the current data
free_storage();
// copy data from p
polynomial_degree=p.get_degree();
create_storage(polynomial_degree+1);
for(int i=0;i<polynomial_degree+1;++i)
mpz_init_set(coef(i),p.coef(i));
return *this;
}
polynomial_degree=p.get_degree();
for(int i=0;i<polynomial_degree+1;++i)
mpz_init_set(coef(i),p.coef(i));
for(int i=polynomial_degree+1;i<capacity;++i)
mpz_init_set_ui(coef(i),0);
return *this;
}
// other constructors
inline
RS_polynomial_1::RS_polynomial_1():
polynomial_degree(0),
is_square_free(false),
square_free_part(polyptr()),
square_free_factorization(sqfrptr()){
fetch_gmp_functions();
create_storage(1);
mpz_init(coef(0));
}
inline
RS_polynomial_1::RS_polynomial_1(unsigned int d):
is_square_free(false),
square_free_part(polyptr()),
square_free_factorization(sqfrptr()){
fetch_gmp_functions();
polynomial_degree=(int)d;
create_storage(polynomial_degree+1);
for(int i=0;i<polynomial_degree+1;++i)
mpz_init(coef(i));
}
inline
RS_polynomial_1::RS_polynomial_1(int d):
is_square_free(false),
square_free_part(polyptr()),
square_free_factorization(sqfrptr()){
fetch_gmp_functions();
polynomial_degree=d<0?0:d;
create_storage(polynomial_degree+1);
for(int i=0;i<polynomial_degree+1;++i)
mpz_init(coef(i));
}
// constructor from a fuckin' c++ string
inline
RS_polynomial_1::RS_polynomial_1(std::string &tstr):
is_square_free(false),
square_free_part(polyptr()),
square_free_factorization(sqfrptr()){
std::vector<Gmpz> Coeff=parse_string(tstr);
std::vector<Gmpz>::iterator it,first,last;
int elements;
first=Coeff.begin();
last=Coeff.end();
elements=Coeff.size();
polynomial_degree=elements-1;
fetch_gmp_functions();
create_storage(elements);
for(it=first;it!=last;++it)
mpz_init_set(coef(polynomial_degree-(--elements)),(*it).mpz());
}
// construct a polynomial whose root is the given rational
inline
RS_polynomial_1::RS_polynomial_1(mpq_srcptr r):
is_square_free(false),
square_free_part(polyptr()),
square_free_factorization(sqfrptr()){
polynomial_degree=1;
fetch_gmp_functions();
create_storage(2);
mpz_init(coef(0));
mpz_neg(coef(0),mpq_numref(r));
mpz_init(coef(1));
mpz_set(coef(1),mpq_denref(r));
}
// p has the address of the coefficients array and d is the degree
inline
RS_polynomial_1::RS_polynomial_1(mpz_t** p,int d):
is_square_free(false),
square_free_part(polyptr()),
square_free_factorization(sqfrptr()){
// at this point, GMP memory functions must be the same
// used to initialize the coefficients and array of the
// polynomial; otherwise, some problem may arise at
// destruction time
fetch_gmp_functions();
capacity=d+1;
polynomial_degree=d;
coefficient_array=*p;
}
inline
RS_polynomial_1::~RS_polynomial_1(){
free_storage();
}
} // namespace CGAL
#endif // CGAL_RS_POLYNOMIAL_1_CONSTRUCTORS_H

126
Algebraic_kernel_d/include/CGAL/RS/polynomial_1_eval.h
View File

@ -1,126 +0,0 @@
// Copyright (c) 2006-2010 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_POLYNOMIAL_1_EVAL_H
#define CGAL_RS_POLYNOMIAL_1_EVAL_H
namespace CGAL{
// This function uses the Horner's method to evaluate the polynomial.
inline
void RS_polynomial_1::eval_dyadic(CGALRS_dyadic_ptr h,
CGALRS_dyadic_srcptr x)const{
int d=get_degree();
CGALRS_dyadic_set_z(h,leading_coefficient());
for(int i=1;i<d+1;++i){
CGALRS_dyadic_mul(h,h,x);
CGALRS_dyadic_add_z(h,h,coef(d-i));
}
}
inline
void RS_polynomial_1::eval_mpfr(mpfr_ptr result,mpfr_srcptr xcoord)const{
// mpfr and dyadic are now the same struct
this->eval_dyadic(result,xcoord);
}
inline
void RS_polynomial_1::inexact_eval_mpfr(mpfr_ptr h,mpfr_srcptr x)const{
int d=get_degree();
mpfr_set_z(h,leading_coefficient(),GMP_RNDN);
for(int i=1;i<d+1;++i){
mpfr_mul(h,h,x,GMP_RNDN);
mpfr_add_z(h,h,coef(d-i),GMP_RNDN);
}
}
// Evaluation of a polynomial in a given interval using Horner's rule.
// y=p(x), y must be already initialized.
inline
void RS_polynomial_1::eval_mpfi(mpfi_ptr y,mpfi_srcptr x)const{
int d=get_degree();
mpfi_set_z(y,leading_coefficient());
for(int i=1;i<d+1;++i){
mpfi_mul(y,y,x);
mpfi_add_z(y,y,coef(d-i));
}
}
// Sign determination is essentially the same that evaluation. As we only
// want to know the sign, we don't finish the Horner's evaluation.
inline
Sign RS_polynomial_1::sign_dyadic(CGALRS_dyadic_srcptr x)const{
int d,s;
CGALRS_dyadic_t h;
if(!(d=get_degree())){
s=mpz_sgn(coef(0));
return(s?(s>0?CGAL::POSITIVE:CGAL::NEGATIVE):CGAL::ZERO);
}
CGALRS_dyadic_init_set_z(h,leading_coefficient());
for(int i=1;i<d;++i){
CGALRS_dyadic_mul(h,h,x);
CGALRS_dyadic_add_z(h,h,coef(d-i));
}
CGALRS_dyadic_mul(h,h,x);
CGALRS_dyadic_neg(h,h);
s=mpfr_cmp_z(h,coef(0));
return(s?(s>0?CGAL::NEGATIVE:CGAL::POSITIVE):CGAL::ZERO);
}
inline
Sign RS_polynomial_1::sign_mpfr(mpfr_srcptr x)const{
// no need to convert mpfr to dyadic anymore:
// they are now the same struct
return sign_dyadic(x);
}
// calculates the sign of the evaluation of a polynomial in a given interval
inline
RS::rs_sign RS_polynomial_1::sign_mpfi(mpfi_srcptr x)const{
mpfi_t y;
int l,r,d;
if(!(d=get_degree())){
l=mpz_sgn(coef(0));
return(l?(l>0?RS::RS_POSITIVE:RS::RS_NEGATIVE):RS::RS_ZERO);
}
// TODO: check if the precision is ok (Fabrice said that we lose
// one bit per operation, there are 2*d operations and with three
// we are on the safe side, but with two it works better)
mpfi_init2(y,mpfi_get_prec(x)+2*d);
mpfi_set_z(y,leading_coefficient());
for(int i=d-1;i>=0;--i){
mpfi_mul(y,y,x);
mpfi_add_z(y,y,coef(i));
}
l=mpfr_sgn(&y->left);
if(l>0){
mpfi_clear(y);
return RS::RS_POSITIVE;
}
r=mpfr_sgn(&y->right);
mpfi_clear(y);
if(r<0)
return RS::RS_NEGATIVE;
if(!l&&!r)
return RS::RS_ZERO;
return RS::RS_UNKNOWN;
}
} // namespace CGAL
#endif // CGAL_RS_POLYNOMIAL_1_EVAL_H

347
Algebraic_kernel_d/include/CGAL/RS/polynomial_1_functors.h
View File

@ -1,347 +0,0 @@
// Copyright (c) 2007-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_POLYNOMIAL_FUNCTORS
#define CGAL_RS_POLYNOMIAL_FUNCTORS
#include <CGAL/basic.h>
#include <CGAL/RS/polynomial_1.h>
#include <CGAL/RS/polynomial_1_utils.h>
#include <CGAL/Exponent_vector.h>
#include <algorithm>
#include <CGAL/assertions.h>
namespace CGAL{
struct RSPolynomialFunctors{
typedef Gmpz Coefficient;
typedef Gmpz Innermost_coefficient;
typedef RS_polynomial_1 Polynomial_1;
typedef RS_polynomial_1 Type;
static const int d=1; // TODO: do this in a clean way...
struct Construct_polynomial{
Polynomial_1 operator()(){
Polynomial_1 p(2);
p.set_coef_ui(0,0);
p.force_degree(0);
return p;
};
Polynomial_1 operator()(int i){
Polynomial_1 p(2);
p.set_coef_si(0,i);
p.force_degree(0);
return p;
};
template <class InputIterator>
inline Polynomial_1 operator()(
InputIterator begin,InputIterator end)const{
return Polynomial_1(begin,end);
};
//--------------------------------------------------
// template <class InputIterator1,class InputIterator2>
// Polynomial_1 operator()
// (InputIterator1 fst_coef,InputIterator1 lst_coef,
// InputIterator2 deg)const{
// // the degree of the polynomial will be
// // the greatest degree
// unsigned greater=0;
// InputIterator1 c;
// InputIterator2 d=deg;
// for(c=fst_coef;c!=lst_coef;++c)
// if(d>greater)
// greater=d++;
// // now, construct the polynomial of degree d
// Polynomial_1 p(d);
// for(c=fst_coef;c!=lst_coef;++c)
// p.set_coef(deg++,*c);
// return p;
// };
//--------------------------------------------------
template <class InputIterator>
Polynomial_1 operator()
(InputIterator begin,InputIterator end,bool is_sorted)const{
// the degree of the polynomial will be
// the greatest degree
int degree;
if(is_sorted)
degree=*((end-1)->first.begin());
else{
degree=0;
for(InputIterator i=begin;i!=end;++i)
if(*(i->first.begin())>degree)
degree=*(i->first.begin());
}
Polynomial_1 p(degree);
for(InputIterator i=begin;i!=end;++i)
p.set_coef(*(i->first.begin()),i->second);
return p;
};
}; // Construct_polynomial
struct Get_coefficient:
public std::binary_function<Polynomial_1,int,Coefficient>{
inline Coefficient operator()
(const Polynomial_1 &p,int i)const{
return Coefficient(p.coef(i));
};
};
struct Get_innermost_coefficient:
public std::binary_function<Polynomial_1,Exponent_vector,Coefficient>{
inline Coefficient operator()(
const Polynomial_1 &p,Exponent_vector v)const{
CGAL_precondition(v.size()==1);
return Coefficient(p.coef(v[0]));
};
};
// this is an in-place swap
struct Swap{
inline Polynomial_1 operator()(Polynomial_1 &p,int i,int j)const{
mpz_t temp;
CGAL_precondition(i<=d&&j<=d);
mpz_init_set(temp,p.coef(j));
p.set_coef(j,p.coef(i));
p.set_coef(i,temp);
mpz_clear(temp);
return p;
};
};
struct Move{
inline Polynomial_1& operator()(Polynomial_1 &p,int i,int j)const{
CGAL_precondition(i<=d&&j<=d);
p.set_coef(j,p.coef(i));
p.set_coef_ui(i,0);
return p;
};
};
struct Degree{
inline int operator()(const Polynomial_1 &p)const{
return p.get_degree();
};
inline int operator()(const Polynomial_1 &p,int i)const{
CGAL_precondition(!d);
return p.get_degree();
};
};
struct Total_degree:
public std::unary_function<Polynomial_1,int>{
inline int operator()(const Polynomial_1 &p)const{
return p.get_degree();
};
};
struct Degree_vector:
public std::unary_function<Polynomial_1,Exponent_vector>{
inline Exponent_vector operator()(const Polynomial_1 &p){
Exponent_vector v(1);
v[1]=p.get_degree();
return v;
// why the hell this doesn't work?
//return Exponent_vector(1,p.get_degree());
};
};
struct Leading_coefficient{
inline Coefficient operator()(const Polynomial_1 &p)const{
return Coefficient(p.leading_coefficient());
};
inline Coefficient operator()
(const Polynomial_1 &p,int i)const{
CGAL_precondition(!i);
return Coefficient(p.leading_coefficient());
};
};
struct Innermost_leading_coefficient:
public std::unary_function<Polynomial_1,Innermost_coefficient>{
inline Innermost_coefficient operator()(const Polynomial_1 &p)const{
return Innermost_coefficient(p.leading_coefficient());
};
};
struct Canonicalize:
public std::unary_function<Polynomial_1,Polynomial_1>{
Polynomial_1& operator()(Polynomial_1 &f)const{
mpz_t temp;
mpz_init(temp);
Cont()(temp,f);
f/=temp;
mpz_clear(temp);
return f;
};
};
struct Derive{
inline Polynomial_1 operator()(const Polynomial_1 &p)const{
return p.derive();
};
inline Polynomial_1 operator()(const Polynomial_1 &p,int i)const{
CGAL_precondition(!i);
return p.derive();
};
};
struct Evaluate{
inline Coefficient operator()(
const Polynomial_1 &p,
const Innermost_coefficient x)const{
return p(x);
};
inline Coefficient operator()(
const Polynomial_1 &p,
Innermost_coefficient c,
int i)const{
CGAL_precondition(!i);
return p(c);
};
};
typedef Null_functor Evaluate_homogeneous;
struct Is_zero_at{
template <class InputIterator>
inline bool operator()(
const Polynomial_1 &p,
InputIterator begin,
InputIterator end)const{
CGAL_assertion(end=begin+1);
return(p(*begin)==0);
};
};
typedef Null_functor Is_zero_at_homogeneous;
struct Sign_at{
template <class InputIterator>
inline Sign operator()(
const Polynomial_1 &p,
InputIterator begin,
InputIterator end)const{
return p.sign_at(*begin);
};
};
typedef Null_functor Sign_at_homogeneous;
struct Compare:
public std::binary_function<Polynomial_1,Polynomial_1,Comparison_result>{
Comparison_result operator()(
const Polynomial_1 &f,
const Polynomial_1 &g){
if(f.get_degree()>g.get_degree())
return LARGER;
if(f.get_degree_static()<g.get_degree_static())
return SMALLER;
int i,c;
for(i=f.get_degree_static();i>=0;--i){
if((c=mpz_cmp(f.coef(i),g.coef(i)))>0)
return LARGER;
if(c<0)
return SMALLER;
}
return EQUAL;
};
}; // Compare
// now, we define five functors which are needed by the algebraic kernel,
// but they are related to polynomials (I think this is the best place
// to do it)
template <class GcdPolicy>
struct Is_square_free_1:
public std::unary_function<Polynomial_1,bool>{
typedef GcdPolicy Gcd;
inline bool operator()(const Polynomial_1 &p){
return(!(Gcd()(p,p.derive()).get_degree_static()));
};
}; // Is_square_free_1
template <class GcdPolicy>
struct Make_square_free_1:
public std::unary_function<Polynomial_1,Polynomial_1>{
typedef GcdPolicy Gcd;
inline Polynomial_1 operator()(const Polynomial_1 &p){
return sfpart_1<Gcd>(p);
};
}; // Make_square_free_1
// this function is not well defined in algebraic kernel concepts; we
// don't care, we do it correctly
template <class GcdPolicy>
struct Square_free_factorize_1{
template <class OutputIterator>
int operator()(const Polynomial_1 &p,OutputIterator oi){
typedef GcdPolicy Gcd;
sqfrvec factorization(sqfr_1<Gcd>(p));
std::copy(factorization.begin(),factorization.end(),oi);
return factorization.size();
};
}; // Square_free_factorize_1
template <class GcdPolicy>
struct Is_coprime_1:
public std::binary_function<Polynomial_1,Polynomial_1,bool>{
inline bool operator()
(const Polynomial_1 &p1,const Polynomial_1 &p2){
typedef GcdPolicy Gcd;
return(!Gcd()(p1,p2).get_degree_static());
};
}; // Is_coprime_1
template <class GcdPolicy>
struct Make_coprime_1{
typedef bool result_type;
typedef Polynomial_1 P;
bool operator()(const P &p1,const P &p2,P &g,P &q1,P &q2){
typedef GcdPolicy Gcd;
g=Gcd()(p1,p2);
// we don't calculate q1 and q2 when g==1, shall we?
if(!(g.get_degree_static()))
return true;
q1=*Ediv_1()(p1,g);
q2=*Ediv_1()(p2,g);
return false;
};
}; // Make_coprime_1
struct IntegralDivision:
public std::binary_function<Polynomial_1,Polynomial_1,Polynomial_1>{
inline Polynomial_1 operator()(
const Polynomial_1 &p1,
const Polynomial_1 &p2){
return *Ediv_1()(p1,p2);
};
}; // IntegralDivision
typedef IntegralDivision IntegralDivisionUpToConstantFactor;
}; // RSPolynomialFunctors
} // namespace CGAL
#endif // CGAL_RS_POLYNOMIAL_FUNCTORS

99
Algebraic_kernel_d/include/CGAL/RS/polynomial_1_impl.h
View File

@ -1,99 +0,0 @@
// 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_POLYNOMIAL_1_IMPL_H
#define CGAL_RS_POLYNOMIAL_1_IMPL_H
namespace CGAL{
// constructor from a container of sorted coefficients
template <class InputIterator>
RS_polynomial_1::RS_polynomial_1(InputIterator first,InputIterator last){
// count the number of elements in the container
int elements=0;
InputIterator it;
for(it=first;it!=last;++it)
++elements;
// now we create the polynomial
if(elements){
polynomial_degree=elements-1;
create_storage(elements);
int i=0;
for(it=first;it!=last;++it)
mpz_init_set(coef(i++),Gmpz(*it).mpz());
}else{
polynomial_degree=0;
create_storage(1);
mpz_init_set_ui(coef(0),0);
}
}
template <class T>
RS_polynomial_1 RS_polynomial_1::operator*(const T &n)const{
RS_polynomial_1 r(*this);
return (r*=n);
}
template<class T>
RS_polynomial_1& RS_polynomial_1::operator*=(const T &s){
Gmpz z(s);
return (*this*=z);
}
template <class T>
RS_polynomial_1 RS_polynomial_1::operator/(const T &d)const{
RS_polynomial_1 r(*this);
return (r/=d);
}
template<class T>
RS_polynomial_1& RS_polynomial_1::operator/=(const T &d){
Gmpz z(d);
return (*this/=z);
}
template<class T>
Sign RS_polynomial_1::sign_at(const T &x)const{
T y(leading_coefficient());
int d=get_degree_static();
for(int i=1;i<d+1;++i)
y=y*x+T(coef(d-i));
if(y>0)
return POSITIVE;
if(y<0)
return NEGATIVE;
return ZERO;
}
template<>
inline Sign RS_polynomial_1::sign_at(const Gmpfr &x)const{
return sign_mpfr(x.fr());
}
template<class T>
T RS_polynomial_1::operator()(const T &x)const{
T y(leading_coefficient());
int d=get_degree_static();
for(int i=1;i<d+1;++i)
y=y*x+T(coef(d-i));
return y;
}
} // namespace CGAL
#endif // CGAL_RS_POLYNOMIAL_1_IMPL_H

134
Algebraic_kernel_d/include/CGAL/RS/polynomial_1_io.h
View File

@ -1,134 +0,0 @@
// Copyright (c) 2006-2009 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_POLYNOMIAL_1_IO_H
#define CGAL_RS_POLYNOMIAL_1_IO_H
namespace CGAL{
inline
std::ostream& operator<<(std::ostream &os,const RS_polynomial_1 &p){
if(is_pretty(os)){
bool printed = false;
int degree=p.get_degree();
mpz_t *coef=p.get_coefs();
if(!degree)
return(os<<Gmpz(coef[0]));
for(int i=degree;i>=0;--i){
if(mpz_sgn(coef[i])){
if(printed&&(mpz_sgn(coef[i])==1))
os<<"+";
printed = true;
bool flag = false;
if((!mpz_cmp_si(coef[i],-1))&&i)
os<<"-";
else
if((mpz_cmp_ui(coef[i],1))||(!i)){
flag = true;
os<<Gmpz(coef[i]);
}
if(i){
if (flag)
os<<"*";
os<<"x";
if(i-1)
os<<"^"<<i;
}
}
}
if (!printed)
os<<'0';
return os;
}else{
int d=p.get_degree();
os<<"[ d="<<d<<" ";
os<<"[ ";
for(int i=0;i<d+1;++i)
os<<Gmpz(p.coef(i))<<" ";
os<<"] ]";
return os;
}
}
inline
std::istream& operator>>(std::istream &is,RS_polynomial_1 &pol){
std::istream::int_type c;
std::ios::fmtflags old_flags=is.flags();
is.unsetf(std::ios::skipws);
gmpz_eat_white_space(is);
if(is_pretty(is)){
std::string s;
while(((c=is.get())>='0'&&c<='9')||
c=='+'||c=='-'||c=='*'||c==' '||c=='x'||c=='^')
s.push_back(c);
is.putback(c);
pol=RS_polynomial_1(s);
goto ret_is;
}else{
int deg;
Gmpz z;
c=is.get(); // eat a '['
if(c!='[')
goto ret_fail_is;
gmpz_eat_white_space(is);
c=is.get(); // c='d'
if(c!='d')
goto ret_fail_is;
gmpz_eat_white_space(is);
c=is.get(); // c='='
if(c!='=')
goto ret_fail_is;
gmpz_eat_white_space(is);
c=is.peek();
deg=0;
while(c>='0'&&c<='9'){
c=is.get();
deg=10*deg+(c-'0');
c=is.peek();
}
pol=RS_polynomial_1(0);
pol.set_degree(deg);
gmpz_eat_white_space(is);
c=is.get();
if(c!='[')
goto ret_fail_is;
gmpz_eat_white_space(is);
for(int k=0;k<deg+1;++k){
is>>z;
pol.set_coef(k,z);
gmpz_eat_white_space(is);
}
c=is.get(); // ']'
if(c!=']')
goto ret_fail_is;
gmpz_eat_white_space(is);
c=is.get(); // ']'
if(c!=']')
goto ret_fail_is;
goto ret_is;
}
ret_fail_is:
is.setstate(std::ios_base::failbit);
ret_is:
is.flags(old_flags);
return is;
}
} // namespace CGAL
#endif // CGAL_RS_POLYNOMIAL_1_IO_H

279
Algebraic_kernel_d/include/CGAL/RS/polynomial_1_member.h
View File

@ -1,279 +0,0 @@
// Copyright (c) 2006-2009 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_POLYNOMIAL_1_MEMBER_H
#define CGAL_RS_POLYNOMIAL_1_MEMBER_H
namespace CGAL{
//////////////////////
// private functions
inline
void RS_polynomial_1::create_storage(int size){
coefficient_array=(mpz_t*)(*alloc_function)(sizeof(mpz_t)*size);
capacity=size;
}
inline
void RS_polynomial_1::free_storage(){
void *(*af)(size_t);
void *(*rf)(void*,size_t,size_t);
void (*ff)(void*,size_t);
mp_get_memory_functions(&af,&rf,&ff);
mp_set_memory_functions(alloc_function,realloc_function,free_function);
for(int i=0;i<polynomial_degree+1;++i)
mpz_clear(coef(i));
mp_set_memory_functions(af,rf,ff);
(*free_function)(coefficient_array,sizeof(mpz_t)*capacity);
capacity=0;
}
inline
void RS_polynomial_1::fetch_gmp_functions(){
mp_get_memory_functions(&alloc_function,
&realloc_function,
&free_function);
}
//////////////////////
// private functions
// dangerous function! this will erase the coefficients if the new coefficients
// don't fit in the previously allocated space; otherwise, it will set to zero
// the new coefficients
inline
void RS_polynomial_1::set_degree(int d){
if(d+1>capacity){
free_storage();
polynomial_degree=d;
++d;
create_storage(d);
for(int i=0;i<=polynomial_degree;++i)
mpz_init(coef(i));
}else{
polynomial_degree=d;
for(int i=polynomial_degree+1;i<=d;++i)
mpz_init(coef(i));
}
return;
}
inline
void RS_polynomial_1::force_degree(int n){
polynomial_degree=n;
}
// to change the storage capacity of a polynomial object
inline
int RS_polynomial_1::resize(int newcap){
if(newcap<=capacity)
return -1;
int i;
mpz_t *newcoef=(mpz_t*)(*alloc_function)(newcap*sizeof(mpz_t));
for(i=0;i<=polynomial_degree;++i){
mpz_init_set(newcoef[i],coefficient_array[i]);
mpz_clear(coefficient_array[i]);
}
for(i=polynomial_degree+1;i<newcap;++i)
mpz_init(newcoef[i]);
(*free_function)(coefficient_array,sizeof(mpz_t)*capacity);
coefficient_array=newcoef;
capacity=newcap;
return newcap;
}
inline
void RS_polynomial_1::set_coef(int pow_x,mpz_srcptr z){
mpz_set(coefficient_array[pow_x],z);
}
inline
void RS_polynomial_1::set_coef(int pow_x,const CGAL::Gmpz &z){
mpz_set(coefficient_array[pow_x],z.mpz());
}
inline
void RS_polynomial_1::set_coef_ui(int pow_x,unsigned long z){
mpz_set_ui(coefficient_array[pow_x],z);
}
inline
void RS_polynomial_1::set_coef_si(int pow_x,long z){
mpz_set_si(coefficient_array[pow_x],z);
}
inline
int RS_polynomial_1::get_degree()const{
while(!mpz_sgn(coef(polynomial_degree))&&polynomial_degree)
--polynomial_degree;
return polynomial_degree;
}
inline
int RS_polynomial_1::get_degree_static()const{
return polynomial_degree;
}
inline
bool RS_polynomial_1::has_sfpart()const{
return (is_square_free?true:(square_free_part.get()!=NULL?true:false));
}
inline
const RS_polynomial_1& RS_polynomial_1::sfpart()const{
if(is_square_free)
return *this;
else
return *square_free_part;
}
inline
void RS_polynomial_1::set_sfpart(RS_polynomial_1 *s)const{
is_square_free=false;
square_free_part=polyptr(s);
}
inline
void RS_polynomial_1::set_sfpart(const polyptr &s)const{
is_square_free=false;
square_free_part=s;
}
inline
void RS_polynomial_1::set_sf()const{
is_square_free=true;
}
inline
bool RS_polynomial_1::has_sqfr()const{
return (square_free_factorization.get()!=NULL?true:false);
}
inline
sqfrvec& RS_polynomial_1::sqfr()const{
return *square_free_factorization;
}
inline
void RS_polynomial_1::set_sqfr(sqfrvec *s)const{
square_free_factorization=sqfrptr(s);
}
inline
void RS_polynomial_1::set_sqfr(const sqfrptr &s)const{
square_free_factorization=s;
}
inline
mpz_ptr RS_polynomial_1::leading_coefficient()const{
return(coefficient_array[get_degree()]);
}
// gets the power of the lowest coefficient non-zero monomial
inline
int RS_polynomial_1::first_non_zero()const{
int i=0;
while(i<=polynomial_degree)
if(mpz_sgn(coef(i)))
return i;
else
++i;
return -1; // if all the coefficients are zero
}
inline
mpz_t* RS_polynomial_1::get_coefs()const{
return coefficient_array;
}
inline
mpz_ptr RS_polynomial_1::coef(int i)const{
return coefficient_array[i];
}
inline
RS_polynomial_1& RS_polynomial_1::derive()const{
int d=get_degree();
RS_polynomial_1 *derivative=new RS_polynomial_1(d-1);
for(int x=0;x<d;++x)
mpz_mul_si(derivative->coef(x),coef(x+1),x+1);
return *derivative;
}
inline
RS_polynomial_1& RS_polynomial_1::minusx()const{
int d=get_degree();
RS_polynomial_1 *mx=new RS_polynomial_1(d);
for(int x=0;x<=d;++x){
if(x%2==1)
mpz_set(mx->coef(x),coef(x));
else
mpz_neg(mx->coef(x),coef(x));
}
return *mx;
}
// This function copies into lb a number which is less than all roots. It just
// looks for the coefficient with the greatest absolute value. This is useful
// when evaluating curves in the infinite. We can refine this, but it makes no
// sense since it will be slower.
inline
void RS_polynomial_1::get_lower_bound(mpfr_ptr lb)const{
int d=get_degree();
mpz_t temp;
mpz_init(temp);
mpz_abs(temp,leading_coefficient());
for(int i=0;i<d;++i)
if(mpz_cmpabs(coef(i),temp)>0)
mpz_abs(temp,coef(i));
mpz_add_ui(temp,temp,1);
mpz_neg(temp,temp);
mpfr_set_z(lb,temp,GMP_RNDD);
mpz_clear(temp);
}
// the same that above, but the copied value is the upper bound of the roots
inline
void RS_polynomial_1::get_upper_bound(mpfr_ptr ub)const{
int d=get_degree();
mpz_t temp;
mpz_init(temp);
mpz_abs(temp,leading_coefficient());
for(int i=0;i<d;++i)
if(mpz_cmpabs(coef(i),temp)>0)
mpz_abs(temp,coef(i));
mpz_add_ui(temp,temp,1);
mpfr_set_z(ub,temp,GMP_RNDU);
mpz_clear(temp);
}
// returns this polynomial times c*x^p
inline
RS_polynomial_1
RS_polynomial_1::times_monomial(mpz_srcptr c,int p)const{
int dr=get_degree()+p;
RS_polynomial_1 r(dr);
for(int i=p;i<=dr;++i)
mpz_mul(r.coef(i),coef(i-p),c);
return r;
}
} // namespace CGAL
#endif // CGAL_RS_POLYNOMIAL_1_MEMBER_H

178
Algebraic_kernel_d/include/CGAL/RS/polynomial_1_operators.h
View File

@ -1,178 +0,0 @@
// 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_POLYNOMIAL_1_OPERATORS_H
#define CGAL_RS_POLYNOMIAL_1_OPERATORS_H
namespace CGAL{
inline
double RS_polynomial_1::operator()(double x)const{
double result=mpz_get_d(leading_coefficient());
int d=get_degree_static();
for(int i=1;i<d+1;++i)
result=x*result+mpz_get_d(coef(d-i));
return result;
}
inline
Gmpz RS_polynomial_1::operator()(int x)const{
Gmpz result(mpz_get_d(leading_coefficient()));
int d=get_degree_static();
for(int i=1;i<d+1;++i)
result=x*result+Gmpz(coef(d-i));
return result;
}
inline
RS_polynomial_1 RS_polynomial_1::operator-()const{
RS_polynomial_1 opposite(polynomial_degree);
int d=get_degree()+1;
for(int i=0;i<d;++i)
mpz_neg(opposite.coef(i),coef(i));
return opposite;
}
inline
RS_polynomial_1& RS_polynomial_1::operator+=(const RS_polynomial_1 &s){
is_square_free=false;
square_free_part.reset();
square_free_factorization.reset();
int sd;
if(polynomial_degree<(sd=s.get_degree())){
mpz_t *coef_sum=(mpz_t*)(*alloc_function)(sizeof(mpz_t)*(sd+1));
for(int i=0;i<polynomial_degree+1;++i){
mpz_init(coef_sum[i]);
mpz_add(coef_sum[i],s.coef(i),coef(i));
mpz_clear(coef(i));
}
for(int i=polynomial_degree+1;i<sd+1;++i)
mpz_init_set(coef_sum[i],s.coef(i));
(*free_function)(coefficient_array,sizeof(mpz_t)*capacity);
coefficient_array=coef_sum;
polynomial_degree=sd;
}else
for(int i=0;i<sd+1;++i)
mpz_add(coef(i),s.coef(i),coef(i));
return *this;
}
inline
RS_polynomial_1& RS_polynomial_1::operator-=(const RS_polynomial_1 &s){
is_square_free=false;
square_free_part.reset();
square_free_factorization.reset();
int sd;
if(polynomial_degree<(sd=s.get_degree())){
mpz_t *coef_sum=(mpz_t*)(*alloc_function)(sizeof(mpz_t)*(sd+1));
for(int i=0;i<polynomial_degree+1;++i){
mpz_init(coef_sum[i]);
mpz_sub(coef_sum[i],coef(i),s.coef(i));
mpz_clear(coef(i));
}
for(int i=polynomial_degree+1;i<sd+1;++i){
mpz_init(coef_sum[i]);
mpz_neg(coef_sum[i],s.coef(i));
}
(*free_function)(coefficient_array,sizeof(mpz_t)*capacity);
coefficient_array=coef_sum;
polynomial_degree=sd;
}else
for(int i=0;i<sd+1;++i)
mpz_sub(coef(i),coef(i),s.coef(i));
return *this;
}
inline
RS_polynomial_1 RS_polynomial_1::operator*(const RS_polynomial_1 &f)const{
int d=get_degree();
int degree_f=f.get_degree();
int degree_p=d+degree_f;
RS_polynomial_1 product(degree_p);
for(int c=0;c<degree_p+1;++c){
int max=(c<d?c:d)+1;
for(int i=0;i<max;++i)
if(c-i<=degree_f)
mpz_addmul(product.coef(c),coef(i),f.coef(c-i));
}
return product;
}
inline
RS_polynomial_1& RS_polynomial_1::operator*=(const RS_polynomial_1 &f){
is_square_free=false;
square_free_part.reset();
square_free_factorization.reset();
RS_polynomial_1 aux(*this);
*this=aux*f;
return *this;
}
inline
RS_polynomial_1& RS_polynomial_1::operator*=(mpz_srcptr s){
for(int i=0;i<polynomial_degree+1;++i)
mpz_mul(coef(i),coef(i),s);
return *this;
}
inline
RS_polynomial_1& RS_polynomial_1::operator*=(const CGAL::Gmpz &s){
for(int i=0;i<polynomial_degree+1;++i)
mpz_mul(coef(i),coef(i),s.mpz());
return *this;
}
inline
RS_polynomial_1& RS_polynomial_1::operator/=(mpz_srcptr d){
for(int i=0;i<polynomial_degree+1;++i)
mpz_divexact(coef(i),coef(i),d);
return *this;
}
inline
RS_polynomial_1& RS_polynomial_1::operator/=(const CGAL::Gmpz &d){
for(int i=0;i<polynomial_degree+1;++i)
mpz_divexact(coef(i),coef(i),d.mpz());
return *this;
}
inline
bool RS_polynomial_1::operator==(const RS_polynomial_1 &p)const{
int d=get_degree();
int degree_p=p.get_degree();
if(degree_p<d){
for(int i=degree_p+1;i<d+1;++i)
if(mpz_sgn(coef(i)))
return false;
for(int i=0;i<degree_p+1;++i)
if(mpz_cmp(coef(i),p.coef(i)))
return false;
}else{
for (int i=d+1;i<degree_p+1;++i)
if(mpz_sgn(p.coef(i)))
return false;
for(int i=0;i<degree_p+1;++i)
if(mpz_cmp(coef(i),p.coef(i)))
return false;
}
return true;
}
} // namespace CGAL
#endif // CGAL_RS_POLYNOMIAL_1_OPERATORS_H

361
Algebraic_kernel_d/include/CGAL/RS/polynomial_1_parser.h
View File

@ -1,361 +0,0 @@
// Copyright (c) 2007-2010 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>
//
// (the first version of this file was written by Elias Tsigaridas)
#ifndef CGAL_RS_PARSER_1_H
#define CGAL_RS_PARSER_1_H
#include <iostream>
#include <string>
#include <algorithm>
// if boost version is 1.38 or newer, we use the new version of spirit
#include <boost/version.hpp>
#if BOOST_VERSION >= 103800
#include <boost/spirit/include/classic_core.hpp>
#define CGAL_BOOST_SPIRIT boost::spirit::classic
#else
#include <boost/spirit/core.hpp>
#define CGAL_BOOST_SPIRIT boost::spirit
#endif
namespace CGAL{
// The semantics.
// In this class we store the current pair and the whole result;
struct Semantics
{
typedef std::pair< int, std::string > pair_t;
std::string variable;
pair_t current;
std::vector< pair_t > result;
std::string sign;
Semantics( const std::string& var) : variable( var) { }
void set_variable( const std::string& var)
{
variable = var;
}
std::string get_variable() const
{
return variable;
}
bool is_first() const
{
return current == pair_t(-1, "F");
}
void add()
{
result.push_back( current);
current = pair_t(-1, "F");
}
void init()
{
result.clear();
current = pair_t(-1, "F");
sign = "";
}
};
//------------- Semantic actions -----------------------------------
// Add a new monomial at the result
struct AddMonomial
{
void operator()( char const* first, char const* last) const
{
if (!sem.is_first()) {
sem.add();
}
}
AddMonomial( Semantics& sem_) : sem( sem_) { }
Semantics& sem;
};
// Set the coefficient of the current monomial
struct SetCoeff
{
void operator() ( char const* first, char const* last) const
{
sem.current.second = std::string( first, last);
if (sem.current.first == -1) {
sem.current.first = 0;
}
}
SetCoeff( Semantics& sem_) : sem( sem_) { }
Semantics& sem;
};
// Set the exponent of the current monomial
struct SetExp
{
void
operator()( const char* first, const char* last) const
{
sem.current.first = 1;
if (sem.current.second == "F") sem.current.second = "1";
}
void operator()( unsigned num) const
{
sem.current.first = num;
if (sem.current.second == "F") sem.current.second = "1";
}
SetExp( Semantics& sem_) : sem( sem_) { }
Semantics& sem;
};
// Set the sign of the current monomial
struct SetSign
{
void
operator()( char const* first, char const* last) const
{
std::string str( first, last);
sem.sign = str;
}
SetSign( Semantics& sem_) : sem( sem_) { }
Semantics& sem;
};
// Adjust the sign of the coefficient of the current monomial
struct AdjustCoeff
{
void
operator()( char const* first, char const* last) const
{
if (sem.sign == "-") {
if (sem.current.second[0] == '-') {
sem.current.second[0] = '+';
} else if (sem.current.second[0] == '+') {
sem.current.second[0] = '-';
} else {
sem.current.second = '-' + sem.current.second;
}
}
}
AdjustCoeff( Semantics& sem_) : sem( sem_) { }
Semantics& sem;
};
//---------------- Grammar -------------------------------------
// The grammar of the polynomial parser.
// The coefficient of the polynomial can also be rationals.
struct polynomial_p : public CGAL_BOOST_SPIRIT::grammar<polynomial_p>
{
template < typename ScannerT >
struct definition
{
definition(polynomial_p const& self)
{
CGAL_BOOST_SPIRIT::chlit<> lpar('(');
CGAL_BOOST_SPIRIT::chlit<> rpar(')');
CGAL_BOOST_SPIRIT::chlit<> plus('+');
CGAL_BOOST_SPIRIT::chlit<> minus('-');
CGAL_BOOST_SPIRIT::chlit<> mul('*');
CGAL_BOOST_SPIRIT::strlit<>
varX(self.sem.get_variable().c_str());
poly = (
(monomial | smonomial) [ AddMonomial( self.sem) ] >>
*( smonomial [ AddMonomial( self.sem) ] )
);
monomial
= ( unumber [ SetCoeff( self.sem) ] >> !( mul >> X ) )
| ( lpar >> number[ SetCoeff( self.sem) ] >> rpar >> !( mul >> X) )
| ( X >> mul >> number [ SetCoeff( self.sem) ] )
| ( X >> mul >> lpar >> number [ SetCoeff( self.sem) ] >> rpar )
| X;
smonomial = ( ( plus | minus ) [ SetSign( self.sem) ] >>
monomial) [ AdjustCoeff( self.sem) ];
X = ( (varX)[ SetExp( self.sem) ] >>
!(
CGAL_BOOST_SPIRIT::ch_p('^') >>
( (CGAL_BOOST_SPIRIT::uint_p) [ SetExp( self.sem) ]
| (lpar >>
CGAL_BOOST_SPIRIT::uint_p [ SetExp( self.sem) ] >>
rpar))
)
);
unumber = +(CGAL_BOOST_SPIRIT::digit_p) >>
!(CGAL_BOOST_SPIRIT::ch_p("/") >>
+(CGAL_BOOST_SPIRIT::digit_p));
number = ( (plus | minus) >> unumber );
}
CGAL_BOOST_SPIRIT::rule<ScannerT>
poly, monomial, smonomial, X, unumber, number;
CGAL_BOOST_SPIRIT::rule<ScannerT> const& start() const
{
return poly;
}
};
polynomial_p( Semantics& sem_) : sem(sem_) { }
Semantics& sem;
};
// The polynomial parser.
// The constructor takes the name of the variable.
struct Polynomial_parser_1
{
protected:
// internal::Semantics sem;
Semantics sem;
// internal::polynomial_p poly_p;
polynomial_p poly_p;
CGAL_BOOST_SPIRIT::parse_info<> info;
// The default conversion is to int
struct Converter
{
int operator()( const std::string str) const
{
return atoi( str.c_str());
}
};
public:
Polynomial_parser_1( const std::string& var = "x")
: sem( var), poly_p( sem)
{
}
void parse( const std::string& in_str)
{
// Eliminate spaces
std::string str(in_str);
while (str.find(" ") < str.size() )
{
size_t pos = str.find(" ");
str.erase( pos, 1);
}
sem.init();
info = CGAL_BOOST_SPIRIT::parse(str.c_str(),
poly_p,
CGAL_BOOST_SPIRIT::space_p);
std::sort( sem.result.begin(), sem.result.end());
}
bool is_correct() const
{
return info.full;
}
std::string error() const
{
if (!is_correct()) {
return info.stop;
}
return "";
}
std::string get_variable() const
{
return sem.get_variable();
}
void set_variable( const std::string& var)
{
sem.set_variable( var);
}
template < typename OutputIterator,
typename ConverterFunction>
OutputIterator
result( OutputIterator oi, ConverterFunction conv) const
{
//for (int i = 0; i < sem.result.size(); ++i) {
// std::cout << "(" << sem.result[i].first << ", " <<
// sem.result[i].second << ") ";
//}
//std::cout << std::endl;
for (unsigned i = 0, j = 0; j < sem.result.size(); ++i) {
if (sem.result[j].first == static_cast<int>( i)) {
*oi++ = conv( sem.result[j].second);
++j;
} else {
*oi++ = conv( "0");
}
}
return oi;
}
template < typename OutputIterator >
OutputIterator
result( OutputIterator oi) const
{
return this->result( oi, Converter());
}
};
template < typename T >
struct Convert_to
{
T
operator()( const std::string str) const {
return static_cast<T>( atoi( str.c_str()));
}
};
struct Convert_to_Gmpz
{
Gmpz
operator()( const std::string str) const {
return Gmpz(str.c_str());
}
};
} // namespace CGAL
#undef CGAL_BOOST_SPIRIT
#endif // CGAL_RS_PARSER_1_H

408
Algebraic_kernel_d/include/CGAL/RS/polynomial_1_utils.h
View File

@ -1,408 +0,0 @@
// Copyright (c) 2007-2009 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_POLYNOMIAL_1_UTILS_H
#define CGAL_RS_POLYNOMIAL_1_UTILS_H
#include <CGAL/basic.h>
#include <gmp.h>
#include <CGAL/RS/polynomial_1.h>
#include <CGAL/RS/polynomial_converter.h>
#include <CGAL/RS/solve_1.h>
#ifdef CGAL_RS_USE_UGCD
#include <CGAL/RS/ugcd/ugcd.h>
#endif
#include <rs_exports.h>
#ifdef CGAL_USE_RS3
#include <rs3_fncts.h>
#endif
namespace CGAL{
#ifdef CGAL_USE_RS3
// Fabrice's modular gcd algorithm
struct Rsgcd_1:
public std::binary_function<RS_polynomial_1,RS_polynomial_1,RS_polynomial_1>{
RS_polynomial_1& operator()(
const RS_polynomial_1 &p1,
const RS_polynomial_1 &p2)const{
int dr,d1,d2;
mpz_t * r_z;
d1=p1.get_degree();
d2=p2.get_degree();
dr=rs3_up_mz_gcd(
&r_z,
(const mpz_t*)p1.get_coefs(),
d1,
(const mpz_t*)p2.get_coefs(),
d2);
RS_polynomial_1 *result=new RS_polynomial_1(&r_z,dr);
return *result;
}
};
#endif
struct Cgalgcd_1:
public std::binary_function<RS_polynomial_1,RS_polynomial_1,RS_polynomial_1>{
RS_polynomial_1& operator()(
const RS_polynomial_1 &p1,
const RS_polynomial_1 &p2)const{
typedef Polynomial<Gmpz> P;
typedef from_rs_poly<P> convert;
typedef to_rs_poly<P> convertback;
return convertback()(CGAL::gcd(convert()(p1),convert()(p2)));
}
};
#ifdef CGAL_RS_USE_UGCD
// my modular gcd algorithm
struct Modgcd_1:
public std::binary_function<RS_polynomial_1,RS_polynomial_1,RS_polynomial_1>{
RS_polynomial_1& operator()(
const RS_polynomial_1 &u,
const RS_polynomial_1 &v)const{
RS_polynomial_1 *result=new RS_polynomial_1(
u.get_degree()<v.get_degree()?
u.get_degree_static():
v.get_degree_static());
int dr=
RS_MGCD::Ugcd::ugcd(
result->get_coefs(),
u.get_coefs(),
u.get_degree_static(),
v.get_coefs(),
v.get_degree_static());
result->force_degree(dr);
return *result;
}
};
#endif // CGAL_RS_USE_UGCD
// Cont()(c,u) stores in c the gcd of the coefficients of u
struct Cont:
public std::binary_function<mpz_ptr,RS_polynomial_1,void>{
void operator()(mpz_ptr c,const RS_polynomial_1 &u){
mpz_set(c,u.coef(0));
for(int i=1;i<u.get_degree_static()+1;++i)
mpz_gcd(c,c,u.coef(i));
}
};
// Pp()(u)=u/Cont()(u)
struct Pp:
public std::unary_function<RS_polynomial_1,RS_polynomial_1>{
RS_polynomial_1& operator()(const RS_polynomial_1 &u){
mpz_t c;
mpz_init(c);
Cont()(c,u);
RS_polynomial_1 *res=new RS_polynomial_1(u/c);
mpz_clear(c);
return *res;
}
};
// pseudo division; returns <q,r>, where:
// d^l*f=qg+r, d=leadingcoef(f), l=max(deg(f)-deg(g)+1,0)
struct Pdiv_1:
public std::binary_function<
RS_polynomial_1,
RS_polynomial_1,
std::pair<RS_polynomial_1,RS_polynomial_1> >{
std::pair<RS_polynomial_1,RS_polynomial_1>
operator()(const RS_polynomial_1 &f,const RS_polynomial_1 &g){
int degf,degg,lambda,i;
mpz_t d;
mpz_ptr lcg;
degf=f.get_degree();
degg=g.get_degree();
RS_polynomial_1 q(degf);
RS_polynomial_1 r(f);
lcg=g.leading_coefficient();
lambda=degf-degg+1;
if(lambda>0){
mpz_init(d);
mpz_pow_ui(d,g.leading_coefficient(),lambda);
r*=d;
mpz_clear(d);
}
if(!degg){
for(int i=0;i<=f.get_degree_static();++i)
mpz_divexact(q.coef(i),r.coef(i),lcg);
return std::make_pair(q,RS_polynomial_1());
}
// don't use get_degree_static
while((i=r.get_degree()-degg)>=0){
mpz_divexact(q.coef(i),r.leading_coefficient(),lcg);
r-=g.times_monomial(q.coef(i),i);
}
return std::make_pair(q,r);
}
};
struct Pdivrem_1:
public std::binary_function<RS_polynomial_1,RS_polynomial_1,RS_polynomial_1>{
RS_polynomial_1
operator()(const RS_polynomial_1 &f,const RS_polynomial_1 &g){
return Pdiv_1()(f,g).second;
}
};
struct Pdivquo_1:
public std::binary_function<RS_polynomial_1,RS_polynomial_1,RS_polynomial_1>{
RS_polynomial_1
operator()(const RS_polynomial_1 &f,const RS_polynomial_1 &g){
return Pdiv_1()(f,g).first;
}
};
// subresultant GCD algorithm:
// Knuth, TAOCP vol. 2, 2nd edition, section 4.6.1, algorithm C
struct Gcd_1:
public std::binary_function<RS_polynomial_1,RS_polynomial_1,RS_polynomial_1>{
RS_polynomial_1& operator()(
const RS_polynomial_1 &uu,
const RS_polynomial_1 &vv)const{
RS_polynomial_1 u,r;
RS_polynomial_1 *v=new RS_polynomial_1(vv.get_degree());
mpz_t contu,g,d,h,den_h;
int delta;
// C1.
mpz_init(contu);
Cont()(contu,uu);
u=uu/contu;
mpz_init(g);
Cont()(g,vv); // g will store cont(v)
*v=vv/g;
mpz_init(d);
mpz_gcd(d,contu,g);
// we don't need cont(u) and cont(v) anymore; so we will use the
// variables which stored them (we will use contu as a temp and g to
// store the value g)
mpz_set_ui(g,1); // now, g will store the g from the algorithm
mpz_init_set_ui(h,1);
mpz_init(den_h);
// C2.
while((r=Pdivrem_1()(u,*v)).get_degree_static()){
delta=u.get_degree()-v->get_degree();
// C3.
u=*v;
if(delta<0){
// v := r(x) * (gh^(-delta))
mpz_pow_ui(contu,h,(unsigned)(-delta));
mpz_mul(contu,g,contu);
*v=r*contu;
}else{
// v := r(x) / (gh^delta)
mpz_pow_ui(contu,h,(unsigned)delta);
mpz_mul(contu,g,contu);
*v=r/contu;
}
mpz_set(g,u.leading_coefficient());
if(delta){ // if delta=0, h remains unchanged
if(delta==1){ // if delta=1, h:=g
mpz_set(h,g);
}else{ // otherwise, h:=h^(1-delta)*g^delta
if(delta>1){
mpz_pow_ui(den_h,h,(unsigned)(delta-1));
mpz_pow_ui(h,g,(unsigned)delta);
mpz_divexact(h,h,den_h);
}else{ // delta<0
mpz_pow_ui(h,h,(unsigned)(1-delta));
mpz_pow_ui(den_h,g,(unsigned)(-delta));
mpz_divexact(h,h,den_h);
}
}
}
}
mpz_clear(contu);
mpz_clear(g);
mpz_clear(h);
mpz_clear(den_h);
if(mpz_sgn(r.coef(0))){
v->force_degree(0);
// v(x)=1, so d*Pp()(v(x))=d
v->set_coef(0,d);
}else{
*v=Pp()(*v);
*v*=d;
}
mpz_clear(d);
return *v;
}
};
// exact division: the result is undefined when the division is not exact
struct Ediv_1:
public std::binary_function<RS_polynomial_1,RS_polynomial_1,RS_polynomial_1*>{
RS_polynomial_1*
operator()(const RS_polynomial_1 &f,const RS_polynomial_1 &g){
//int degf,degg,i;
//mpz_ptr lcg;
//mpz_t r;
//degf=f.get_degree();
//degg=g.get_degree();
//RS_polynomial_1 *q=new RS_polynomial_1(degf-degg);
//lcg=g.leading_coefficient();
//if(!degg){
// for(int i=0;i<=degf;++i)
// mpz_divexact(q->coef(i),f.coef(i),lcg);
// return q;
//}
//mpz_init(r);
//mpz_set(r,f.leading_coefficient());
//std::cout<<"f="<<f<<std::endl;
//std::cout<<"g="<<g<<std::endl;
//for(i=degf-degg;i>0;--i){
// //std::cout<<"\ni="<<i;
// mpz_divexact(q->coef(i),r/ *f.coef(i+degg)* /,lcg);
// //--------------------------------------------------
// // mpz_mul(r,q->coef(i),lcg);
// // mpz_sub(r,f.coef(i+degf-degg-1),r);
// //--------------------------------------------------
// mpz_mul(r,q->coef(i),g.coef(degg-1));
// mpz_sub(r,f.coef(i+degf-degg-1),r);
// std::cout<<"i="<<i<<", q[i]="<<q->coef(i)<<", r="<<Gmpz(r)<<std::endl;
//}
//mpz_divexact(q->coef(0),r/ *f.coef(degg)* /,lcg);
//mpz_clear(r);
////std::cout<<"\nediv("<<f<<","<<g<<") = "<<(*q)<<std::endl;
//return q;
// naive implementation
RS_polynomial_1 *ret=new RS_polynomial_1(Pdivquo_1()(f,g));
return ret;
}
};
struct Constantpoly_1:
public std::unary_function<long,RS_polynomial_1>{
RS_polynomial_1&
operator()(long i){
RS_polynomial_1 *cpoly=new RS_polynomial_1(0);
cpoly->set_coef_si(0,i);
return *cpoly;
}
};
// squarefree factorization, Yun's algorithm (1976)
// if P=sum(P_i^i), this function returns a vector with pairs (P_i,i)
// it gives for free the square-free part of P, which is C_1
template<class GcdPolicy>
struct do_sqfr_1:
public std::unary_function<RS_polynomial_1,sqfrvec*>{
typedef GcdPolicy Gcd;
sqfrvec* operator()(const RS_polynomial_1 &P)const{
sqfrvec *res=new sqfrvec();
if(!P.get_degree()){
res->push_back(
std::make_pair(Constantpoly_1()(mpz_sgn(P.coef(0))?1:0),1));
return res;
}
RS_polynomial_1 dP(P.derive());
RS_polynomial_1 G(Pp()(Gcd()(P,dP)));
if(!G.get_degree()){
res->push_back(std::make_pair(P,1));
return res;
}
if(!P.has_sfpart()){
polyptr pp(Ediv_1()(P,G));
P.set_sfpart(pp);
}
RS_polynomial_1 C_i(P.sfpart()); // C_1 is P/G, i.e. the sf part
RS_polynomial_1 D_i(*Ediv_1()(dP,G)-(C_i.derive()));//D_1=dP/G-dC_1
RS_polynomial_1 P_i;
for(int i=1;;++i){
if(D_i.get_degree_static()||mpz_sgn(D_i.coef(0)))
P_i=Pp()(Gcd()(C_i,D_i));
else
P_i=C_i;
if(P_i.get_degree_static())
res->push_back(std::make_pair(P_i,i));
C_i=*Ediv_1()(C_i,P_i); // C_{i+1}=C_i/P_i
if(!C_i.get_degree_static())
return res;
D_i=*Ediv_1()(D_i,P_i)-(C_i.derive()); // D_{i+1}=D_i/P_i-dC_{i+1}
}
}
};
template<class GcdPolicy>
struct sqfr_1:
public std::unary_function<RS_polynomial_1,sqfrvec>{
typedef GcdPolicy Gcd;
sqfrvec operator()(const RS_polynomial_1 &P){
if(!P.has_sqfr()){
sqfrptr sp(do_sqfr_1<Gcd>()(P));
P.set_sqfr(sp);
}
return P.sqfr();
}
};
// the square-free part of P is P/gcd(P,dP), but it can be also calculated
// as the product of its sf factors
template<class GcdPolicy>
struct do_sfpart_1:
public std::unary_function<RS_polynomial_1,RS_polynomial_1*>{
typedef GcdPolicy Gcd;
RS_polynomial_1* operator()(const RS_polynomial_1 &P)const{
if(P.get_degree()){
// TODO: not optimal
return &Pp()(*Ediv_1()(P,Pp()(Gcd()(P,P.derive()))));
}else
return (new RS_polynomial_1(
Constantpoly_1()(mpz_sgn(P.coef(0))?1:0)));
}
};
#ifdef CGAL_USE_RS3
template<>
struct do_sfpart_1<Rsgcd_1>:
public std::unary_function<RS_polynomial_1,RS_polynomial_1*>{
RS_polynomial_1* operator()(const RS_polynomial_1 &P)const{
if(P.get_degree()){
int d_sfp;
mpz_t* sfp_z;
d_sfp=rs3_up_mz_sfp(
&sfp_z,
(const mpz_t*)P.get_coefs(),
P.get_degree());
RS_polynomial_1 *result=new RS_polynomial_1(&sfp_z,d_sfp);
return result;
}else
return (new RS_polynomial_1(
Constantpoly_1()(mpz_sgn(P.coef(0))?1:0)));
}
};
#endif
template<class GcdPolicy>
struct sfpart_1:
public std::unary_function<RS_polynomial_1,RS_polynomial_1*>{
typedef GcdPolicy Gcd;
const RS_polynomial_1& operator()(const RS_polynomial_1 &P){
if(!P.has_sfpart()){
polyptr pp(do_sfpart_1<Gcd>()(P));
P.set_sfpart(pp);
}
return P.sfpart();
}
};
} // namespace CGAL
#endif // CGAL_RS_POLYNOMIAL_1_UTILS_H

178
Algebraic_kernel_d/include/CGAL/RS/polynomial_converter.h
View File

@ -1,178 +0,0 @@
// Copyright (c) 2009 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_POLYNOMIAL_CONVERTER
#define CGAL_RS_POLYNOMIAL_CONVERTER
#include <CGAL/Polynomial.h>
#include <CGAL/RS/polynomial_1.h>
namespace CGAL{
template<class Polynomial_>
struct to_rs_poly:
public std::unary_function<Polynomial_,RS_polynomial_1>{
typedef Polynomial_ Polynomial;
RS_polynomial_1 operator()(const Polynomial &p)const{
std::cerr<<"can't convert to integer polynomial"<<std::endl;
exit(-1);
}
};
template<>
struct to_rs_poly<RS_polynomial_1>:
public std::unary_function<RS_polynomial_1,RS_polynomial_1>{
const RS_polynomial_1& operator()(const RS_polynomial_1 &p)const{
return p;
}
};
// The conversions using this macro are efficient, since this construct an
// RS polynomial only by copying pointers. The only implementation detail
// is that it should not free the occuped by the pointed coefficients.
#define CGALRS_POLYNOMIAL_CONVERTER_REF(TYPE,CONVERT_FUNCTION) \
template<> \
struct to_rs_poly<Polynomial<TYPE> >: \
public std::unary_function<Polynomial<TYPE>,RS_polynomial_1>{ \
RS_polynomial_1& operator()(const Polynomial<TYPE> &p)const{ \
void *(*af)(size_t); \
void *(*rf)(void*,size_t,size_t); \
void (*ff)(void*,size_t); \
int d=p.degree(); \
mpz_t* c=(mpz_t*)malloc((d+1)*sizeof(mpz_t)); \
for(int i=0;i<=d;++i) \
CONVERT_FUNCTION; \
mp_get_memory_functions(&af,&rf,&ff); \
mp_set_memory_functions(af,rf,cgalrs_dummy_free); \
RS_polynomial_1 *r=new RS_polynomial_1(&c,d); \
mp_set_memory_functions(af,rf,ff); \
return *r; \
} \
}
// The conversions using this macro are not intended to be efficient, since
// there is no direct way to convert from these types to mpz_t and we need
// thus to create new mpz_t's.
#define CGALRS_POLYNOMIAL_CONVERTER_COPY(TYPE,CONVERT_FUNCTION) \
template<> \
struct to_rs_poly<Polynomial<TYPE> >: \
public std::unary_function<Polynomial<TYPE>,RS_polynomial_1>{ \
RS_polynomial_1& operator()(const Polynomial<TYPE> &p)const{ \
int d=p.degree(); \
mpz_t* c=(mpz_t*)malloc((d+1)*sizeof(mpz_t)); \
for(int i=0;i<=d;++i){ \
mpz_init(c[i]); \
CONVERT_FUNCTION; \
} \
return *(new RS_polynomial_1(&c,d)); \
} \
}
//CGALRS_POLYNOMIAL_CONVERTER_REF(Gmpz,c[i][0]=*(p[i].mpz()));
CGALRS_POLYNOMIAL_CONVERTER_COPY(Gmpz,mpz_set(c[i],p[i].mpz()));
CGALRS_POLYNOMIAL_CONVERTER_COPY(int,mpz_set_si(c[i],(long)p[i]));
CGALRS_POLYNOMIAL_CONVERTER_COPY(long,mpz_set_si(c[i],p[i]));
CGALRS_POLYNOMIAL_CONVERTER_COPY(unsigned,mpz_set_ui(c[i],(unsigned long)p[i]));
CGALRS_POLYNOMIAL_CONVERTER_COPY(unsigned long,mpz_set_ui(c[i],p[i]));
#undef CGALRS_POLYNOMIAL_CONVERTER_REF
#undef CGALRS_POLYNOMIAL_CONVERTER_COPY
// convert a Gmpz rational polynomial to an integer one
template<>
struct to_rs_poly<Polynomial<Gmpq> >:
public std::unary_function<Polynomial<Gmpq>,RS_polynomial_1>{
RS_polynomial_1& operator()(const Polynomial<Gmpq> &p)const{
int d=p.degree();
mpz_t denominator;
mpz_init(denominator);
mpz_lcm(denominator,
mpq_denref(p[0].mpq()),
mpq_denref(p[d].mpq()));
for(int j=1;j<d;++j)
mpz_lcm(denominator,
denominator,
mpq_denref(p[j].mpq()));
mpz_t* c=(mpz_t*)malloc((d+1)*sizeof(mpz_t));
for(int i=0;i<=d;++i){
mpz_init(c[i]);
mpz_div(c[i],denominator,mpq_denref(p[i].mpq()));
mpz_mul(c[i],c[i],mpq_numref(p[i].mpq()));
}
mpz_clear(denominator);
return *(new RS_polynomial_1(&c,d));
}
};
template<class Polynomial_>
struct from_rs_poly:
public std::unary_function<RS_polynomial_1,Polynomial_>{
typedef Polynomial_ Polynomial;
Polynomial operator()(const RS_polynomial_1 &p)const{
std::cerr<<"can't convert to integer polynomial"<<std::endl;
exit(-1);
}
};
template<>
struct from_rs_poly<RS_polynomial_1>:
public std::unary_function<RS_polynomial_1,RS_polynomial_1>{
const RS_polynomial_1& operator()(const RS_polynomial_1 &p)const{
return p;
}
};
template<>
struct from_rs_poly<Polynomial<Gmpz> >:
public std::unary_function<RS_polynomial_1,Polynomial<Gmpz> >{
Polynomial<Gmpz> operator()(const RS_polynomial_1 &p)const{
typedef Polynomial_traits_d<Polynomial<Gmpz> > PT;
mpz_t* pcoef=p.get_coefs();
std::vector<Gmpz> coeffs;
int d=p.get_degree();
for(int i=0;i<=d;++i){
// Gmpz c(1);
// *c.mpz()=*pcoef[i];
Gmpz c(pcoef[i]);
coeffs.push_back(c);
}
return PT::Construct_polynomial()(coeffs.begin(),coeffs.end());
}
};
template<>
struct from_rs_poly<Polynomial<Gmpq> >:
public std::unary_function<RS_polynomial_1,Polynomial<Gmpq> >{
Polynomial<Gmpq> operator()(const RS_polynomial_1 &p)const{
typedef Polynomial_traits_d<Polynomial<Gmpq> > PT;
mpz_t* pcoef=p.get_coefs();
std::vector<Gmpq> coeffs;
int d=p.get_degree();
for(int i=0;i<=d;++i){
// Gmpq c(1);
// *mpq_numref(c.mpq())=*pcoef[i];
Gmpq c(pcoef[i]);
coeffs.push_back(c);
}
return PT::Construct_polynomial()(coeffs.begin(),coeffs.end());
}
};
} // namespace CGAL
#endif // CGAL_RS_POLYNOMIAL_CONVERTER

485
Algebraic_kernel_d/include/CGAL/RS/refine_1.h
View File

@ -1,485 +0,0 @@
// 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 GcdPolicy>
Comparison_result bisect_at(const Algebraic_1 &a,mpfr_srcptr p){
typedef GcdPolicy 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 GcdPolicy>
int bisect_at_endpoints(const Algebraic_1 &a,Algebraic_1 &b){
typedef GcdPolicy 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 GcdPolicy>
int bisect_n(const Algebraic_1 &a,unsigned long n=1){
typedef GcdPolicy 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 GcdPolicy>
int bisect_n_dyadic(const Algebraic_1 &a,unsigned long n=1){
typedef GcdPolicy 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 GcdPolicy>
int bisect(const Algebraic_1 &a,mp_exp_t s){
typedef GcdPolicy 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 GcdPolicy>
int qir_n(const Algebraic_1 &a,unsigned t=1){
typedef GcdPolicy 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 GcdPolicy>
int qir(const Algebraic_1 &a,mp_exp_t t){
typedef GcdPolicy 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

75
Algebraic_kernel_d/include/CGAL/RS/refine_1_rs.h
View File

@ -1,75 +0,0 @@
// Copyright (c) 2009 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_RS_H
#define CGAL_RS_REFINE_1_RS_H
#include <gmp.h>
#include <mpfr.h>
#include <mpfi.h>
#include <CGAL/RS/algebraic_1.h>
#include <CGAL/assertions.h>
namespace RS3{
inline void refine_1(const CGAL::Algebraic_1 &a,unsigned int s=10000){
CGAL_precondition(a.inf()<=a.sup());
// If the algebraic endpoints have room for a refinement bigger
// than desired, refine as much is possible. This would ensure that
// small refinements are never made. Other possibility to avoid
// small refinements is to refine the number just after it was
// isolated (this early refinement would also have the advantage
// that one can choose to refine until reaching a certain criterion
// and assume it is true for later refinements, but this can slow
// down the Solve_1 functor). And a last option would be to
// determine a lower bound on the precisions RS likes to refine.
if(s<mpfr_get_prec(&((mpfi_ptr)(a.mpfi()))->left)-2)
s=mpfr_get_prec(&((mpfi_ptr)(a.mpfi()))->left)-2;
if(s<mpfr_get_prec(&((mpfi_ptr)(a.mpfi()))->right)-2)
s=mpfr_get_prec(&((mpfi_ptr)(a.mpfi()))->right)-2;
// If the precision of the endpoints is not enough for the desired
// refinement, allocate a new mpfi and swap later the result.
if(mpfr_get_prec(&((mpfi_ptr)(a.mpfi()))->left)<s+2||
mpfr_get_prec(&((mpfi_ptr)(a.mpfi()))->right)<s+2){
mpfi_t n;
mpfi_init2(n,s+2);
mpfi_set(n,a.mpfi());
rs3_refine_u_root(n,
a.pol().get_coefs(),
a.pol().get_degree(),
s+1,
0,
0);
mpfi_swap(n,(mpfi_ptr)a.mpfi());
// It is not possible to clear the original mpfi because it
// could be allocated inside RS memory.
// TODO: clear the mpfi except the first time it is refined
}else{
rs3_refine_u_root((mpfi_ptr)a.mpfi(),
a.pol().get_coefs(),
a.pol().get_degree(),
s+1,
0,
0);
}
CGAL_postcondition(a.inf()<=a.sup());
}
} // namespace RS3
#endif // CGAL_RS_REFINE_1_RS_H

175
Algebraic_kernel_d/include/CGAL/RS/rs_calls_1.h
View File

@ -1,175 +0,0 @@
// Copyright (c) 2006-2009 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_RS_CALLS_1_H
#define CGAL_RS_RS_CALLS_1_H
#include <CGAL/RS/basic.h>
#include <gmp.h>
#include <mpfr.h>
#include <mpfi.h>
#include <rs_exports.h>
namespace CGAL{
#define CGALRS_CSTR(S) ((char*)(S))
#ifdef CGAL_RS_OLD_INCLUDES
#define CGALRS_PTR(a) long int a
#else
#define CGALRS_PTR(a) void *a
#endif
// initialize RS solver
inline void init_solver(){
static bool first=true;
if(first){
first=false;
//rs_version(stderr);
rs_init_rs();
}else{
rs_reset_all();
}
}
// reset RS memory
inline void reset_solver(){
rs_reset_all();
}
inline int affiche_sols_eqs(mpfi_ptr *x){
CGALRS_PTR(ident_sols_eqs);
CGALRS_PTR(ident_node);
CGALRS_PTR(ident_vect);
CGALRS_PTR(ident_elt);
int nb_elts;
ident_sols_eqs=rs_get_default_sols_eqs();
nb_elts=rs_export_list_vect_ibfr_nb(ident_sols_eqs);
ident_node=rs_export_list_vect_ibfr_firstnode(ident_sols_eqs);
//x=(mpfi_ptr*)malloc(nb_elts*sizeof(mpfi_ptr));
for(int i=0;i<nb_elts;++i){
ident_vect=rs_export_list_vect_ibfr_monnode(ident_node);
CGAL_assertion_msg(rs_export_dim_vect_ibfr(ident_vect)==1,
"the dimension of vector must be 1");
ident_elt=rs_export_elt_vect_ibfr(ident_vect,0);
x[i]=(mpfi_ptr)rs_export_ibfr_mpfi(ident_elt);
ident_node=rs_export_list_vect_ibfr_nextnode(ident_node);
}
return nb_elts;
}
inline void affiche_sols_constr(int nr,mpfi_ptr p){
CGALRS_PTR(ident_sols_eqs);
CGALRS_PTR(ident_node);
CGALRS_PTR(ident_vect);
CGALRS_PTR(ident_elt);
int nb_elts;
ident_sols_eqs=rs_get_default_sols_ineqs();
nb_elts=rs_export_list_vect_ibfr_nb(ident_sols_eqs);
ident_node=rs_export_list_vect_ibfr_firstnode(ident_sols_eqs);
for(int i=0;i<nb_elts;++i){
ident_vect=rs_export_list_vect_ibfr_monnode(ident_node);
CGAL_assertion_code(int nb=)
rs_export_dim_vect_ibfr(ident_vect);
CGAL_assertion_msg((nb==1),
"the vector must contain one element");
ident_elt=rs_export_elt_vect_ibfr(ident_vect,0);
if(i==nr){
mpfi_set(p,(mpfi_ptr)rs_export_ibfr_mpfi(ident_elt));
//break;
}
ident_node=rs_export_list_vect_ibfr_nextnode(ident_node);
}
}
// nroot is the number of root of the algebraic number
inline Sign affiche_signs_constr(int nroot){
CGALRS_PTR(ident_sols_eqs);
CGALRS_PTR(ident_node);
CGALRS_PTR(ident_vect);
CGALRS_PTR(ident_elt);
int nb_elts;
mpfi_t tmp;
mpfi_init(tmp);
ident_sols_eqs=rs_get_default_sols_ineqs();
nb_elts=rs_export_list_vect_ibfr_nb(ident_sols_eqs);
ident_node=rs_export_list_vect_ibfr_firstnode(ident_sols_eqs);
for(int i=1;i<nb_elts+1;++i){
ident_vect=rs_export_list_vect_ibfr_monnode(ident_node);
CGAL_assertion_code(int nb=)
rs_export_dim_vect_ibfr(ident_vect);
CGAL_assertion_msg((nb==1),
"the vector must contain one element");
ident_elt=rs_export_elt_vect_ibfr(ident_vect,0);
if(i==nroot+1){
mpfi_set(tmp,(mpfi_ptr)rs_export_ibfr_mpfi(ident_elt));
break;
}
ident_node=rs_export_list_vect_ibfr_nextnode(ident_node);
}
//std::cerr << "\nreturned value: ";
//mpfi_out_str(stderr,10,0,tmp);
//std::cerr << std::endl;
// mpfi_is_zero(tmp) doesn't work. The reason is that MPFR_SIGN in
// the mpfi code returns 1 when applied to the left and right zeros.
// This is not surprising because zero is signed in IEEE 754, and MPFR
// adopts it. Nevertheless, mpfr_sgn returns 0, but mpfi doesn't use
// it to implement mpfi_is_zero.
// Here is the difference (from MPFR source code):
// define mpfr_sgn(_x) (mpfr_zero_p(_x) ? 0 : MPFR_SIGN(_x))
//
if(mpfr_zero_p(&(tmp->right))&&mpfr_zero_p(&(tmp->left)))
return ZERO;
// the same holds for mpfi_is_pos and mpfi_is_neg
if((mpfr_sgn(&(tmp->left))>=0)&&(mpfr_sgn(&(tmp->right)))>0)
return POSITIVE;
if((mpfr_sgn(&(tmp->left))<0)&&(mpfr_sgn(&(tmp->right))<=0))
return NEGATIVE;
// if we arrive here, it is because the signs of the endpoints are -
// and +, and (I think) RS guarantees that this never happens
CGAL_assertion_msg(false,"error in sign calculation");
return ZERO;
}
inline void create_rs_upoly(mpz_t *poly,const int deg,CGALRS_PTR(ident_pol)){
CGALRS_PTR(ident_mon);
CGALRS_PTR(ident_coeff);
rs_import_uppring(CGALRS_CSTR("T"));
for(int i=0;i<=deg;++i)
if(mpz_sgn(poly[i])){ // don't add if == 0
ident_mon=rs_export_new_mon_upp_bz();
ident_coeff=rs_export_new_gmp();
rs_import_bz_gmp
(ident_coeff,TO_RSPTR_IN(&(poly[i])));
rs_dset_mon_upp_bz(ident_mon,ident_coeff,i);
rs_dappend_list_mon_upp_bz(ident_pol,ident_mon);
}
}
inline void create_rs_uconstr(mpz_t **list_constr,
const int *list_degs,
CGALRS_PTR(ident_list)){
CGALRS_PTR(ident_poly);
ident_poly=rs_export_new_list_mon_upp_bz();
create_rs_upoly(*list_constr,*list_degs,ident_poly);
rs_dappend_list_sup_bz(ident_list,ident_poly);
}
} // namespace CGAL
#endif // CGAL_RS_RS_CALLS_1_H

89
Algebraic_kernel_d/include/CGAL/RS/sign_1.h
View File

@ -1,89 +0,0 @@
// Copyright (c) 2006-2010 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_SIGN_1_H
#define CGAL_RS_SIGN_1_H
#include <gmp.h>
#include <mpfr.h>
#include <mpfi.h>
#include <CGAL/RS/basic.h>
#include <CGAL/RS/dyadic.h>
#include <CGAL/RS/polynomial_1.h>
#include <CGAL/RS/algebraic_1.h>
#include <CGAL/RS/solve_1.h>
#include <CGAL/assertions.h>
namespace CGAL{
struct RSSign{
// This function calculates the sign of the evaluation of a polynomial
// at a given algebraic number. If it is impossible to know the sign
// evaluating the interval, it calls sign_1_rs, which uses RS to do it.
static CGAL::Sign sign_1(const RS_polynomial_1 &p,const Algebraic_1 &x){
RS::rs_sign s=p.sign_mpfi(x.mpfi());
if(s!=RS::RS_UNKNOWN)
return RS::convert_rs_sign(s);
return sign_1_rs(p,x);
}
// compute the sign of the polynomial at a given dyadic
static inline CGAL::Sign
exactsignatdyadic(const RS_polynomial_1 &p,CGALRS_dyadic_srcptr d){
return p.sign_dyadic(d);
}
// compute the sign of the polynomial at a given mpfr
static inline CGAL::Sign
exactsignat(const RS_polynomial_1 &p,mpfr_srcptr m){
return p.sign_mpfr(m);
}
// This function uses interval arithmetic to calculate the sign of the
// evaluation. First, it converts the given mpfr number to an mpfi
// interval with the given precision. Later, it evaluates the
// polynomial in the interval an looks for the sign. If it is not able
// to determine it, it calls the exact signat function.
static CGAL::Sign quicksignat
(const RS_polynomial_1 &p,mpfr_srcptr xcoord,mp_prec_t prec){
mpfi_t x_approx;
RS::rs_sign s;
mpfi_init2(x_approx,prec);
mpfi_set_fr(x_approx,xcoord);
s=p.sign_mpfi(x_approx);
mpfi_clear(x_approx);
if(s!=RS::RS_UNKNOWN)
return RS::convert_rs_sign(s);
return p.sign_mpfr(xcoord);
}
// This function is supposed to return the signat calculated in
// the best way.
static CGAL::Sign signat(const RS_polynomial_1 &p,mpfr_srcptr xcoord){
mp_prec_t xprec=mpfr_get_prec(xcoord);
if(xprec>>(mp_prec_t)12) // switch to exact if xprec>=8192
return p.sign_mpfr(xcoord);
return quicksignat(p,xcoord,xprec);
}
}; // struct RSSign
} // namespace CGAL
#endif // CGAL_RS_SIGN_1_H

88
Algebraic_kernel_d/include/CGAL/RS/sign_1_no_rs.h
View File

@ -1,88 +0,0 @@
// Copyright (c) 2006-2010 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_SIGN_1_NO_RS_H
#define CGAL_RS_SIGN_1_NO_RS_H
#include <mpfi.h>
#include <CGAL/RS/basic.h>
#include <CGAL/RS/sign_1_no_rs.h>
#include <CGAL/RS/polynomial_1.h>
#include <CGAL/RS/polynomial_1_utils.h>
#include <CGAL/RS/algebraic_1.h>
#include <CGAL/RS/sign_1.h>
#include <CGAL/RS/refine_1_rs.h>
namespace CGAL{
// compute the sign of the polynomial at a given algebraic number
template <class GcdPolicy>
CGAL::Sign sign_1_no_rs(const RS_polynomial_1 &p,const Algebraic_1 &x){
typedef GcdPolicy Gcd;
//unsigned bisect_steps=4;
unsigned bisect_steps=1000;
rs_sign s;
CGAL::Sign sleft,sright;
RS_polynomial_1 *gcd,*deriv;
// we check first by evaluating intervals
s=p.sign_mpfi(x.mpfi());
if(s!=RS::RS_UNKNOWN)
return convert_rs_sign(s);
// we are not sure, we calculate the gcd in order to know if both
// polynomials have common roots
gcd=&(Gcd()(sfpart_1<Gcd>()(p),sfpart_1<Gcd>()(x.pol())));
if(!gcd->get_degree_static())
goto refineandreturn;
// at this point, gcd is not 1; we proceed as follows:
// -we refine x until having p monotonic in x's interval (to be
// sure that p has at most one root on that interval)
// -if the gcd has a root on this interval, both roots are equal
// (we return 0), otherwise, we refine until having a result
// how to assure that p is monotonic: when its derivative is never zero
deriv=&(sfpart_1<Gcd>()(p).derive());
while(deriv->sign_mpfi(x.mpfi())==RS::RS_UNKNOWN)
RS3::refine_1(x,(bisect_steps*=2));
//bisect_n<Gcd>(x,(bisect_steps*=2));
// how to know that the gcd has a root: just evaluating endpoints
if((sleft=RSSign::exactsignat(*gcd,x.left()))==CGAL::ZERO||
(sright=RSSign::exactsignat(*gcd,x.right()))==CGAL::ZERO||
(sleft!=sright)){
return CGAL::ZERO;
}
refineandreturn:
// the sign is not zero, we refine until having a result
for(;;){
RS3::refine_1(x,bisect_steps);
s=p.sign_mpfi(x.mpfi());
if(s==RS::RS_POSITIVE)
return CGAL::POSITIVE;
if(s==RS::RS_NEGATIVE)
return CGAL::NEGATIVE;
bisect_steps*=2;
}
}
} // namespace CGAL
#endif // CGAL_RS_SIGN_1_NO_RS_H

59
Algebraic_kernel_d/include/CGAL/RS/sign_1_rs.h
View File

@ -1,59 +0,0 @@
// Copyright (c) 2009-2010 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_SIGN_1_RS_H
#define CGAL_RS_SIGN_1_RS_H
#include <gmp.h>
#include <mpfr.h>
#include <mpfi.h>
#include <CGAL/RS/algebraic_1.h>
#include <CGAL/assertions.h>
namespace RS3{
inline CGAL::Sign sign_1(const CGAL::RS_polynomial_1 &p,
const CGAL::Algebraic_1 &a){
mpz_t* list_constr[]={p.get_coefs()};
int list_degs[]={p.get_degree()};
mpfi_t* sols_u=(mpfi_t*)malloc(sizeof(mpfi_t));
*(sols_u[0])=*(a.mpfi());
mpfi_t **vals_constr=(mpfi_t**)malloc(sizeof(mpfi_t*));
vals_constr[0]=(mpfi_t*)malloc(sizeof(mpfi_t));
mpfi_init(vals_constr[0][0]);
rs3_refine_eval_u(a.pol().get_coefs(),a.pol().get_degree(),NULL,0,
(const mpz_t **)list_constr,list_degs,
1,sols_u,1,
vals_constr,NULL,MPFR_PREC_MIN,1,1,1);
if(mpfr_zero_p(&vals_constr[0][0]->left)!=0 &&
mpfr_zero_p(&vals_constr[0][0]->right)!=0){
return CGAL::ZERO;
}
CGAL_assertion(mpfi_has_zero(vals_constr[0][0])<=0);
if(mpfi_is_pos(vals_constr[0][0])){
return CGAL::POSITIVE;
}
return CGAL::NEGATIVE;
}
} // namespace RS3
#endif // CGAL_RS_SIGN_1_RS_H

73
Algebraic_kernel_d/include/CGAL/RS/solve_1.h
View File

@ -1,73 +0,0 @@
// Copyright (c) 2006-2009 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_SOLVE_1_H
#define CGAL_RS_SOLVE_1_H
#include <CGAL/basic.h>
#include <CGAL/RS/basic.h>
#include <CGAL/RS/dyadic.h>
#include <CGAL/RS/polynomial_1.h>
#include <CGAL/RS/algebraic_1.h>
#include <CGAL/RS/rs_calls_1.h>
#include <CGAL/Gmpfi.h>
#include <vector>
namespace CGAL{
class RS_polynomial_1;
class Algebraic_1;
// solve given the precision, returns the number of roots
inline int solve_1(mpfi_ptr *x,
const RS_polynomial_1 &p1,
unsigned int prec=CGAL_RS_DEF_PREC){
rs_reset_all();
create_rs_upoly(p1.get_coefs(),p1.get_degree(),rs_get_default_up());
set_rs_precisol(prec);
set_rs_verbose(CGAL_RS_VERB);
rs_run_algo(CGALRS_CSTR("UISOLE"));
return affiche_sols_eqs(x);
}
// calculate the sign of a polynomial evaluated at the root of another
inline Sign sign_1_rs(const RS_polynomial_1 &p1,
const Algebraic_1 &a,
unsigned int prec=CGAL_RS_MIN_PREC){
mpz_t **constr;
int *degs;
CGAL_assertion(a.is_consistent());
rs_reset_all ();
// tell RS to find the roots of this polynomial
create_rs_upoly (a.pol().get_coefs (), a.pol().get_degree (),
rs_get_default_up ());
// the constraint vector will have just one element
constr = (mpz_t**)malloc(sizeof(mpz_t*));
*constr = p1.get_coefs ();
degs = (int*)malloc(sizeof(int));
*degs = p1.get_degree ();
create_rs_uconstr (constr, degs, rs_get_default_ineqs_u ());
set_rs_precisol (prec);
set_rs_verbose (CGAL_RS_VERB);
rs_run_algo(CGALRS_CSTR("UISOLES"));
return affiche_signs_constr(a.nr());
}
} // namespace CGAL
#endif // CGAL_RS_SOLVE_1_H

114
Algebraic_kernel_d/include/CGAL/RS/ugcd/crt.h
View File

@ -1,114 +0,0 @@
// Copyright (c) 2007 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__CRT_H
#define CGAL_RS__CRT_H
#include "pp.h"
#include <gmp.h>
#include <vector>
#include <boost/multi_array.hpp>
namespace CGAL{
namespace RS_MGCD{
class Crt:public Prime_polynomial{
protected:
// chinese remainder torture (GCL, page 180)
static
void cra(mpz_ptr r,CGALRS_PN * m,CGALRS_PN *u,int n){
int k,i;
CGALRS_PN product,temp;
std::vector<CGALRS_PN> v(n);
v[0]=u[0];
for(k=1;k<n;++k){
p_set_prime(m[k]);
// step 1, gamma_k is p_inv(product)
product=p_convert(m[0]);
for(i=1;i<k;++i)
//product=p_mul(product,p_convert(m[i]));
product=p_mulc(product,m[i]);
// step 2
temp=p_convert(v[k-1]);
for(i=k-2;i>=0;--i)
//temp=p_add(p_convert(v[i]),p_mul(temp,p_convert(m[i])));
temp=p_mulcaddc(temp,m[i],v[i]);
//v[k]=p_mul(p_sub(p_convert(u[k]),temp),p_inv(product));
v[k]=p_convsubdiv(u[k],temp,product);
}
// step 3: operations are done in Zm, not in Z
CGALRS_mpz_set_spn(r,p_pntospn(v[n-1]));
for(k=n-2;k>=0;--k){
CGALRS_mpz_mul_pn(r,r,m[k]);
CGALRS_mpz_add_pn(r,r,v[k]);
}
return;
};
// polynomial chinese remainder algorithm, it is the same:
// m are the modules, and m has size size_y, p is the residue
// vector and also has size size_y, but every one of its elements
// has size size_x (which will be de degree of the output
// polynomial);
// size_y is what is called n in the book
static
void pcra(mpz_t *r,
CGALRS_PN *m,
std::vector<CGALRS_PN* > p,
int size_x,
int size_y){
typedef boost::multi_array<CGALRS_PN,2> pn_matrix;
typedef pn_matrix::index pn_matrix_index;
pn_matrix v(boost::extents[size_x+1][size_y]);
pn_matrix_index i,j,k;
CGALRS_PN product,temp;
for(j=0;j<=size_x;++j){
v[j][0]=p[0][j];
for(k=1;k<size_y;++k){
p_set_prime(m[k]);
// step 1, gamma_k is p_inv(product)
product=p_convert(m[0]);
for(i=1;i<k;++i)
product=p_mulc(product,m[i]);
// step 2
temp=p_convert(v[j][k-1]);
for(i=k-2;i>=0;--i)
temp=p_mulcaddc(temp,m[i],v[j][i]);
v[j][k]=p_convsubdiv(p[k][j],temp,product);
}
}
// step 3
// be careful: operations are done in Zm, not in Z
for(j=0;j<=size_x;++j){
CGALRS_mpz_set_spn(r[j],p_pntospn(v[j][size_y-1]));
for(k=size_y-2;k>=0;--k){
CGALRS_mpz_mul_pn(r[j],r[j],m[k]);
CGALRS_mpz_add_pn(r[j],r[j],v[j][k]);
}
}
};
}; // class Crt
} // namespace RS_MGCD
} // namespace CGAL
#endif // CGAL_RS__CRT_H

75
Algebraic_kernel_d/include/CGAL/RS/ugcd/inverse.h
View File

@ -1,75 +0,0 @@
// Copyright (c) 2007-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__INVERSE_H
#define CGAL_RS__INVERSE_H
#ifdef _MSC_VER
# define CGALRS_S64 __int64
# define CGALRS_U64 unsigned __int64
# define CGALRS_U32 unsigned __int32
#else
# include <stdint.h>
# define CGALRS_S64 int64_t
# define CGALRS_U64 uint64_t
# define CGALRS_U32 uint32_t
#endif
namespace CGAL{
namespace RS_MGCD{
#define CGALRS_N(A) (A<0?-A:A)
//#define CGALRS_U(A) (A<0?-1:1)
#define CGALRS_U(A,C) (A<0?(C<0?1:-1):(C<0?-1:1))
class Inverse{
protected:
// given a and b, returns s such that gcd(a,b)=s*a+t*b (GCL, page 36)
// s*a+t*q=1 => s is the inverse of a, mod q (pafe 173)
static CGALRS_S64 eea_s(CGALRS_U32 a,CGALRS_U32 b){
CGALRS_S64 c1,d1,r1;//,c2,d2,r2,t,s;
CGALRS_U32 r,c,d;//,q;
// a and b are positive, so a=CGALRS_N(a) and b=CGALRS_N(b)
//c=CGALRS_N(a); d=CGALRS_N(b);
c=a; d=b;
c1=1; d1=0;
//c2=0; d2=1;
while(d){
//q=c/d;
r=c%d;
r1=c1-d1*(c/d); //r2=c2-d2*q;
c=d; c1=d1; //c2=d2;
d=r; d1=r1; //d2=r2;
}
// gcd(a,b) is CGALRS_N(c)
//t=c2/(CGALRS_U(b)*CGALRS_U(c));
// a and c are always positive, so s=c1/CGALRS_U(a,c) equals c1
//s=c1/CGALRS_U(a,c);
//return s;
return c1;
};
}; // class Inverse
#undef CGALRS_N
#undef CGALRS_U
} // namespace RS_MGCD
} // namespace CGAL
#endif // CGAL_RS__INVERSE_H

148
Algebraic_kernel_d/include/CGAL/RS/ugcd/p.h
View File

@ -1,148 +0,0 @@
// Copyright (c) 2007-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__P_H
#define CGAL_RS__P_H
#include <CGAL/RS/basic.h>
#include "inverse.h"
namespace CGAL{
namespace RS_MGCD{
// CGALRS_PN size is 32 bits, the sizes of CGALRS_LPN and CGALRS_SPN must be,
// at least, as twice as the size of CGALRS_PN
#define CGALRS_PN_BITS 32
#define CGALRS_PN CGALRS_U32 // unsigned
#define CGALRS_LPN CGALRS_U64 // unsigned long long
#define CGALRS_SPN CGALRS_S64 // long long
#define CGALRS_mpz_set_pn(A,PN) mpz_set_ui(A,(unsigned long)(PN))
#define CGALRS_mpz_mul_pn(A,B,PN) mpz_mul_ui(A,B,(unsigned long)(PN))
#define CGALRS_mpz_add_pn(A,B,PN) mpz_add_ui(A,B,(unsigned long)(PN))
#define CGALRS_mpz_sub_pn(A,B,PN) mpz_sub_ui(A,B,(unsigned long)(PN))
#define CGALRS_mpz_set_spn(A,SPN) mpz_set_si(A,(long)(SPN))
#define CGALRS_mpz_add_spn(A,B,SPN) \
(SPN<0?CGALRS_mpz_sub_pn(A,B,-(SPN)):CGALRS_mpz_add_pn(A,B,(SPN)))
#define CGALRS_mpz_mul_spn(A,B,SPN) mpz_mul_si(A,B,SPN)
CGALRS_THREAD_ATTR CGALRS_PN prime; class Prime:public Inverse{
protected:
static CGALRS_SPN p_pntospn(CGALRS_PN p){
if(p>(prime-1)/2)
return (CGALRS_SPN)p-prime;
return (CGALRS_SPN)p;
};
static void p_set_prime(CGALRS_PN p){prime=p;};
static CGALRS_PN p_prime(){return prime;};
static CGALRS_PN p_add(CGALRS_PN a,CGALRS_PN b){
CGALRS_LPN c=(CGALRS_LPN)a+b;
return (c<prime?(CGALRS_PN)c:(CGALRS_PN)(c-prime));
};
static CGALRS_PN p_sub(CGALRS_PN a,CGALRS_PN b){
return (a<b?prime-(b-a):a-b);
};
// returns a*b-c*d
static CGALRS_PN p_submuls(CGALRS_PN a,CGALRS_PN b,CGALRS_PN c,CGALRS_PN d){
CGALRS_LPN mul;
CGALRS_PN pnm1,pnm2;
mul=(CGALRS_LPN)a*b;
pnm1=(mul<prime?(CGALRS_PN)mul:(CGALRS_PN)(mul%prime));
mul=(CGALRS_LPN)c*d;
pnm2=(mul<prime?(CGALRS_PN)mul:(CGALRS_PN)(mul%prime));
return (pnm1<pnm2?prime-(pnm2-pnm1):pnm1-pnm2);
};
static CGALRS_PN p_inv(CGALRS_PN a){
CGALRS_SPN c=eea_s(a,prime);
return(c<0?(CGALRS_PN)(c+prime):(CGALRS_PN)c);
};
static CGALRS_PN p_convert(CGALRS_PN a){
return (a<prime?a:a%prime);
};
// define p_mul(A,B) ((CGALRS_PN)(((CGALRS_LPN)A*B)%p_prime()))
static CGALRS_PN p_mul(CGALRS_PN a,CGALRS_PN b){
CGALRS_LPN c=(CGALRS_LPN)a*b;
return (CGALRS_PN)(c<prime?c:c%prime);
};
static CGALRS_PN p_mul3(CGALRS_PN a,CGALRS_PN b,CGALRS_PN c){
CGALRS_LPN d;
if((d=(CGALRS_LPN)a*b)<prime)
d*=c;
else
d=(d%prime)*c;
return (CGALRS_PN)(d<prime?d:d%prime);
};
// returns a*conv(b)
static CGALRS_PN p_mulc(CGALRS_PN a,CGALRS_PN b){
CGALRS_LPN temp=(CGALRS_LPN)a*(b<prime?b:b%prime);
return (temp<prime?(CGALRS_PN)temp:(CGALRS_PN)(temp%prime));
};
// returns a*conv(b)+conv(c)
static CGALRS_PN p_mulcaddc(CGALRS_PN a,CGALRS_PN b,CGALRS_PN c){
CGALRS_LPN temp=(CGALRS_LPN)a*(b<prime?b:b%prime);
temp=(temp<prime?temp:temp%prime)+(c<prime?c:c%prime);
return (temp<prime?(CGALRS_PN)temp:(CGALRS_PN)(temp-prime));
};
// returns (conv(a)-b)*inv(c)
static CGALRS_PN p_convsubdiv(CGALRS_PN a,CGALRS_PN b,CGALRS_PN c){
CGALRS_SPN inv_c=eea_s(c,prime);
CGALRS_PN pninv_c=(inv_c<0?(CGALRS_PN)(inv_c+prime):(CGALRS_PN)inv_c);
CGALRS_PN conv_a=(a<prime?a:a%prime);
CGALRS_LPN mult=
(CGALRS_LPN)(conv_a<b?prime-(b-conv_a):conv_a-b)*pninv_c;
return (CGALRS_PN)(mult<prime?mult:mult%prime);
};
// vzGG, p. 73
static CGALRS_PN p_pow(CGALRS_PN a,CGALRS_PN n){
CGALRS_PN b,i;
i=1<<(CGALRS_PN_BITS-1);
while(!(i&n))
i>>=1;
b=a;
while(i>>=1)
b=(i&n?p_mul3(b,b,a):p_mul(b,b));
return b;
};
#define CGALRS_P_DIV(A,B) (p_mul(A,p_inv(B)))
static CGALRS_PN p_gcd(CGALRS_PN a,CGALRS_PN b){
if(!b)
return a;
return p_gcd(b,a%b);
};
}; // class Prime
} // namespace RS_MGCD
} // namespace CGAL
#endif // CGAL_RS__P_H

169
Algebraic_kernel_d/include/CGAL/RS/ugcd/pagealloc.h
View File

@ -1,169 +0,0 @@
// Copyright (c) 2007 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__PAGEALLOC_H
#define CGAL_RS__PAGEALLOC_H
#include <cstdlib>
#include <cstring>
namespace CGAL{
namespace RS_MGCD{
#define CGALRS_PAGESIZE 4194304
#define CGALRS_TABLESIZE 2048
#define CGALRS_PAGES 8
#define CGALRS_VOIDSCAST unsigned long
struct pinfo{
void *start;
size_t size;
};
CGALRS_THREAD_ATTR void** pages_startptr;
CGALRS_THREAD_ATTR size_t pages_max;//=CGALRS_PAGES;
CGALRS_THREAD_ATTR size_t pages_allocated;
CGALRS_THREAD_ATTR size_t pages_current;
CGALRS_THREAD_ATTR size_t page_remainingbytes;
CGALRS_THREAD_ATTR void* page_currentptr;
CGALRS_THREAD_ATTR struct pinfo *nodes_allocated;
CGALRS_THREAD_ATTR size_t nodes_total;
CGALRS_THREAD_ATTR size_t nodes_assigned;
class Page_alloc{
protected:
static
void* meminit(){
pages_startptr=(void**)malloc(CGALRS_PAGES*sizeof(void*));
pages_startptr[0]=malloc(CGALRS_PAGESIZE);
pages_allocated=1;
pages_current=0;
page_remainingbytes=CGALRS_PAGESIZE;
page_currentptr=pages_startptr[0];
nodes_total=CGALRS_TABLESIZE;
nodes_allocated=
(struct pinfo*)malloc(CGALRS_TABLESIZE*sizeof(struct pinfo));
nodes_assigned=0;
return page_currentptr;
};
static
void* newpage(){
void *r;
if(pages_allocated>pages_current+1){
++pages_current;
r=pages_startptr[pages_current];
page_currentptr=r;
page_remainingbytes=CGALRS_PAGESIZE;
return r;
}
// iso c++ forbids to initialize a static member (pages_max),
// so we have to start using pages_max when the amount of
// allocated pages reaches the value CGALRS_PAGES (this is not of
// course the cleanest way to do it)
if(pages_allocated==CGALRS_PAGES)
pages_max=2*CGALRS_PAGES;
else
pages_max=0;
if(pages_allocated==pages_max){
pages_max*=2;
pages_startptr=
(void**)realloc(pages_startptr,pages_max*sizeof(void*));
}
r=malloc(CGALRS_PAGESIZE);
pages_startptr[pages_allocated]=r;
page_currentptr=r;
++pages_allocated;
page_remainingbytes=CGALRS_PAGESIZE;
return r;
};
// if they ask us to allocate a memory size bigger than the page,
// we are lost (we could in that case make a bigger page)
static
void* palloc(size_t size){
void* r;
if(size>page_remainingbytes){
newpage();
return palloc(size);
}
if(nodes_assigned==nodes_total){
nodes_total*=2;
nodes_allocated=(struct pinfo*)realloc
(nodes_allocated,nodes_total*sizeof(struct pinfo));
}
page_remainingbytes-=size;
r=page_currentptr;
page_currentptr=(void*)((CGALRS_VOIDSCAST)page_currentptr+size);
// c++ does not support nodes_allocated[nodes_assigned]={r,s}
nodes_allocated[nodes_assigned].start=r;
nodes_allocated[nodes_assigned].size=size;
++nodes_assigned;
return r;
};
static
void* prealloc(void *ptr,size_t size){
void *dest;
size_t i=0;
while(nodes_allocated[i].start!=ptr)
++i;
if(nodes_allocated[i].size<size){
dest=palloc(size);
nodes_allocated[i].start=dest;
nodes_allocated[i].size=size;
return memcpy(dest,ptr,nodes_allocated[i].size);
}
return ptr;
};
#define CGALRS_PFREE(X) {}
//void CGALRS_PFREE(void* ptr){
// size_t i=0;
// while(nodes_allocated[i].start!=ptr)
// ++i;
// nodes_allocated[i].start=0;
// return;
//};
static
void* memclear(){
pages_current=0;
page_remainingbytes=CGALRS_PAGESIZE;
page_currentptr=pages_startptr[0];
nodes_assigned=0;
return page_currentptr;
}
static
void memrelease(){
size_t i;
for(i=0;i<pages_allocated;++i)
free(pages_startptr[i]);
return;
};
}; // class Page_alloc
} // namespace RS_MGCD
} // namespace CGAL
#endif // CGAL_RS__PAGEALLOC_H

135
Algebraic_kernel_d/include/CGAL/RS/ugcd/pp.h
View File

@ -1,135 +0,0 @@
// Copyright (c) 2007 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__PP_H
#define CGAL_RS__PP_H
#include <gmp.h>
#include "p.h"
#include "pagealloc.h"
#include <cstdio>
namespace CGAL{
namespace RS_MGCD{
class Prime_polynomial:public Prime,public Page_alloc{
protected:
static
int pp_from_poly(CGALRS_PN *pp,mpz_t *poly,int n){
int i;
if(!(pp[n]=mpz_fdiv_ui(poly[n],p_prime())))
return -1;
for(i=0;i<n;++i)
pp[i]=mpz_fdiv_ui(poly[i],p_prime());
return n;
};
static
void pp_out_str(FILE *stream,CGALRS_PN* pp,int n){
int i;
for(i=n;i;--i)
fprintf(stream,"%d*x^%d+",pp[i],i);
fprintf(stream,"%d",pp[0]);
return;
};
// Knuth 2; m>=n
static
int pp_pdivrem(CGALRS_PN *r,CGALRS_PN *u,int m,CGALRS_PN *v,int n){
int k,j;
for(k=0;k<=m;++k)
r[k]=u[k];
// division p. 402
//for(k=m-n;k>=0;--k)
// for(j=n+k-1;j>=k;--j)
// r[j]=p_sub(r[j],p_mul(qk,v[j-k]));
// pseudo-division, p. 407
for(k=m-n;k>=0;--k)
for(j=n+k-1;j>=0;--j)
//r[j]=p_sub(
// p_mul(v[n],r[j]),
// p_mul(r[n+k],(j<k?0:v[j-k])));
r[j]=p_submuls(v[n],r[j],r[n+k],(j<k?0:v[j-k]));
--n;
while(!r[n]&&n)
--n;
return n;
};
static
CGALRS_PN pp_pp(CGALRS_PN *pp,CGALRS_PN *p,int dp){
int i;
CGALRS_PN inv,cont=p[dp];
for(i=0;i<dp;++i)
cont=p_gcd(cont,p[i]);
inv=p_inv(cont);
for(i=0;i<=dp;++i)
pp[i]=p_mul(p[i],inv);
return cont;
};
// GCL, page 280; da>=db
static
int pp_gcd(CGALRS_PN *g,CGALRS_PN *a,int da,CGALRS_PN *b,int db){
CGALRS_PN *r0,*r1,*r2;
int i,d0,d1,d2;
d0=da;
r0=(CGALRS_PN*)palloc((1+da)*sizeof(CGALRS_PN));
//for(i=0;i<=da;++i)
// r0[i]=a[i];
pp_pp(r0,a,da);
d1=db;
r1=(CGALRS_PN*)palloc((1+da)*sizeof(CGALRS_PN));
//for(i=0;i<=db;++i)
// r1[i]=b[i];
pp_pp(r1,b,db);
r2=(CGALRS_PN*)palloc((1+da)*sizeof(CGALRS_PN));
d2=pp_pdivrem(r2,r0,d0,r1,d1);
while(d2){
for(i=0;i<=d1;++i)
r0[i]=r1[i];
d0=d1;
//for(i=0;i<=d2;++i)
// r1[i]=r2[i];
pp_pp(r1,r2,d2);
d1=d2;
d2=pp_pdivrem(r2,r0,d0,r1,d1);
}
if(!r2[0]){
//CGALRS_PN inv=p_inv(r1[d1]);
//g[d1]=1;
//for(i=0;i<d1;++i)
// g[i]=p_mul(inv,r1[i]);
for(i=0;i<=d1;++i)
g[i]=r1[i];
}else{
d1=0;
g[0]=1;
}
CGALRS_PFREE(r0);
CGALRS_PFREE(r1);
CGALRS_PFREE(r2);
return d1;
};
}; // class Prime_polynomial
} // namespace RS_MGCD
} // namespace CGAL
#endif // CGAL_RS__PP_H

142
Algebraic_kernel_d/include/CGAL/RS/ugcd/primes.h
View File

@ -1,142 +0,0 @@
// Copyright (c) 2007-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__PRIMES_H
#define CGAL_RS__PRIMES_H
#include "crt.h"
// I borrowed these numbers from Fabrice, this leaves us around 250000 primes
#define CGALRS_PR_MIN 2145338339
#define CGALRS_PR_MAX 2147483647
//#define CGALRS_PR_IS_PRIME(N) (pr_fermat(N)?pr_is_prime_bruteforce(N):0)
#define CGALRS_PR_IS_PRIME(N) pr_mrj(N)
#ifdef _MSC_VER
# define CGAL_RS_RANDOM rand
#else
# define CGAL_RS_RANDOM random
#endif
namespace CGAL{
namespace RS_MGCD{
CGALRS_THREAD_ATTR CGALRS_PN currentprime;
class Primes:public Crt{
private:
static int pr_is_prime_bruteforce(CGALRS_PN n){
int i,sqrtn;
sqrtn=(CGALRS_PN)(sqrt((double)n));
for(i=3;i<=sqrtn;++i)
if(!(n%i))
return 0;
return 1;
}
// vzGG, p. 507; returns 0 if n is composite
static int pr_fermat(CGALRS_PN n){
p_set_prime(n);
return(p_pow(2+((CGALRS_PN)CGAL_RS_RANDOM())%(n-4),n-1)==1);
}
// Solovay-Strassen
static int pr_ss(CGALRS_PN n){
CGALRS_PN a,x;
p_set_prime(n);
a=1+(CGALRS_PN)CGAL_RS_RANDOM()%(n-2);
x=CGALRS_P_DIV(a,n);
return(!x||p_pow(a,n>>1)!=x);
}
// Miller-Rabin
static int pr_mr(CGALRS_PN n){
CGALRS_PN s,d,a,x,r;
s=1;
d=(n-1)>>1;
while(!(d&1)){
++s;
d=d>>1;
}
p_set_prime(n);
a=2+(CGALRS_PN)CGAL_RS_RANDOM()%(n-4);
x=p_pow(a,d);
if(x==1||x==n-1)
return 1; // pobably prime
for(r=1;r<s;++r){
x=p_mul(x,x);
if(x==1)
return 0; // composite
if(x==n-1)
return 1; // probably
}
return 0; // composite
}
// Miller-Rabin-Jaeschke
static int pr_mrj(CGALRS_PN n){
CGALRS_PN s,d,a[3],r;//,x;
int index;
a[0]=2;
a[1]=7;
a[2]=61;
s=1;
d=(n-1)>>1;
while(!(d&1)){
++s;
d=d>>1;
}
p_set_prime(n);
index=-1;
start_test:
++index;
//if(index=3)
if(index==3)
return 1; // prime
if(p_pow(a[index],d)==1)
goto start_test;
for(r=0;r<s;++r)
if(p_pow(a[index],(d<<r))==n-1)
goto start_test;
return 0; // composite
}
static CGALRS_PN pr_actual(){
return currentprime;
}
protected:
static int pr_init(){
currentprime=CGALRS_PR_MIN;
return 0;
}
static CGALRS_PN pr_next(){
do{
currentprime+=2;
}while(!CGALRS_PR_IS_PRIME(currentprime));
return currentprime;
}
}; // class Primes
} // namespace RS_MGCD
} // namespace CGAL
#endif // CGAL_RS__PRIMES_H

197
Algebraic_kernel_d/include/CGAL/RS/ugcd/ugcd.h
View File

@ -1,197 +0,0 @@
// Copyright (c) 2007 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__UGCD_H
#define CGAL_RS__UGCD_H
#include <gmp.h>
#include "primes.h"
// let's assume that 300 is enough for degree 500 gcds
#define CGALRS_MOD_QTY 300
namespace CGAL{
namespace RS_MGCD{
class Ugcd:public Primes{
public:
static int ugcd (mpz_t *gcd,mpz_t *Anp,int degA,mpz_t *Bnp,int degB){
mpz_t *A,*B;
mpz_t lcgcd,cA,cB;
mpz_ptr m,bound;
// dG is initialized to zero only to avoid compiler complaints
int dA,dB,dG=0,maxd,i,maxA,maxB;
size_t modsize,modalloc;
std::vector<CGALRS_PN* > p;
CGALRS_PN *mA,*mB,*mG,*mod;
CGALRS_PN lc=0,scaleG;
if(degB>degA){
if(!degA){
mpz_set_ui(gcd[0],1);
return 0;
}else
return ugcd(gcd,Bnp,degB,Anp,degA);
}
if(!degB){
mpz_set_ui(gcd[0],1);
return 0;
}
// initialize the memory
meminit();
A=(mpz_t*)malloc((1+degA)*sizeof(mpz_t));
B=(mpz_t*)malloc((1+degB)*sizeof(mpz_t));
mpz_init_set(cA,Anp[degA]);
for(i=0;i<degA;++i)
mpz_gcd(cA,cA,Anp[i]);
mpz_init_set(cB,Bnp[degB]);
for(i=0;i<degB;++i)
mpz_gcd(cB,cB,Bnp[i]);
for(i=0;i<=degA;++i){
mpz_init(A[i]);
mpz_divexact(A[i],Anp[i],cA);
}
for(i=0;i<=degB;++i){
mpz_init(B[i]);
mpz_divexact(B[i],Bnp[i],cB);
}
// calculate the gcd of the principal coefficients
mpz_init(lcgcd);
mpz_gcd(lcgcd,A[degA],B[degB]);
// find the limit of modular image computation
maxA=degA;
for(i=0;i<degA;++i)
if(mpz_cmpabs(A[i],A[maxA])>0)
maxA=i;
maxB=degB;
for(i=0;i<degB;++i)
if(mpz_cmpabs(B[i],B[maxB])>0)
maxB=i;
mpz_pow_ui(cA,A[maxA],2);
mpz_mul_ui(cA,cA,degA+1);
mpz_mul_2exp(cA,cA,2*degA);
mpz_sqrt(cA,cA);
mpz_pow_ui(cB,B[maxB],2);
mpz_mul_ui(cB,cB,degB+1);
mpz_mul_2exp(cB,cB,2*degB);
mpz_sqrt(cB,cB);
if(mpz_cmp(cA,cB)<0){
bound=(mpz_ptr)cA;
m=(mpz_ptr)cB;
}else{
bound=(mpz_ptr)cB;
m=(mpz_ptr)cA;
}
mpz_mul(bound,bound,lcgcd);
mpz_mul_2exp(bound,bound,1);
mpz_setbit(bound,0);
mA=(CGALRS_PN*)palloc((1+degA)*sizeof(CGALRS_PN));
mB=(CGALRS_PN*)palloc((1+degB)*sizeof(CGALRS_PN));
maxd=degA; // we know that degA>=degB
mG=(CGALRS_PN*)palloc((1+maxd)*sizeof(CGALRS_PN));
pr_init();
mpz_set_ui(m,1);
mod=(CGALRS_PN*)palloc(CGALRS_MOD_QTY*sizeof(CGALRS_PN));
modalloc=CGALRS_MOD_QTY;
modsize=0;
while(mpz_cmp(m,bound)<=0){
do{
p_set_prime(pr_next());
dA=pp_from_poly(mA,A,degA);
if(dA!=-1){
dB=pp_from_poly(mB,B,degB);
if(dB!=-1)
lc=mpz_fdiv_ui(lcgcd,p_prime());
// lc is the image of the principal coefficient
}
}while(dA==-1||dB==-1||!lc
||mpz_divisible_ui_p(A[degA],p_prime())
||mpz_divisible_ui_p(B[degB],p_prime()));
// now we calculate the gcd mod p_prime
dG=pp_gcd(mG,mA,degA,mB,degB);
scaleG=CGALRS_P_DIV(lc,mG[dG]);
mG[dG]=lc;
for(i=0;i<dG;++i)
mG[i]=p_mul(mG[i],scaleG);
if(!dG){ // done, we know that the gcd is constant
mpz_set_ui(gcd[0],1);
dG=0;
goto cleanandexit;
}
if(dG<maxd){
CGALRS_mpz_set_pn(m,p_prime());
maxd=dG;
p.clear();
p.push_back(mG);
mod[0]=p_prime();
modsize=1;
mG=(CGALRS_PN*)palloc((1+maxd)*sizeof(CGALRS_PN));
// TODO: clean the CGALRS_PN* that are in p
}else{
if(dG==maxd){
CGALRS_mpz_mul_pn(m,m,p_prime());
if(modsize==modalloc){
modalloc*=2;
mod=(CGALRS_PN*)
prealloc(mod,modalloc*sizeof(CGALRS_PN));
}
mod[modsize]=p_prime();
++modsize;
p.push_back(mG);
mG=(CGALRS_PN*)palloc((1+maxd)*sizeof(CGALRS_PN));
}
}
}
pcra(gcd,mod,p,maxd,modsize);
cleanandexit:
CGALRS_PFREE(mA);
CGALRS_PFREE(mB);
CGALRS_PFREE(mG);
CGALRS_PFREE(mod);
// TODO: clean the CGALRS_PN* that are in p
for(i=0;i<=degA;++i)
mpz_clear(A[i]);
for(i=0;i<=degB;++i)
mpz_clear(B[i]);
mpz_clear(m);
mpz_clear(bound);
mpz_clear(lcgcd);
free(A);
free(B);
memrelease();
return dG;
};
}; // class Ugcd
} // namespace RS_MGCD
} // namespace CGAL
#endif // CGAL_RS__UGCD_H