My Project
Loading...
Searching...
No Matches
linearAlgebra.cc
Go to the documentation of this file.
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"
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"
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 * preferred, 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 **/
50int pivotScore(number n, const ring r)
51{
52 int s = n_Size(n, r->cf);
53 if (rField_is_long_C(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 **/
68bool 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
95bool 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
103void 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 **/
200bool 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 */
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
230int 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
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
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 **/
352bool 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
377bool 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 */
506void 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 */
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 **/
541number 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 **/
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 */
575void 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
601bool 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
626int 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
705number 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)));
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. */
725number absValue(poly p)
726{
727 if (p == NULL) return nInit(0);
728 number result = nCopy(pGetCoeff(p));
730 return result;
731}
732
733bool 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
748bool 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
781void 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
793void 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
805void 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 **/
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
909void 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 **/
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 */
1090bool 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*)&currentMat); 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*)&currentMat);
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 **/
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. */
1219void 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*)&notUsedMat);
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
1343void 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
1461bool 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
int degree(const CanonicalForm &f)
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
int p
Definition: cfModGcd.cc:4078
g
Definition: cfModGcd.cc:4090
CanonicalForm b
Definition: cfModGcd.cc:4103
FILE * f
Definition: checklibs.c:9
gmp_complex numbers based on
Definition: mpr_complex.h:179
int & cols()
Definition: matpol.h:24
int & rows()
Definition: matpol.h:23
Coefficient rings, fields and other domains suitable for Singular polynomials.
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition: coeffs.h:448
static FORCE_INLINE number n_Gcd(number a, number b, const coeffs r)
in Z: return the gcd of 'a' and 'b' in Z/nZ, Z/2^kZ: computed as in the case Z in Z/pZ,...
Definition: coeffs.h:661
static FORCE_INLINE number n_Invers(number a, const coeffs r)
return the multiplicative inverse of 'a'; raise an error if 'a' is not invertible
Definition: coeffs.h:561
static FORCE_INLINE number n_InpNeg(number n, const coeffs r)
in-place negation of n MUST BE USED: n = n_InpNeg(n) (no copy is returned)
Definition: coeffs.h:554
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
Definition: coeffs.h:612
static FORCE_INLINE int n_Size(number n, const coeffs r)
return a non-negative measure for the complexity of n; return 0 only when n represents zero; (used fo...
Definition: coeffs.h:567
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:452
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition: coeffs.h:575
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int s
Definition: facAbsFact.cc:51
const CanonicalForm & w
Definition: facAbsFact.cc:51
CanonicalForm H
Definition: facAbsFact.cc:60
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
CFList tmp1
Definition: facFqBivar.cc:72
CFList tmp2
Definition: facFqBivar.cc:72
int j
Definition: facHensel.cc:110
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
#define exponent
STATIC_VAR Poly * h
Definition: janet.cc:971
bool 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)
Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LDU-decomposit...
number absValue(poly p)
void hessenberg(const matrix aMat, matrix &pMat, matrix &hessenbergMat, const number tolerance, const ring R)
Computes the Hessenberg form of a given square matrix.
number hessenbergStep(const matrix vVec, matrix &uVec, matrix &pMat, const number tolerance)
Computes information related to one householder transformation step for constructing the Hessenberg f...
bool realSqrt(const number n, const number tolerance, number &root)
Computes the square root of a non-negative real number and returns it as a new number.
void mpTrafo(matrix &H, int it, const number tolerance, const ring R)
Performs one transformation step on the given matrix H as part of the gouverning QR double shift algo...
number complexNumber(const double r, const double i)
Creates a new complex number from real and imaginary parts given by doubles.
void swapColumns(int column1, int column2, matrix &aMat)
Swaps two columns of a given matrix and thereby modifies it.
void printNumber(const number z)
bool qrDS(const int, matrix *queue, int &queueL, number *eigenValues, int &eigenValuesL, const number tol1, const number tol2, const ring R)
int luRank(const matrix aMat, const bool isRowEchelon, const ring R)
Computes the rank of a given (m x n)-matrix.
int rankFromRowEchelonForm(const matrix aMat)
int quadraticSolve(const poly p, number &s1, number &s2, const number tolerance)
Returns all solutions of a quadratic polynomial relation with real coefficients.
void printMatrix(const matrix m)
bool luInverseFromLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, matrix &iMat, const ring R)
This code computes the inverse by inverting lMat and uMat, and then performing two matrix multiplicat...
void matrixBlock(const matrix aMat, const matrix bMat, matrix &block)
Creates a new matrix which contains the first argument in the top left corner, and the second in the ...
void henselFactors(const int xIndex, const int yIndex, const poly h, const poly f0, const poly g0, const int d, poly &f, poly &g)
Computes a factorization of a polynomial h(x, y) in K[[x]][y] up to a certain degree in x,...
bool upperRightTriangleInverse(const matrix uMat, matrix &iMat, bool diagonalIsOne, const ring R)
Computes the inverse of a given (n x n)-matrix in upper right triangular form.
bool unitMatrix(const int n, matrix &unitMat, const ring R)
Creates a new matrix which is the (nxn) unit matrix, and returns true in case of success.
int pivotScore(number n, const ring r)
The returned score is based on the implementation of 'nSize' for numbers (, see numbers....
bool pivot(const matrix aMat, const int r1, const int r2, const int c1, const int c2, int *bestR, int *bestC, const ring R)
This code computes a score for each non-zero matrix entry in aMat[r1..r2, c1..c2].
int similar(const number *nn, const int nnLength, const number n, const number tolerance)
Tries to find the number n in the array nn[0..nnLength-1].
bool charPoly(const matrix aMat, poly &charPoly)
Computes the characteristic polynomial from a quadratic (2x2) matrix and returns true in case of succ...
void lduDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &dMat, matrix &uMat, poly &l, poly &u, poly &lTimesU)
LU-decomposition of a given (m x n)-matrix with performing only those divisions that yield zero remai...
bool luInverse(const matrix aMat, matrix &iMat, const ring R)
This code first computes the LU-decomposition of aMat, and then calls the method for inverting a matr...
void swapRows(int row1, int row2, matrix &aMat)
Swaps two rows of a given matrix and thereby modifies it.
number tenToTheMinus(const int exponent)
Returns one over the exponent-th power of ten.
bool lowerLeftTriangleInverse(const matrix lMat, matrix &iMat, bool diagonalIsOne)
Computes the inverse of a given (n x n)-matrix in lower left triangular form.
number euclideanNormSquared(const matrix aMat)
bool subMatrix(const matrix aMat, const int rowIndex1, const int rowIndex2, const int colIndex1, const int colIndex2, matrix &subMat)
Creates a new matrix which is a submatrix of the first argument, and returns true in case of success.
void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat, const ring R)
LU-decomposition of a given (m x n)-matrix.
void printSolutions(const int a, const int b, const int c)
bool luSolveViaLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, const matrix bVec, matrix &xVec, matrix &H)
Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LU-decompositi...
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:37
matrix mp_Mult(matrix a, matrix b, const ring R)
Definition: matpol.cc:206
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition: matpol.cc:57
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
#define MATROWS(i)
Definition: matpol.h:26
#define MATCOLS(i)
Definition: matpol.h:27
#define assume(x)
Definition: mod2.h:389
#define pNext(p)
Definition: monomials.h:36
#define p_GetCoeff(p, r)
Definition: monomials.h:50
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
#define nDiv(a, b)
Definition: numbers.h:32
#define nDelete(n)
Definition: numbers.h:16
#define nInpNeg(n)
Definition: numbers.h:21
#define nIsZero(n)
Definition: numbers.h:19
#define nSub(n1, n2)
Definition: numbers.h:22
#define nCopy(n)
Definition: numbers.h:15
#define nGreater(a, b)
Definition: numbers.h:28
#define nAdd(n1, n2)
Definition: numbers.h:18
#define nGreaterZero(n)
Definition: numbers.h:27
#define nInvers(a)
Definition: numbers.h:33
#define nIsOne(n)
Definition: numbers.h:25
#define nNormalize(n)
Definition: numbers.h:30
#define nInit(i)
Definition: numbers.h:24
#define nMult(n1, n2)
Definition: numbers.h:17
#define NULL
Definition: omList.c:12
void p_Normalize(poly p, const ring r)
Definition: p_polys.cc:3813
poly p_One(const ring r)
Definition: p_polys.cc:1313
poly p_NSet(number n, const ring r)
returns the poly representing the number n, destroys n
Definition: p_polys.cc:1473
long p_Deg(poly a, const ring r)
Definition: p_polys.cc:587
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1105
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:934
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1112
#define __pp_Mult_nn(p, n, r)
Definition: p_polys.h:1000
static poly pp_Mult_nn(poly p, number n, const ring r)
Definition: p_polys.h:990
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:899
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition: p_polys.h:1149
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:844
#define __p_Mult_nn(p, n, r)
Definition: p_polys.h:969
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
#define pAdd(p, q)
Definition: polys.h:203
#define pDelete(p_ptr)
Definition: polys.h:186
#define pSetm(p)
Definition: polys.h:271
#define pNeg(p)
Definition: polys.h:198
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
#define pNSet(n)
Definition: polys.h:313
#define ppMult_qq(p, q)
Definition: polys.h:208
#define pLmDelete(p)
assume p != NULL, deletes Lm(p)->coef and Lm(p)
Definition: polys.h:76
#define pMult(p, q)
Definition: polys.h:207
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
#define pNormalize(p)
Definition: polys.h:317
#define pSetExp(p, i, v)
Definition: polys.h:42
char * pString(poly p)
Definition: polys.h:306
#define pCopy(p)
return a copy of the poly
Definition: polys.h:185
#define pOne()
Definition: polys.h:315
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:518
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:545
static BOOLEAN rField_is_long_R(const ring r)
Definition: ring.h:542
#define block
Definition: scanner.cc:646
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
#define R
Definition: sirandom.c:27
int dim(ideal I, ring r)