source: git/kernel/linearAlgebra.cc @ 9e269a0

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