1 // Copyright (c) 2008 Max-Planck-Institute Saarbruecken (Germany).
2 // All rights reserved.
3 //
4 // This file is part of CGAL (www.cgal.org)
5 //
6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Polynomial/include/CGAL/Polynomial/resultant.h $
7 // $Id: resultant.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot
8 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Michael Hemmer <hemmer@mpi-inf.mpg.de>
12
13
14 #ifndef CGAL_POLYNOMIAL_RESULTANT_H
15 #define CGAL_POLYNOMIAL_RESULTANT_H
16
17 // Modular arithmetic is slower, hence the default is 0
18 #ifndef CGAL_RESULTANT_USE_MODULAR_ARITHMETIC
19 #define CGAL_RESULTANT_USE_MODULAR_ARITHMETIC 0
20 #endif
21
22 #ifndef CGAL_RESULTANT_USE_DECOMPOSE
23 #define CGAL_RESULTANT_USE_DECOMPOSE 1
24 #endif
25
26
27 #include <CGAL/basic.h>
28 #include <CGAL/Polynomial.h>
29
30 #include <CGAL/Polynomial_traits_d.h>
31 #include <CGAL/Polynomial/Interpolator.h>
32 #include <CGAL/Polynomial/prs_resultant.h>
33 #include <CGAL/Polynomial/bezout_matrix.h>
34
35 #include <CGAL/Residue.h>
36 #include <CGAL/Modular_traits.h>
37 #include <CGAL/Chinese_remainder_traits.h>
38 #include <CGAL/primes.h>
39 #include <CGAL/Polynomial/Cached_extended_euclidean_algorithm.h>
40
41 namespace CGAL {
42
43
44 // The main function provided within this file is CGAL::internal::resultant(F,G),
45 // all other functions are used for dispatching.
46 // The implementation uses interpolatation for multivariate polynomials
47 // Due to the recursive structuture of CGAL::Polynomial<Coeff> it is better
48 // to write the function such that the inner most variabel is eliminated.
49 // However, CGAL::internal::resultant(F,G) eliminates the outer most variabel.
50 // This is due to backward compatibility issues with code base on EXACUS.
51 // In turn CGAL::internal::resultant_(F,G) eliminates the innermost variable.
52
53 // Dispatching
54 // CGAL::internal::resultant_decompose applies if Coeff is a Fraction
55 // CGAL::internal::resultant_modularize applies if Coeff is Modularizable
56 // CGAL::internal::resultant_interpolate applies for multivairate polynomials
57 // CGAL::internal::resultant_univariate selects the proper algorithm for IC
58
59 // CGAL_RESULTANT_USE_DECOMPOSE ( default = 1 )
60 // CGAL_RESULTANT_USE_MODULAR_ARITHMETIC (default = 0 )
61
62 namespace internal{
63
64 template <class Coeff>
65 inline Coeff resultant_interpolate(
66 const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>& );
67 template <class Coeff>
68 inline Coeff resultant_modularize(
69 const CGAL::Polynomial<Coeff>&,
70 const CGAL::Polynomial<Coeff>&, CGAL::Tag_true);
71 template <class Coeff>
72 inline Coeff resultant_modularize(
73 const CGAL::Polynomial<Coeff>&,
74 const CGAL::Polynomial<Coeff>&, CGAL::Tag_false);
75 template <class Coeff>
76 inline Coeff resultant_decompose(
77 const CGAL::Polynomial<Coeff>&,
78 const CGAL::Polynomial<Coeff>&, CGAL::Tag_true);
79 template <class Coeff>
80 inline Coeff resultant_decompose(
81 const CGAL::Polynomial<Coeff>&,
82 const CGAL::Polynomial<Coeff>&, CGAL::Tag_false);
83 template <class Coeff>
84 inline Coeff resultant_(
85 const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&);
86
87 template <class Coeff>
resultant_univariate(const CGAL::Polynomial<Coeff> & A,const CGAL::Polynomial<Coeff> & B,CGAL::Integral_domain_without_division_tag)88 inline Coeff resultant_univariate(
89 const CGAL::Polynomial<Coeff>& A,
90 const CGAL::Polynomial<Coeff>& B,
91 CGAL::Integral_domain_without_division_tag){
92 return hybrid_bezout_subresultant(A,B,0);
93 }
94 template <class Coeff>
resultant_univariate(const CGAL::Polynomial<Coeff> & A,const CGAL::Polynomial<Coeff> & B,CGAL::Integral_domain_tag)95 inline Coeff resultant_univariate(
96 const CGAL::Polynomial<Coeff>& A,
97 const CGAL::Polynomial<Coeff>& B, CGAL::Integral_domain_tag){
98 // this seems to help for for large polynomials
99 return prs_resultant_integral_domain(A,B);
100 }
101 template <class Coeff>
resultant_univariate(const CGAL::Polynomial<Coeff> & A,const CGAL::Polynomial<Coeff> & B,CGAL::Unique_factorization_domain_tag)102 inline Coeff resultant_univariate(
103 const CGAL::Polynomial<Coeff>& A,
104 const CGAL::Polynomial<Coeff>& B, CGAL::Unique_factorization_domain_tag){
105 return prs_resultant_ufd(A,B);
106 }
107
108 template <class Coeff>
resultant_univariate(const CGAL::Polynomial<Coeff> & A,const CGAL::Polynomial<Coeff> & B,CGAL::Field_tag)109 inline Coeff resultant_univariate(
110 const CGAL::Polynomial<Coeff>& A,
111 const CGAL::Polynomial<Coeff>& B, CGAL::Field_tag){
112 return prs_resultant_field(A,B);
113 }
114
115 } // namespace internal
116
117 namespace internal{
118
119
120 template <class IC>
121 inline IC
resultant_interpolate(const CGAL::Polynomial<IC> & F,const CGAL::Polynomial<IC> & G)122 resultant_interpolate(
123 const CGAL::Polynomial<IC>& F,
124 const CGAL::Polynomial<IC>& G){
125 CGAL_precondition(CGAL::Polynomial_traits_d<CGAL::Polynomial<IC> >::d == 1);
126 typedef CGAL::Algebraic_structure_traits<IC> AST_IC;
127 typedef typename AST_IC::Algebraic_category Algebraic_category;
128 return internal::resultant_univariate(F,G,Algebraic_category());
129 }
130
131 template <class Coeff_2>
132 inline
resultant_interpolate(const CGAL::Polynomial<CGAL::Polynomial<Coeff_2>> & F,const CGAL::Polynomial<CGAL::Polynomial<Coeff_2>> & G)133 CGAL::Polynomial<Coeff_2> resultant_interpolate(
134 const CGAL::Polynomial<CGAL::Polynomial<Coeff_2> >& F,
135 const CGAL::Polynomial<CGAL::Polynomial<Coeff_2> >& G){
136
137 typedef CGAL::Polynomial<Coeff_2> Coeff_1;
138 typedef CGAL::Polynomial<Coeff_1> POLY;
139 typedef CGAL::Polynomial_traits_d<POLY> PT;
140 typedef typename PT::Innermost_coefficient_type IC;
141
142 CGAL_precondition(PT::d >= 2);
143
144 typename PT::Degree degree;
145 int maxdegree = degree(F,0)*degree(G,PT::d-1) + degree(F,PT::d-1)*degree(G,0);
146
147 typedef std::pair<IC,Coeff_2> Point;
148 std::vector<Point> points; // interpolation points
149
150
151 typename CGAL::Polynomial_traits_d<Coeff_1>::Degree coeff_degree;
152 int i(-maxdegree/2);
153 int deg_f(0);
154 int deg_g(0);
155
156
157 while((int) points.size() <= maxdegree + 1){
158 i++;
159 // timer1.start();
160 Coeff_1 c_i(i);
161 Coeff_1 Fat_i(typename PT::Evaluate()(F,c_i));
162 Coeff_1 Gat_i(typename PT::Evaluate()(G,c_i));
163 // timer1.stop();
164
165 int deg_f_at_i = coeff_degree(Fat_i,0);
166 int deg_g_at_i = coeff_degree(Gat_i,0);
167
168 // std::cout << F << std::endl;
169 // std::cout << Fat_i << std::endl;
170 // std::cout << deg_f_at_i << " vs. " << deg_f << std::endl;
171 if(deg_f_at_i > deg_f ){
172 points.clear();
173 deg_f = deg_f_at_i;
174 CGAL_postcondition(points.size() == 0);
175 }
176
177 if(deg_g_at_i > deg_g){
178 points.clear();
179 deg_g = deg_g_at_i;
180 CGAL_postcondition(points.size() == 0);
181 }
182
183 if(deg_f_at_i == deg_f && deg_g_at_i == deg_g){
184 // timer2.start();
185 Coeff_2 res_at_i = resultant_interpolate(Fat_i, Gat_i);
186 // timer2.stop();
187 points.push_back(Point(IC(i),res_at_i));
188
189 // std::cout << typename Polynomial_traits_d<Coeff_2>::Degree()(res_at_i) << std::endl ;
190 }
191 }
192
193 // timer3.start();
194 CGAL::internal::Interpolator<Coeff_1> interpolator(points.begin(),points.end());
195 Coeff_1 result = interpolator.get_interpolant();
196 // timer3.stop();
197
198 #ifndef CGAL_NDEBUG
199 while((int) points.size() <= maxdegree + 3){
200 i++;
201
202 Coeff_1 c_i(i);
203 Coeff_1 Fat_i(typename PT::Evaluate()(F,c_i));
204 Coeff_1 Gat_i(typename PT::Evaluate()(G,c_i));
205
206 CGAL_assertion(coeff_degree(Fat_i,0) <= deg_f);
207 CGAL_assertion(coeff_degree(Gat_i,0) <= deg_g);
208
209 if(coeff_degree( Fat_i , 0) == deg_f && coeff_degree( Gat_i , 0 ) == deg_g){
210 Coeff_2 res_at_i = resultant_interpolate(Fat_i, Gat_i);
211 points.push_back(Point(IC(i), res_at_i));
212 }
213 }
214 CGAL::internal::Interpolator<Coeff_1>
215 interpolator_(points.begin(),points.end());
216 Coeff_1 result_= interpolator_.get_interpolant();
217
218 // the interpolate polynomial has to be stable !
219 CGAL_assertion(result_ == result);
220 #endif
221 return result;
222 }
223
224 template <class Coeff>
225 inline
resultant_modularize(const CGAL::Polynomial<Coeff> & F,const CGAL::Polynomial<Coeff> & G,CGAL::Tag_false)226 Coeff resultant_modularize(
227 const CGAL::Polynomial<Coeff>& F,
228 const CGAL::Polynomial<Coeff>& G,
229 CGAL::Tag_false){
230 return resultant_interpolate(F,G);
231 }
232
233 template <class Coeff>
234 inline
resultant_modularize(const CGAL::Polynomial<Coeff> & F,const CGAL::Polynomial<Coeff> & G,CGAL::Tag_true)235 Coeff resultant_modularize(
236 const CGAL::Polynomial<Coeff>& F,
237 const CGAL::Polynomial<Coeff>& G,
238 CGAL::Tag_true){
239
240 // Enforce IEEE double precision and to nearest before using modular arithmetic
241 CGAL::Protect_FPU_rounding<true> pfr(CGAL_FE_TONEAREST);
242
243
244 typedef Polynomial_traits_d<CGAL::Polynomial<Coeff> > PT;
245 typedef typename PT::Polynomial_d Polynomial;
246
247 typedef Chinese_remainder_traits<Coeff> CRT;
248 typedef typename CRT::Scalar_type Scalar;
249
250
251 typedef typename CGAL::Modular_traits<Polynomial>::Residue_type MPolynomial;
252 typedef typename CGAL::Modular_traits<Coeff>::Residue_type MCoeff;
253
254 typename CRT::Chinese_remainder chinese_remainder;
255 typename CGAL::Modular_traits<Coeff>::Modular_image_representative inv_map;
256
257
258 typename PT::Degree_vector degree_vector;
259 typename CGAL::Polynomial_traits_d<MPolynomial>::Degree_vector mdegree_vector;
260
261 bool solved = false;
262 int prime_index = 0;
263 int n = 0;
264 Scalar p,q,pq,s,t;
265 Coeff R, R_old;
266
267 // CGAL::Timer timer_evaluate, timer_resultant, timer_cr;
268
269 do{
270 MPolynomial mF, mG;
271 MCoeff mR;
272 //timer_evaluate.start();
273 do{
274 // select a prime number
275 int current_prime = -1;
276 prime_index++;
277 if(prime_index >= 2000){
278 std::cerr<<"primes in the array exhausted"<<std::endl;
279 CGAL_assertion(false);
280 current_prime = internal::get_next_lower_prime(current_prime);
281 } else{
282 current_prime = internal::primes[prime_index];
283 }
284 CGAL::Residue::set_current_prime(current_prime);
285
286 mF = CGAL::modular_image(F);
287 mG = CGAL::modular_image(G);
288
289 }while( degree_vector(F) != mdegree_vector(mF) ||
290 degree_vector(G) != mdegree_vector(mG));
291 //timer_evaluate.stop();
292
293 //timer_resultant.start();
294 n++;
295 mR = resultant_interpolate(mF,mG);
296 //timer_resultant.stop();
297 //timer_cr.start();
298 if(n == 1){
299 // init chinese remainder
300 q = CGAL::Residue::get_current_prime(); // implicit !
301 R = inv_map(mR);
302 }else{
303 // continue chinese remainder
304 p = CGAL::Residue::get_current_prime(); // implicit!
305 R_old = R ;
306 // chinese_remainder(q,Gs ,p,inv_map(mG_),pq,Gs);
307 // cached_extended_euclidean_algorithm(q,p,s,t);
308 internal::Cached_extended_euclidean_algorithm
309 <typename CRT::Scalar_type> ceea;
310 ceea(q,p,s,t);
311 pq =p*q;
312 chinese_remainder(q,p,pq,s,t,R_old,inv_map(mR),R);
313 q=pq;
314 }
315 solved = (R==R_old);
316 //timer_cr.stop();
317 } while(!solved);
318
319 //std::cout << "Time Evaluate : " << timer_evaluate.time() << std::endl;
320 //std::cout << "Time Resultant : " << timer_resultant.time() << std::endl;
321 //std::cout << "Time Chinese R : " << timer_cr.time() << std::endl;
322 // CGAL_postcondition(R == resultant_interpolate(F,G));
323 return R;
324 // return resultant_interpolate(F,G);
325 }
326
327
328 template <class Coeff>
329 inline
resultant_decompose(const CGAL::Polynomial<Coeff> & F,const CGAL::Polynomial<Coeff> & G,CGAL::Tag_false)330 Coeff resultant_decompose(
331 const CGAL::Polynomial<Coeff>& F,
332 const CGAL::Polynomial<Coeff>& G,
333 CGAL::Tag_false){
334 #if CGAL_RESULTANT_USE_MODULAR_ARITHMETIC
335 typedef CGAL::Polynomial<Coeff> Polynomial;
336 typedef typename Modular_traits<Polynomial>::Is_modularizable Is_modularizable;
337 return resultant_modularize(F,G,Is_modularizable());
338 #else
339 return resultant_modularize(F,G,CGAL::Tag_false());
340 #endif
341 }
342
343 template <class Coeff>
344 inline
resultant_decompose(const CGAL::Polynomial<Coeff> & F,const CGAL::Polynomial<Coeff> & G,CGAL::Tag_true)345 Coeff resultant_decompose(
346 const CGAL::Polynomial<Coeff>& F,
347 const CGAL::Polynomial<Coeff>& G,
348 CGAL::Tag_true){
349
350 typedef Polynomial<Coeff> POLY;
351 typedef typename Fraction_traits<POLY>::Numerator_type Numerator;
352 typedef typename Fraction_traits<POLY>::Denominator_type Denominator;
353 typename Fraction_traits<POLY>::Decompose decompose;
354 typedef typename Numerator::NT RES;
355
356 Denominator a, b;
357 // F.simplify_coefficients(); not const
358 // G.simplify_coefficients(); not const
359 Numerator F0; decompose(F,F0,a);
360 Numerator G0; decompose(G,G0,b);
361 Denominator c = CGAL::ipower(a, G.degree()) * CGAL::ipower(b, F.degree());
362
363 RES res0 = CGAL::internal::resultant_(F0, G0);
364 typename Fraction_traits<Coeff>::Compose comp_frac;
365 Coeff res = comp_frac(res0, c);
366 typename Algebraic_structure_traits<Coeff>::Simplify simplify;
367 simplify(res);
368 return res;
369 }
370
371
372 template <class Coeff>
373 inline
resultant_(const CGAL::Polynomial<Coeff> & F,const CGAL::Polynomial<Coeff> & G)374 Coeff resultant_(
375 const CGAL::Polynomial<Coeff>& F,
376 const CGAL::Polynomial<Coeff>& G){
377 #if CGAL_RESULTANT_USE_DECOMPOSE
378 typedef CGAL::Fraction_traits<Polynomial<Coeff > > FT;
379 typedef typename FT::Is_fraction Is_fraction;
380 return resultant_decompose(F,G,Is_fraction());
381 #else
382 return resultant_decompose(F,G,CGAL::Tag_false());
383 #endif
384 }
385
386
387
388 template <class Coeff>
389 inline
resultant(const CGAL::Polynomial<Coeff> & F_,const CGAL::Polynomial<Coeff> & G_)390 Coeff resultant(
391 const CGAL::Polynomial<Coeff>& F_,
392 const CGAL::Polynomial<Coeff>& G_){
393 // make the variable to be elimnated the innermost one.
394 typedef CGAL::Polynomial_traits_d<CGAL::Polynomial<Coeff> > PT;
395 CGAL::Polynomial<Coeff> F = typename PT::Move()(F_, PT::d-1, 0);
396 CGAL::Polynomial<Coeff> G = typename PT::Move()(G_, PT::d-1, 0);
397 return internal::resultant_(F,G);
398 }
399
400 } // namespace internal
401 } //namespace CGAL
402
403
404
405 #endif // CGAL_POLYNOMIAL_RESULTANT_H
406