diff --git a/Polynomial/include/CGAL/Polynomial_traits_d.h b/Polynomial/include/CGAL/Polynomial_traits_d.h index 02ee317289a..20f1045d8ea 100644 --- a/Polynomial/include/CGAL/Polynomial_traits_d.h +++ b/Polynomial/include/CGAL/Polynomial_traits_d.h @@ -16,11 +16,12 @@ #include - #include #include #include +#include + #define CGAL_POLYNOMIAL_TRAITS_D_BASE_TYPEDEFS \ typedef Polynomial_traits_d< Polynomial< Coefficient_ > > PT; \ typedef Polynomial_traits_d< Coefficient_ > PTC; \ @@ -77,6 +78,7 @@ template< class Coefficient_, class ICoeffAlgebraicCategory > class Polynomial_traits_d_base_icoeff_algebraic_category { public: typedef Null_functor Multivariate_content; + typedef Null_functor Interpolate; }; // Specializations @@ -117,6 +119,7 @@ class Polynomial_traits_d_base_icoeff_algebraic_category< return content; } }; + }; template< class Coefficient_ > @@ -146,6 +149,52 @@ class Polynomial_traits_d_base_icoeff_algebraic_category< return Innermost_coefficient(1); } }; + + struct Interpolate{ + typedef Polynomial Polynomial_1; + + void operator() ( + const Polynomial_1& m1, const Polynomial_d& u1, + const Polynomial_1& m2, const Polynomial_d& u2, + Polynomial_1& m, Polynomial_d& u) const { + Polynomial_1 s,t; + + CGAL::extended_euclidean_algorithm(m1,m2,s,t); + m = m1 * m2; + this->operator()(m1,m2,m,s,t,u1,u2,u); + } + + void operator() ( + const Polynomial_1& m1, const Polynomial_1& m2, const Polynomial_1& m, + const Polynomial_1& s, const Polynomial_1& t, + Polynomial_d u1, Polynomial_d u2, + Polynomial_d& u) const { +#ifndef NDEBUG + Polynomial_1 tmp,s_,t_; + tmp = CGAL::extended_euclidean_algorithm(m1,m2,s_,t_); + + CGAL_precondition(tmp == Polynomial_1(1)); + CGAL_precondition(s_ == s); + CGAL_precondition(t_ == t); +#endif + + typename CGAL::Coercion_traits::Cast cast; + typename Polynomial_traits_d::Canonicalize canonicalize; + typename Polynomial_traits_d::Pseudo_division_remainder + pseudo_remainder; + + CGAL_precondition(u1.degree() < m1.degree() || u1.is_zero()); + CGAL_precondition(u2.degree() < m2.degree() || u2.is_zero()); + if(m1.degree() < m2.degree()){ + Polynomial_d v = pseudo_remainder(cast(s)*(u2-u1),cast(canonicalize(m2))); + u = cast(m1)*v + u1; + } + else{ + Polynomial_d v = pseudo_remainder(cast(t)*(u1-u2),cast(canonicalize(m1))); + u = cast(m2)*v + u2; + } + } + }; }; template< class Coefficient_ > diff --git a/Polynomial/test/Polynomial/Polynomial_traits_d.cpp b/Polynomial/test/Polynomial/Polynomial_traits_d.cpp index d253f0d2e4e..b91931e3639 100644 --- a/Polynomial/test/Polynomial/Polynomial_traits_d.cpp +++ b/Polynomial/test/Polynomial/Polynomial_traits_d.cpp @@ -13,7 +13,6 @@ #include - #include static CGAL::Random my_rnd(346); // some seed @@ -84,40 +83,6 @@ generate_sparse_random_polynomial(int max_degree = 10){ return result; } -void print_monom(std::vector ev){ - if(ev.size() >= 1) - if (ev[0] != 0) - if (ev[0] == 1) - std::cout << "x"; - else - std::cout << "x^" << ev[0]; - - if(ev.size() >= 2) - if (ev[0] != 0) - if (ev[0] == 1) - std::cout << "y"; - else - std::cout << "y^" << ev[1]; - if(ev.size() >= 3) std::cout << "z^" << ev[2]; - if(ev.size() >= 4) std::cout << "x^" << ev[3]; - -} - -template -void print_monom_rep(Monom_rep monom_rep){ - std::sort(monom_rep.begin(), monom_rep.end()); - for(typename Monom_rep::iterator it = monom_rep.begin(); - it != monom_rep.end(); ++it){ - //if(it->second != 0){ - std::cout << it->second <<"*" ; - print_monom(it->first); - std::cout <<" + "; - //} - } - std::cout << std::endl; -} - - template void test_construct_polynomial(){ @@ -487,6 +452,51 @@ void test_multivariate_content(){ std::cerr << " ok "<< std::endl; } +// // Multivariate_content; +template +void test_interpolate(){ + std::cerr << "start test_interpolate "; std::cerr.flush(); + typedef Polynomial_traits_d PT_d; + typedef typename PT_d::Innermost_coefficient ICoeff; + typedef typename PT_d::Polynomial_d Polynomial_d; + typedef typename PT_d:: template Rebind::Other PT_1; + typedef typename PT_1::Polynomial_d Polynomial_1; + + typename PT_d::Interpolate interpolate; + typename PT_d::Evaluate eval; + + for(int i = 0; i < 5; i++){ + Polynomial_d p = generate_sparse_random_polynomial(i); + + Polynomial_1 m(1); + Polynomial_d u(0); + for (int j = 0; j <= i; j++){ + Polynomial_1 m1 = m; + Polynomial_d u1 = u; + Polynomial_1 m2 = Polynomial_1(ICoeff(-j),ICoeff(2)); + Polynomial_d u2 = eval(p,ICoeff(j)/ICoeff(2)); + interpolate(m1,u1,m2,u2,m,u); + } + CGAL_test_assert(u == p); + } + + for(int i = 0; i < 5; i++){ + Polynomial_d p = generate_sparse_random_polynomial(i); + + Polynomial_1 m(1); + Polynomial_d u(0); + for (int j = 0; j <= i; j++){ + Polynomial_1 m2 = m; + Polynomial_d u2 = u; + Polynomial_1 m1 = Polynomial_1(ICoeff(-j),ICoeff(2)); + Polynomial_d u1 = eval(p,ICoeff(j)/ICoeff(2)); + interpolate(m1,u1,m2,u2,m,u); + } + CGAL_test_assert(u == p); + } + std::cerr << " ok "<< std::endl; +} + // // Shift; template void test_shift(){ @@ -1361,6 +1371,8 @@ Polynomial_traits_d, CGAL::Integral_domain_tag > { test_univariate_content(); // Multivariate_content; test_multivariate_content(); + // Interpolate; + test_interpolate(); // Square_free_factorization; test_square_free_factorization(); } @@ -1551,9 +1563,10 @@ void test_AT(){ CGAL::set_pretty_mode(std::cout); CGAL::set_pretty_mode(std::cerr); + typedef typename AT::Integer Integer; - typedef typename AT::Rational Rational; - + typedef typename AT::Rational Rational; + std::cerr << std::endl; std::cerr << "Test for coefficient type Integer" << std::endl; std::cerr << "--------------------------------------" << std::endl; @@ -1596,6 +1609,9 @@ void test_AT(){ int main(){ +#if 1 + test_AT(); +#else #ifdef CGAL_USE_LEDA { typedef CGAL::LEDA_arithmetic_kernel AT; @@ -1607,6 +1623,7 @@ int main(){ typedef CGAL::CORE_arithmetic_kernel AT; test_AT(); } +#endif #endif return 0; }