1 // Copyright (c) 2008 Max-Planck-Institute Saarbruecken (Germany)
2 //
3 // This file is part of CGAL (www.cgal.org)
4 //
5 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Polynomial/include/CGAL/Polynomial/bezout_matrix.h $
6 // $Id: bezout_matrix.h e9d41d7 2020-04-21T10:03:00+02:00 Maxime Gimeno
7 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
8 //
9 //
10 // Author(s) : Michael Hemmer
11 //
12 // ============================================================================
13
14 // TODO: The comments are all original EXACUS comments and aren't adapted. So
15 // they may be wrong now.
16
17 #ifndef CGAL_POLYNOMIAL_BEZOUT_MATRIX_H
18 #define CGAL_POLYNOMIAL_BEZOUT_MATRIX_H
19
20 #include <algorithm>
21
22 #include <CGAL/basic.h>
23 #include <CGAL/Polynomial_traits_d.h>
24 #include <CGAL/Polynomial/determinant.h>
25 #include <CGAL/use.h>
26
27 #include <vector>
28
29 namespace CGAL {
30
31 namespace internal {
32
33 /*! \ingroup CGAL_resultant_matrix
34 * \brief construct hybrid Bezout matrix of two polynomials
35 *
36 * If \c sub=0 , this function returns the hybrid Bezout matrix
37 * of \c f and \c g.
38 * The hybrid Bezout matrix of two polynomials \e f and \e g
39 * (seen as polynomials in one variable) is a
40 * square matrix of size max(deg(<I>f</I>), deg(<I>g</I>)) whose entries
41 * are expressions in the polynomials' coefficients.
42 * Its determinant is the resultant of \e f and \e g, maybe up to sign.
43 * The function computes the same matrix as the Maple command
44 * <I>BezoutMatrix</I>.
45 *
46 * If \c sub>0 , this function returns the matrix obtained by chopping
47 * off the \c sub topmost rows and the \c sub rightmost columns.
48 * Its determinant is the <I>sub</I>-th (scalar) subresultant
49 * of \e f and \e g, maybe up to sign.
50 *
51 * If specified, \c sub must be less than the degrees of both \e f and \e g.
52 *
53 * See also \c CGAL::hybrid_bezout_subresultant() and \c CGAL::sylvester_matrix() .
54 *
55 * A formal definition of the hybrid Bezout matrix and a proof for the
56 * subresultant property can be found in:
57 *
58 * Gema M.Diaz-Toca, Laureano Gonzalez-Vega: Various New Expressions for
59 * Subresultants and Their Applications. AAECC 15, 233-266 (2004)
60 *
61 */
62 template <typename Polynomial_traits_d>
63 typename internal::Simple_matrix< typename Polynomial_traits_d::Coefficient_type >
64 hybrid_bezout_matrix(typename Polynomial_traits_d::Polynomial_d f,
65 typename Polynomial_traits_d::Polynomial_d g,
66 int sub = 0)
67 {
68
69 typedef typename Polynomial_traits_d::Polynomial_d Polynomial;
70 typedef typename Polynomial_traits_d::Coefficient_type NT;
71 typename Polynomial_traits_d::Degree degree;
72 typename CGAL::Algebraic_structure_traits<Polynomial>::Is_zero is_zero;
73 typename Polynomial_traits_d::Get_coefficient coeff;
74
75 typedef typename internal::Simple_matrix<NT> Matrix;
76
77 int n = degree(f);
78 int m = degree(g);
79 CGAL_precondition((n >= 0) && !is_zero(f));
80 CGAL_precondition((m >= 0) && !is_zero(g));
81 CGAL_precondition(n > sub || sub == 0);
82 CGAL_precondition(m > sub || sub == 0);
83
84 int i, j, k, l;
85 NT s;
86
87 if (m > n) {
88 std::swap(f, g);
89 std::swap(m, n);
90 }
91
92 Matrix B(n-sub);
93
94 for (i = 1+sub; i <= m; i++) {
95 for (j = 1; j <= n-sub; j++) {
96 s = NT(0);
97 for (k = 0; k <= i-1; k++) {
98 l = n+i-j-k;
99 if ((l <= n) and (l >= n-(m-i))) {
100 s += coeff(f,l) * coeff(g,k);
101 }
102 }
103 for (k = 0; k <= n-(m-i+1); k++) {
104 l = n+i-j-k;
105 if ((l <= m) and (l >= i)) {
106 s -= coeff(f,k) * coeff(g,l);
107 }
108 }
109 B[i-sub-1][j-1] = s;
110 }
111 }
112 for (i = (std::max)(m+1, 1+sub); i <= n; i++) {
113 for (j = i-m; j <= (std::min)(i, n-sub); j++) {
114 B[i-sub-1][j-1] = coeff(g,i-j);
115 }
116 }
117
118 return B; // g++ 3.1+ has NRVO, so this should not be too expensive
119 }
120
121 /*! \ingroup CGAL_resultant_matrix
122 * \brief construct the symmetric Bezout matrix of two polynomials
123 *
124 * This function returns the (symmetric) Bezout matrix of \c f and \c g.
125 * The Bezout matrix of two polynomials \e f and \e g
126 * (seen as polynomials in one variable) is a
127 * square matrix of size max(deg(<I>f</I>), deg(<I>g</I>)) whose entries
128 * are expressions in the polynomials' coefficients.
129 * Its determinant is the resultant of \e f and \e g, maybe up to sign.
130 *
131 */
132 template <typename Polynomial_traits_d>
133 typename internal::Simple_matrix<typename Polynomial_traits_d::Coefficient_type>
134 symmetric_bezout_matrix
135 (typename Polynomial_traits_d::Polynomial_d f,
136 typename Polynomial_traits_d::Polynomial_d g,
137 int sub = 0)
138 {
139
140
141
142 // Note: The algorithm is taken from:
143 // Chionh, Zhang, Goldman: Fast Computation of the Bezout and Dixon Resultant
144 // Matrices. J.Symbolic Computation 33, 13-29 (2002)
145
146 typedef typename Polynomial_traits_d::Polynomial_d Polynomial;
147 typedef typename Polynomial_traits_d::Coefficient_type NT;
148 typename Polynomial_traits_d::Degree degree;
149 CGAL_assertion_code(typename CGAL::Algebraic_structure_traits<Polynomial>::Is_zero is_zero;)
150 CGAL_USE_TYPE(Polynomial);
151 typename Polynomial_traits_d::Get_coefficient coeff;
152
153 typedef typename internal::Simple_matrix<NT> Matrix;
154
155 int n = degree(f);
156 int m = degree(g);
157 CGAL_precondition((n >= 0) && !is_zero(f));
158 CGAL_precondition((m >= 0) && !is_zero(g));
159
160 int i,j,stop;
161
162 NT sum1,sum2;
163
164 if (m > n) {
165 std::swap(f, g);
166 std::swap(m, n);
167 }
168
169 CGAL_precondition((sub>=0) && sub < n);
170
171 int d = n - sub;
172
173 Matrix B(d);
174
175 // 1st step: Initialisation
176 for(i=0;i<d;i++) {
177 for(j=i;j<d;j++) {
178 sum1 = ((j+sub)+1>m) ? NT(0) : -coeff(f,i+sub)*coeff(g,(j+sub)+1);
179 sum2 = ((i+sub)>m) ? NT(0) : coeff(g,i+sub)*coeff(f,(j+sub)+1);
180 B[i][j]=sum1+sum2;
181 }
182 }
183
184 // 2nd Step: Recursion adding
185
186 // First, set up the first line correctly
187 for(i=0;i<d-1;i++) {
188 stop = (sub<d-1-i) ? sub : d-i-1;
189 for(j=1;j<=stop;j++) {
190 sum1 = ((i+sub+j)+1>m) ? NT(0)
191 : -coeff(f,sub-j)*coeff(g,(i+sub+j)+1);
192 sum2 = ((sub-j)>m) ? NT(0)
193 : coeff(g,sub-j)*coeff(f,(i+sub+j)+1);
194
195 B[0][i]+=sum1+sum2;
196 }
197 }
198 // Now, compute the rest
199 for(i=1;i<d-1;i++) {
200 for(j=i;j<d-1;j++) {
201 B[i][j]+=B[i-1][j+1];
202 }
203 }
204
205
206 //3rd Step: Exploit symmetry
207 for(i=1;i<d;i++) {
208 for(j=0;j<i;j++) {
209 B[i][j]=B[j][i];
210 }
211 }
212
213 return B;
214 }
215
216
217
218 /*! \ingroup CGAL_resultant_matrix
219 * \brief compute (sub)resultant as Bezout matrix determinant
220 *
221 * This function returns the determinant of the matrix returned
222 * by <TT>hybrid_bezout_matrix(f, g, sub)</TT> which is the
223 * resultant of \c f and \c g, maybe up to sign;
224 * or rather the <I>sub</I>-th (scalar) subresultant, if a non-zero third
225 * argument is given.
226 *
227 * If specified, \c sub must be less than the degrees of both \e f and \e g.
228 *
229 * This function might be faster than \c CGAL::Polynomial<..>::resultant() ,
230 * which computes the resultant from a subresultant remainder sequence.
231 * See also \c CGAL::sylvester_subresultant().
232 */
233 template <class Polynomial_traits_d>
234 typename Polynomial_traits_d::Coefficient_type hybrid_bezout_subresultant(
235 typename Polynomial_traits_d::Polynomial_d f,
236 typename Polynomial_traits_d::Polynomial_d g,
237 int sub = 0
238 ) {
239
240 typedef typename Polynomial_traits_d::Polynomial_d Polynomial;
241 typedef typename Polynomial_traits_d::Coefficient_type NT;
242 typename Polynomial_traits_d::Degree degree;
243 typename CGAL::Algebraic_structure_traits<Polynomial>::Is_zero is_zero;
244
245 typedef internal::Simple_matrix<NT> Matrix;
246
247 CGAL_precondition((degree(f) >= 0));
248 CGAL_precondition((degree(g) >= 0));
249
250 if (is_zero(f) || is_zero(g)) return NT(0);
251
252 Matrix S = hybrid_bezout_matrix<Polynomial_traits_d>(f, g, sub);
253 CGAL_assertion(S.row_dimension() == S.column_dimension());
254 if (S.row_dimension() == 0) {
255 return NT(1);
256 } else {
257 return internal::determinant(S);
258 }
259 }
260
261 // Transforms the minors of the symmetric bezout matrix into the subresultant.
262 // Needs the appropriate power of the leading coedfficient of f and the
263 // degrees of f and g
264 template<class InputIterator,class OutputIterator,class NT>
symmetric_minors_to_subresultants(InputIterator in,OutputIterator out,NT divisor,int n,int m,bool swapped)265 void symmetric_minors_to_subresultants(InputIterator in,
266 OutputIterator out,
267 NT divisor,
268 int n,
269 int m,
270 bool swapped) {
271
272 typename CGAL::Algebraic_structure_traits<NT>::Integral_division idiv;
273
274 for(int i=0;i<m;i++) {
275 bool negate = ((n-m+i+1) & 2)>>1; // (n-m+i+1)==2 or 3 mod 4
276 negate=negate ^ (swapped & ((n-m+i+1)*(i+1)));
277 //...XOR (swapped AND (n-m+i+1)* (i+1) is odd)
278
279 *out = idiv(*in, negate ? -divisor : divisor);
280 in++;
281 out++;
282 }
283 }
284
285
286 /*! \ingroup CGAL_resultant_matrix
287 * \brief compute the principal subresultant coefficients as minors
288 * of the symmetric Bezout matrix.
289 *
290 * Returns the sequence sres<sub>0</sub>,..,sres<sub>m</sub>, where
291 * sres<sub>i</sub> denotes the ith principal subresultant coefficient
292 *
293 * The function uses an extension of the Berkowitz method to compute the
294 * determinant
295 * See also \c CGAL::minors_berkowitz
296 */
297 template<class Polynomial_traits_d,class OutputIterator>
symmetric_bezout_subresultants(typename Polynomial_traits_d::Polynomial_d f,typename Polynomial_traits_d::Polynomial_d g,OutputIterator sres)298 OutputIterator symmetric_bezout_subresultants(
299 typename Polynomial_traits_d::Polynomial_d f,
300 typename Polynomial_traits_d::Polynomial_d g,
301 OutputIterator sres)
302 {
303
304 typedef typename Polynomial_traits_d::Polynomial_d Polynomial;
305 typedef typename Polynomial_traits_d::Coefficient_type NT;
306 typename Polynomial_traits_d::Degree degree;
307 typename CGAL::Algebraic_structure_traits<Polynomial>::Is_zero is_zero;
308 typename Polynomial_traits_d::Leading_coefficient lcoeff;
309
310 typedef typename internal::Simple_matrix<NT> Matrix;
311
312 int n = degree(f);
313 int m = degree(g);
314
315 bool swapped=false;
316
317 if(n < m) {
318 std::swap(f,g);
319 std::swap(n,m);
320 swapped=true;
321
322 }
323
324 Matrix B = symmetric_bezout_matrix<Polynomial_traits_d>(f,g);
325
326 // Compute a_0^{n-m}
327
328 NT divisor=ipower(lcoeff(f),n-m);
329
330 std::vector<NT> minors;
331 minors_berkowitz(B,std::back_inserter(minors),n,m);
332 CGAL::internal::symmetric_minors_to_subresultants(minors.begin(),sres,
333 divisor,n,m,swapped);
334
335 return sres;
336 }
337
338 /*
339 * Return a modified version of the hybrid bezout matrix such that the minors
340 * from the last k rows and columns give the subresultants
341 */
342 template<class Polynomial_traits_d>
343 typename internal::Simple_matrix<typename Polynomial_traits_d::Coefficient_type>
modified_hybrid_bezout_matrix(typename Polynomial_traits_d::Polynomial_d f,typename Polynomial_traits_d::Polynomial_d g)344 modified_hybrid_bezout_matrix
345 (typename Polynomial_traits_d::Polynomial_d f,
346 typename Polynomial_traits_d::Polynomial_d g) {
347
348 typedef typename Polynomial_traits_d::Coefficient_type NT;
349
350 typedef typename internal::Simple_matrix<NT> Matrix;
351
352 typename Polynomial_traits_d::Degree degree;
353
354 int n = degree(f);
355 int m = degree(g);
356
357 int i,j;
358
359 bool negate, swapped=false;
360
361 if(n < m) {
362 std::swap(f,g); //(*)
363 std::swap(n,m);
364 swapped=true;
365 }
366
367 Matrix B = CGAL::internal::hybrid_bezout_matrix<Polynomial_traits_d>(f,g);
368
369
370 // swap columns
371 i=0;
372 while(i<n-i-1) {
373 B.swap_columns(i,n-i-1); // (**)
374 i++;
375 }
376 for(i=0;i<n;i++) {
377 negate=(n-i-1) & 1; // Negate every second column because of (**)
378 negate=negate ^ (swapped & (n-m+1)); // XOR negate everything because of(*)
379 if(negate) {
380 for(j=0;j<n;j++) {
381 B[j][i] *= -1;
382 }
383 }
384 }
385 return B;
386 }
387
388 /*! \ingroup CGAL_resultant_matrix
389 * \brief compute the principal subresultant coefficients as minors
390 * of the hybrid Bezout matrix.
391 *
392 * Returns the sequence sres<sub>0</sub>,...,sres<sub>m</sub>$, where
393 * sres<sub>i</sub> denotes the ith principal subresultant coefficient
394 *
395 * The function uses an extension of the Berkowitz method to compute the
396 * determinant
397 * See also \c CGAL::minors_berkowitz
398 */
399 template<class Polynomial_traits_d,class OutputIterator>
hybrid_bezout_subresultants(typename Polynomial_traits_d::Polynomial_d f,typename Polynomial_traits_d::Polynomial_d g,OutputIterator sres)400 OutputIterator hybrid_bezout_subresultants(
401 typename Polynomial_traits_d::Polynomial_d f,
402 typename Polynomial_traits_d::Polynomial_d g,
403 OutputIterator sres)
404 {
405
406 typedef typename Polynomial_traits_d::Coefficient_type NT;
407 typename Polynomial_traits_d::Degree degree;
408
409 typedef typename internal::Simple_matrix<NT> Matrix;
410
411 int n = degree(f);
412 int m = degree(g);
413
414 Matrix B = CGAL::internal::modified_hybrid_bezout_matrix<Polynomial_traits_d>
415 (f,g);
416
417 if(n<m) {
418 std::swap(n,m);
419 }
420
421 return minors_berkowitz(B,sres,n,m);
422 }
423
424
425 // Swap entry A_ij with A_(n-i)(n-j) for square matrix A of dimension n
426 template<class NT>
swap_entries(typename internal::Simple_matrix<NT> & A)427 void swap_entries(typename internal::Simple_matrix<NT> & A) {
428 CGAL_precondition(A.row_dimension()==A.column_dimension());
429 int n = A.row_dimension();
430 int i=0;
431 while(i<n-i-1) {
432 A.swap_rows(i,n-i-1);
433 A.swap_columns(i,n-i-1);
434 i++;
435 }
436 }
437
438 // Produce S-matrix with the given matrix and integers.
439 template<class NT,class InputIterator>
s_matrix(const typename internal::Simple_matrix<NT> & B,InputIterator num,int size)440 typename internal::Simple_matrix<NT> s_matrix(
441 const typename internal::Simple_matrix<NT>& B,
442 InputIterator num,int size)
443 {
444 typename internal::Simple_matrix<NT> S(size);
445 CGAL_precondition_code(int n = B.row_dimension();)
446 CGAL_precondition(n==(int)B.column_dimension());
447 int curr_num;
448 bool negate;
449
450 for(int i=0;i<size;i++) {
451 curr_num=(*num);
452 num++;
453 negate = curr_num<0;
454 if(curr_num<0) {
455 curr_num=-curr_num;
456 }
457 for(int j=0;j<size;j++) {
458
459 S[j][i]=negate ? -B[j][curr_num-1] : B[j][curr_num-1];
460
461 }
462 }
463 return S;
464 }
465
466 // Produces the integer sequence for the S-matrix, where c is the first entry
467 // of the sequence, s the number of desired diagonals and n the dimension
468 // of the base matrix
469 template<class OutputIterator>
s_matrix_integer_sequence(OutputIterator it,int c,int s,int n)470 OutputIterator s_matrix_integer_sequence(OutputIterator it,
471 int c,int s,int n) {
472 CGAL_precondition(0<s);
473 CGAL_precondition(s<=n);
474 // c is interpreted modulo s wrt to the representants {1,..,s}
475 c=c%s;
476 if(c==0) {
477 c=s;
478 }
479
480 int i, p=0, q=c;
481 while(q<=n) {
482 *it = q;
483 it++;
484 for(i=p+1;i<q;i++) {
485 *it = -i;
486 it++;
487 }
488 p = q;
489 q = q+s;
490 }
491 return it;
492 }
493
494 /*! \ingroup CGAL_resultant_matrix
495 * \brief computes the coefficients of the polynomial subresultant sequence
496 *
497 * Returns an upper triangular matrix <I>A</I> such that A<sub>i,j</sub> is
498 * the coefficient of <I>x<sup>j-1</sup></I> in the <I>i</I>th polynomial
499 * subresultant. In particular, the main diagonal contains the scalar
500 * subresultants.
501 *
502 * If \c d > 0 is specified, only the first \c d diagonals of <I>A</I> are
503 * computed. In particular, setting \c d to one yields exactly the same
504 * result as applying \c hybrid_subresultants or \c symmetric_subresultants
505 * (except the different output format).
506 *
507 * These coefficients are computed as special minors of the hybrid Bezout matrix.
508 * See also \c CGAL::minors_berkowitz
509 */
510 template<typename Polynomial_traits_d>
511 typename internal::Simple_matrix<typename Polynomial_traits_d::Coefficient_type>
512 polynomial_subresultant_matrix(typename Polynomial_traits_d::Polynomial_d f,
513 typename Polynomial_traits_d::Polynomial_d g,
514 int d=0) {
515
516 typedef typename Polynomial_traits_d::Coefficient_type NT;
517 typename Polynomial_traits_d::Degree degree;
518 typename Polynomial_traits_d::Leading_coefficient lcoeff;
519
520 int n = degree(f);
521 int m = degree(g);
522
523 CGAL_precondition(n>=0);
524 CGAL_precondition(m>=0);
525 CGAL_precondition(d>=0);
526
527 typedef internal::Simple_matrix<NT> Matrix;
528
529 bool swapped=false;
530
531 if(n < m) {
532 std::swap(f,g);
533 std::swap(n,m);
534 swapped=true;
535 }
536
537 if(d==0) {
538 d=m;
539 };
540
541
542 Matrix B = CGAL::internal::symmetric_bezout_matrix<Polynomial_traits_d>(f,g);
543
544 // For easier notation, we swap all entries:
545 internal::swap_entries<NT>(B);
546
547 // Compute the S-matrices and collect the minors
548 std::vector<Matrix> s_mat(m);
549 std::vector<std::vector<NT> > coeffs(d);
550 for(int i = 1; i<=d;i++) {
551 std::vector<int> intseq;
552 internal::s_matrix_integer_sequence(std::back_inserter(intseq),i,d,n);
553
554 Matrix S = internal::s_matrix<NT>(B,intseq.begin(),(int)intseq.size());
555 internal::swap_entries<NT>(S);
556 //std::cout << S << std::endl;
557 int Sdim = S.row_dimension();
558 int number_of_minors=(Sdim < m) ? Sdim : Sdim;
559
560 internal::minors_berkowitz(S,std::back_inserter(coeffs[i-1]),
561 Sdim,number_of_minors);
562
563 }
564
565 // Now, rearrange the minors in the matrix
566
567 Matrix Ret(m,m,NT(0));
568
569 for(int i = 0; i < d; i++) {
570 for(int j = 0;j < m-i ; j++) {
571 int s_index=(n-m+j+i+1)%d;
572 if(s_index==0) {
573 s_index=d;
574 }
575 s_index--;
576 Ret[j][j+i]=coeffs[s_index][n-m+j];
577
578 }
579 }
580
581 typename CGAL::Algebraic_structure_traits<NT>::Integral_division idiv;
582
583 NT divisor = ipower(lcoeff(f),n-m);
584
585 int bit_mask = swapped ? 1 : 0;
586 // Divide through the divisor and set the correct sign
587 for(int i=0;i<m;i++) {
588 for(int j = i;j<m;j++) {
589 int negate = ((n-m+i+1) & 2)>>1; // (n-m+i+1)==2 or 3 mod 4
590 negate^=(bit_mask & ((n-m+i+1)*(i+1)));
591 //...XOR (swapped AND (n-m+i+1)* (i+1) is odd)
592 Ret[i][j] = idiv(Ret[i][j], negate>0 ? -divisor : divisor);
593 }
594 }
595
596 return Ret;
597 }
598
599 }
600
601 } //namespace CGAL
602
603
604
605 #endif // CGAL_POLYNOMIAL_BEZOUT_MATRIX_H
606 // EOF
607