source: git/kernel/linearAlgebra.cc @ fbc7cb

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