My Project  UNKNOWN_GIT_VERSION
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  * 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  **/
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  **/
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 
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 
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  **/
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 */
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 
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 
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 
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  **/
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 
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 */
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 */
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  **/
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  **/
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 */
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 
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 
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 
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)));
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. */
725 number absValue(poly p)
726 {
727  if (p == NULL) return nInit(0);
728  number result = nCopy(pGetCoeff(p));
730  return result;
731 }
732 
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 
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 
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 
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 
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  **/
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 
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 algorith.
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  **/
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 */
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*)&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. */
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*)&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 
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 
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 
dim
int dim(ideal I, ring r)
Definition: tropicalStrategy.cc:22
exponent
int exponent(const CanonicalForm &f, int q)
int exponent ( const CanonicalForm & f, int q )
Definition: gengftables-conway.cc:92
rField_is_long_R
static BOOLEAN rField_is_long_R(const ring r)
Definition: ring.h:532
mpr_complex.h
p_GetCoeff
#define p_GetCoeff(p, r)
Definition: monomials.h:48
ip_smatrix
Definition: matpol.h:13
nNormalize
#define nNormalize(n)
Definition: numbers.h:30
j
int j
Definition: facHensel.cc:105
f
FILE * f
Definition: checklibs.c:9
p_Normalize
void p_Normalize(poly p, const ring r)
Definition: p_polys.cc:3720
k
int k
Definition: cfEzgcd.cc:92
idDelete
#define idDelete(H)
delete an ideal
Definition: ideals.h:28
MATELEM
#define MATELEM(mat, i, j)
Definition: matpol.h:28
result
return result
Definition: facAbsBiFact.cc:76
luInverse
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...
Definition: linearAlgebra.cc:200
nAdd
#define nAdd(n1, n2)
Definition: numbers.h:18
unitMatrix
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.
Definition: linearAlgebra.cc:95
nGreater
#define nGreater(a, b)
Definition: numbers.h:28
pGetExp
#define pGetExp(p, i)
Exponent.
Definition: polys.h:40
realSqrt
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.
Definition: linearAlgebra.cc:601
qrDS
bool qrDS(const int, matrix *queue, int &queueL, number *eigenValues, int &eigenValuesL, const number tol1, const number tol2, const ring R)
Definition: linearAlgebra.cc:1090
p_Neg
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1029
simpleideals.h
ip_smatrix::cols
int & cols()
Definition: matpol.h:24
luSolveViaLUDecomp
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...
Definition: linearAlgebra.cc:377
g
g
Definition: cfModGcd.cc:4031
n_Delete
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:454
pNeg
#define pNeg(p)
Definition: polys.h:185
pDelete
#define pDelete(p_ptr)
Definition: polys.h:174
N
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:48
hessenberg
void hessenberg(const matrix aMat, matrix &pMat, matrix &hessenbergMat, const number tolerance, const ring R)
Computes the Hessenberg form of a given square matrix.
Definition: linearAlgebra.cc:909
pMult
#define pMult(p, q)
Definition: polys.h:194
quadraticSolve
int quadraticSolve(const poly p, number &s1, number &s2, const number tolerance)
Returns all solutions of a quadratic polynomial relation with real coefficients.
Definition: linearAlgebra.cc:626
w
const CanonicalForm & w
Definition: facAbsFact.cc:55
pString
char * pString(poly p)
Definition: polys.h:288
b
CanonicalForm b
Definition: cfModGcd.cc:4044
__p_Mult_nn
#define __p_Mult_nn(p, n, r)
Definition: p_polys.h:913
charPoly
bool charPoly(const matrix aMat, poly &charPoly)
Computes the characteristic polynomial from a quadratic (2x2) matrix and returns true in case of succ...
Definition: linearAlgebra.cc:748
n_Normalize
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition: coeffs.h:577
ppMult_qq
#define ppMult_qq(p, q)
Definition: polys.h:195
swapColumns
void swapColumns(int column1, int column2, matrix &aMat)
Swaps two columns of a given matrix and thereby modifies it.
Definition: linearAlgebra.cc:793
nGreaterZero
#define nGreaterZero(n)
Definition: numbers.h:27
similar
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].
Definition: linearAlgebra.cc:1188
nInpNeg
#define nInpNeg(n)
Definition: numbers.h:21
p_Copy
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:798
currRing
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
hessenbergStep
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...
Definition: linearAlgebra.cc:837
i
int i
Definition: cfEzgcd.cc:125
id_Delete
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
Definition: simpleideals.cc:113
ip_smatrix::rows
int & rows()
Definition: matpol.h:23
nIsOne
#define nIsOne(n)
Definition: numbers.h:25
complexNumber
number complexNumber(const double r, const double i)
Creates a new complex number from real and imaginary parts given by doubles.
Definition: linearAlgebra.cc:541
block
#define block
Definition: scanner.cc:664
matpol.h
tenToTheMinus
number tenToTheMinus(const int exponent)
Returns one over the exponent-th power of ten.
Definition: linearAlgebra.cc:554
lowerLeftTriangleInverse
bool lowerLeftTriangleInverse(const matrix lMat, matrix &iMat, bool diagonalIsOne)
Computes the inverse of a given (n x n)-matrix in lower left triangular form.
Definition: linearAlgebra.cc:300
luSolveViaLDUDecomp
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...
Definition: linearAlgebra.cc:1461
euclideanNormSquared
number euclideanNormSquared(const matrix aMat)
Definition: linearAlgebra.cc:705
printNumber
void printNumber(const number z)
Definition: linearAlgebra.cc:506
printMatrix
void printMatrix(const matrix m)
Definition: linearAlgebra.cc:522
lduDecomp
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...
Definition: linearAlgebra.cc:1343
structs.h
henselFactors
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,...
Definition: linearAlgebra.cc:1219
h
static Poly * h
Definition: janet.cc:972
matrixBlock
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 ...
Definition: linearAlgebra.cc:805
mod2.h
pOne
#define pOne()
Definition: polys.h:297
pp_Mult_qq
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition: p_polys.h:1073
pivot
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].
Definition: linearAlgebra.cc:68
p_polys.h
nDiv
#define nDiv(a, b)
Definition: numbers.h:32
rField_is_R
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:508
pNSet
#define pNSet(n)
Definition: polys.h:295
__pp_Mult_nn
#define __pp_Mult_nn(p, n, r)
Definition: p_polys.h:944
tmp1
CFList tmp1
Definition: facFqBivar.cc:70
n_InpNeg
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:556
mpNew
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:36
pAdd
#define pAdd(p, q)
Definition: polys.h:190
luInverseFromLUDecomp
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...
Definition: linearAlgebra.cc:352
pSetCoeff
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:30
mp_Mult
matrix mp_Mult(matrix a, matrix b, const ring R)
Definition: matpol.cc:212
nIsZero
#define nIsZero(n)
Definition: numbers.h:19
p_Deg
long p_Deg(poly a, const ring r)
Definition: p_polys.cc:578
nMult
#define nMult(n1, n2)
Definition: numbers.h:17
p_Delete
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:843
p_Add_q
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:878
p_One
poly p_One(const ring r)
Definition: p_polys.cc:1302
n_Invers
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:563
pivotScore
int pivotScore(number n, const ring r)
The returned score is based on the implementation of 'nSize' for numbers (, see numbers....
Definition: linearAlgebra.cc:50
nSub
#define nSub(n1, n2)
Definition: numbers.h:22
nInvers
#define nInvers(a)
Definition: numbers.h:33
p_NSet
poly p_NSet(number n, const ring r)
returns the poly representing the number n, destroys n
Definition: p_polys.cc:1432
n_Copy
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition: coeffs.h:450
n_Gcd
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:685
rankFromRowEchelonForm
int rankFromRowEchelonForm(const matrix aMat)
Definition: linearAlgebra.cc:217
m
int m
Definition: cfEzgcd.cc:121
MATCOLS
#define MATCOLS(i)
Definition: matpol.h:27
assume
#define assume(x)
Definition: mod2.h:384
NULL
#define NULL
Definition: omList.c:9
pLmDelete
#define pLmDelete(p)
assume p != NULL, deletes Lm(p)->coef and Lm(p)
Definition: polys.h:74
swapRows
void swapRows(int row1, int row2, matrix &aMat)
Swaps two rows of a given matrix and thereby modifies it.
Definition: linearAlgebra.cc:781
pSetm
#define pSetm(p)
Definition: polys.h:254
ideals.h
l
int l
Definition: cfEzgcd.cc:93
nDelete
#define nDelete(n)
Definition: numbers.h:16
R
#define R
Definition: sirandom.c:26
printSolutions
void printSolutions(const int a, const int b, const int c)
Definition: linearAlgebra.cc:575
pp_Mult_nn
static poly pp_Mult_nn(poly p, number n, const ring r)
Definition: p_polys.h:934
luRank
int luRank(const matrix aMat, const bool isRowEchelon, const ring R)
Computes the rank of a given (m x n)-matrix.
Definition: linearAlgebra.cc:230
pSetExp
#define pSetExp(p, i, v)
Definition: polys.h:41
absValue
number absValue(poly p)
Definition: linearAlgebra.cc:725
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
p
int p
Definition: cfModGcd.cc:4019
s
const CanonicalForm int s
Definition: facAbsFact.cc:55
subMatrix
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.
Definition: linearAlgebra.cc:733
nInit
#define nInit(i)
Definition: numbers.h:24
H
CanonicalForm H
Definition: facAbsFact.cc:64
upperRightTriangleInverse
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.
Definition: linearAlgebra.cc:251
n_Div
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:614
pCopy
#define pCopy(p)
return a copy of the poly
Definition: polys.h:173
pNormalize
#define pNormalize(p)
Definition: polys.h:299
p_Mult_q
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1036
pGetCoeff
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:42
tmp2
CFList tmp2
Definition: facFqBivar.cc:70
n_Size
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:569
mp_Copy
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition: matpol.cc:63
rField_is_long_C
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:535
MATROWS
#define MATROWS(i)
Definition: matpol.h:26
nCopy
#define nCopy(n)
Definition: numbers.h:15
numbers.h
pNext
#define pNext(p)
Definition: monomials.h:34
degree
int degree(const CanonicalForm &f)
Definition: canonicalform.h:309
luDecomp
void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat, const ring R)
LU-decomposition of a given (m x n)-matrix.
Definition: linearAlgebra.cc:103
gmp_complex
gmp_complex numbers based on
Definition: mpr_complex.h:177
linearAlgebra.h
coeffs.h
mpTrafo
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...
Definition: linearAlgebra.cc:979