mirror of https://github.com/CGAL/cgal
new multivariate resultant using interpolation
This commit is contained in:
Michael Hemmer
parent
147b11c193
commit
5808b8076e
|
|
@ -10,21 +10,21 @@
|
|||
CGAL_BEGIN_NAMESPACE
|
||||
|
||||
template <class Coeff>
|
||||
inline Coeff new_resultant( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff>, int);
|
||||
inline Coeff new_resultant( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&, int);
|
||||
namespace CGALi{
|
||||
|
||||
template <class Coeff>
|
||||
inline Coeff new_resultant_interpolate( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff> );
|
||||
inline Coeff new_resultant_interpolate( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>& );
|
||||
template <class Coeff>
|
||||
inline Coeff new_resultant_modularize( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff>, CGAL::Tag_true);
|
||||
inline Coeff new_resultant_modularize( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&, CGAL::Tag_true);
|
||||
template <class Coeff>
|
||||
inline Coeff new_resultant_modularize( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff>, CGAL::Tag_false);
|
||||
inline Coeff new_resultant_modularize( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&, CGAL::Tag_false);
|
||||
template <class Coeff>
|
||||
inline Coeff new_resultant_decompose( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff>, CGAL::Tag_true);
|
||||
inline Coeff new_resultant_decompose( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&, CGAL::Tag_true);
|
||||
template <class Coeff>
|
||||
inline Coeff new_resultant_decompose( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff>, CGAL::Tag_false);
|
||||
inline Coeff new_resultant_decompose( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&, CGAL::Tag_false);
|
||||
template <class Coeff>
|
||||
inline Coeff new_resultant_( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff>);
|
||||
inline Coeff new_resultant_( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&);
|
||||
|
||||
} // namespace CGALi
|
||||
|
||||
|
|
@ -32,7 +32,7 @@ namespace CGALi{
|
|||
|
||||
template <class IC>
|
||||
inline IC
|
||||
new_resultant_interpolate( CGAL::Polynomial<IC> F, CGAL::Polynomial<IC> G){
|
||||
new_resultant_interpolate( const CGAL::Polynomial<IC>& F, const CGAL::Polynomial<IC>& G){
|
||||
CGAL_precondition(CGAL::Polynomial_traits_d<CGAL::Polynomial<IC> >::d == 1);
|
||||
return CGALi::resultant(F,G);
|
||||
}
|
||||
|
|
@ -40,9 +40,8 @@ new_resultant_interpolate( CGAL::Polynomial<IC> F, CGAL::Polynomial<IC> G){
|
|||
template <class Coeff_2>
|
||||
inline
|
||||
CGAL::Polynomial<Coeff_2> new_resultant_interpolate(
|
||||
CGAL::Polynomial<CGAL::Polynomial<Coeff_2> > F,
|
||||
CGAL::Polynomial<CGAL::Polynomial<Coeff_2> > G,
|
||||
int index = 0 ){
|
||||
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;
|
||||
|
|
@ -50,59 +49,65 @@ CGAL::Polynomial<Coeff_2> new_resultant_interpolate(
|
|||
typedef typename PT::Innermost_coefficient IC;
|
||||
|
||||
CGAL_precondition(PT::d >= 2);
|
||||
|
||||
CGAL_precondition(index >= 0);
|
||||
CGAL_precondition(index < PT::d);
|
||||
|
||||
POLY F_ = typename PT::Move()(F, index, 0);
|
||||
POLY G_ = typename PT::Move()(G, index, 0);
|
||||
CGAL_postcondition(index != 0 || F_ == F);
|
||||
|
||||
typename PT::Degree degree;
|
||||
int maxdegree = degree(F_,0)*degree(G_,PT::d-1) + degree(F_,PT::d-1)*degree(G_,0);
|
||||
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
|
||||
|
||||
typename CGAL::Polynomial_traits_d<Coeff_1>::Degree_vector degree_vector;
|
||||
|
||||
|
||||
|
||||
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++;
|
||||
Coeff_1 F_at_i = typename PT::Evaluate()(F_,IC(i));
|
||||
Coeff_1 G_at_i = typename PT::Evaluate()(G_,IC(i));
|
||||
if(degree_vector(F_at_i) > ev_f || degree_vector(G_at_i) > ev_g){
|
||||
// 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(F_at_i);
|
||||
ev_g = degree_vector(G_at_i);
|
||||
ev_f = degree_vector(Fat_i);
|
||||
ev_g = degree_vector(Gat_i);
|
||||
CGAL_postcondition(points.size() == 0);
|
||||
}
|
||||
if(degree_vector(F_at_i) == ev_f && degree_vector(G_at_i) == ev_g){
|
||||
Coeff_2 res_at_i = new_resultant(F_at_i, G_at_i, 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));
|
||||
}
|
||||
}
|
||||
}
|
||||
Coeff_1 result = CGAL::Interpolator<Coeff_1>(points.begin(),points.end()).get_interpolant();
|
||||
|
||||
|
||||
// timer3.start();
|
||||
CGAL::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 F_at_i = typename PT::Evaluate()(F_,IC(i));
|
||||
Coeff_1 G_at_i = typename PT::Evaluate()(G_,IC(i));
|
||||
Coeff_1 Fat_i = typename PT::Evaluate()(F,IC(i));
|
||||
Coeff_1 Gat_i = typename PT::Evaluate()(G,IC(i));
|
||||
|
||||
assert(degree_vector(F_at_i) <= ev_f);
|
||||
assert(degree_vector(G_at_i) <= ev_g);
|
||||
assert(degree_vector(Fat_i) <= ev_f);
|
||||
assert(degree_vector(Gat_i) <= ev_g);
|
||||
|
||||
if(degree_vector(F_at_i) == ev_f && degree_vector(G_at_i) == ev_g){
|
||||
Coeff_2 res_at_i = new_resultant(F_at_i, G_at_i, 0);
|
||||
if(degree_vector(Fat_i) == ev_f && degree_vector(Gat_i) == ev_g){
|
||||
Coeff_2 res_at_i = new_resultant(Fat_i, Gat_i, 0);
|
||||
points.push_back(Point(IC(i), res_at_i));
|
||||
}
|
||||
}
|
||||
CGAL::Interpolator<Coeff_1> interpolator_(points.begin(),points.end());
|
||||
Coeff_1 result_= interpolator_.get_interpolant();
|
||||
|
||||
Coeff_1 result_ = CGAL::Interpolator<Coeff_1>(points.begin(), points.end()).get_interpolant();
|
||||
// the interpolate polynomial has to be stable !
|
||||
// the interpolate polynomial has to be stable !
|
||||
assert(result_ == result);
|
||||
#endif
|
||||
return result;
|
||||
|
|
@ -112,31 +117,119 @@ CGAL::Polynomial<Coeff_2> new_resultant_interpolate(
|
|||
template <class Coeff>
|
||||
inline
|
||||
Coeff new_resultant_modularize(
|
||||
CGAL::Polynomial<Coeff> F, CGAL::Polynomial<Coeff> G, CGAL::Tag_false){
|
||||
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(
|
||||
CGAL::Polynomial<Coeff> F, CGAL::Polynomial<Coeff> G, CGAL::Tag_true){
|
||||
return new_resultant_interpolate(F,G);
|
||||
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_algorithm_M::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(
|
||||
CGAL::Polynomial<Coeff> F, CGAL::Polynomial<Coeff> G, CGAL::Tag_false){
|
||||
const CGAL::Polynomial<Coeff>& F, const CGAL::Polynomial<Coeff>& G, CGAL::Tag_false){
|
||||
#ifdef CGAL_RESULTANT_NUSE_MODULAR_ARITHMETIC
|
||||
return new_resultant_modularize(F,G,CGAL::Tag_false());
|
||||
#else
|
||||
typedef CGAL::Polynomial<Coeff> Polynomial;
|
||||
typedef typename Modular_traits<Polynomial>::Is_modularizable Is_modularizable;
|
||||
return new_resultant_modularize(F,G,Is_modularizable());
|
||||
#endif
|
||||
}
|
||||
|
||||
template <class Coeff>
|
||||
inline
|
||||
Coeff new_resultant_decompose(
|
||||
CGAL::Polynomial<Coeff> F, CGAL::Polynomial<Coeff> G, CGAL::Tag_true){
|
||||
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;
|
||||
|
|
@ -163,7 +256,7 @@ Coeff new_resultant_decompose(
|
|||
template <class Coeff>
|
||||
inline
|
||||
Coeff new_resultant_(
|
||||
CGAL::Polynomial<Coeff> F, CGAL::Polynomial<Coeff> G){
|
||||
const CGAL::Polynomial<Coeff>& F, const CGAL::Polynomial<Coeff>& G){
|
||||
typedef CGAL::Fraction_traits<Polynomial<Coeff > > FT;
|
||||
typedef typename FT::Is_fraction Is_fraction;
|
||||
return new_resultant_decompose(F,G,Is_fraction());
|
||||
|
|
@ -175,13 +268,13 @@ Coeff new_resultant_(
|
|||
template <class Coeff>
|
||||
inline
|
||||
Coeff new_resultant(
|
||||
CGAL::Polynomial<Coeff> F,
|
||||
CGAL::Polynomial<Coeff> G,
|
||||
const CGAL::Polynomial<Coeff>& F_,
|
||||
const CGAL::Polynomial<Coeff>& G_,
|
||||
int index = CGAL::Polynomial_traits_d< CGAL::Polynomial<Coeff> >::d-1){
|
||||
typedef CGAL::Polynomial_traits_d<CGAL::Polynomial<Coeff> > PT;
|
||||
CGAL::Polynomial<Coeff> F_ = typename PT::Move()(F, index, 0);
|
||||
CGAL::Polynomial<Coeff> G_ = typename PT::Move()(G, index, 0);
|
||||
return CGALi::new_resultant_(F_,G_);
|
||||
CGAL::Polynomial<Coeff> F = typename PT::Move()(F_, index, 0);
|
||||
CGAL::Polynomial<Coeff> G = typename PT::Move()(G_, index, 0);
|
||||
return CGALi::new_resultant_(F,G);
|
||||
}
|
||||
|
||||
CGAL_END_NAMESPACE
|
||||
|
|
|
|||
|
|
@ -0,0 +1,130 @@
|
|||
|
||||
// #define CGAL_RESULTANT_NUSE_MODULAR_ARITHMETIC 1
|
||||
|
||||
#include <CGAL/basic.h>
|
||||
|
||||
#include <CGAL/Timer.h>
|
||||
CGAL::Timer timer1, timer2, timer3;
|
||||
|
||||
#include <cassert>
|
||||
#include <iostream>
|
||||
|
||||
#include <CGAL/Arithmetic_kernel.h>
|
||||
#include <CGAL/Modular.h>
|
||||
#include <CGAL/Polynomial.h>
|
||||
#include <CGAL/Polynomial_traits_d.h>
|
||||
#include <CGAL/Random.h>
|
||||
|
||||
#include <CGAL/gen_sparse_polynomial.h>
|
||||
#include <CGAL/new_resultant.h>
|
||||
|
||||
|
||||
static CGAL::Random my_rnd(346); // some seed
|
||||
|
||||
template<class Polynomial_d>
|
||||
void test_resultant(){
|
||||
|
||||
typedef typename CGAL::Polynomial_traits_d<Polynomial_d> PT;
|
||||
typedef typename PT::Coefficient Coeff;
|
||||
for (int k = 0; k < 1; k++){
|
||||
Polynomial_d F2 =
|
||||
CGAL::generate_sparse_random_polynomial<Polynomial_d>(my_rnd,12);
|
||||
Polynomial_d F1 =
|
||||
CGAL::generate_sparse_random_polynomial<Polynomial_d>(my_rnd,12);
|
||||
Polynomial_d G1 = typename PT::Move()(F1,PT::d-1,0);
|
||||
Polynomial_d G2 = typename PT::Move()(F2,PT::d-1,0);
|
||||
|
||||
CGAL::Timer timer_new;
|
||||
timer_new.start();
|
||||
Coeff F = CGAL::new_resultant(F1,F2,PT::d-1);
|
||||
// Coeff G = CGAL::new_resultant(G1,G2,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>
|
||||
void test_multi_resultant(){
|
||||
|
||||
typedef CGAL::Polynomial<Coeff> Polynomial_1;
|
||||
typedef CGAL::Polynomial<Polynomial_1> Polynomial_2;
|
||||
typedef CGAL::Polynomial<Polynomial_2> Polynomial_3;
|
||||
|
||||
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;
|
||||
timer_new.start();
|
||||
Polynomial_2 R1 = CGAL::new_resultant(Q,Polynomial_3(B1),0);
|
||||
Polynomial_1 R2 = CGAL::new_resultant(R1,B2,0);
|
||||
std::cout << std::endl;
|
||||
std::cout << R1.degree() << std::endl;
|
||||
std::cout << R2.degree() << std::endl;
|
||||
timer_new.stop();
|
||||
|
||||
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 multi resultant: " << std::endl;
|
||||
}
|
||||
|
||||
|
||||
int main(){
|
||||
|
||||
//CGAL::force_ieee_double_precision();
|
||||
CGAL::set_pretty_mode(std::cout);
|
||||
|
||||
typedef CGAL::Arithmetic_kernel AK;
|
||||
typedef AK::Integer Integer;
|
||||
typedef CGAL::Polynomial<Integer> Polynomial_1;
|
||||
typedef CGAL::Polynomial<Polynomial_1> Polynomial_2;
|
||||
typedef CGAL::Polynomial<Polynomial_2> Polynomial_3;
|
||||
|
||||
test_multi_resultant<Integer>();
|
||||
|
||||
// test_resultant<Polynomial_1>();
|
||||
// test_resultant<Polynomial_2>();
|
||||
//test_resultant<Polynomial_3>();
|
||||
|
||||
|
||||
typedef CGAL::Polynomial<CGAL::Modular> MPolynomial_1;
|
||||
typedef CGAL::Polynomial<MPolynomial_1> MPolynomial_2;
|
||||
typedef CGAL::Polynomial<MPolynomial_2> MPolynomial_3;
|
||||
|
||||
//test_resultant<MPolynomial_1>();
|
||||
//test_resultant<MPolynomial_2>();
|
||||
//test_resultant<MPolynomial_3>();
|
||||
|
||||
typedef CGAL::Polynomial<AK::Rational> RPolynomial_1;
|
||||
typedef CGAL::Polynomial<RPolynomial_1> RPolynomial_2;
|
||||
typedef CGAL::Polynomial<RPolynomial_2> RPolynomial_3;
|
||||
|
||||
//test_resultant<RPolynomial_1>();
|
||||
//test_resultant<RPolynomial_2>();
|
||||
//test_resultant<RPolynomial_3>();
|
||||
}
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue