1 /*****************************************************************************\
2 * Computer Algebra System SINGULAR
3 \*****************************************************************************/
4 /** @file lineareAlgebra.cc
5 *
6 * This file implements basic linear algebra functionality.
7 *
8 * For more general information, see the documentation in
9 * lineareAlgebra.h.
10 *
11 * @author Frank Seelisch
12 *
13 *
14 **/
15 /*****************************************************************************/
16
17 // include header files
18
19
20
21 #include "kernel/mod2.h"
22
23 #include "coeffs/coeffs.h"
24 #include "coeffs/numbers.h"
25
26 #include "coeffs/mpr_complex.h"
27 #include "polys/monomials/p_polys.h"
28 #include "polys/simpleideals.h"
29 #include "polys/matpol.h"
30
31 // #include "kernel/polys.h"
32
33 #include "kernel/structs.h"
34 #include "kernel/ideals.h"
35 #include "kernel/linear_algebra/linearAlgebra.h"
36
37 /**
38 * The returned score is based on the implementation of 'nSize' for
39 * numbers (, see numbers.h): nSize(n) provides a measure for the
40 * complexity of n. Thus, less complex pivot elements will be
41 * prefered, and get therefore a smaller pivot score. Consequently,
42 * we simply return the value of nSize.
43 * An exception to this rule are the ground fields R, long R, and
44 * long C: In these, the result of nSize relates to |n|. And since
45 * a larger modulus of the pivot element ensures a numerically more
46 * stable Gauss elimination, we would like to have a smaller score
47 * for larger values of |n|. Therefore, in R, long R, and long C,
48 * the negative of nSize will be returned.
49 **/
pivotScore(number n,const ring r)50 int pivotScore(number n, const ring r)
51 {
52 int s = n_Size(n, r->cf);
53 if (rField_is_long_C(r) ||
54 rField_is_long_R(r) ||
55 rField_is_R(r))
56 return -s;
57 else
58 return s;
59 }
60
61 /**
62 * This code computes a score for each non-zero matrix entry in
63 * aMat[r1..r2, c1..c2]. If all entries are zero, false is returned,
64 * otherwise true. In the latter case, the minimum of all scores
65 * is sought, and the row and column indices of the corresponding
66 * matrix entry are stored in bestR and bestC.
67 **/
pivot(const matrix aMat,const int r1,const int r2,const int c1,const int c2,int * bestR,int * bestC,const ring R)68 bool pivot(const matrix aMat, const int r1, const int r2, const int c1,
69 const int c2, int* bestR, int* bestC, const ring R)
70 {
71 int bestScore; int score; bool foundBestScore = false; poly matEntry;
72
73 for (int c = c1; c <= c2; c++)
74 {
75 for (int r = r1; r <= r2; r++)
76 {
77 matEntry = MATELEM(aMat, r, c);
78 if (matEntry != NULL)
79 {
80 score = pivotScore(pGetCoeff(matEntry), R);
81 if ((!foundBestScore) || (score < bestScore))
82 {
83 bestScore = score;
84 *bestR = r;
85 *bestC = c;
86 }
87 foundBestScore = true;
88 }
89 }
90 }
91
92 return foundBestScore;
93 }
94
unitMatrix(const int n,matrix & unitMat,const ring R)95 bool unitMatrix(const int n, matrix &unitMat, const ring R)
96 {
97 if (n < 1) return false;
98 unitMat = mpNew(n, n);
99 for (int r = 1; r <= n; r++) MATELEM(unitMat, r, r) = p_One(R);
100 return true;
101 }
102
luDecomp(const matrix aMat,matrix & pMat,matrix & lMat,matrix & uMat,const ring R)103 void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat,
104 const ring R)
105 {
106 int rr = aMat->rows();
107 int cc = aMat->cols();
108 pMat = mpNew(rr, rr);
109 uMat = mp_Copy(aMat,R); /* copy aMat into uMat: */
110
111 /* we use an int array to store all row permutations;
112 note that we only make use of the entries [1..rr] */
113 int* permut = new int[rr + 1];
114 for (int i = 1; i <= rr; i++) permut[i] = i;
115
116 /* fill lMat with the (rr x rr) unit matrix */
117 unitMatrix(rr, lMat,R);
118
119 int bestR; int bestC; int intSwap; poly pSwap; int cOffset = 0;
120 for (int r = 1; r < rr; r++)
121 {
122 if (r > cc) break;
123 while ((r + cOffset <= cc) &&
124 (!pivot(uMat, r, rr, r + cOffset, r + cOffset, &bestR, &bestC, R)))
125 cOffset++;
126 if (r + cOffset <= cc)
127 {
128 /* swap rows with indices r and bestR in permut */
129 intSwap = permut[r];
130 permut[r] = permut[bestR];
131 permut[bestR] = intSwap;
132
133 /* swap rows with indices r and bestR in uMat;
134 it is sufficient to do this for columns >= r + cOffset*/
135 for (int c = r + cOffset; c <= cc; c++)
136 {
137 pSwap = MATELEM(uMat, r, c);
138 MATELEM(uMat, r, c) = MATELEM(uMat, bestR, c);
139 MATELEM(uMat, bestR, c) = pSwap;
140 }
141
142 /* swap rows with indices r and bestR in lMat;
143 we must do this only for columns < r */
144 for (int c = 1; c < r; c++)
145 {
146 pSwap = MATELEM(lMat, r, c);
147 MATELEM(lMat, r, c) = MATELEM(lMat, bestR, c);
148 MATELEM(lMat, bestR, c) = pSwap;
149 }
150
151 /* perform next Gauss elimination step, i.e., below the
152 row with index r;
153 we need to adjust lMat and uMat;
154 we are certain that the matrix entry at [r, r + cOffset]
155 is non-zero: */
156 number pivotElement = pGetCoeff(MATELEM(uMat, r, r + cOffset));
157 poly p;
158 for (int rGauss = r + 1; rGauss <= rr; rGauss++)
159 {
160 p = MATELEM(uMat, rGauss, r + cOffset);
161 if (p != NULL)
162 {
163 number n = n_Div(pGetCoeff(p), pivotElement, R->cf);
164 n_Normalize(n,R->cf);
165
166 /* filling lMat;
167 old entry was zero, so no need to call pDelete(.) */
168 MATELEM(lMat, rGauss, r) = p_NSet(n_Copy(n,R->cf),R);
169
170 /* adjusting uMat: */
171 MATELEM(uMat, rGauss, r + cOffset) = NULL; p_Delete(&p,R);
172 n = n_InpNeg(n,R->cf);
173 for (int cGauss = r + cOffset + 1; cGauss <= cc; cGauss++)
174 {
175 MATELEM(uMat, rGauss, cGauss)
176 = p_Add_q(MATELEM(uMat, rGauss, cGauss),
177 pp_Mult_nn(MATELEM(uMat, r, cGauss), n, R), R);
178 p_Normalize(MATELEM(uMat, rGauss, cGauss),R);
179 }
180
181 n_Delete(&n,R->cf); /* clean up n */
182 }
183 }
184 }
185 }
186
187 /* building the permutation matrix from 'permut' */
188 for (int r = 1; r <= rr; r++)
189 MATELEM(pMat, r, permut[r]) = p_One(R);
190 delete[] permut;
191
192 return;
193 }
194
195 /**
196 * This code first computes the LU-decomposition of aMat,
197 * and then calls the method for inverting a matrix which
198 * is given by its LU-decomposition.
199 **/
luInverse(const matrix aMat,matrix & iMat,const ring R)200 bool luInverse(const matrix aMat, matrix &iMat, const ring R)
201 { /* aMat is guaranteed to be an (n x n)-matrix */
202 matrix pMat;
203 matrix lMat;
204 matrix uMat;
205 luDecomp(aMat, pMat, lMat, uMat, R);
206 bool result = luInverseFromLUDecomp(pMat, lMat, uMat, iMat, R);
207
208 /* clean-up */
209 id_Delete((ideal*)&pMat,R);
210 id_Delete((ideal*)&lMat,R);
211 id_Delete((ideal*)&uMat,R);
212
213 return result;
214 }
215
216 /* Assumes that aMat is already in row echelon form */
rankFromRowEchelonForm(const matrix aMat)217 int rankFromRowEchelonForm(const matrix aMat)
218 {
219 int rank = 0;
220 int rr = aMat->rows(); int cc = aMat->cols();
221 int r = 1; int c = 1;
222 while ((r <= rr) && (c <= cc))
223 {
224 if (MATELEM(aMat, r, c) == NULL) c++;
225 else { rank++; r++; }
226 }
227 return rank;
228 }
229
luRank(const matrix aMat,const bool isRowEchelon,const ring R)230 int luRank(const matrix aMat, const bool isRowEchelon, const ring R)
231 {
232 if (isRowEchelon) return rankFromRowEchelonForm(aMat);
233 else
234 { /* compute the LU-decomposition and read off the rank from
235 the upper triangular matrix of that decomposition */
236 matrix pMat;
237 matrix lMat;
238 matrix uMat;
239 luDecomp(aMat, pMat, lMat, uMat,R);
240 int result = rankFromRowEchelonForm(uMat);
241
242 /* clean-up */
243 id_Delete((ideal*)&pMat,R);
244 id_Delete((ideal*)&lMat,R);
245 id_Delete((ideal*)&uMat,R);
246
247 return result;
248 }
249 }
250
upperRightTriangleInverse(const matrix uMat,matrix & iMat,bool diagonalIsOne,const ring R)251 bool upperRightTriangleInverse(const matrix uMat, matrix &iMat,
252 bool diagonalIsOne, const ring R)
253 {
254 int d = uMat->rows();
255 poly p; poly q;
256
257 /* check whether uMat is invertible */
258 bool invertible = diagonalIsOne;
259 if (!invertible)
260 {
261 invertible = true;
262 for (int r = 1; r <= d; r++)
263 {
264 if (MATELEM(uMat, r, r) == NULL)
265 {
266 invertible = false;
267 break;
268 }
269 }
270 }
271
272 if (invertible)
273 {
274 iMat = mpNew(d, d);
275 for (int r = d; r >= 1; r--)
276 {
277 if (diagonalIsOne)
278 MATELEM(iMat, r, r) = p_One(R);
279 else
280 MATELEM(iMat, r, r) = p_NSet(n_Invers(p_GetCoeff(MATELEM(uMat, r, r),R),R->cf),R);
281 for (int c = r + 1; c <= d; c++)
282 {
283 p = NULL;
284 for (int k = r + 1; k <= c; k++)
285 {
286 q = pp_Mult_qq(MATELEM(uMat, r, k), MATELEM(iMat, k, c),R);
287 p = p_Add_q(p, q,R);
288 }
289 p = p_Neg(p,R);
290 p = p_Mult_q(p, p_Copy(MATELEM(iMat, r, r),R),R);
291 p_Normalize(p,R);
292 MATELEM(iMat, r, c) = p;
293 }
294 }
295 }
296
297 return invertible;
298 }
299
lowerLeftTriangleInverse(const matrix lMat,matrix & iMat,bool diagonalIsOne)300 bool lowerLeftTriangleInverse(const matrix lMat, matrix &iMat,
301 bool diagonalIsOne)
302 {
303 int d = lMat->rows(); poly p; poly q;
304
305 /* check whether lMat is invertible */
306 bool invertible = diagonalIsOne;
307 if (!invertible)
308 {
309 invertible = true;
310 for (int r = 1; r <= d; r++)
311 {
312 if (MATELEM(lMat, r, r) == NULL)
313 {
314 invertible = false;
315 break;
316 }
317 }
318 }
319
320 if (invertible)
321 {
322 iMat = mpNew(d, d);
323 for (int c = d; c >= 1; c--)
324 {
325 if (diagonalIsOne)
326 MATELEM(iMat, c, c) = pOne();
327 else
328 MATELEM(iMat, c, c) = pNSet(nInvers(pGetCoeff(MATELEM(lMat, c, c))));
329 for (int r = c + 1; r <= d; r++)
330 {
331 p = NULL;
332 for (int k = c; k <= r - 1; k++)
333 {
334 q = ppMult_qq(MATELEM(lMat, r, k), MATELEM(iMat, k, c));
335 p = pAdd(p, q);
336 }
337 p = pNeg(p);
338 p = pMult(p, pCopy(MATELEM(iMat, c, c)));
339 pNormalize(p);
340 MATELEM(iMat, r, c) = p;
341 }
342 }
343 }
344
345 return invertible;
346 }
347
348 /**
349 * This code computes the inverse by inverting lMat and uMat, and
350 * then performing two matrix multiplications.
351 **/
luInverseFromLUDecomp(const matrix pMat,const matrix lMat,const matrix uMat,matrix & iMat,const ring R)352 bool luInverseFromLUDecomp(const matrix pMat, const matrix lMat,
353 const matrix uMat, matrix &iMat, const ring R)
354 { /* uMat is guaranteed to be quadratic */
355 //int d = uMat->rows();
356
357 matrix lMatInverse; /* for storing the inverse of lMat;
358 this inversion will always work */
359 matrix uMatInverse; /* for storing the inverse of uMat, if invertible */
360
361 bool result = upperRightTriangleInverse(uMat, uMatInverse, false);
362 if (result)
363 {
364 /* next will always work, since lMat is known to have all diagonal
365 entries equal to 1 */
366 lowerLeftTriangleInverse(lMat, lMatInverse, true);
367 iMat = mp_Mult(mp_Mult(uMatInverse, lMatInverse,R), pMat,R);
368
369 /* clean-up */
370 idDelete((ideal*)&lMatInverse);
371 idDelete((ideal*)&uMatInverse);
372 }
373
374 return result;
375 }
376
luSolveViaLUDecomp(const matrix pMat,const matrix lMat,const matrix uMat,const matrix bVec,matrix & xVec,matrix & H)377 bool luSolveViaLUDecomp(const matrix pMat, const matrix lMat,
378 const matrix uMat, const matrix bVec,
379 matrix &xVec, matrix &H)
380 {
381 int m = uMat->rows(); int n = uMat->cols();
382 matrix cVec = mpNew(m, 1); /* for storing pMat * bVec */
383 matrix yVec = mpNew(m, 1); /* for storing the unique solution of
384 lMat * yVec = cVec */
385
386 /* compute cVec = pMat * bVec but without actual multiplications */
387 for (int r = 1; r <= m; r++)
388 {
389 for (int c = 1; c <= m; c++)
390 {
391 if (MATELEM(pMat, r, c) != NULL)
392 { MATELEM(cVec, r, 1) = pCopy(MATELEM(bVec, c, 1)); break; }
393 }
394 }
395
396 /* solve lMat * yVec = cVec; this will always work as lMat is invertible;
397 moreover, no divisions are needed, since lMat[i, i] = 1, for all i */
398 for (int r = 1; r <= m; r++)
399 {
400 poly p = pNeg(pCopy(MATELEM(cVec, r, 1)));
401 for (int c = 1; c < r; c++)
402 p = pAdd(p, ppMult_qq(MATELEM(lMat, r, c), MATELEM(yVec, c, 1) ));
403 MATELEM(yVec, r, 1) = pNeg(p);
404 pNormalize(MATELEM(yVec, r, 1));
405 }
406
407 /* determine whether uMat * xVec = yVec is solvable */
408 bool isSolvable = true;
409 bool isZeroRow;
410 int nonZeroRowIndex = 0 ; // handle case that the matrix is zero
411 for (int r = m; r >= 1; r--)
412 {
413 isZeroRow = true;
414 for (int c = 1; c <= n; c++)
415 if (MATELEM(uMat, r, c) != NULL) { isZeroRow = false; break; }
416 if (isZeroRow)
417 {
418 if (MATELEM(yVec, r, 1) != NULL) { isSolvable = false; break; }
419 }
420 else { nonZeroRowIndex = r; break; }
421 }
422
423 if (isSolvable)
424 {
425 xVec = mpNew(n, 1);
426 matrix N = mpNew(n, n); int dim = 0;
427 poly p; poly q;
428 /* solve uMat * xVec = yVec and determine a basis of the solution
429 space of the homogeneous system uMat * xVec = 0;
430 We do not know in advance what the dimension (dim) of the latter
431 solution space will be. Thus, we start with the possibly too wide
432 matrix N and later copy the relevant columns of N into H. */
433 int nonZeroC = 0 ;
434 int lastNonZeroC = n + 1;
435
436 for (int r = nonZeroRowIndex; r >= 1; r--)
437 {
438 for (nonZeroC = 1; nonZeroC <= n; nonZeroC++)
439 if (MATELEM(uMat, r, nonZeroC) != NULL) break;
440
441 for (int w = lastNonZeroC - 1; w >= nonZeroC + 1; w--)
442 {
443 /* this loop will only be done when the given linear system has
444 more than one, i.e., infinitely many solutions */
445 dim++;
446 /* now we fill two entries of the dim-th column of N */
447 MATELEM(N, w, dim) = pNeg(pCopy(MATELEM(uMat, r, nonZeroC)));
448 MATELEM(N, nonZeroC, dim) = pCopy(MATELEM(uMat, r, w));
449 }
450 for (int d = 1; d <= dim; d++)
451 {
452 /* here we fill the entry of N at [r, d], for each valid vector
453 that we already have in N;
454 again, this will only be done when the given linear system has
455 more than one, i.e., infinitely many solutions */
456 p = NULL;
457 for (int c = nonZeroC + 1; c <= n; c++)
458 if (MATELEM(N, c, d) != NULL)
459 p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(N, c, d)));
460 q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
461 MATELEM(N, nonZeroC, d) = pMult(pNeg(p), q);
462 pNormalize(MATELEM(N, nonZeroC, d));
463 }
464 p = pNeg(pCopy(MATELEM(yVec, r, 1)));
465 for (int c = nonZeroC + 1; c <= n; c++)
466 if (MATELEM(xVec, c, 1) != NULL)
467 p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(xVec, c, 1)));
468 q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
469 MATELEM(xVec, nonZeroC, 1) = pMult(pNeg(p), q);
470 pNormalize(MATELEM(xVec, nonZeroC, 1));
471 lastNonZeroC = nonZeroC;
472 }
473 for (int w = lastNonZeroC - 1; w >= 1; w--)
474 {
475 // remaining variables are free
476 dim++;
477 MATELEM(N, w, dim) = pOne();
478 }
479
480 if (dim == 0)
481 {
482 /* that means the given linear system has exactly one solution;
483 we return as H the 1x1 matrix with entry zero */
484 H = mpNew(1, 1);
485 }
486 else
487 {
488 /* copy the first 'dim' columns of N into H */
489 H = mpNew(n, dim);
490 for (int r = 1; r <= n; r++)
491 for (int c = 1; c <= dim; c++)
492 MATELEM(H, r, c) = pCopy(MATELEM(N, r, c));
493 }
494 /* clean up N */
495 idDelete((ideal*)&N);
496 }
497 idDelete((ideal*)&cVec);
498 idDelete((ideal*)&yVec);
499
500 return isSolvable;
501 }
502
503 /* for debugging:
504 for printing numbers to the console
505 DELETE LATER */
printNumber(const number z)506 void printNumber(const number z)
507 {
508 if (nIsZero(z)) printf("number = 0\n");
509 else
510 {
511 poly p = pOne();
512 pSetCoeff(p, nCopy(z));
513 pSetm(p);
514 printf("number = %s\n", pString(p));
515 pDelete(&p);
516 }
517 }
518
519 /* for debugging:
520 for printing matrices to the console
521 DELETE LATER */
printMatrix(const matrix m)522 void printMatrix(const matrix m)
523 {
524 int rr = MATROWS(m); int cc = MATCOLS(m);
525 printf("\n-------------\n");
526 for (int r = 1; r <= rr; r++)
527 {
528 for (int c = 1; c <= cc; c++)
529 printf("%s ", pString(MATELEM(m, r, c)));
530 printf("\n");
531 }
532 printf("-------------\n");
533 }
534
535 /**
536 * Creates a new complex number from real and imaginary parts given
537 * by doubles.
538 *
539 * @return the new complex number
540 **/
complexNumber(const double r,const double i)541 number complexNumber(const double r, const double i)
542 {
543 gmp_complex* n= new gmp_complex(r, i);
544 return (number)n;
545 }
546
547 /**
548 * Returns one over the exponent-th power of ten.
549 *
550 * The method assumes that the base ring is the complex numbers.
551 *
552 * return one over the exponent-th power of 10
553 **/
tenToTheMinus(const int exponent)554 number tenToTheMinus(
555 const int exponent /**< [in] the exponent for 1/10 */
556 )
557 {
558 number ten = complexNumber(10.0, 0.0);
559 number result = complexNumber(1.0, 0.0);
560 number tmp;
561 /* compute 10^{-exponent} inside result by subsequent divisions by 10 */
562 for (int i = 1; i <= exponent; i++)
563 {
564 tmp = nDiv(result, ten);
565 nDelete(&result);
566 result = tmp;
567 }
568 nDelete(&ten);
569 return result;
570 }
571
572 /* for debugging:
573 for printing numbers to the console
574 DELETE LATER */
printSolutions(const int a,const int b,const int c)575 void printSolutions(const int a, const int b, const int c)
576 {
577 printf("\n------\n");
578 /* build the polynomial a*x^2 + b*x + c: */
579 poly p = NULL; poly q = NULL; poly r = NULL;
580 if (a != 0)
581 { p = pOne(); pSetExp(p, 1, 2); pSetm(p); pSetCoeff(p, nInit(a)); }
582 if (b != 0)
583 { q = pOne(); pSetExp(q, 1, 1); pSetm(q); pSetCoeff(q, nInit(b)); }
584 if (c != 0)
585 { r = pOne(); pSetCoeff(r, nInit(c)); }
586 p = pAdd(p, q); p = pAdd(p, r);
587 printf("poly = %s\n", pString(p));
588 number tol = tenToTheMinus(20);
589 number s1; number s2; int nSol = quadraticSolve(p, s1, s2, tol);
590 nDelete(&tol);
591 printf("solution code = %d\n", nSol);
592 if ((1 <= nSol) && (nSol <= 3))
593 {
594 if (nSol != 3) { printNumber(s1); nDelete(&s1); }
595 else { printNumber(s1); nDelete(&s1); printNumber(s2); nDelete(&s2); }
596 }
597 printf("------\n");
598 pDelete(&p);
599 }
600
realSqrt(const number n,const number tolerance,number & root)601 bool realSqrt(const number n, const number tolerance, number &root)
602 {
603 if (!nGreaterZero(n)) return false;
604 if (nIsZero(n)) return nInit(0);
605
606 number oneHalf = complexNumber(0.5, 0.0);
607 number nHalf = nMult(n, oneHalf);
608 root = nCopy(n);
609 number nOld = complexNumber(10.0, 0.0);
610 number nDiff = nCopy(nOld);
611
612 while (nGreater(nDiff, tolerance))
613 {
614 nDelete(&nOld);
615 nOld = root;
616 root = nAdd(nMult(oneHalf, nOld), nDiv(nHalf, nOld));
617 nDelete(&nDiff);
618 nDiff = nSub(nOld, root);
619 if (!nGreaterZero(nDiff)) nDiff = nInpNeg(nDiff);
620 }
621
622 nDelete(&nOld); nDelete(&nDiff); nDelete(&oneHalf); nDelete(&nHalf);
623 return true;
624 }
625
quadraticSolve(const poly p,number & s1,number & s2,const number tolerance)626 int quadraticSolve(const poly p, number &s1, number &s2,
627 const number tolerance)
628 {
629 poly q = pCopy(p);
630 int result;
631
632 if (q == NULL) result = -1;
633 else
634 {
635 int degree = pGetExp(q, 1);
636 if (degree == 0) result = 0; /* constant polynomial <> 0 */
637 else
638 {
639 number c2 = nInit(0); /* coefficient of var(1)^2 */
640 number c1 = nInit(0); /* coefficient of var(1)^1 */
641 number c0 = nInit(0); /* coefficient of var(1)^0 */
642 if (pGetExp(q, 1) == 2)
643 { nDelete(&c2); c2 = nCopy(pGetCoeff(q)); q = q->next; }
644 if ((q != NULL) && (pGetExp(q, 1) == 1))
645 { nDelete(&c1); c1 = nCopy(pGetCoeff(q)); q = q->next; }
646 if ((q != NULL) && (pGetExp(q, 1) == 0))
647 { nDelete(&c0); c0 = nCopy(pGetCoeff(q)); q = q->next; }
648
649 if (degree == 1)
650 {
651 c0 = nInpNeg(c0);
652 s1 = nDiv(c0, c1);
653 result = 1;
654 }
655 else
656 {
657 number tmp = nMult(c0, c2);
658 number tmp2 = nAdd(tmp, tmp); nDelete(&tmp);
659 number tmp4 = nAdd(tmp2, tmp2); nDelete(&tmp2);
660 number discr = nSub(nMult(c1, c1), tmp4); nDelete(&tmp4);
661 if (nIsZero(discr))
662 {
663 tmp = nAdd(c2, c2);
664 s1 = nDiv(c1, tmp); nDelete(&tmp);
665 s1 = nInpNeg(s1);
666 result = 2;
667 }
668 else if (nGreaterZero(discr))
669 {
670 realSqrt(discr, tolerance, tmp); /* sqrt of the discriminant */
671 tmp2 = nSub(tmp, c1);
672 tmp4 = nAdd(c2, c2);
673 s1 = nDiv(tmp2, tmp4); nDelete(&tmp2);
674 tmp = nInpNeg(tmp);
675 tmp2 = nSub(tmp, c1); nDelete(&tmp);
676 s2 = nDiv(tmp2, tmp4); nDelete(&tmp2); nDelete(&tmp4);
677 result = 3;
678 }
679 else
680 {
681 discr = nInpNeg(discr);
682 realSqrt(discr, tolerance, tmp); /* sqrt of |discriminant| */
683 tmp2 = nAdd(c2, c2);
684 tmp4 = nDiv(tmp, tmp2); nDelete(&tmp);
685 tmp = nDiv(c1, tmp2); nDelete(&tmp2);
686 tmp = nInpNeg(tmp);
687 s1 = (number)new gmp_complex(((gmp_complex*)tmp)->real(),
688 ((gmp_complex*)tmp4)->real());
689 tmp4 = nInpNeg(tmp4);
690 s2 = (number)new gmp_complex(((gmp_complex*)tmp)->real(),
691 ((gmp_complex*)tmp4)->real());
692 nDelete(&tmp); nDelete(&tmp4);
693 result = 3;
694 }
695 nDelete(&discr);
696 }
697 nDelete(&c0); nDelete(&c1); nDelete(&c2);
698 }
699 }
700 pDelete(&q);
701
702 return result;
703 }
704
euclideanNormSquared(const matrix aMat)705 number euclideanNormSquared(const matrix aMat)
706 {
707 int rr = MATROWS(aMat);
708 number result = nInit(0);
709 number tmp1; number tmp2;
710 for (int r = 1; r <= rr; r++)
711 if (MATELEM(aMat, r, 1) != NULL)
712 {
713 tmp1 = nMult(pGetCoeff(MATELEM(aMat, r, 1)),
714 pGetCoeff(MATELEM(aMat, r, 1)));
715 tmp2 = nAdd(result, tmp1); nDelete(&result); nDelete(&tmp1);
716 result = tmp2;
717 }
718 return result;
719 }
720
721 /* Returns a new number which is the absolute value of the coefficient
722 of the given polynomial.
723 *
724 * The method assumes that the coefficient has imaginary part zero. */
absValue(poly p)725 number absValue(poly p)
726 {
727 if (p == NULL) return nInit(0);
728 number result = nCopy(pGetCoeff(p));
729 if (!nGreaterZero(result)) result = nInpNeg(result);
730 return result;
731 }
732
subMatrix(const matrix aMat,const int rowIndex1,const int rowIndex2,const int colIndex1,const int colIndex2,matrix & subMat)733 bool subMatrix(const matrix aMat, const int rowIndex1, const int rowIndex2,
734 const int colIndex1, const int colIndex2, matrix &subMat)
735 {
736 if (rowIndex1 > rowIndex2) return false;
737 if (colIndex1 > colIndex2) return false;
738 int rr = rowIndex2 - rowIndex1 + 1;
739 int cc = colIndex2 - colIndex1 + 1;
740 subMat = mpNew(rr, cc);
741 for (int r = 1; r <= rr; r++)
742 for (int c = 1; c <= cc; c++)
743 MATELEM(subMat, r, c) =
744 pCopy(MATELEM(aMat, rowIndex1 + r - 1, colIndex1 + c - 1));
745 return true;
746 }
747
charPoly(const matrix aMat,poly & charPoly)748 bool charPoly(const matrix aMat, poly &charPoly)
749 {
750 if (MATROWS(aMat) != 2) return false;
751 if (MATCOLS(aMat) != 2) return false;
752 number b = nInit(0); number t;
753 if (MATELEM(aMat, 1, 1) != NULL)
754 { t = nAdd(b, pGetCoeff(MATELEM(aMat, 1, 1))); nDelete(&b); b = t;}
755 if (MATELEM(aMat, 2, 2) != NULL)
756 { t = nAdd(b, pGetCoeff(MATELEM(aMat, 2, 2))); nDelete(&b); b = t;}
757 b = nInpNeg(b);
758 number t1;
759 if ((MATELEM(aMat, 1, 1) != NULL) && (MATELEM(aMat, 2, 2) != NULL))
760 t1 = nMult(pGetCoeff(MATELEM(aMat, 1, 1)),
761 pGetCoeff(MATELEM(aMat, 2, 2)));
762 else t1 = nInit(0);
763 number t2;
764 if ((MATELEM(aMat, 1, 2) != NULL) && (MATELEM(aMat, 2, 1) != NULL))
765 t2 = nMult(pGetCoeff(MATELEM(aMat, 1, 2)),
766 pGetCoeff(MATELEM(aMat, 2, 1)));
767 else t2 = nInit(0);
768 number c = nSub(t1, t2); nDelete(&t1); nDelete(&t2);
769 poly p = pOne(); pSetExp(p, 1, 2); pSetm(p);
770 poly q = NULL;
771 if (!nIsZero(b))
772 { q = pOne(); pSetExp(q, 1, 1); pSetm(q); pSetCoeff(q, b); }
773 poly r = NULL;
774 if (!nIsZero(c))
775 { r = pOne(); pSetCoeff(r, c); }
776 p = pAdd(p, q); p = pAdd(p, r);
777 charPoly = p;
778 return true;
779 }
780
swapRows(int row1,int row2,matrix & aMat)781 void swapRows(int row1, int row2, matrix& aMat)
782 {
783 poly p;
784 int cc = MATCOLS(aMat);
785 for (int c = 1; c <= cc; c++)
786 {
787 p = MATELEM(aMat, row1, c);
788 MATELEM(aMat, row1, c) = MATELEM(aMat, row2, c);
789 MATELEM(aMat, row2, c) = p;
790 }
791 }
792
swapColumns(int column1,int column2,matrix & aMat)793 void swapColumns(int column1, int column2, matrix& aMat)
794 {
795 poly p;
796 int rr = MATROWS(aMat);
797 for (int r = 1; r <= rr; r++)
798 {
799 p = MATELEM(aMat, r, column1);
800 MATELEM(aMat, r, column1) = MATELEM(aMat, r, column2);
801 MATELEM(aMat, r, column2) = p;
802 }
803 }
804
matrixBlock(const matrix aMat,const matrix bMat,matrix & block)805 void matrixBlock(const matrix aMat, const matrix bMat, matrix &block)
806 {
807 int rowsA = MATROWS(aMat);
808 int rowsB = MATROWS(bMat);
809 int n = rowsA + rowsB;
810 block = mpNew(n, n);
811 for (int i = 1; i <= rowsA; i++)
812 for (int j = 1; j <= rowsA; j++)
813 MATELEM(block, i, j) = pCopy(MATELEM(aMat, i, j));
814 for (int i = 1; i <= rowsB; i++)
815 for (int j = 1; j <= rowsB; j++)
816 MATELEM(block, i + rowsA, j + rowsA) = pCopy(MATELEM(bMat, i, j));
817 }
818
819 /**
820 * Computes information related to one householder transformation step for
821 * constructing the Hessenberg form of a given non-derogative matrix.
822 *
823 * The method assumes that all matrix entries are numbers coming from some
824 * subfield of the reals. And that v has a non-zero first entry v_1 and a
825 * second non-zero entry somewhere else.
826 * Given such a vector v, it computes a number r (which will be the return
827 * value of the method), a vector u and a matrix P such that:
828 * 1) P * v = r * e_1,
829 * 2) P = E - u * u^T,
830 * 3) P = P^{-1}.
831 * Note that enforcing norm(u) = sqrt(2), which is done in the algorithm,
832 * guarantees property 3). This can be checked by expanding P^2 using
833 * property 2).
834 *
835 * @return the number r
836 **/
hessenbergStep(const matrix vVec,matrix & uVec,matrix & pMat,const number tolerance)837 number hessenbergStep(
838 const matrix vVec, /**< [in] the input vector v */
839 matrix &uVec, /**< [out] the output vector u */
840 matrix &pMat, /**< [out] the output matrix P */
841 const number tolerance /**< [in] accuracy for square roots */
842 )
843 {
844 int rr = MATROWS(vVec);
845 number vNormSquared = euclideanNormSquared(vVec);
846 number vNorm; realSqrt(vNormSquared, tolerance, vNorm);
847 /* v1 is guaranteed to be non-zero: */
848 number v1 = pGetCoeff(MATELEM(vVec, 1, 1));
849 bool v1Sign = true; if (nGreaterZero(v1)) v1Sign = false;
850
851 number v1Abs = nCopy(v1); if (v1Sign) v1Abs = nInpNeg(v1Abs);
852 number t1 = nDiv(v1Abs, vNorm);
853 number one = nInit(1);
854 number t2 = nAdd(t1, one); nDelete(&t1);
855 number denominator; realSqrt(t2, tolerance, denominator); nDelete(&t2);
856 uVec = mpNew(rr, 1);
857 t1 = nDiv(v1Abs, vNorm);
858 t2 = nAdd(t1, one); nDelete(&t1);
859 t1 = nDiv(t2, denominator); nDelete(&t2);
860 MATELEM(uVec, 1, 1) = pOne();
861 pSetCoeff(MATELEM(uVec, 1, 1), t1); /* we know that t1 != 0 */
862 for (int r = 2; r <= rr; r++)
863 {
864 if (MATELEM(vVec, r, 1) != NULL)
865 t1 = nCopy(pGetCoeff(MATELEM(vVec, r, 1)));
866 else t1 = nInit(0);
867 if (v1Sign) t1 = nInpNeg(t1);
868 t2 = nDiv(t1, vNorm); nDelete(&t1);
869 t1 = nDiv(t2, denominator); nDelete(&t2);
870 if (!nIsZero(t1))
871 {
872 MATELEM(uVec, r, 1) = pOne();
873 pSetCoeff(MATELEM(uVec, r, 1), t1);
874 }
875 else nDelete(&t1);
876 }
877 nDelete(&denominator);
878
879 /* finished building vector u; now turn to pMat */
880 pMat = mpNew(rr, rr);
881 /* we set P := E - u * u^T, as desired */
882 for (int r = 1; r <= rr; r++)
883 for (int c = 1; c <= rr; c++)
884 {
885 if ((MATELEM(uVec, r, 1) != NULL) && (MATELEM(uVec, c, 1) != NULL))
886 t1 = nMult(pGetCoeff(MATELEM(uVec, r, 1)),
887 pGetCoeff(MATELEM(uVec, c, 1)));
888 else t1 = nInit(0);
889 if (r == c) { t2 = nSub(one, t1); nDelete(&t1); }
890 else t2 = nInpNeg(t1);
891 if (!nIsZero(t2))
892 {
893 MATELEM(pMat, r, c) = pOne();
894 pSetCoeff(MATELEM(pMat, r, c), t2);
895 }
896 else nDelete(&t2);
897 }
898 nDelete(&one);
899 /* finished building pMat; now compute the return value */
900 t1 = vNormSquared; if (v1Sign) t1 = nInpNeg(t1);
901 t2 = nMult(v1, vNorm);
902 number t3 = nAdd(t1, t2); nDelete(&t1); nDelete(&t2);
903 t1 = nAdd(v1Abs, vNorm); nDelete(&v1Abs); nDelete(&vNorm);
904 t2 = nDiv(t3, t1); nDelete(&t1); nDelete(&t3);
905 t2 = nInpNeg(t2);
906 return t2;
907 }
908
hessenberg(const matrix aMat,matrix & pMat,matrix & hessenbergMat,const number tolerance,const ring R)909 void hessenberg(const matrix aMat, matrix &pMat, matrix &hessenbergMat,
910 const number tolerance, const ring R)
911 {
912 int n = MATROWS(aMat);
913 unitMatrix(n, pMat);
914 subMatrix(aMat, 1, n, 1, n, hessenbergMat);
915 for (int c = 1; c <= n; c++)
916 {
917 /* find one or two non-zero entries in the current column */
918 int r1 = 0; int r2 = 0;
919 for (int r = c + 1; r <= n; r++)
920 if (MATELEM(hessenbergMat, r, c) != NULL)
921 {
922 if (r1 == 0) r1 = r;
923 else if (r2 == 0) { r2 = r; break; }
924 }
925 if (r1 != 0)
926 { /* At least one entry in the current column is non-zero. */
927 if (r1 != c + 1)
928 { /* swap rows to bring non-zero element to row with index c + 1 */
929 swapRows(r1, c + 1, hessenbergMat);
930 /* now also swap columns to reflect action of permutation
931 from the right-hand side */
932 swapColumns(r1, c + 1, hessenbergMat);
933 /* include action of permutation also in pMat */
934 swapRows(r1, c + 1, pMat);
935 }
936 if (r2 != 0)
937 { /* There is at least one more non-zero element in the current
938 column. So let us perform a hessenberg step in order to make
939 all additional non-zero elements zero. */
940 matrix v; subMatrix(hessenbergMat, c + 1, n, c, c, v);
941 matrix u; matrix pTmp;
942 number r = hessenbergStep(v, u, pTmp, tolerance);
943 idDelete((ideal*)&v); idDelete((ideal*)&u); nDelete(&r);
944 /* pTmp just acts on the lower right block of hessenbergMat;
945 i.e., it needs to be extended by a unit matrix block at the top
946 left in order to define a whole transformation matrix;
947 this code may be optimized */
948 unitMatrix(c, u);
949 matrix pTmpFull; matrixBlock(u, pTmp, pTmpFull);
950 idDelete((ideal*)&u); idDelete((ideal*)&pTmp);
951 /* now include pTmpFull in pMat (by letting it act from the left) */
952 pTmp = mp_Mult(pTmpFull, pMat,R); idDelete((ideal*)&pMat);
953 pMat = pTmp;
954 /* now let pTmpFull act on hessenbergMat from the left and from the
955 right (note that pTmpFull is self-inverse) */
956 pTmp = mp_Mult(pTmpFull, hessenbergMat,R);
957 idDelete((ideal*)&hessenbergMat);
958 hessenbergMat = mp_Mult(pTmp, pTmpFull, R);
959 idDelete((ideal*)&pTmp); idDelete((ideal*)&pTmpFull);
960 /* as there may be inaccuracy, we erase those entries of hessenbergMat
961 which must have become zero by the last transformation */
962 for (int r = c + 2; r <= n; r++)
963 pDelete(&MATELEM(hessenbergMat, r, c));
964 }
965 }
966 }
967 }
968
969 /**
970 * Performs one transformation step on the given matrix H as part of
971 * the gouverning QR double shift algorithm.
972 * The method will change the given matrix H side-effect-wise. The resulting
973 * matrix H' will be in Hessenberg form.
974 * The iteration index is needed, since for the 11th and 21st iteration,
975 * the transformation step is different from the usual step, to avoid
976 * convergence problems of the gouverning QR double shift process (who is
977 * also the only caller of this method).
978 **/
mpTrafo(matrix & H,int it,const number tolerance,const ring R)979 void mpTrafo(
980 matrix &H, /**< [in/out] the matrix to be transformed */
981 int it, /**< [in] iteration index */
982 const number tolerance,/**< [in] accuracy for square roots */
983 const ring R
984 )
985 {
986 int n = MATROWS(H);
987 number trace; number det; number tmp1; number tmp2; number tmp3;
988
989 if ((it != 11) && (it != 21)) /* the standard case */
990 {
991 /* in this case 'trace' will really be the trace of the lowermost
992 (2x2) block of hMat */
993 trace = nInit(0);
994 det = nInit(0);
995 if (MATELEM(H, n - 1, n - 1) != NULL)
996 {
997 tmp1 = nAdd(trace, pGetCoeff(MATELEM(H, n - 1, n - 1)));
998 nDelete(&trace);
999 trace = tmp1;
1000 }
1001 if (MATELEM(H, n, n) != NULL)
1002 {
1003 tmp1 = nAdd(trace, pGetCoeff(MATELEM(H, n, n)));
1004 nDelete(&trace);
1005 trace = tmp1;
1006 }
1007 /* likewise 'det' will really be the determinante of the lowermost
1008 (2x2) block of hMat */
1009 if ((MATELEM(H, n - 1, n - 1 ) != NULL) && (MATELEM(H, n, n) != NULL))
1010 {
1011 tmp1 = nMult(pGetCoeff(MATELEM(H, n - 1, n - 1)),
1012 pGetCoeff(MATELEM(H, n, n)));
1013 tmp2 = nAdd(tmp1, det); nDelete(&tmp1); nDelete(&det);
1014 det = tmp2;
1015 }
1016 if ((MATELEM(H, n - 1, n) != NULL) && (MATELEM(H, n, n - 1) != NULL))
1017 {
1018 tmp1 = nMult(pGetCoeff(MATELEM(H, n - 1, n)),
1019 pGetCoeff(MATELEM(H, n, n - 1)));
1020 tmp2 = nSub(det, tmp1); nDelete(&tmp1); nDelete(&det);
1021 det = tmp2;
1022 }
1023 }
1024 else
1025 {
1026 /* for it = 11 or it = 21, we use special formulae to avoid convergence
1027 problems of the gouverning QR double shift algorithm (who is the only
1028 caller of this method) */
1029 /* trace = 3/2 * (|hMat[n, n-1]| + |hMat[n-1, n-2]|) */
1030 tmp1 = nInit(0);
1031 if (MATELEM(H, n, n - 1) != NULL)
1032 { nDelete(&tmp1); tmp1 = nCopy(pGetCoeff(MATELEM(H, n, n - 1))); }
1033 if (!nGreaterZero(tmp1)) tmp1 = nInpNeg(tmp1);
1034 tmp2 = nInit(0);
1035 if (MATELEM(H, n - 1, n - 2) != NULL)
1036 { nDelete(&tmp2); tmp2 = nCopy(pGetCoeff(MATELEM(H, n - 1, n - 2))); }
1037 if (!nGreaterZero(tmp2)) tmp2 = nInpNeg(tmp2);
1038 tmp3 = nAdd(tmp1, tmp2); nDelete(&tmp1); nDelete(&tmp2);
1039 tmp1 = nInit(3); tmp2 = nInit(2);
1040 trace = nDiv(tmp1, tmp2); nDelete(&tmp1); nDelete(&tmp2);
1041 tmp1 = nMult(tmp3, trace); nDelete(&trace);
1042 trace = tmp1;
1043 /* det = (|hMat[n, n-1]| + |hMat[n-1, n-2]|)^2 */
1044 det = nMult(tmp3, tmp3); nDelete(&tmp3);
1045 }
1046 matrix c = mpNew(n, 1);
1047 trace = nInpNeg(trace);
1048 MATELEM(c,1,1) = pAdd(pAdd(pAdd(ppMult_qq(MATELEM(H,1,1), MATELEM(H,1,1)),
1049 ppMult_qq(MATELEM(H,1,2), MATELEM(H,2,1))),
1050 __pp_Mult_nn(MATELEM(H,1,1), trace, currRing)),
1051 __p_Mult_nn(pOne(), det,currRing));
1052 MATELEM(c,2,1) = pAdd(pMult(pCopy(MATELEM(H,2,1)),
1053 pAdd(pCopy(MATELEM(H,1,1)),
1054 pCopy(MATELEM(H,2,2)))),
1055 __pp_Mult_nn(MATELEM(H,2,1), trace,currRing));
1056 MATELEM(c,3,1) = ppMult_qq(MATELEM(H,2,1), MATELEM(H,3,2));
1057 nDelete(&trace); nDelete(&det);
1058
1059 /* for applying hessenbergStep, we need to make sure that c[1, 1] is
1060 not zero */
1061 if ((MATELEM(c,1,1) != NULL) &&
1062 ((MATELEM(c,2,1) != NULL) || (MATELEM(c,3,1) != NULL)))
1063 {
1064 matrix uVec; matrix hMat;
1065 tmp1 = hessenbergStep(c, uVec, hMat, tolerance); nDelete(&tmp1);
1066 /* now replace H by hMat * H * hMat: */
1067 matrix wMat = mp_Mult(hMat, H,R); idDelete((ideal*)&H);
1068 matrix H1 = mp_Mult(wMat, hMat,R);
1069 idDelete((ideal*)&wMat); idDelete((ideal*)&hMat);
1070 /* now need to re-establish Hessenberg form of H1 and put it in H */
1071 hessenberg(H1, wMat, H, tolerance,R);
1072 idDelete((ideal*)&wMat); idDelete((ideal*)&H1);
1073 }
1074 else if ((MATELEM(c,1,1) == NULL) && (MATELEM(c,2,1) != NULL))
1075 {
1076 swapRows(1, 2, H);
1077 swapColumns(1, 2, H);
1078 }
1079 else if ((MATELEM(c,1,1) == NULL) && (MATELEM(c,3,1) != NULL))
1080 {
1081 swapRows(1, 3, H);
1082 swapColumns(1, 3, H);
1083 }
1084 else
1085 { /* c is the zero vector or a multiple of e_1;
1086 no hessenbergStep needed */ }
1087 }
1088
1089 /* helper for qrDoubleShift */
qrDS(const int,matrix * queue,int & queueL,number * eigenValues,int & eigenValuesL,const number tol1,const number tol2,const ring R)1090 bool qrDS(
1091 const int /*n*/,
1092 matrix* queue,
1093 int& queueL,
1094 number* eigenValues,
1095 int& eigenValuesL,
1096 const number tol1,
1097 const number tol2,
1098 const ring R
1099 )
1100 {
1101 bool deflationFound = true;
1102 /* we loop until the working queue is empty,
1103 provided we always find deflation */
1104 while (deflationFound && (queueL > 0))
1105 {
1106 /* take out last queue entry */
1107 matrix currentMat = queue[queueL - 1]; queueL--;
1108 int m = MATROWS(currentMat);
1109 if (m == 1)
1110 {
1111 number newEigenvalue;
1112 /* the entry at [1, 1] is the eigenvalue */
1113 if (MATELEM(currentMat, 1, 1) == NULL) newEigenvalue = nInit(0);
1114 else newEigenvalue = nCopy(pGetCoeff(MATELEM(currentMat, 1, 1)));
1115 eigenValues[eigenValuesL++] = newEigenvalue;
1116 }
1117 else if (m == 2)
1118 {
1119 /* there are two eigenvalues which come as zeros of the characteristic
1120 polynomial */
1121 poly p; charPoly(currentMat, p);
1122 number s1; number s2;
1123 int nSol = quadraticSolve(p, s1, s2, tol2); pDelete(&p);
1124 assume(nSol >= 2);
1125 eigenValues[eigenValuesL++] = s1;
1126 /* if nSol = 2, then s1 is a double zero, and s2 is invalid: */
1127 if (nSol == 2) s2 = nCopy(s1);
1128 eigenValues[eigenValuesL++] = s2;
1129 }
1130 else /* m > 2 */
1131 {
1132 /* bring currentMat into Hessenberg form to fasten computations: */
1133 matrix mm1; matrix mm2;
1134 hessenberg(currentMat, mm1, mm2, tol2,R);
1135 idDelete((ideal*)¤tMat); idDelete((ideal*)&mm1);
1136 currentMat = mm2;
1137 int it = 1; bool doLoop = true;
1138 while (doLoop && (it <= 30 * m))
1139 {
1140 /* search for deflation */
1141 number w1; number w2;
1142 number test1; number test2; bool stopCriterion = false; int k;
1143 for (k = 1; k < m; k++)
1144 {
1145 test1 = absValue(MATELEM(currentMat, k + 1, k));
1146 w1 = absValue(MATELEM(currentMat, k, k));
1147 w2 = absValue(MATELEM(currentMat, k + 1, k + 1));
1148 test2 = nMult(tol1, nAdd(w1, w2));
1149 nDelete(&w1); nDelete(&w2);
1150 if (!nGreater(test1, test2)) stopCriterion = true;
1151 nDelete(&test1); nDelete(&test2);
1152 if (stopCriterion) break;
1153 }
1154 if (k < m) /* found deflation at position (k + 1, k) */
1155 {
1156 pDelete(&MATELEM(currentMat, k + 1, k)); /* make this entry zero */
1157 subMatrix(currentMat, 1, k, 1, k, queue[queueL++]);
1158 subMatrix(currentMat, k + 1, m, k + 1, m, queue[queueL++]);
1159 doLoop = false;
1160 }
1161 else /* no deflation found yet */
1162 {
1163 mpTrafo(currentMat, it, tol2,R);
1164 it++; /* try again */
1165 }
1166 }
1167 if (doLoop) /* could not find deflation for currentMat */
1168 {
1169 deflationFound = false;
1170 }
1171 idDelete((ideal*)¤tMat);
1172 }
1173 }
1174 return deflationFound;
1175 }
1176
1177 /**
1178 * Tries to find the number n in the array nn[0..nnLength-1].
1179 *
1180 * The method assumes that the ground field is the complex numbers.
1181 * n and an entry of nn will be regarded equal when the absolute
1182 * value of their difference is not larger than the given tolerance.
1183 * In this case, the index of the respective entry of nn is returned,
1184 * otherwise -1.
1185 *
1186 * @return the index of n in nn (up to tolerance) or -1
1187 **/
similar(const number * nn,const int nnLength,const number n,const number tolerance)1188 int similar(
1189 const number* nn, /**< [in] array of numbers to look-up */
1190 const int nnLength, /**< [in] length of nn */
1191 const number n, /**< [in] number to loop-up in nn */
1192 const number tolerance /**< [in] tolerance for comparison */
1193 )
1194 {
1195 int result = -1;
1196 number tt = nMult(tolerance, tolerance);
1197 number nr = (number)new gmp_complex(((gmp_complex*)n)->real());
1198 number ni = (number)new gmp_complex(((gmp_complex*)n)->imag());
1199 number rr; number ii;
1200 number w1; number w2; number w3; number w4; number w5;
1201 for (int i = 0; i < nnLength; i++)
1202 {
1203 rr = (number)new gmp_complex(((gmp_complex*)nn[i])->real());
1204 ii = (number)new gmp_complex(((gmp_complex*)nn[i])->imag());
1205 w1 = nSub(nr, rr); w2 = nMult(w1, w1);
1206 w3 = nSub(ni, ii); w4 = nMult(w3, w3);
1207 w5 = nAdd(w2, w4);
1208 if (!nGreater(w5, tt)) result = i;
1209 nDelete(&w1); nDelete(&w2); nDelete(&w3); nDelete(&w4);
1210 nDelete(&w5); nDelete(&rr); nDelete(&ii);
1211 if (result != -1) break;
1212 }
1213 nDelete(&tt); nDelete(&nr); nDelete(&ni);
1214 return result;
1215 }
1216
1217 /* This codes assumes that there are at least two variables in the current
1218 base ring. No assumption is made regarding the monomial ordering. */
henselFactors(const int xIndex,const int yIndex,const poly h,const poly f0,const poly g0,const int d,poly & f,poly & g)1219 void henselFactors(const int xIndex, const int yIndex, const poly h,
1220 const poly f0, const poly g0, const int d, poly &f, poly &g)
1221 {
1222 int n = (int)p_Deg(f0,currRing);
1223 int m = (int)p_Deg(g0,currRing);
1224 matrix aMat = mpNew(n + m, n + m); /* matrix A for linear system */
1225 matrix pMat; matrix lMat; matrix uMat; /* for the decomposition of A */
1226 f = pCopy(f0); g = pCopy(g0); /* initially: h = f*g mod <x^1> */
1227
1228 /* initial step: read off coefficients of f0, and g0 */
1229 poly p = f0; poly matEntry; number c;
1230 while (p != NULL)
1231 {
1232 c = nCopy(pGetCoeff(p));
1233 matEntry = pOne(); pSetCoeff(matEntry, c);
1234 MATELEM(aMat, pGetExp(p, yIndex) + 1, 1) = matEntry;
1235 p = pNext(p);
1236 }
1237 p = g0;
1238 while (p != NULL)
1239 {
1240 c = nCopy(pGetCoeff(p));
1241 matEntry = pOne(); pSetCoeff(matEntry, c);
1242 MATELEM(aMat, pGetExp(p, yIndex) + 1, m + 1) = matEntry;
1243 p = pNext(p);
1244 }
1245 /* fill the rest of A */
1246 for (int row = 2; row <= n + 1; row++)
1247 for (int col = 2; col <= m; col++)
1248 {
1249 if (col > row) break;
1250 MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1));
1251 }
1252 for (int row = n + 2; row <= n + m; row++)
1253 for (int col = row - n; col <= m; col++)
1254 MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1));
1255 for (int row = 2; row <= m + 1; row++)
1256 for (int col = m + 2; col <= m + n; col++)
1257 {
1258 if (col - m > row) break;
1259 MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1));
1260 }
1261 for (int row = m + 2; row <= n + m; row++)
1262 for (int col = row; col <= m + n; col++)
1263 MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1));
1264
1265 /* constructing the LU-decomposition of A */
1266 luDecomp(aMat, pMat, lMat, uMat);
1267
1268 /* Before the xExp-th loop, we know that h = f*g mod <x^xExp>.
1269 Afterwards the algorithm ensures h = f*g mod <x^(xExp + 1)>.
1270 Hence in the end we obtain f and g as required, i.e.,
1271 h = f*g mod <x^(d+1)>.
1272 The algorithm works by solving a (m+n)x(m+n) linear equation system
1273 A*x = b with constant matrix A (as decomposed above). By theory, the
1274 system is guaranteed to have a unique solution. */
1275 poly fg = ppMult_qq(f, g); /* for storing the product of f and g */
1276 for (int xExp = 1; xExp <= d; xExp++)
1277 {
1278 matrix bVec = mpNew(n + m, 1); /* b */
1279 matrix xVec = mpNew(n + m, 1); /* x */
1280
1281 p = pCopy(fg);
1282 p = pAdd(pCopy(h), pNeg(p)); /* p = h - f*g */
1283 /* we collect all terms in p with x-exponent = xExp and use their
1284 coefficients to build the vector b */
1285 int bIsZeroVector = true;
1286 while (p != NULL)
1287 {
1288 if (pGetExp(p, xIndex) == xExp)
1289 {
1290 number c = nCopy(pGetCoeff(p));
1291 poly matEntry = pOne(); pSetCoeff(matEntry, c);
1292 MATELEM(bVec, pGetExp(p, yIndex) + 1, 1) = matEntry;
1293 if (matEntry != NULL) bIsZeroVector = false;
1294 }
1295 pLmDelete(&p); /* destruct leading term of p and move to next term */
1296 }
1297 /* solve the linear equation system */
1298 if (!bIsZeroVector) /* otherwise x = 0 and f, g do not change */
1299 {
1300 matrix notUsedMat;
1301 luSolveViaLUDecomp(pMat, lMat, uMat, bVec, xVec, notUsedMat);
1302 idDelete((ideal*)¬UsedMat);
1303 /* update f and g by newly computed terms, and update f*g */
1304 poly fNew = NULL; poly gNew = NULL;
1305 for (int row = 1; row <= m; row++)
1306 {
1307 if (MATELEM(xVec, row, 1) != NULL)
1308 {
1309 p = pCopy(MATELEM(xVec, row, 1)); /* p = c */
1310 pSetExp(p, xIndex, xExp); /* p = c * x^xExp */
1311 pSetExp(p, yIndex, row - 1); /* p = c * x^xExp * y^i */
1312 pSetm(p);
1313 gNew = pAdd(gNew, p);
1314 }
1315 }
1316 for (int row = m + 1; row <= m + n; row++)
1317 {
1318 if (MATELEM(xVec, row, 1) != NULL)
1319 {
1320 p = pCopy(MATELEM(xVec, row, 1)); /* p = c */
1321 pSetExp(p, xIndex, xExp); /* p = c * x^xExp */
1322 pSetExp(p, yIndex, row - m - 1); /* p = c * x^xExp * y^i */
1323 pSetm(p);
1324 fNew = pAdd(fNew, p);
1325 }
1326 }
1327 fg = pAdd(fg, ppMult_qq(f, gNew));
1328 fg = pAdd(fg, ppMult_qq(g, fNew));
1329 fg = pAdd(fg, ppMult_qq(fNew, gNew));
1330 f = pAdd(f, fNew);
1331 g = pAdd(g, gNew);
1332 }
1333 /* clean-up loop-dependent vectors */
1334 idDelete((ideal*)&bVec); idDelete((ideal*)&xVec);
1335 }
1336
1337 /* clean-up matrices A, P, L and U, and polynomial fg */
1338 idDelete((ideal*)&aMat); idDelete((ideal*)&pMat);
1339 idDelete((ideal*)&lMat); idDelete((ideal*)&uMat);
1340 pDelete(&fg);
1341 }
1342
lduDecomp(const matrix aMat,matrix & pMat,matrix & lMat,matrix & dMat,matrix & uMat,poly & l,poly & u,poly & lTimesU)1343 void lduDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &dMat,
1344 matrix &uMat, poly &l, poly &u, poly &lTimesU)
1345 {
1346 int rr = aMat->rows();
1347 int cc = aMat->cols();
1348 /* we use an int array to store all row permutations;
1349 note that we only make use of the entries [1..rr] */
1350 int* permut = new int[rr + 1];
1351 for (int i = 1; i <= rr; i++) permut[i] = i;
1352 /* fill lMat and dMat with the (rr x rr) unit matrix */
1353 unitMatrix(rr, lMat);
1354 unitMatrix(rr, dMat);
1355 uMat = mpNew(rr, cc);
1356 /* copy aMat into uMat: */
1357 for (int r = 1; r <= rr; r++)
1358 for (int c = 1; c <= cc; c++)
1359 MATELEM(uMat, r, c) = pCopy(MATELEM(aMat, r, c));
1360 u = pOne(); l = pOne();
1361
1362 int col = 1; int row = 1;
1363 while ((col <= cc) & (row < rr))
1364 {
1365 int pivotR; int pivotC; bool pivotValid = false;
1366 while (col <= cc)
1367 {
1368 pivotValid = pivot(uMat, row, rr, col, col, &pivotR, &pivotC);
1369 if (pivotValid) break;
1370 col++;
1371 }
1372 if (pivotValid)
1373 {
1374 if (pivotR != row)
1375 {
1376 swapRows(row, pivotR, uMat);
1377 poly p = MATELEM(dMat, row, row);
1378 MATELEM(dMat, row, row) = MATELEM(dMat, pivotR, pivotR);
1379 MATELEM(dMat, pivotR, pivotR) = p;
1380 swapColumns(row, pivotR, lMat);
1381 swapRows(row, pivotR, lMat);
1382 int temp = permut[row];
1383 permut[row] = permut[pivotR]; permut[pivotR] = temp;
1384 }
1385 /* in gg, we compute the gcd of all non-zero elements in
1386 uMat[row..rr, col];
1387 the next number is the pivot and thus guaranteed to be different
1388 from zero: */
1389 number gg = nCopy(pGetCoeff(MATELEM(uMat, row, col))); number t;
1390 for (int r = row + 1; r <= rr; r++)
1391 {
1392 if (MATELEM(uMat, r, col) != NULL)
1393 {
1394 t = gg;
1395 gg = n_Gcd(t, pGetCoeff(MATELEM(uMat, r, col)),currRing->cf);
1396 nDelete(&t);
1397 }
1398 }
1399 t = nDiv(pGetCoeff(MATELEM(uMat, row, col)), gg);
1400 nNormalize(t); /* this division works without remainder */
1401 if (!nIsOne(t))
1402 {
1403 for (int r = row; r <= rr; r++)
1404 MATELEM(dMat, r, r)=__p_Mult_nn(MATELEM(dMat, r, r), t,currRing);
1405 MATELEM(lMat, row, row)=__p_Mult_nn(MATELEM(lMat, row, row), t,currRing);
1406 }
1407 l = pMult(l, pCopy(MATELEM(lMat, row, row)));
1408 u = pMult(u, pCopy(MATELEM(uMat, row, col)));
1409
1410 for (int r = row + 1; r <= rr; r++)
1411 {
1412 if (MATELEM(uMat, r, col) != NULL)
1413 {
1414 number g = n_Gcd(pGetCoeff(MATELEM(uMat, row, col)),
1415 pGetCoeff(MATELEM(uMat, r, col)),
1416 currRing->cf);
1417 number f1 = nDiv(pGetCoeff(MATELEM(uMat, r, col)), g);
1418 nNormalize(f1); /* this division works without remainder */
1419 number f2 = nDiv(pGetCoeff(MATELEM(uMat, row, col)), g);
1420 nNormalize(f2); /* this division works without remainder */
1421 pDelete(&MATELEM(uMat, r, col)); MATELEM(uMat, r, col) = NULL;
1422 for (int c = col + 1; c <= cc; c++)
1423 {
1424 poly p = MATELEM(uMat, r, c);
1425 p=__p_Mult_nn(p, f2,currRing);
1426 poly q = pCopy(MATELEM(uMat, row, c));
1427 q=__p_Mult_nn(q, f1,currRing); q = pNeg(q);
1428 MATELEM(uMat, r, c) = pAdd(p, q);
1429 }
1430 number tt = nDiv(g, gg);
1431 nNormalize(tt); /* this division works without remainder */
1432 MATELEM(lMat, r, r)=__p_Mult_nn(MATELEM(lMat, r, r), tt, currRing);
1433 nDelete(&tt);
1434 MATELEM(lMat, r, row) = pCopy(MATELEM(lMat, r, r));
1435 MATELEM(lMat, r, row)=__p_Mult_nn(MATELEM(lMat, r, row), f1,currRing);
1436 nDelete(&f1); nDelete(&f2); nDelete(&g);
1437 }
1438 else
1439 MATELEM(lMat, r, r)=__p_Mult_nn(MATELEM(lMat, r, r), t, currRing);
1440 }
1441 nDelete(&t); nDelete(&gg);
1442 col++; row++;
1443 }
1444 }
1445 /* one factor in the product u might be missing: */
1446 if (row == rr)
1447 {
1448 while ((col <= cc) && (MATELEM(uMat, row, col) == NULL)) col++;
1449 if (col <= cc) u = pMult(u, pCopy(MATELEM(uMat, row, col)));
1450 }
1451
1452 /* building the permutation matrix from 'permut' and computing l */
1453 pMat = mpNew(rr, rr);
1454 for (int r = 1; r <= rr; r++)
1455 MATELEM(pMat, r, permut[r]) = pOne();
1456 delete[] permut;
1457
1458 lTimesU = ppMult_qq(l, u);
1459 }
1460
luSolveViaLDUDecomp(const matrix pMat,const matrix lMat,const matrix dMat,const matrix uMat,const poly l,const poly u,const poly lTimesU,const matrix bVec,matrix & xVec,matrix & H)1461 bool luSolveViaLDUDecomp(const matrix pMat, const matrix lMat,
1462 const matrix dMat, const matrix uMat,
1463 const poly l, const poly u, const poly lTimesU,
1464 const matrix bVec, matrix &xVec, matrix &H)
1465 {
1466 int m = uMat->rows(); int n = uMat->cols();
1467 matrix cVec = mpNew(m, 1); /* for storing l * pMat * bVec */
1468 matrix yVec = mpNew(m, 1); /* for storing the unique solution of
1469 lMat * yVec = cVec */
1470
1471 /* compute cVec = l * pMat * bVec but without actual matrix mult. */
1472 for (int r = 1; r <= m; r++)
1473 {
1474 for (int c = 1; c <= m; c++)
1475 {
1476 if (MATELEM(pMat, r, c) != NULL)
1477 { MATELEM(cVec, r, 1) = ppMult_qq(l, MATELEM(bVec, c, 1)); break; }
1478 }
1479 }
1480
1481 /* solve lMat * yVec = cVec; this will always work as lMat is invertible;
1482 moreover, all divisions are guaranteed to be without remainder */
1483 poly p; poly q;
1484 for (int r = 1; r <= m; r++)
1485 {
1486 p = pNeg(pCopy(MATELEM(cVec, r, 1)));
1487 for (int c = 1; c < r; c++)
1488 p = pAdd(p, ppMult_qq(MATELEM(lMat, r, c), MATELEM(yVec, c, 1) ));
1489 /* The following division is without remainder. */
1490 q = pNSet(nInvers(pGetCoeff(MATELEM(lMat, r, r))));
1491 MATELEM(yVec, r, 1) = pMult(pNeg(p), q);
1492 pNormalize(MATELEM(yVec, r, 1));
1493 }
1494
1495 /* compute u * dMat * yVec and put result into yVec */
1496 for (int r = 1; r <= m; r++)
1497 {
1498 p = ppMult_qq(u, MATELEM(dMat, r, r));
1499 MATELEM(yVec, r, 1) = pMult(p, MATELEM(yVec, r, 1));
1500 }
1501
1502 /* determine whether uMat * xVec = yVec is solvable */
1503 bool isSolvable = true;
1504 bool isZeroRow; int nonZeroRowIndex;
1505 for (int r = m; r >= 1; r--)
1506 {
1507 isZeroRow = true;
1508 for (int c = 1; c <= n; c++)
1509 if (MATELEM(uMat, r, c) != NULL) { isZeroRow = false; break; }
1510 if (isZeroRow)
1511 {
1512 if (MATELEM(yVec, r, 1) != NULL) { isSolvable = false; break; }
1513 }
1514 else { nonZeroRowIndex = r; break; }
1515 }
1516
1517 if (isSolvable)
1518 {
1519 xVec = mpNew(n, 1);
1520 matrix N = mpNew(n, n); int dim = 0;
1521 /* solve uMat * xVec = yVec and determine a basis of the solution
1522 space of the homogeneous system uMat * xVec = 0;
1523 We do not know in advance what the dimension (dim) of the latter
1524 solution space will be. Thus, we start with the possibly too wide
1525 matrix N and later copy the relevant columns of N into H. */
1526 int nonZeroC; int lastNonZeroC = n + 1;
1527 for (int r = nonZeroRowIndex; r >= 1; r--)
1528 {
1529 for (nonZeroC = 1; nonZeroC <= n; nonZeroC++)
1530 if (MATELEM(uMat, r, nonZeroC) != NULL) break;
1531 for (int w = lastNonZeroC - 1; w >= nonZeroC + 1; w--)
1532 {
1533 /* this loop will only be done when the given linear system has
1534 more than one, i.e., infinitely many solutions */
1535 dim++;
1536 /* now we fill two entries of the dim-th column of N */
1537 MATELEM(N, w, dim) = pNeg(pCopy(MATELEM(uMat, r, nonZeroC)));
1538 MATELEM(N, nonZeroC, dim) = pCopy(MATELEM(uMat, r, w));
1539 }
1540 for (int d = 1; d <= dim; d++)
1541 {
1542 /* here we fill the entry of N at [r, d], for each valid vector
1543 that we already have in N;
1544 again, this will only be done when the given linear system has
1545 more than one, i.e., infinitely many solutions */
1546 p = NULL;
1547 for (int c = nonZeroC + 1; c <= n; c++)
1548 if (MATELEM(N, c, d) != NULL)
1549 p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(N, c, d)));
1550 /* The following division may be with remainder but only takes place
1551 when dim > 0. */
1552 q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
1553 MATELEM(N, nonZeroC, d) = pMult(pNeg(p), q);
1554 pNormalize(MATELEM(N, nonZeroC, d));
1555 }
1556 p = pNeg(pCopy(MATELEM(yVec, r, 1)));
1557 for (int c = nonZeroC + 1; c <= n; c++)
1558 if (MATELEM(xVec, c, 1) != NULL)
1559 p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(xVec, c, 1)));
1560 /* The following division is without remainder. */
1561 q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
1562 MATELEM(xVec, nonZeroC, 1) = pMult(pNeg(p), q);
1563 pNormalize(MATELEM(xVec, nonZeroC, 1));
1564 lastNonZeroC = nonZeroC;
1565 }
1566
1567 /* divide xVec by l*u = lTimesU and put result in xVec */
1568 number z = nInvers(pGetCoeff(lTimesU));
1569 for (int c = 1; c <= n; c++)
1570 {
1571 MATELEM(xVec, c, 1)=__p_Mult_nn(MATELEM(xVec, c, 1), z,currRing);
1572 pNormalize(MATELEM(xVec, c, 1));
1573 }
1574 nDelete(&z);
1575
1576 if (dim == 0)
1577 {
1578 /* that means the given linear system has exactly one solution;
1579 we return as H the 1x1 matrix with entry zero */
1580 H = mpNew(1, 1);
1581 }
1582 else
1583 {
1584 /* copy the first 'dim' columns of N into H */
1585 H = mpNew(n, dim);
1586 for (int r = 1; r <= n; r++)
1587 for (int c = 1; c <= dim; c++)
1588 MATELEM(H, r, c) = pCopy(MATELEM(N, r, c));
1589 }
1590 /* clean up N */
1591 idDelete((ideal*)&N);
1592 }
1593
1594 idDelete((ideal*)&cVec);
1595 idDelete((ideal*)&yVec);
1596
1597 return isSolvable;
1598 }
1599
1600