# source:git/kernel/linearAlgebra.cc@51386e9

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