mirror of https://github.com/CGAL/cgal
new_resultant implementation using interpolation.
This commit is contained in:
Michael Hemmer
parent
b6d0b72e78
commit
acb79b6084
|
|
@ -117,6 +117,15 @@ template<typename OutputIterator, typename NT> OutputIterator
|
||||||
OutputIterator2 out_f,
|
OutputIterator2 out_f,
|
||||||
OutputIterator3 out_fx);
|
OutputIterator3 out_fx);
|
||||||
|
|
||||||
|
// eliminates outermost variable
|
||||||
|
template <class Coeff>
|
||||||
|
inline Coeff new_resultant(
|
||||||
|
const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&);
|
||||||
|
// eliminates innermost variable
|
||||||
|
template <class Coeff>
|
||||||
|
inline Coeff new_resultant_(
|
||||||
|
const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&);
|
||||||
|
|
||||||
|
|
||||||
} // namespace CGALi
|
} // namespace CGALi
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -1,109 +1,396 @@
|
||||||
// TODO: Add licence
|
// Copyright (c) 1997-2000 Max-Planck-Institute Saarbruecken (Germany).
|
||||||
|
// All rights reserved.
|
||||||
|
//
|
||||||
|
// This file is part of CGAL (www.cgal.org); you may redistribute it under
|
||||||
|
// the terms of the Q Public License version 1.0.
|
||||||
|
// See the file LICENSE.QPL distributed with CGAL.
|
||||||
|
//
|
||||||
|
// 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
|
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
||||||
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||||
//
|
//
|
||||||
// $URL:$
|
// $URL:$
|
||||||
// $Id: $
|
// $Id:$
|
||||||
//
|
//
|
||||||
//
|
//
|
||||||
// Author(s) :
|
// Author(s) : Michael Hemmer <hemmer@mpi-inf.mpg.de>
|
||||||
//
|
|
||||||
// ============================================================================
|
|
||||||
|
#ifndef CGAL_NEW_RESULTANT_H
|
||||||
|
#define CGAL_NEW_RESULTANT_H
|
||||||
|
|
||||||
|
// Modular arithmetic is slower, hence the default is 0
|
||||||
|
#ifndef CGAL_RESULTANT_USE_MODULAR_ARITHMETIC
|
||||||
|
#define CGAL_RESULTANT_USE_MODULAR_ARITHMETIC 0
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifndef CGAL_RESULTANT_USE_DECOMPOSE
|
||||||
|
#define CGAL_RESULTANT_USE_DECOMPOSE 1
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifndef CGAL_RESULTANT_USE_INTERPOLATION
|
||||||
|
#define CGAL_RESULTANT_USE_INTERPOLATION 1
|
||||||
|
#endif
|
||||||
|
|
||||||
// TODO: The comments are all original EXACUS comments and aren't adapted. So
|
|
||||||
// they may be wrong now.
|
|
||||||
|
|
||||||
#ifndef CGAL_POLYNOMIAL_RESULTANT_H
|
|
||||||
#define CGAL_POLYNOMIAL_RESULTANT_H
|
|
||||||
|
|
||||||
#include <CGAL/basic.h>
|
#include <CGAL/basic.h>
|
||||||
#include <CGAL/Polynomial.h>
|
#include <CGAL/Polynomial.h>
|
||||||
|
|
||||||
|
#include <CGAL/Polynomial_traits_d.h>
|
||||||
|
#include <CGAL/Polynomial/Interpolator.h>
|
||||||
#include <CGAL/Polynomial/prs_resultant.h>
|
#include <CGAL/Polynomial/prs_resultant.h>
|
||||||
#include <CGAL/Polynomial/bezout_matrix.h>
|
#include <CGAL/Polynomial/bezout_matrix.h>
|
||||||
|
#include <CGAL/Polynomial/subresultants.h>
|
||||||
|
|
||||||
#include <boost/type_traits.hpp>
|
#include <CGAL/Modular.h>
|
||||||
#include <boost/mpl/if.hpp>
|
#include <CGAL/Modular_traits.h>
|
||||||
|
#include <CGAL/Chinese_remainder_traits.h>
|
||||||
|
#include <CGAL/primes.h>
|
||||||
|
#include <CGAL/Polynomial/Cached_extended_euclidean_algorithm.h>
|
||||||
|
|
||||||
CGAL_BEGIN_NAMESPACE
|
CGAL_BEGIN_NAMESPACE
|
||||||
|
|
||||||
namespace CGALi {
|
|
||||||
|
|
||||||
template <class NT> inline
|
// The main function provided within this file is CGAL::CGALi::resultant(F,G),
|
||||||
NT resultant_(Polynomial<NT> A, Polynomial<NT> B) {
|
// all other functions are used for dispatching.
|
||||||
typedef typename Algebraic_structure_traits<NT>::Algebraic_category Algebraic_category;
|
// The implementation uses interpolatation for multivariate polynomials
|
||||||
return resultant_(A,B,Algebraic_category());
|
// Due to the recursive structuture of CGAL::Polynomial<Coeff> it is better
|
||||||
}
|
// to write the function such that the inner most variabel is eliminated.
|
||||||
|
// However, CGAL::CGALi::resultant(F,G) eliminates the outer most variabel.
|
||||||
|
// This is due to backward compatibility issues with code base on EXACUS.
|
||||||
|
// In turn CGAL::CGALi::resultant_(F,G) eliminates the innermost variable.
|
||||||
|
|
||||||
// the general function for CGAL::Integral_domain_without_div_tag
|
// Dispatching
|
||||||
template < class NT> inline
|
// CGAL::CGALi::new_resultant_decompose applies if Coeff is a Fraction
|
||||||
NT resultant_(Polynomial<NT> A,
|
// CGAL::CGALi::new_resultant_modularize applies if Coeff is Modularizable
|
||||||
Polynomial<NT> B,
|
// CGAL::CGALi::new_resultant_interpolate applies for multivairate polynomials
|
||||||
Integral_domain_without_division_tag){
|
// CGAL::CGALi::new_resultant_univariate selects the proper algorithm for IC
|
||||||
return hybrid_bezout_subresultant(A,B,0);
|
|
||||||
}
|
|
||||||
|
|
||||||
// the specialization for CGAL::Unique_factorization_domain_tag
|
namespace CGALi{
|
||||||
template < class NT> inline
|
|
||||||
NT resultant_(Polynomial<NT> A,
|
|
||||||
Polynomial<NT> B,
|
|
||||||
Unique_factorization_domain_tag ){
|
|
||||||
return prs_resultant_ufd(A,B);
|
|
||||||
}
|
|
||||||
|
|
||||||
template <class NT> inline
|
template <class Coeff>
|
||||||
NT resultant_field(Polynomial<NT> A, Polynomial<NT> B, ::CGAL::Tag_false) {
|
inline Coeff new_resultant_interpolate(
|
||||||
return prs_resultant_field(A,B);
|
const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>& );
|
||||||
}
|
template <class Coeff>
|
||||||
template <class NT> inline
|
inline Coeff new_resultant_modularize(
|
||||||
NT resultant_field(Polynomial<NT> A, Polynomial<NT> B, ::CGAL::Tag_true) {
|
const CGAL::Polynomial<Coeff>&,
|
||||||
return prs_resultant_decompose(A,B);
|
const CGAL::Polynomial<Coeff>&, CGAL::Tag_true);
|
||||||
}
|
template <class Coeff>
|
||||||
|
inline Coeff new_resultant_modularize(
|
||||||
|
const CGAL::Polynomial<Coeff>&,
|
||||||
|
const CGAL::Polynomial<Coeff>&, CGAL::Tag_false);
|
||||||
|
template <class Coeff>
|
||||||
|
inline Coeff new_resultant_decompose(
|
||||||
|
const CGAL::Polynomial<Coeff>&,
|
||||||
|
const CGAL::Polynomial<Coeff>&, CGAL::Tag_true);
|
||||||
|
template <class Coeff>
|
||||||
|
inline Coeff new_resultant_decompose(
|
||||||
|
const CGAL::Polynomial<Coeff>&,
|
||||||
|
const CGAL::Polynomial<Coeff>&, CGAL::Tag_false);
|
||||||
|
template <class Coeff>
|
||||||
|
inline Coeff new_resultant_(
|
||||||
|
const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&);
|
||||||
|
|
||||||
// the specialization for CGAL::Field_tag
|
template <class Coeff>
|
||||||
template < class NT> inline
|
inline Coeff new_resultant_univariate(
|
||||||
NT resultant_(Polynomial<NT> A, Polynomial<NT> B, Field_tag){
|
const CGAL::Polynomial<Coeff>& A,
|
||||||
// prs_resultant_decompose can only be used if
|
const CGAL::Polynomial<Coeff>& B, CGAL::Integral_domain_without_division_tag){
|
||||||
// NT is decomposable && the resulting type is at least a UFDomain
|
return hybrid_bezout_subresultant(A,B,0);
|
||||||
typedef typename Fraction_traits<NT>::Is_fraction Is_decomposable;
|
}
|
||||||
const bool is_decomposable
|
template <class Coeff>
|
||||||
= ::boost::is_same<Is_decomposable, ::CGAL::Tag_true>::value;
|
inline Coeff new_resultant_univariate(
|
||||||
|
const CGAL::Polynomial<Coeff>& A,
|
||||||
|
const CGAL::Polynomial<Coeff>& B, CGAL::Unique_factorization_domain_tag){
|
||||||
|
return prs_resultant_ufd(A,B);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class Coeff>
|
||||||
typedef typename Fraction_traits<NT>::Numerator_type NUM;
|
inline Coeff new_resultant_univariate(
|
||||||
typedef typename Algebraic_structure_traits<NUM>::Algebraic_category ALG_TYPE;
|
const CGAL::Polynomial<Coeff>& A,
|
||||||
const bool is_inheritance
|
const CGAL::Polynomial<Coeff>& B, CGAL::Field_tag){
|
||||||
= ::boost::is_base_and_derived<Unique_factorization_domain_tag,ALG_TYPE>::value ||
|
return prs_resultant_field(A,B);
|
||||||
::boost::is_same<Unique_factorization_domain_tag,ALG_TYPE>::value;
|
}
|
||||||
|
|
||||||
return resultant_field(
|
|
||||||
A, B,
|
|
||||||
typename ::boost::mpl::if_c< is_decomposable && is_inheritance,
|
|
||||||
::CGAL::Tag_true,
|
|
||||||
::CGAL::Tag_false >::type() );
|
|
||||||
}
|
|
||||||
|
|
||||||
/*! \ingroup CGAL_Polynomial
|
|
||||||
* \relates CGAL::Polynomial
|
|
||||||
* \brief compute the resultant of the polynomials \c A and \c B
|
|
||||||
*
|
|
||||||
* The way the resultant is computed depends on the Algebraic_category.
|
|
||||||
* In general the resultant will be computed by the function
|
|
||||||
* CGAL::hybrid_bezout_subresultant, but if possible the function
|
|
||||||
* CGAL::prs_resultant_ufd or CGAL::prs_resultant_field are used.
|
|
||||||
*
|
|
||||||
* Up to now it is not clear, that the functions based on the polynomial
|
|
||||||
* remainder sequence are faster than the one based on the bezoutian.
|
|
||||||
* Thus you can use CGAL::hybrid_bezout_subresultant instead, which will
|
|
||||||
* work for any Algebraic_category
|
|
||||||
*/
|
|
||||||
template <class NT> inline
|
|
||||||
NT resultant(Polynomial<NT> A, Polynomial<NT> B) {
|
|
||||||
return CGALi::resultant_(A, B);
|
|
||||||
}
|
|
||||||
|
|
||||||
} // namespace CGALi
|
} // namespace CGALi
|
||||||
|
|
||||||
|
namespace CGALi{
|
||||||
|
|
||||||
|
|
||||||
|
template <class IC>
|
||||||
|
inline IC
|
||||||
|
new_resultant_interpolate(
|
||||||
|
const CGAL::Polynomial<IC>& F,
|
||||||
|
const CGAL::Polynomial<IC>& G){
|
||||||
|
CGAL_precondition(CGAL::Polynomial_traits_d<CGAL::Polynomial<IC> >::d == 1);
|
||||||
|
typedef CGAL::Algebraic_structure_traits<IC> AST_IC;
|
||||||
|
typedef typename AST_IC::Algebraic_category Algebraic_category;
|
||||||
|
return CGALi::new_resultant_univariate(F,G,Algebraic_category());
|
||||||
|
}
|
||||||
|
|
||||||
|
#if CGAL_RESULTANT_USE_INTERPOLATION
|
||||||
|
template <class Coeff_2>
|
||||||
|
inline
|
||||||
|
CGAL::Polynomial<Coeff_2> new_resultant_interpolate(
|
||||||
|
const CGAL::Polynomial<CGAL::Polynomial<Coeff_2> >& F,
|
||||||
|
const CGAL::Polynomial<CGAL::Polynomial<Coeff_2> >& G){
|
||||||
|
|
||||||
|
typedef CGAL::Polynomial<Coeff_2> Coeff_1;
|
||||||
|
typedef CGAL::Polynomial<Coeff_1> POLY;
|
||||||
|
typedef CGAL::Polynomial_traits_d<POLY> PT;
|
||||||
|
typedef typename PT::Innermost_coefficient IC;
|
||||||
|
|
||||||
|
CGAL_precondition(PT::d >= 2);
|
||||||
|
|
||||||
|
typename PT::Degree degree;
|
||||||
|
typename CGAL::Polynomial_traits_d<Coeff_1>::Degree_vector degree_vector;
|
||||||
|
|
||||||
|
int maxdegree = degree(F,0)*degree(G,PT::d-1) + degree(F,PT::d-1)*degree(G,0);
|
||||||
|
|
||||||
|
typedef std::pair<IC,Coeff_2> Point;
|
||||||
|
std::vector<Point> points; // interpolation points
|
||||||
|
|
||||||
|
|
||||||
|
int i(0);
|
||||||
|
CGAL::Exponent_vector ev_f(degree_vector(Coeff_1()));
|
||||||
|
CGAL::Exponent_vector ev_g(degree_vector(Coeff_1()));
|
||||||
|
|
||||||
|
|
||||||
|
while((int) points.size() <= maxdegree + 1){
|
||||||
|
i++;
|
||||||
|
// timer1.start();
|
||||||
|
Coeff_1 Fat_i = F.evaluate(Coeff_1(i));
|
||||||
|
Coeff_1 Gat_i = G.evaluate(Coeff_1(i));
|
||||||
|
// timer1.stop();
|
||||||
|
|
||||||
|
if(degree_vector(Fat_i) > ev_f || degree_vector(Gat_i) > ev_g){
|
||||||
|
points.clear();
|
||||||
|
ev_f = degree_vector(Fat_i);
|
||||||
|
ev_g = degree_vector(Gat_i);
|
||||||
|
CGAL_postcondition(points.size() == 0);
|
||||||
|
}
|
||||||
|
if(degree_vector(Fat_i) == ev_f && degree_vector(Gat_i) == ev_g){
|
||||||
|
// timer2.start();
|
||||||
|
Coeff_2 res_at_i = new_resultant_interpolate(Fat_i, Gat_i);
|
||||||
|
// timer2.stop();
|
||||||
|
points.push_back(Point(IC(i),res_at_i));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// timer3.start();
|
||||||
|
CGAL::CGALi::Interpolator<Coeff_1> interpolator(points.begin(),points.end());
|
||||||
|
Coeff_1 result = interpolator.get_interpolant();
|
||||||
|
// timer3.stop();
|
||||||
|
|
||||||
|
#ifndef CGAL_NDEBUG
|
||||||
|
while((int) points.size() <= maxdegree + 3){
|
||||||
|
i++;
|
||||||
|
Coeff_1 Fat_i = typename PT::Evaluate()(F,IC(i));
|
||||||
|
Coeff_1 Gat_i = typename PT::Evaluate()(G,IC(i));
|
||||||
|
|
||||||
|
assert(degree_vector(Fat_i) <= ev_f);
|
||||||
|
assert(degree_vector(Gat_i) <= ev_g);
|
||||||
|
|
||||||
|
if(degree_vector(Fat_i) == ev_f && degree_vector(Gat_i) == ev_g){
|
||||||
|
Coeff_2 res_at_i = new_resultant_interpolate(Fat_i, Gat_i);
|
||||||
|
points.push_back(Point(IC(i), res_at_i));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
CGAL::CGALi::Interpolator<Coeff_1>
|
||||||
|
interpolator_(points.begin(),points.end());
|
||||||
|
Coeff_1 result_= interpolator_.get_interpolant();
|
||||||
|
|
||||||
|
// the interpolate polynomial has to be stable !
|
||||||
|
assert(result_ == result);
|
||||||
|
#endif
|
||||||
|
return result;
|
||||||
|
}
|
||||||
|
#endif
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
template <class Coeff>
|
||||||
|
inline
|
||||||
|
Coeff new_resultant_modularize(
|
||||||
|
const CGAL::Polynomial<Coeff>& F,
|
||||||
|
const CGAL::Polynomial<Coeff>& G,
|
||||||
|
CGAL::Tag_false){
|
||||||
|
return new_resultant_interpolate(F,G);
|
||||||
|
};
|
||||||
|
|
||||||
|
template <class Coeff>
|
||||||
|
inline
|
||||||
|
Coeff new_resultant_modularize(
|
||||||
|
const CGAL::Polynomial<Coeff>& F,
|
||||||
|
const CGAL::Polynomial<Coeff>& G,
|
||||||
|
CGAL::Tag_true){
|
||||||
|
|
||||||
|
typedef Polynomial_traits_d<CGAL::Polynomial<Coeff> > PT;
|
||||||
|
typedef typename PT::Polynomial_d Polynomial;
|
||||||
|
typedef typename PT::Innermost_coefficient IC;
|
||||||
|
|
||||||
|
|
||||||
|
typedef Chinese_remainder_traits<Coeff> CRT;
|
||||||
|
typedef typename CRT::Scalar_type Scalar;
|
||||||
|
|
||||||
|
|
||||||
|
typedef typename CGAL::Modular_traits<Polynomial>::Modular_NT MPolynomial;
|
||||||
|
typedef typename CGAL::Modular_traits<Coeff>::Modular_NT MCoeff;
|
||||||
|
typedef typename CGAL::Modular_traits<Scalar>::Modular_NT MScalar;
|
||||||
|
|
||||||
|
typename CRT::Chinese_remainder chinese_remainder;
|
||||||
|
typename CGAL::Modular_traits<Coeff>::Modular_image_inv inv_map;
|
||||||
|
|
||||||
|
|
||||||
|
typename PT::Degree_vector degree_vector;
|
||||||
|
typename CGAL::Polynomial_traits_d<MPolynomial>::Degree_vector mdegree_vector;
|
||||||
|
|
||||||
|
bool solved = false;
|
||||||
|
int prime_index = 0;
|
||||||
|
int n = 0;
|
||||||
|
Scalar p,q,pq,s,t;
|
||||||
|
Coeff R, R_old;
|
||||||
|
|
||||||
|
// CGAL::Timer timer_evaluate, timer_resultant, timer_cr;
|
||||||
|
|
||||||
|
do{
|
||||||
|
MPolynomial mF, mG;
|
||||||
|
MCoeff mR;
|
||||||
|
//timer_evaluate.start();
|
||||||
|
do{
|
||||||
|
// select a prime number
|
||||||
|
int current_prime = -1;
|
||||||
|
prime_index++;
|
||||||
|
if(prime_index >= 2000){
|
||||||
|
std::cerr<<"primes in the array exhausted"<<std::endl;
|
||||||
|
assert(false);
|
||||||
|
current_prime = CGALi::get_next_lower_prime(current_prime);
|
||||||
|
} else{
|
||||||
|
current_prime = CGALi::primes[prime_index];
|
||||||
|
}
|
||||||
|
CGAL::Modular::set_current_prime(current_prime);
|
||||||
|
|
||||||
|
mF = CGAL::modular_image(F);
|
||||||
|
mG = CGAL::modular_image(G);
|
||||||
|
|
||||||
|
}while( degree_vector(F) != mdegree_vector(mF) ||
|
||||||
|
degree_vector(G) != mdegree_vector(mG));
|
||||||
|
//timer_evaluate.stop();
|
||||||
|
|
||||||
|
//timer_resultant.start();
|
||||||
|
n++;
|
||||||
|
mR = new_resultant_interpolate(mF,mG);
|
||||||
|
//timer_resultant.stop();
|
||||||
|
//timer_cr.start();
|
||||||
|
if(n == 1){
|
||||||
|
// init chinese remainder
|
||||||
|
q = CGAL::Modular::get_current_prime(); // implicit !
|
||||||
|
R = inv_map(mR);
|
||||||
|
}else{
|
||||||
|
// continue chinese remainder
|
||||||
|
p = CGAL::Modular::get_current_prime(); // implicit!
|
||||||
|
R_old = R ;
|
||||||
|
// chinese_remainder(q,Gs ,p,inv_map(mG_),pq,Gs);
|
||||||
|
// cached_extended_euclidean_algorithm(q,p,s,t);
|
||||||
|
CGALi::Cached_extended_euclidean_algorithm
|
||||||
|
<typename CRT::Scalar_type> ceea;
|
||||||
|
ceea(q,p,s,t);
|
||||||
|
pq =p*q;
|
||||||
|
chinese_remainder(q,p,pq,s,t,R_old,inv_map(mR),R);
|
||||||
|
q=pq;
|
||||||
|
}
|
||||||
|
solved = (R==R_old);
|
||||||
|
//timer_cr.stop();
|
||||||
|
} while(!solved);
|
||||||
|
|
||||||
|
//std::cout << "Time Evaluate : " << timer_evaluate.time() << std::endl;
|
||||||
|
//std::cout << "Time Resultant : " << timer_resultant.time() << std::endl;
|
||||||
|
//std::cout << "Time Chinese R : " << timer_cr.time() << std::endl;
|
||||||
|
// CGAL_postcondition(R == new_resultant_interpolate(F,G));
|
||||||
|
return R;
|
||||||
|
// return new_resultant_interpolate(F,G);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
template <class Coeff>
|
||||||
|
inline
|
||||||
|
Coeff new_resultant_decompose(
|
||||||
|
const CGAL::Polynomial<Coeff>& F,
|
||||||
|
const CGAL::Polynomial<Coeff>& G,
|
||||||
|
CGAL::Tag_false){
|
||||||
|
#if CGAL_RESULTANT_USE_MODULAR_ARITHMETIC
|
||||||
|
typedef CGAL::Polynomial<Coeff> Polynomial;
|
||||||
|
typedef typename Modular_traits<Polynomial>::Is_modularizable Is_modularizable;
|
||||||
|
return new_resultant_modularize(F,G,Is_modularizable());
|
||||||
|
#else
|
||||||
|
return new_resultant_modularize(F,G,CGAL::Tag_false());
|
||||||
|
#endif
|
||||||
|
}
|
||||||
|
|
||||||
|
template <class Coeff>
|
||||||
|
inline
|
||||||
|
Coeff new_resultant_decompose(
|
||||||
|
const CGAL::Polynomial<Coeff>& F,
|
||||||
|
const CGAL::Polynomial<Coeff>& G,
|
||||||
|
CGAL::Tag_true){
|
||||||
|
|
||||||
|
typedef Polynomial<Coeff> POLY;
|
||||||
|
typedef typename Fraction_traits<POLY>::Numerator_type Numerator;
|
||||||
|
typedef typename Fraction_traits<POLY>::Denominator_type Denominator;
|
||||||
|
typename Fraction_traits<POLY>::Decompose decompose;
|
||||||
|
typedef typename Numerator::NT RES;
|
||||||
|
|
||||||
|
Denominator a, b;
|
||||||
|
// F.simplify_coefficients(); not const
|
||||||
|
// G.simplify_coefficients(); not const
|
||||||
|
Numerator F0; decompose(F,F0,a);
|
||||||
|
Numerator G0; decompose(G,G0,b);
|
||||||
|
Denominator c = CGAL::ipower(a, G.degree()) * CGAL::ipower(b, F.degree());
|
||||||
|
typedef Algebraic_structure_traits<RES> AST_RES;
|
||||||
|
typedef typename AST_RES::Algebraic_category Algebraic_category;
|
||||||
|
RES res0 = CGAL::CGALi::new_resultant_(F0, G0);
|
||||||
|
typename Fraction_traits<Coeff>::Compose comp_frac;
|
||||||
|
Coeff res = comp_frac(res0, c);
|
||||||
|
typename Algebraic_structure_traits<Coeff>::Simplify simplify;
|
||||||
|
simplify(res);
|
||||||
|
return res;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
template <class Coeff>
|
||||||
|
inline
|
||||||
|
Coeff new_resultant_(
|
||||||
|
const CGAL::Polynomial<Coeff>& F,
|
||||||
|
const CGAL::Polynomial<Coeff>& G){
|
||||||
|
#if CGAL_RESULTANT_USE_DECOMPOSE
|
||||||
|
typedef CGAL::Fraction_traits<Polynomial<Coeff > > FT;
|
||||||
|
typedef typename FT::Is_fraction Is_fraction;
|
||||||
|
return new_resultant_decompose(F,G,Is_fraction());
|
||||||
|
#else
|
||||||
|
return new_resultant_decompose(F,G,CGAL::Tag_false());
|
||||||
|
#endif
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
template <class Coeff>
|
||||||
|
inline
|
||||||
|
Coeff new_resultant(
|
||||||
|
const CGAL::Polynomial<Coeff>& F_,
|
||||||
|
const CGAL::Polynomial<Coeff>& G_){
|
||||||
|
// make the variable to be elimnated the innermost one.
|
||||||
|
typedef CGAL::Polynomial_traits_d<CGAL::Polynomial<Coeff> > PT;
|
||||||
|
CGAL::Polynomial<Coeff> F = typename PT::Move()(F_, PT::d-1, 0);
|
||||||
|
CGAL::Polynomial<Coeff> G = typename PT::Move()(G_, PT::d-1, 0);
|
||||||
|
return CGALi::new_resultant_(F,G);
|
||||||
|
}
|
||||||
|
|
||||||
|
} // namespace CGALi
|
||||||
CGAL_END_NAMESPACE
|
CGAL_END_NAMESPACE
|
||||||
|
|
||||||
#endif// CGAL_POLYNOMIAL_RESULTANT_H
|
|
||||||
|
|
||||||
|
#endif // CGAL_NEW_RESULTANT_H
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -1145,7 +1145,7 @@ void test_compare(const Polynomial_traits_d&) {
|
||||||
|
|
||||||
CGAL_SNAP_CGALi_TRAITS_D(Polynomial_traits_d);
|
CGAL_SNAP_CGALi_TRAITS_D(Polynomial_traits_d);
|
||||||
|
|
||||||
typename PT::Compare compare;
|
typename PT::Compare compare; (void) compare;
|
||||||
|
|
||||||
Polynomial_d p0(Coeff(0));
|
Polynomial_d p0(Coeff(0));
|
||||||
Polynomial_d pp2(Coeff(2));
|
Polynomial_d pp2(Coeff(2));
|
||||||
|
|
|
||||||
|
|
@ -1,129 +1,123 @@
|
||||||
|
|
||||||
// #define CGAL_RESULTANT_NUSE_MODULAR_ARITHMETIC 1
|
// currently using modular arithmetic is far to slow,
|
||||||
|
// that is, we use interpolation for multivariate resultants
|
||||||
|
// but no modular arithmetic.
|
||||||
|
|
||||||
|
// These are the defaults
|
||||||
|
//#define CGAL_RESULTANT_USE_MODULAR_ARITHMETIC 0
|
||||||
|
//#define CGAL_RESULTANT_USE_DECOMPOSE 1
|
||||||
|
//#define CGAL_RESULTANT_USE_INTERPOLATE 1
|
||||||
|
|
||||||
#include <CGAL/basic.h>
|
#include <CGAL/basic.h>
|
||||||
|
|
||||||
#include <CGAL/Timer.h>
|
|
||||||
CGAL::Timer timer1, timer2, timer3;
|
|
||||||
|
|
||||||
#include <cassert>
|
#include <cassert>
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
|
|
||||||
|
#include <CGAL/Polynomial/resultant.h>
|
||||||
|
|
||||||
#include <CGAL/Arithmetic_kernel.h>
|
#include <CGAL/Arithmetic_kernel.h>
|
||||||
#include <CGAL/Modular.h>
|
#include <CGAL/Modular.h>
|
||||||
|
#include <CGAL/Sqrt_extension.h>
|
||||||
#include <CGAL/Polynomial.h>
|
#include <CGAL/Polynomial.h>
|
||||||
#include <CGAL/Polynomial_traits_d.h>
|
#include <CGAL/Polynomial_traits_d.h>
|
||||||
|
|
||||||
#include <CGAL/Random.h>
|
#include <CGAL/Random.h>
|
||||||
|
#include <CGAL/Timer.h>
|
||||||
|
|
||||||
#include <CGAL/gen_sparse_polynomial.h>
|
#include <CGAL/gen_sparse_polynomial.h>
|
||||||
#include <CGAL/new_resultant.h>
|
|
||||||
|
|
||||||
|
|
||||||
static CGAL::Random my_rnd(346); // some seed
|
static CGAL::Random my_rnd(346); // some seed
|
||||||
|
|
||||||
template<class Polynomial_d>
|
template<class Polynomial_d>
|
||||||
void test_resultant(){
|
void test_resultant(){
|
||||||
|
typedef typename CGAL::Polynomial_traits_d<Polynomial_d> PT;
|
||||||
|
typedef typename PT::Coefficient Coeff;
|
||||||
|
typename PT::Resultant resultant;
|
||||||
|
|
||||||
typedef typename CGAL::Polynomial_traits_d<Polynomial_d> PT;
|
{
|
||||||
typedef typename PT::Coefficient Coeff;
|
Polynomial_d A(0);
|
||||||
for (int k = 0; k < 1; k++){
|
Polynomial_d B(0);
|
||||||
Polynomial_d F2 =
|
assert(resultant(A,B)==Coeff(0));
|
||||||
CGAL::generate_sparse_random_polynomial<Polynomial_d>(my_rnd,12);
|
}{
|
||||||
Polynomial_d F1 =
|
Polynomial_d A(4);
|
||||||
CGAL::generate_sparse_random_polynomial<Polynomial_d>(my_rnd,12);
|
Polynomial_d B(8);
|
||||||
Polynomial_d G1 = typename PT::Move()(F1,PT::d-1,0);
|
assert(resultant(A,B)==Coeff(1));
|
||||||
Polynomial_d G2 = typename PT::Move()(F2,PT::d-1,0);
|
}{
|
||||||
|
Polynomial_d f(Coeff(2),Coeff(7),Coeff(1),Coeff(8),Coeff(1),Coeff(8));
|
||||||
|
Polynomial_d g(Coeff(3),Coeff(1),Coeff(4),Coeff(1),Coeff(5),Coeff(9));
|
||||||
|
assert(resultant(f,g) == Coeff(230664271L)); // Maple
|
||||||
|
|
||||||
CGAL::Timer timer_new;
|
Polynomial_d h(Coeff(3),Coeff(4),Coeff(7),Coeff(7));
|
||||||
timer_new.start();
|
Polynomial_d fh(f*h);
|
||||||
Coeff F = CGAL::new_resultant(F1,F2,PT::d-1);
|
Polynomial_d gh(g*h);
|
||||||
// Coeff G = CGAL::new_resultant(G1,G2,0);
|
assert(resultant(fh,gh) == Coeff(0) );
|
||||||
timer_new.stop();
|
}
|
||||||
|
|
||||||
assert(F==F);
|
|
||||||
std::cout <<"timer1: "<< timer1.time() << std::endl;
|
|
||||||
std::cout <<"timer2: "<< timer2.time() << std::endl;
|
|
||||||
std::cout <<"timer3: "<< timer3.time() << std::endl;
|
|
||||||
timer1.reset();
|
|
||||||
timer2.reset();
|
|
||||||
timer3.reset();
|
|
||||||
std::cout <<"new time: "<< timer_new.time() << std::endl;
|
|
||||||
}
|
|
||||||
std::cout << "end resultant: " << PT::d << std::endl;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
template<class Coeff>
|
for (int k = 0; k < 1; k++){
|
||||||
void test_multi_resultant(){
|
Polynomial_d F2 =
|
||||||
|
CGAL::generate_sparse_random_polynomial<Polynomial_d>(my_rnd,5);
|
||||||
typedef CGAL::Polynomial<Coeff> Polynomial_1;
|
Polynomial_d F1 =
|
||||||
typedef CGAL::Polynomial<Polynomial_1> Polynomial_2;
|
CGAL::generate_sparse_random_polynomial<Polynomial_d>(my_rnd,5);
|
||||||
typedef CGAL::Polynomial<Polynomial_2> Polynomial_3;
|
Polynomial_d G1 = typename PT::Move()(F1,PT::d-1,0);
|
||||||
|
Polynomial_d G2 = typename PT::Move()(F2,PT::d-1,0);
|
||||||
Polynomial_3 Q =
|
|
||||||
CGAL::generate_polynomial_degree_each<Polynomial_3>(my_rnd,4,100);
|
|
||||||
Polynomial_2 B1 =
|
|
||||||
CGAL::generate_polynomial_degree_each<Polynomial_2>(my_rnd,6,80);
|
|
||||||
Polynomial_2 B2 =
|
|
||||||
CGAL::generate_polynomial_degree_each<Polynomial_2>(my_rnd,6,80);
|
|
||||||
std::cout << "Q: "<< Q << std::endl;
|
|
||||||
std::cout << "B1: "<< B1 << std::endl;
|
|
||||||
std::cout << "B2: "<< B2 <<std::endl;
|
|
||||||
|
|
||||||
CGAL::Timer timer_new;
|
CGAL::Timer timer_new;
|
||||||
timer_new.start();
|
timer_new.start();
|
||||||
Polynomial_2 R1 = CGAL::new_resultant(Q,Polynomial_3(B1),0);
|
Coeff F = resultant(F1,F2,PT::d-1);
|
||||||
Polynomial_1 R2 = CGAL::new_resultant(R1,B2,0);
|
Coeff G = resultant(G1,G2,0);
|
||||||
std::cout << std::endl;
|
|
||||||
std::cout << R1.degree() << std::endl;
|
|
||||||
std::cout << R2.degree() << std::endl;
|
|
||||||
timer_new.stop();
|
timer_new.stop();
|
||||||
|
|
||||||
std::cout <<"timer1: "<< timer1.time() << std::endl;
|
assert(F==G);
|
||||||
std::cout <<"timer2: "<< timer2.time() << std::endl;
|
}
|
||||||
std::cout <<"timer3: "<< timer3.time() << std::endl;
|
std::cout << "end resultant: " << PT::d << std::endl;
|
||||||
timer1.reset();
|
|
||||||
timer2.reset();
|
|
||||||
timer3.reset();
|
|
||||||
|
|
||||||
std::cout <<"new time: "<< timer_new.time() << std::endl;
|
|
||||||
std::cout << "end multi resultant: " << std::endl;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
int main(){
|
int main(){
|
||||||
|
|
||||||
//CGAL::force_ieee_double_precision();
|
CGAL::Timer timer;
|
||||||
CGAL::set_pretty_mode(std::cout);
|
timer.start();
|
||||||
|
//CGAL::force_ieee_double_precision();
|
||||||
|
CGAL::set_pretty_mode(std::cout);
|
||||||
|
|
||||||
typedef CGAL::Arithmetic_kernel AK;
|
typedef CGAL::Arithmetic_kernel AK;
|
||||||
typedef AK::Integer Integer;
|
typedef AK::Integer Integer;
|
||||||
typedef CGAL::Polynomial<Integer> Polynomial_1;
|
typedef CGAL::Polynomial<Integer> Polynomial_1;
|
||||||
typedef CGAL::Polynomial<Polynomial_1> Polynomial_2;
|
typedef CGAL::Polynomial<Polynomial_1> Polynomial_2;
|
||||||
typedef CGAL::Polynomial<Polynomial_2> Polynomial_3;
|
typedef CGAL::Polynomial<Polynomial_2> Polynomial_3;
|
||||||
|
|
||||||
test_multi_resultant<Integer>();
|
test_resultant<Polynomial_1>();
|
||||||
|
test_resultant<Polynomial_2>();
|
||||||
|
test_resultant<Polynomial_3>();
|
||||||
|
|
||||||
// test_resultant<Polynomial_1>();
|
typedef CGAL::Polynomial<CGAL::Modular> MPolynomial_1;
|
||||||
// test_resultant<Polynomial_2>();
|
typedef CGAL::Polynomial<MPolynomial_1> MPolynomial_2;
|
||||||
//test_resultant<Polynomial_3>();
|
typedef CGAL::Polynomial<MPolynomial_2> MPolynomial_3;
|
||||||
|
|
||||||
|
test_resultant<MPolynomial_1>();
|
||||||
|
test_resultant<MPolynomial_2>();
|
||||||
|
test_resultant<MPolynomial_3>();
|
||||||
|
|
||||||
typedef CGAL::Polynomial<CGAL::Modular> MPolynomial_1;
|
typedef CGAL::Polynomial<AK::Rational> RPolynomial_1;
|
||||||
typedef CGAL::Polynomial<MPolynomial_1> MPolynomial_2;
|
typedef CGAL::Polynomial<RPolynomial_1> RPolynomial_2;
|
||||||
typedef CGAL::Polynomial<MPolynomial_2> MPolynomial_3;
|
typedef CGAL::Polynomial<RPolynomial_2> RPolynomial_3;
|
||||||
|
|
||||||
//test_resultant<MPolynomial_1>();
|
test_resultant<RPolynomial_1>();
|
||||||
//test_resultant<MPolynomial_2>();
|
test_resultant<RPolynomial_2>();
|
||||||
//test_resultant<MPolynomial_3>();
|
test_resultant<RPolynomial_3>();
|
||||||
|
|
||||||
typedef CGAL::Polynomial<AK::Rational> RPolynomial_1;
|
typedef CGAL::Sqrt_extension<AK::Integer, AK::Integer> EXT;
|
||||||
typedef CGAL::Polynomial<RPolynomial_1> RPolynomial_2;
|
typedef CGAL::Polynomial<EXT> EPolynomial_1;
|
||||||
typedef CGAL::Polynomial<RPolynomial_2> RPolynomial_3;
|
typedef CGAL::Polynomial<EPolynomial_1> EPolynomial_2;
|
||||||
|
typedef CGAL::Polynomial<EPolynomial_2> EPolynomial_3;
|
||||||
|
|
||||||
//test_resultant<RPolynomial_1>();
|
test_resultant<EPolynomial_1>();
|
||||||
//test_resultant<RPolynomial_2>();
|
test_resultant<EPolynomial_2>();
|
||||||
//test_resultant<RPolynomial_3>();
|
test_resultant<EPolynomial_3>();
|
||||||
|
timer.stop();
|
||||||
|
|
||||||
|
std::cout <<" TOTAL TIME: " << timer.time();
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue