diff --git a/Algebraic_kernel_d/include/CGAL/RS/Algebraic_kernel_rs_1.h b/Algebraic_kernel_d/include/CGAL/RS/Algebraic_kernel_rs_1.h deleted file mode 100644 index 15e54a38013..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/Algebraic_kernel_rs_1.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_ALGEBRAIC_KERNEL_RS_1 -#define CGAL_RS_ALGEBRAIC_KERNEL_RS_1 - -#include -#include -#include -#include - -template -struct Algebraic_kernel_rs_1{ - - typedef CoefficientType Coefficient; - typedef GcdFunction Gcd; - typedef CGAL::Polynomial 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 - Construct_algebraic_real_1; - typedef RSFunctors::Compute_polynomial_1 - Compute_polynomial_1; - typedef RSFunctors::Isolate_1 Isolate_1; - typedef RSFunctors::Is_square_free_1 - Is_square_free_1; - typedef RSFunctors::Make_square_free_1 - Make_square_free_1; - typedef RSFunctors::Square_free_factorize_1 - Square_free_factorize_1; - typedef RSFunctors::Is_coprime_1 - Is_coprime_1; - typedef RSFunctors::Make_coprime_1 - Make_coprime_1; - typedef RSFunctors::Solve_1 Solve_1; - typedef RSFunctors::Number_of_solutions_1 - Number_of_solutions_1; - typedef RSFunctors::Sign_at_1 Sign_at_1; - typedef RSFunctors::Is_zero_at_1 - Is_zero_at_1; - typedef RSFunctors::Compare_1 Compare_1; - typedef RSFunctors::Bound_between_1 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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_comparisons.h b/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_comparisons.h deleted file mode 100644 index de029659810..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_comparisons.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_ALGEBRAIC_1_COMPARISONS_H -#define CGAL_RS_ALGEBRAIC_1_COMPARISONS_H - -#include -#include -#include - -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(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(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(n1,n2)==EQUAL); -} - -inline -Algebraic_1 min BOOST_PREVENT_MACRO_SUBSTITUTION (const Algebraic_1 &a,const Algebraic_1 &b){ - return (ab?a:b); -} - -} // namespace CGAL - -#endif // CGAL_RS_ALGEBRAIC_1_COMPARISONS_H diff --git a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_constructors.h b/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_constructors.h deleted file mode 100644 index e67b5d364aa..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_constructors.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_ALGEBRAIC_1_CONSTRUCTORS_H -#define CGAL_RS_ALGEBRAIC_1_CONSTRUCTORS_H - -#include -#include -#include - -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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_member.h b/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_member.h deleted file mode 100644 index 1ca45dde492..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_member.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_ALGEBRAIC_1_MEMBER_H -#define CGAL_RS_ALGEBRAIC_1_MEMBER_H - -#include -#include - -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(&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 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 -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(*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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_operators.h b/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_operators.h deleted file mode 100644 index 43877d693dc..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_operators.h +++ /dev/null @@ -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 - -#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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_other.h b/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_other.h deleted file mode 100644 index ae81caf43ba..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_other.h +++ /dev/null @@ -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 - -#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 -inline -double to_double(Algebraic_1 &n){ - //typedef GcdPolicy Gcd; - //return n.to_double(); - return n.to_double(); -} - -inline -std::pairto_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<<'('<>(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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_real_embeddable.h b/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_real_embeddable.h deleted file mode 100644 index 2ee20e75c83..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/algebraic_1_real_embeddable.h +++ /dev/null @@ -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 -// -// ============================================================================ - - -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 > { - public: - std::pair 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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/basic.h b/Algebraic_kernel_d/include/CGAL/RS/basic.h deleted file mode 100644 index d93e990f09c..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/basic.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_BASIC_H -#define CGAL_RS_BASIC_H - -#include -#include - -#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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/compare_1.h b/Algebraic_kernel_d/include/CGAL/RS/compare_1.h deleted file mode 100644 index dd692fcb9ee..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/compare_1.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_COMPARE_1_H -#define CGAL_RS_COMPARE_1_H - -#include -#include -#include -#include -#include -#include - -// 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 -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() - -#ifndef CGAL_RS_DYADIC_H -#define CGAL_RS_DYADIC_H - -#include -#include -#include -#include -#include - -// for c++, compile with -lgmpxx -#ifdef __cplusplus -#include -#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<<"["< -// Michael Hemmer -// Luis Peñaranda -// -// ============================================================================ - -/// \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 -#include -#include -#include - -namespace CGAL { - -namespace internal { - -/// A model of concept RealRootIsolator. -/// -/// Polynomial_ must be Polynomial, and Bound_ must be Gmpfr. -template -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()(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 real_roots; - - //! restricted interval? - // TODO bool interval_given; - -}; - -} // namespace internal - -} //namespace CGAL - -#endif // CGAL_ALGEBRAIC_KERNEL_D_RS_ISOLATOR_H diff --git a/Algebraic_kernel_d/include/CGAL/RS/isole_1.h b/Algebraic_kernel_d/include/CGAL/RS/isole_1.h deleted file mode 100644 index 6681805d753..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/isole_1.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_ISOLE_1_H -#define CGAL_RS_ISOLE_1_H - -#include -#include -#include -#include - -namespace CGAL{ - -namespace RS{ - -// CGAL::RS::isolator

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. - -template -struct isolator{ - typedef PolynomialType PolynomialT; - inline std::vector - operator()(const PolynomialT&,unsigned int prec=CGAL_RS_DEF_PREC); -}; - -template -inline std::vector -isolator::operator()(const PolynomialType &p,unsigned int prec){ - CGAL_error_msg( - "isolator not implemented for this type of polynomials"); - return std::vector(); -} - -template <> -inline std::vector -isolator >::operator()(const Polynomial &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 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 - -#ifndef CGAL_RS_MEMORY_H -#define CGAL_RS_MEMORY_H - -#include -#include -#include -#include -#ifdef CGAL_USE_OLD_RS3 -extern "C"{ -#include -} -#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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1.h b/Algebraic_kernel_d/include/CGAL/RS/polynomial_1.h deleted file mode 100644 index 03c384c4912..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_POLYNOMIAL_1_H -#define CGAL_RS_POLYNOMIAL_1_H - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -namespace CGAL{ - -class Algebraic_1; -class RS_polynomial_1; -typedef boost::shared_ptr polyptr; -typedef std::pair polypow; -typedef std::vector sqfrvec; -typedef boost::shared_ptr sqfrptr; - -class RS_polynomial_1: - boost::addable1 > -{ - 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 - templateRS_polynomial_1(InputIt,InputIt); - templateT operator()(const T&)const; - templateSign sign_at(const T&)const; - templateRS_polynomial_1 operator*(const T&)const; - templateRS_polynomial_1& operator*=(const T&); - templateRS_polynomial_1& operator/=(const T&); - templateRS_polynomial_1 operator/(const T&)const; -}; - -} // namespace CGAL - -#include -#include -#include -#include -#include -#include - -#endif // CGAL_GRBS_POLYNOMIAL_1_H diff --git a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_constructors.h b/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_constructors.h deleted file mode 100644 index 353ac448987..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_constructors.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_POLYNOMIAL_1_CONSTRUCTORS_H -#define CGAL_RS_POLYNOMIAL_1_CONSTRUCTORS_H - -#include -#include - -namespace CGAL{ - -// auxiliar function to parse a c++ string and return a vector -// containing the coefficients as Gmpz's -inline -std::vector parse_string(std::string tstr){ - Polynomial_parser_1 parser; - parser.parse(tstr); - CGAL_assertion(parser.is_correct()); - std::vector 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 Coeff=parse_string(tstr); - std::vector::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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_eval.h b/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_eval.h deleted file mode 100644 index be38c255335..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_eval.h +++ /dev/null @@ -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 - -#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;ieval_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;i0?CGAL::POSITIVE:CGAL::NEGATIVE):CGAL::ZERO); - } - CGALRS_dyadic_init_set_z(h,leading_coefficient()); - for(int i=1;i0?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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_functors.h b/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_functors.h deleted file mode 100644 index 5f1937ed3a9..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_functors.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_POLYNOMIAL_FUNCTORS -#define CGAL_RS_POLYNOMIAL_FUNCTORS - -#include -#include -#include -#include -#include -#include - -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 - inline Polynomial_1 operator()( - InputIterator begin,InputIterator end)const{ - return Polynomial_1(begin,end); - }; - -//-------------------------------------------------- -// template -// 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 - 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{ - inline Coefficient operator() - (const Polynomial_1 &p,int i)const{ - return Coefficient(p.coef(i)); - }; -}; - -struct Get_innermost_coefficient: -public std::binary_function{ - 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{ - inline int operator()(const Polynomial_1 &p)const{ - return p.get_degree(); - }; -}; - -struct Degree_vector: -public std::unary_function{ - 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{ - inline Innermost_coefficient operator()(const Polynomial_1 &p)const{ - return Innermost_coefficient(p.leading_coefficient()); - }; -}; - -struct Canonicalize: -public std::unary_function{ - 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 - 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 - 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{ - 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()=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 -struct Is_square_free_1: - public std::unary_function{ - typedef GcdPolicy Gcd; - inline bool operator()(const Polynomial_1 &p){ - return(!(Gcd()(p,p.derive()).get_degree_static())); - }; - }; // Is_square_free_1 - -template -struct Make_square_free_1: - public std::unary_function{ - typedef GcdPolicy Gcd; - inline Polynomial_1 operator()(const Polynomial_1 &p){ - return sfpart_1(p); - }; - }; // Make_square_free_1 - -// this function is not well defined in algebraic kernel concepts; we -// don't care, we do it correctly -template -struct Square_free_factorize_1{ - template - int operator()(const Polynomial_1 &p,OutputIterator oi){ - typedef GcdPolicy Gcd; - sqfrvec factorization(sqfr_1(p)); - std::copy(factorization.begin(),factorization.end(),oi); - return factorization.size(); - }; -}; // Square_free_factorize_1 - -template -struct Is_coprime_1: - public std::binary_function{ - 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 -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{ - 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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_impl.h b/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_impl.h deleted file mode 100644 index f769f9f045d..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_impl.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_POLYNOMIAL_1_IMPL_H -#define CGAL_RS_POLYNOMIAL_1_IMPL_H - -namespace CGAL{ - -// constructor from a container of sorted coefficients -template -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 -RS_polynomial_1 RS_polynomial_1::operator*(const T &n)const{ - RS_polynomial_1 r(*this); - return (r*=n); -} - -template -RS_polynomial_1& RS_polynomial_1::operator*=(const T &s){ - Gmpz z(s); - return (*this*=z); -} - -template -RS_polynomial_1 RS_polynomial_1::operator/(const T &d)const{ - RS_polynomial_1 r(*this); - return (r/=d); -} - -template -RS_polynomial_1& RS_polynomial_1::operator/=(const T &d){ - Gmpz z(d); - return (*this/=z); -} - -template -Sign RS_polynomial_1::sign_at(const T &x)const{ - T y(leading_coefficient()); - int d=get_degree_static(); - for(int i=1;i0) - 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 -T RS_polynomial_1::operator()(const T &x)const{ - T y(leading_coefficient()); - int d=get_degree_static(); - for(int i=1;i - -#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<=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<>(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>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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_member.h b/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_member.h deleted file mode 100644 index e99b5a84b30..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_member.h +++ /dev/null @@ -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 - -#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;icapacity){ - 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;icoef(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;i0) - 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;i0) - 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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_operators.h b/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_operators.h deleted file mode 100644 index 80543f762dd..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_operators.h +++ /dev/null @@ -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 - -#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 -// -// (the first version of this file was written by Elias Tsigaridas) - -#ifndef CGAL_RS_PARSER_1_H -#define CGAL_RS_PARSER_1_H - -#include -#include -#include - -// if boost version is 1.38 or newer, we use the new version of spirit -#include -#if BOOST_VERSION >= 103800 -#include -#define CGAL_BOOST_SPIRIT boost::spirit::classic -#else -#include -#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 -{ - - 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 - poly, monomial, smonomial, X, unumber, number; - - CGAL_BOOST_SPIRIT::rule 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( 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( 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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_utils.h b/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_utils.h deleted file mode 100644 index a492a9fb867..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/polynomial_1_utils.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_POLYNOMIAL_1_UTILS_H -#define CGAL_RS_POLYNOMIAL_1_UTILS_H - -#include -#include -#include -#include -#include -#ifdef CGAL_RS_USE_UGCD -#include -#endif -#include -#ifdef CGAL_USE_RS3 -#include -#endif - -namespace CGAL{ - -#ifdef CGAL_USE_RS3 -// Fabrice's modular gcd algorithm -struct Rsgcd_1: -public std::binary_function{ - 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& operator()( - const RS_polynomial_1 &p1, - const RS_polynomial_1 &p2)const{ - typedef Polynomial P; - typedef from_rs_poly

convert; - typedef to_rs_poly

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& operator()( - const RS_polynomial_1 &u, - const RS_polynomial_1 &v)const{ - RS_polynomial_1 *result=new RS_polynomial_1( - u.get_degree()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{ - void operator()(mpz_ptr c,const RS_polynomial_1 &u){ - mpz_set(c,u.coef(0)); - for(int i=1;i{ - 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 , 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 >{ - std::pair - 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 - 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 - 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& 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* - 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="< -struct to_rs_poly: -public std::unary_function{ - 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 >: \ - public std::unary_function,RS_polynomial_1>{ \ - RS_polynomial_1& operator()(const Polynomial &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 >: \ - public std::unary_function,RS_polynomial_1>{ \ - RS_polynomial_1& operator()(const Polynomial &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 >: -public std::unary_function,RS_polynomial_1>{ - RS_polynomial_1& operator()(const Polynomial &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 -struct from_rs_poly: -public std::unary_function{ - typedef Polynomial_ Polynomial; - Polynomial operator()(const RS_polynomial_1 &p)const{ - std::cerr<<"can't convert to integer polynomial"< -struct from_rs_poly: -public std::unary_function{ - const RS_polynomial_1& operator()(const RS_polynomial_1 &p)const{ - return p; - } -}; - -template<> -struct from_rs_poly >: -public std::unary_function >{ - Polynomial operator()(const RS_polynomial_1 &p)const{ - typedef Polynomial_traits_d > PT; - mpz_t* pcoef=p.get_coefs(); - std::vector 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 >: -public std::unary_function >{ - Polynomial operator()(const RS_polynomial_1 &p)const{ - typedef Polynomial_traits_d > PT; - mpz_t* pcoef=p.get_coefs(); - std::vector 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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/refine_1.h b/Algebraic_kernel_d/include/CGAL/RS/refine_1.h deleted file mode 100644 index 87d11580f4e..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/refine_1.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_REFINE_1_H -#define CGAL_RS_REFINE_1_H - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -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 -Comparison_result bisect_at(const Algebraic_1 &a,mpfr_srcptr p){ - typedef GcdPolicy Gcd; - Sign sl,sp; - int round; - sp=RSSign::signat(sfpart_1()(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 ab -// precondition: a and b overlap -template -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(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(b,&(a.mpfi()->right))==LARGER) - return -1; - return 0; -} - -// bisect n times an interval; this function returns the number of -// refinements made -template -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;ileft), - &(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()(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 -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()(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 -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(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 -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()(a.pol()).eval_dyadic(f_lo,x_lo); - sfpart_1()(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()(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()(a.pol()),a.left())); - - return count; -} - -// applies qir until having an interval of size less than 2^(-t) -template -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() - (a.pol()).eval_dyadic(f_lo,(CGALRS_dyadic_ptr)(a.left())); - sfpart_1() - (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)()(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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/refine_1_rs.h b/Algebraic_kernel_d/include/CGAL/RS/refine_1_rs.h deleted file mode 100644 index 6c6192c173e..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/refine_1_rs.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_REFINE_1_RS_H -#define CGAL_RS_REFINE_1_RS_H - -#include -#include -#include -#include -#include - -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(sleft)-2) - s=mpfr_get_prec(&((mpfi_ptr)(a.mpfi()))->left)-2; - if(sright)-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)right) - -#ifndef CGAL_RS_RS_CALLS_1_H -#define CGAL_RS_RS_CALLS_1_H - -#include -#include -#include -#include -#include - -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;iright))&&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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/sign_1.h b/Algebraic_kernel_d/include/CGAL/RS/sign_1.h deleted file mode 100644 index 0da9eecf109..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/sign_1.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_SIGN_1_H -#define CGAL_RS_SIGN_1_H - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/sign_1_no_rs.h b/Algebraic_kernel_d/include/CGAL/RS/sign_1_no_rs.h deleted file mode 100644 index de493744a3a..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/sign_1_no_rs.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_SIGN_1_NO_RS_H -#define CGAL_RS_SIGN_1_NO_RS_H - -#include -#include -#include -#include -#include -#include -#include -#include - -namespace CGAL{ - -// compute the sign of the polynomial at a given algebraic number -template -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()(p),sfpart_1()(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()(p).derive()); - while(deriv->sign_mpfi(x.mpfi())==RS::RS_UNKNOWN) - RS3::refine_1(x,(bisect_steps*=2)); - //bisect_n(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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/sign_1_rs.h b/Algebraic_kernel_d/include/CGAL/RS/sign_1_rs.h deleted file mode 100644 index 2441e01e2b5..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/sign_1_rs.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_SIGN_1_RS_H -#define CGAL_RS_SIGN_1_RS_H - -#include -#include -#include -#include -#include - -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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/solve_1.h b/Algebraic_kernel_d/include/CGAL/RS/solve_1.h deleted file mode 100644 index c70dc5bdbad..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/solve_1.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS_SOLVE_1_H -#define CGAL_RS_SOLVE_1_H - -#include -#include -#include -#include -#include -#include -#include -#include - -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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/ugcd/crt.h b/Algebraic_kernel_d/include/CGAL/RS/ugcd/crt.h deleted file mode 100644 index 4b9514605ce..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/ugcd/crt.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS__CRT_H -#define CGAL_RS__CRT_H - -#include "pp.h" -#include -#include -#include - -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 v(n); - - v[0]=u[0]; - for(k=1;k=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 p, - int size_x, - int size_y){ - typedef boost::multi_array 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=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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/ugcd/inverse.h b/Algebraic_kernel_d/include/CGAL/RS/ugcd/inverse.h deleted file mode 100644 index 76dcb3eae68..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/ugcd/inverse.h +++ /dev/null @@ -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 - -#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 -# 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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/ugcd/p.h b/Algebraic_kernel_d/include/CGAL/RS/ugcd/p.h deleted file mode 100644 index 0f6e878d369..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/ugcd/p.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS__P_H -#define CGAL_RS__P_H - -#include -#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>=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 diff --git a/Algebraic_kernel_d/include/CGAL/RS/ugcd/pagealloc.h b/Algebraic_kernel_d/include/CGAL/RS/ugcd/pagealloc.h deleted file mode 100644 index 1e395f4b313..00000000000 --- a/Algebraic_kernel_d/include/CGAL/RS/ugcd/pagealloc.h +++ /dev/null @@ -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 - -#ifndef CGAL_RS__PAGEALLOC_H -#define CGAL_RS__PAGEALLOC_H - -#include -#include - -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 - -#ifndef CGAL_RS__PP_H -#define CGAL_RS__PP_H - -#include -#include "p.h" -#include "pagealloc.h" -#include - -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 - 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=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 - -#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>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 - -#ifndef CGAL_RS__UGCD_H -#define CGAL_RS__UGCD_H - -#include -#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 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;i0) - maxA=i; - maxB=degB; - for(i=0;i0) - 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