source:git/kernel/linearAlgebra.cc@6ed8c4

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