My Project
Functions
linearAlgebra.cc File Reference
#include "kernel/mod2.h"
#include "coeffs/coeffs.h"
#include "coeffs/numbers.h"
#include "coeffs/mpr_complex.h"
#include "polys/monomials/p_polys.h"
#include "polys/simpleideals.h"
#include "polys/matpol.h"
#include "kernel/structs.h"
#include "kernel/ideals.h"
#include "kernel/linear_algebra/linearAlgebra.h"

Go to the source code of this file.

Functions

int pivotScore (number n, const ring r)
 The returned score is based on the implementation of 'nSize' for numbers (, see numbers.h): nSize(n) provides a measure for the complexity of n. More...
 
bool pivot (const matrix aMat, const int r1, const int r2, const int c1, const int c2, int *bestR, int *bestC, const ring R)
 This code computes a score for each non-zero matrix entry in aMat[r1..r2, c1..c2]. More...
 
bool unitMatrix (const int n, matrix &unitMat, const ring R)
 Creates a new matrix which is the (nxn) unit matrix, and returns true in case of success. More...
 
void luDecomp (const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat, const ring R)
 LU-decomposition of a given (m x n)-matrix. More...
 
bool luInverse (const matrix aMat, matrix &iMat, const ring R)
 This code first computes the LU-decomposition of aMat, and then calls the method for inverting a matrix which is given by its LU-decomposition. More...
 
int rankFromRowEchelonForm (const matrix aMat)
 
int luRank (const matrix aMat, const bool isRowEchelon, const ring R)
 Computes the rank of a given (m x n)-matrix. More...
 
bool upperRightTriangleInverse (const matrix uMat, matrix &iMat, bool diagonalIsOne, const ring R)
 Computes the inverse of a given (n x n)-matrix in upper right triangular form. More...
 
bool lowerLeftTriangleInverse (const matrix lMat, matrix &iMat, bool diagonalIsOne)
 Computes the inverse of a given (n x n)-matrix in lower left triangular form. More...
 
bool luInverseFromLUDecomp (const matrix pMat, const matrix lMat, const matrix uMat, matrix &iMat, const ring R)
 This code computes the inverse by inverting lMat and uMat, and then performing two matrix multiplications. More...
 
bool luSolveViaLUDecomp (const matrix pMat, const matrix lMat, const matrix uMat, const matrix bVec, matrix &xVec, matrix &H)
 Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LU-decomposition. More...
 
void printNumber (const number z)
 
void printMatrix (const matrix m)
 
number complexNumber (const double r, const double i)
 Creates a new complex number from real and imaginary parts given by doubles. More...
 
number tenToTheMinus (const int exponent)
 Returns one over the exponent-th power of ten. More...
 
void printSolutions (const int a, const int b, const int c)
 
bool realSqrt (const number n, const number tolerance, number &root)
 Computes the square root of a non-negative real number and returns it as a new number. More...
 
int quadraticSolve (const poly p, number &s1, number &s2, const number tolerance)
 Returns all solutions of a quadratic polynomial relation with real coefficients. More...
 
number euclideanNormSquared (const matrix aMat)
 
number absValue (poly p)
 
bool subMatrix (const matrix aMat, const int rowIndex1, const int rowIndex2, const int colIndex1, const int colIndex2, matrix &subMat)
 Creates a new matrix which is a submatrix of the first argument, and returns true in case of success. More...
 
bool charPoly (const matrix aMat, poly &charPoly)
 Computes the characteristic polynomial from a quadratic (2x2) matrix and returns true in case of success. More...
 
void swapRows (int row1, int row2, matrix &aMat)
 Swaps two rows of a given matrix and thereby modifies it. More...
 
void swapColumns (int column1, int column2, matrix &aMat)
 Swaps two columns of a given matrix and thereby modifies it. More...
 
void matrixBlock (const matrix aMat, const matrix bMat, matrix &block)
 Creates a new matrix which contains the first argument in the top left corner, and the second in the bottom right. More...
 
number hessenbergStep (const matrix vVec, matrix &uVec, matrix &pMat, const number tolerance)
 Computes information related to one householder transformation step for constructing the Hessenberg form of a given non-derogative matrix. More...
 
void hessenberg (const matrix aMat, matrix &pMat, matrix &hessenbergMat, const number tolerance, const ring R)
 Computes the Hessenberg form of a given square matrix. More...
 
void mpTrafo (matrix &H, int it, const number tolerance, const ring R)
 Performs one transformation step on the given matrix H as part of the gouverning QR double shift algorithm. More...
 
bool qrDS (const int, matrix *queue, int &queueL, number *eigenValues, int &eigenValuesL, const number tol1, const number tol2, const ring R)
 
int similar (const number *nn, const int nnLength, const number n, const number tolerance)
 Tries to find the number n in the array nn[0..nnLength-1]. More...
 
void henselFactors (const int xIndex, const int yIndex, const poly h, const poly f0, const poly g0, const int d, poly &f, poly &g)
 Computes a factorization of a polynomial h(x, y) in K[[x]][y] up to a certain degree in x, whenever a factorization of h(0, y) is given. More...
 
void lduDecomp (const matrix aMat, matrix &pMat, matrix &lMat, matrix &dMat, matrix &uMat, poly &l, poly &u, poly &lTimesU)
 LU-decomposition of a given (m x n)-matrix with performing only those divisions that yield zero remainders. More...
 
bool luSolveViaLDUDecomp (const matrix pMat, const matrix lMat, const matrix dMat, const matrix uMat, const poly l, const poly u, const poly lTimesU, const matrix bVec, matrix &xVec, matrix &H)
 Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LDU-decomposition. More...
 

Function Documentation

◆ absValue()

number absValue ( poly  p)

Definition at line 725 of file linearAlgebra.cc.

726 {
727  if (p == NULL) return nInit(0);
728  number result = nCopy(pGetCoeff(p));
730  return result;
731 }
int p
Definition: cfModGcd.cc:4078
return result
Definition: facAbsBiFact.cc:75
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
#define nInpNeg(n)
Definition: numbers.h:21
#define nCopy(n)
Definition: numbers.h:15
#define nGreaterZero(n)
Definition: numbers.h:27
#define nInit(i)
Definition: numbers.h:24
#define NULL
Definition: omList.c:12

◆ charPoly()

bool charPoly ( const matrix  aMat,
poly &  charPoly 
)

Computes the characteristic polynomial from a quadratic (2x2) matrix and returns true in case of success.

The method will be successful whenever the input matrix is a (2x2) matrix. In this case, the resulting polynomial will be a univariate polynomial in the first ring variable of degree 2 and it will reside in the second argument. The method assumes that the matrix entries are all numbers, i.e., elements from the ground field/ring.

Returns
true iff the input matrix was (2x2)
Parameters
[in]aMatthe input matrix
[out]charPolythe characteristic polynomial

Definition at line 748 of file linearAlgebra.cc.

749 {
750  if (MATROWS(aMat) != 2) return false;
751  if (MATCOLS(aMat) != 2) return false;
752  number b = nInit(0); number t;
753  if (MATELEM(aMat, 1, 1) != NULL)
754  { t = nAdd(b, pGetCoeff(MATELEM(aMat, 1, 1))); nDelete(&b); b = t;}
755  if (MATELEM(aMat, 2, 2) != NULL)
756  { t = nAdd(b, pGetCoeff(MATELEM(aMat, 2, 2))); nDelete(&b); b = t;}
757  b = nInpNeg(b);
758  number t1;
759  if ((MATELEM(aMat, 1, 1) != NULL) && (MATELEM(aMat, 2, 2) != NULL))
760  t1 = nMult(pGetCoeff(MATELEM(aMat, 1, 1)),
761  pGetCoeff(MATELEM(aMat, 2, 2)));
762  else t1 = nInit(0);
763  number t2;
764  if ((MATELEM(aMat, 1, 2) != NULL) && (MATELEM(aMat, 2, 1) != NULL))
765  t2 = nMult(pGetCoeff(MATELEM(aMat, 1, 2)),
766  pGetCoeff(MATELEM(aMat, 2, 1)));
767  else t2 = nInit(0);
768  number c = nSub(t1, t2); nDelete(&t1); nDelete(&t2);
769  poly p = pOne(); pSetExp(p, 1, 2); pSetm(p);
770  poly q = NULL;
771  if (!nIsZero(b))
772  { q = pOne(); pSetExp(q, 1, 1); pSetm(q); pSetCoeff(q, b); }
773  poly r = NULL;
774  if (!nIsZero(c))
775  { r = pOne(); pSetCoeff(r, c); }
776  p = pAdd(p, q); p = pAdd(p, r);
777  charPoly = p;
778  return true;
779 }
CanonicalForm b
Definition: cfModGcd.cc:4103
bool charPoly(const matrix aMat, poly &charPoly)
Computes the characteristic polynomial from a quadratic (2x2) matrix and returns true in case of succ...
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
#define MATROWS(i)
Definition: matpol.h:26
#define MATCOLS(i)
Definition: matpol.h:27
#define nDelete(n)
Definition: numbers.h:16
#define nIsZero(n)
Definition: numbers.h:19
#define nSub(n1, n2)
Definition: numbers.h:22
#define nAdd(n1, n2)
Definition: numbers.h:18
#define nMult(n1, n2)
Definition: numbers.h:17
#define pAdd(p, q)
Definition: polys.h:203
#define pSetm(p)
Definition: polys.h:271
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
#define pSetExp(p, i, v)
Definition: polys.h:42
#define pOne()
Definition: polys.h:315

◆ complexNumber()

number complexNumber ( const double  r,
const double  i 
)

Creates a new complex number from real and imaginary parts given by doubles.

Returns
the new complex number

Definition at line 541 of file linearAlgebra.cc.

542 {
543  gmp_complex* n= new gmp_complex(r, i);
544  return (number)n;
545 }
int i
Definition: cfEzgcd.cc:132
gmp_complex numbers based on
Definition: mpr_complex.h:179

◆ euclideanNormSquared()

number euclideanNormSquared ( const matrix  aMat)

Definition at line 705 of file linearAlgebra.cc.

706 {
707  int rr = MATROWS(aMat);
708  number result = nInit(0);
709  number tmp1; number tmp2;
710  for (int r = 1; r <= rr; r++)
711  if (MATELEM(aMat, r, 1) != NULL)
712  {
713  tmp1 = nMult(pGetCoeff(MATELEM(aMat, r, 1)),
714  pGetCoeff(MATELEM(aMat, r, 1)));
716  result = tmp2;
717  }
718  return result;
719 }
CFList tmp1
Definition: facFqBivar.cc:72
CFList tmp2
Definition: facFqBivar.cc:72

◆ henselFactors()

void henselFactors ( const int  xIndex,
const int  yIndex,
const poly  h,
const poly  f0,
const poly  g0,
const int  d,
poly &  f,
poly &  g 
)

Computes a factorization of a polynomial h(x, y) in K[[x]][y] up to a certain degree in x, whenever a factorization of h(0, y) is given.

The algorithm is based on Hensel's lemma: Let h(x, y) denote a monic polynomial in y of degree m + n with coefficients in K[[x]]. Suppose there are two monic factors f_0(y) (of degree n) and g_0(y) of degree (m) such that h(0, y) = f_0(y) * g_0(y) and <f_0, g_0> = K[y]. Fix an integer d >= 0. Then there are monic polynomials in y with coefficients in K[[x]], namely f(x, y) of degree n and g(x, y) of degree m such that h(x, y) = f(x, y) * g(x, y) modulo <x^(d+1)> (*).

This implementation expects h, f0, g0, and d as described and computes the factors f and g. Effectively, h will be given as an element of K[x, y] since all terms of h with x-degree larger than d can be ignored due to (*). The method expects the ground ring to contain at least two variables; then x is the first ring variable, specified by the input index xIndex, and y the second one, specified by yIndex.

This code was placed here since the algorithm works by successively solving d linear equation systems. It is hence an application of other methods defined in this h-file and its corresponding cc-file.

Parameters
[in]xIndexindex of x in {1, ..., nvars(basering)}
[in]yIndexindex of y in {1, ..., nvars(basering)}
[in]hthe polynomial h(x, y) as above
[in]f0the first univariate factor of h(0, y)
[in]g0the second univariate factor of h(0, y)
[in]dthe degree bound, d >= 0
[out]fthe first factor of h(x, y)
[out]gthe second factor of h(x, y)

Definition at line 1219 of file linearAlgebra.cc.

1221 {
1222  int n = (int)p_Deg(f0,currRing);
1223  int m = (int)p_Deg(g0,currRing);
1224  matrix aMat = mpNew(n + m, n + m); /* matrix A for linear system */
1225  matrix pMat; matrix lMat; matrix uMat; /* for the decomposition of A */
1226  f = pCopy(f0); g = pCopy(g0); /* initially: h = f*g mod <x^1> */
1227 
1228  /* initial step: read off coefficients of f0, and g0 */
1229  poly p = f0; poly matEntry; number c;
1230  while (p != NULL)
1231  {
1232  c = nCopy(pGetCoeff(p));
1233  matEntry = pOne(); pSetCoeff(matEntry, c);
1234  MATELEM(aMat, pGetExp(p, yIndex) + 1, 1) = matEntry;
1235  p = pNext(p);
1236  }
1237  p = g0;
1238  while (p != NULL)
1239  {
1240  c = nCopy(pGetCoeff(p));
1241  matEntry = pOne(); pSetCoeff(matEntry, c);
1242  MATELEM(aMat, pGetExp(p, yIndex) + 1, m + 1) = matEntry;
1243  p = pNext(p);
1244  }
1245  /* fill the rest of A */
1246  for (int row = 2; row <= n + 1; row++)
1247  for (int col = 2; col <= m; col++)
1248  {
1249  if (col > row) break;
1250  MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1));
1251  }
1252  for (int row = n + 2; row <= n + m; row++)
1253  for (int col = row - n; col <= m; col++)
1254  MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1));
1255  for (int row = 2; row <= m + 1; row++)
1256  for (int col = m + 2; col <= m + n; col++)
1257  {
1258  if (col - m > row) break;
1259  MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1));
1260  }
1261  for (int row = m + 2; row <= n + m; row++)
1262  for (int col = row; col <= m + n; col++)
1263  MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1));
1264 
1265  /* constructing the LU-decomposition of A */
1266  luDecomp(aMat, pMat, lMat, uMat);
1267 
1268  /* Before the xExp-th loop, we know that h = f*g mod <x^xExp>.
1269  Afterwards the algorithm ensures h = f*g mod <x^(xExp + 1)>.
1270  Hence in the end we obtain f and g as required, i.e.,
1271  h = f*g mod <x^(d+1)>.
1272  The algorithm works by solving a (m+n)x(m+n) linear equation system
1273  A*x = b with constant matrix A (as decomposed above). By theory, the
1274  system is guaranteed to have a unique solution. */
1275  poly fg = ppMult_qq(f, g); /* for storing the product of f and g */
1276  for (int xExp = 1; xExp <= d; xExp++)
1277  {
1278  matrix bVec = mpNew(n + m, 1); /* b */
1279  matrix xVec = mpNew(n + m, 1); /* x */
1280 
1281  p = pCopy(fg);
1282  p = pAdd(pCopy(h), pNeg(p)); /* p = h - f*g */
1283  /* we collect all terms in p with x-exponent = xExp and use their
1284  coefficients to build the vector b */
1285  int bIsZeroVector = true;
1286  while (p != NULL)
1287  {
1288  if (pGetExp(p, xIndex) == xExp)
1289  {
1290  number c = nCopy(pGetCoeff(p));
1291  poly matEntry = pOne(); pSetCoeff(matEntry, c);
1292  MATELEM(bVec, pGetExp(p, yIndex) + 1, 1) = matEntry;
1293  if (matEntry != NULL) bIsZeroVector = false;
1294  }
1295  pLmDelete(&p); /* destruct leading term of p and move to next term */
1296  }
1297  /* solve the linear equation system */
1298  if (!bIsZeroVector) /* otherwise x = 0 and f, g do not change */
1299  {
1300  matrix notUsedMat;
1301  luSolveViaLUDecomp(pMat, lMat, uMat, bVec, xVec, notUsedMat);
1302  idDelete((ideal*)&notUsedMat);
1303  /* update f and g by newly computed terms, and update f*g */
1304  poly fNew = NULL; poly gNew = NULL;
1305  for (int row = 1; row <= m; row++)
1306  {
1307  if (MATELEM(xVec, row, 1) != NULL)
1308  {
1309  p = pCopy(MATELEM(xVec, row, 1)); /* p = c */
1310  pSetExp(p, xIndex, xExp); /* p = c * x^xExp */
1311  pSetExp(p, yIndex, row - 1); /* p = c * x^xExp * y^i */
1312  pSetm(p);
1313  gNew = pAdd(gNew, p);
1314  }
1315  }
1316  for (int row = m + 1; row <= m + n; row++)
1317  {
1318  if (MATELEM(xVec, row, 1) != NULL)
1319  {
1320  p = pCopy(MATELEM(xVec, row, 1)); /* p = c */
1321  pSetExp(p, xIndex, xExp); /* p = c * x^xExp */
1322  pSetExp(p, yIndex, row - m - 1); /* p = c * x^xExp * y^i */
1323  pSetm(p);
1324  fNew = pAdd(fNew, p);
1325  }
1326  }
1327  fg = pAdd(fg, ppMult_qq(f, gNew));
1328  fg = pAdd(fg, ppMult_qq(g, fNew));
1329  fg = pAdd(fg, ppMult_qq(fNew, gNew));
1330  f = pAdd(f, fNew);
1331  g = pAdd(g, gNew);
1332  }
1333  /* clean-up loop-dependent vectors */
1334  idDelete((ideal*)&bVec); idDelete((ideal*)&xVec);
1335  }
1336 
1337  /* clean-up matrices A, P, L and U, and polynomial fg */
1338  idDelete((ideal*)&aMat); idDelete((ideal*)&pMat);
1339  idDelete((ideal*)&lMat); idDelete((ideal*)&uMat);
1340  pDelete(&fg);
1341 }
int m
Definition: cfEzgcd.cc:128
g
Definition: cfModGcd.cc:4090
FILE * f
Definition: checklibs.c:9
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
STATIC_VAR Poly * h
Definition: janet.cc:971
void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat, const ring R)
LU-decomposition of a given (m x n)-matrix.
bool luSolveViaLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, const matrix bVec, matrix &xVec, matrix &H)
Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LU-decompositi...
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:37
#define pNext(p)
Definition: monomials.h:36
long p_Deg(poly a, const ring r)
Definition: p_polys.cc:583
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
#define pDelete(p_ptr)
Definition: polys.h:186
#define pNeg(p)
Definition: polys.h:198
#define ppMult_qq(p, q)
Definition: polys.h:208
#define pLmDelete(p)
assume p != NULL, deletes Lm(p)->coef and Lm(p)
Definition: polys.h:76
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
#define pCopy(p)
return a copy of the poly
Definition: polys.h:185

◆ hessenberg()

void hessenberg ( const matrix  aMat,
matrix pMat,
matrix hessenbergMat,
const number  tolerance,
const ring  r 
)

Computes the Hessenberg form of a given square matrix.

The method assumes that all matrix entries are numbers coming from some subfield of the reals.. Afterwards, the following conditions will hold: 1) hessenbergMat = pMat * aMat * pMat^{-1}, 2) hessenbergMat is in Hessenberg form. During the algorithm, pMat will be constructed as the product of self- inverse matrices. The algorithm relies on computing square roots of floating point numbers. These will be combuted by Heron's iteration formula, with iteration stopping when two successive approximations of the square root differ by no more than the given tolerance, which is assumed to be a positive real number.

Parameters
[in]aMatthe square input matrix
[out]pMatthe transformation matrix
[out]hessenbergMatthe Hessenberg form of aMat
[in]toleranceaccuracy

Definition at line 909 of file linearAlgebra.cc.

911 {
912  int n = MATROWS(aMat);
913  unitMatrix(n, pMat);
914  subMatrix(aMat, 1, n, 1, n, hessenbergMat);
915  for (int c = 1; c <= n; c++)
916  {
917  /* find one or two non-zero entries in the current column */
918  int r1 = 0; int r2 = 0;
919  for (int r = c + 1; r <= n; r++)
920  if (MATELEM(hessenbergMat, r, c) != NULL)
921  {
922  if (r1 == 0) r1 = r;
923  else if (r2 == 0) { r2 = r; break; }
924  }
925  if (r1 != 0)
926  { /* At least one entry in the current column is non-zero. */
927  if (r1 != c + 1)
928  { /* swap rows to bring non-zero element to row with index c + 1 */
929  swapRows(r1, c + 1, hessenbergMat);
930  /* now also swap columns to reflect action of permutation
931  from the right-hand side */
932  swapColumns(r1, c + 1, hessenbergMat);
933  /* include action of permutation also in pMat */
934  swapRows(r1, c + 1, pMat);
935  }
936  if (r2 != 0)
937  { /* There is at least one more non-zero element in the current
938  column. So let us perform a hessenberg step in order to make
939  all additional non-zero elements zero. */
940  matrix v; subMatrix(hessenbergMat, c + 1, n, c, c, v);
941  matrix u; matrix pTmp;
942  number r = hessenbergStep(v, u, pTmp, tolerance);
943  idDelete((ideal*)&v); idDelete((ideal*)&u); nDelete(&r);
944  /* pTmp just acts on the lower right block of hessenbergMat;
945  i.e., it needs to be extended by a unit matrix block at the top
946  left in order to define a whole transformation matrix;
947  this code may be optimized */
948  unitMatrix(c, u);
949  matrix pTmpFull; matrixBlock(u, pTmp, pTmpFull);
950  idDelete((ideal*)&u); idDelete((ideal*)&pTmp);
951  /* now include pTmpFull in pMat (by letting it act from the left) */
952  pTmp = mp_Mult(pTmpFull, pMat,R); idDelete((ideal*)&pMat);
953  pMat = pTmp;
954  /* now let pTmpFull act on hessenbergMat from the left and from the
955  right (note that pTmpFull is self-inverse) */
956  pTmp = mp_Mult(pTmpFull, hessenbergMat,R);
957  idDelete((ideal*)&hessenbergMat);
958  hessenbergMat = mp_Mult(pTmp, pTmpFull, R);
959  idDelete((ideal*)&pTmp); idDelete((ideal*)&pTmpFull);
960  /* as there may be inaccuracy, we erase those entries of hessenbergMat
961  which must have become zero by the last transformation */
962  for (int r = c + 2; r <= n; r++)
963  pDelete(&MATELEM(hessenbergMat, r, c));
964  }
965  }
966  }
967 }
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
number hessenbergStep(const matrix vVec, matrix &uVec, matrix &pMat, const number tolerance)
Computes information related to one householder transformation step for constructing the Hessenberg f...
void swapColumns(int column1, int column2, matrix &aMat)
Swaps two columns of a given matrix and thereby modifies it.
void matrixBlock(const matrix aMat, const matrix bMat, matrix &block)
Creates a new matrix which contains the first argument in the top left corner, and the second in the ...
bool unitMatrix(const int n, matrix &unitMat, const ring R)
Creates a new matrix which is the (nxn) unit matrix, and returns true in case of success.
void swapRows(int row1, int row2, matrix &aMat)
Swaps two rows of a given matrix and thereby modifies it.
bool subMatrix(const matrix aMat, const int rowIndex1, const int rowIndex2, const int colIndex1, const int colIndex2, matrix &subMat)
Creates a new matrix which is a submatrix of the first argument, and returns true in case of success.
matrix mp_Mult(matrix a, matrix b, const ring R)
Definition: matpol.cc:213
#define R
Definition: sirandom.c:27

◆ hessenbergStep()

number hessenbergStep ( const matrix  vVec,
matrix uVec,
matrix pMat,
const number  tolerance 
)

Computes information related to one householder transformation step for constructing the Hessenberg form of a given non-derogative matrix.

The method assumes that all matrix entries are numbers coming from some subfield of the reals. And that v has a non-zero first entry v_1 and a second non-zero entry somewhere else. Given such a vector v, it computes a number r (which will be the return value of the method), a vector u and a matrix P such that: 1) P * v = r * e_1, 2) P = E - u * u^T, 3) P = P^{-1}. Note that enforcing norm(u) = sqrt(2), which is done in the algorithm, guarantees property 3). This can be checked by expanding P^2 using property 2).

Returns
the number r
Parameters
[in]vVecthe input vector v
[out]uVecthe output vector u
[out]pMatthe output matrix P
[in]toleranceaccuracy for square roots

Definition at line 837 of file linearAlgebra.cc.

843 {
844  int rr = MATROWS(vVec);
845  number vNormSquared = euclideanNormSquared(vVec);
846  number vNorm; realSqrt(vNormSquared, tolerance, vNorm);
847  /* v1 is guaranteed to be non-zero: */
848  number v1 = pGetCoeff(MATELEM(vVec, 1, 1));
849  bool v1Sign = true; if (nGreaterZero(v1)) v1Sign = false;
850 
851  number v1Abs = nCopy(v1); if (v1Sign) v1Abs = nInpNeg(v1Abs);
852  number t1 = nDiv(v1Abs, vNorm);
853  number one = nInit(1);
854  number t2 = nAdd(t1, one); nDelete(&t1);
855  number denominator; realSqrt(t2, tolerance, denominator); nDelete(&t2);
856  uVec = mpNew(rr, 1);
857  t1 = nDiv(v1Abs, vNorm);
858  t2 = nAdd(t1, one); nDelete(&t1);
859  t1 = nDiv(t2, denominator); nDelete(&t2);
860  MATELEM(uVec, 1, 1) = pOne();
861  pSetCoeff(MATELEM(uVec, 1, 1), t1); /* we know that t1 != 0 */
862  for (int r = 2; r <= rr; r++)
863  {
864  if (MATELEM(vVec, r, 1) != NULL)
865  t1 = nCopy(pGetCoeff(MATELEM(vVec, r, 1)));
866  else t1 = nInit(0);
867  if (v1Sign) t1 = nInpNeg(t1);
868  t2 = nDiv(t1, vNorm); nDelete(&t1);
869  t1 = nDiv(t2, denominator); nDelete(&t2);
870  if (!nIsZero(t1))
871  {
872  MATELEM(uVec, r, 1) = pOne();
873  pSetCoeff(MATELEM(uVec, r, 1), t1);
874  }
875  else nDelete(&t1);
876  }
877  nDelete(&denominator);
878 
879  /* finished building vector u; now turn to pMat */
880  pMat = mpNew(rr, rr);
881  /* we set P := E - u * u^T, as desired */
882  for (int r = 1; r <= rr; r++)
883  for (int c = 1; c <= rr; c++)
884  {
885  if ((MATELEM(uVec, r, 1) != NULL) && (MATELEM(uVec, c, 1) != NULL))
886  t1 = nMult(pGetCoeff(MATELEM(uVec, r, 1)),
887  pGetCoeff(MATELEM(uVec, c, 1)));
888  else t1 = nInit(0);
889  if (r == c) { t2 = nSub(one, t1); nDelete(&t1); }
890  else t2 = nInpNeg(t1);
891  if (!nIsZero(t2))
892  {
893  MATELEM(pMat, r, c) = pOne();
894  pSetCoeff(MATELEM(pMat, r, c), t2);
895  }
896  else nDelete(&t2);
897  }
898  nDelete(&one);
899  /* finished building pMat; now compute the return value */
900  t1 = vNormSquared; if (v1Sign) t1 = nInpNeg(t1);
901  t2 = nMult(v1, vNorm);
902  number t3 = nAdd(t1, t2); nDelete(&t1); nDelete(&t2);
903  t1 = nAdd(v1Abs, vNorm); nDelete(&v1Abs); nDelete(&vNorm);
904  t2 = nDiv(t3, t1); nDelete(&t1); nDelete(&t3);
905  t2 = nInpNeg(t2);
906  return t2;
907 }
bool realSqrt(const number n, const number tolerance, number &root)
Computes the square root of a non-negative real number and returns it as a new number.
number euclideanNormSquared(const matrix aMat)
#define nDiv(a, b)
Definition: numbers.h:32

◆ lduDecomp()

void lduDecomp ( const matrix  aMat,
matrix pMat,
matrix lMat,
matrix dMat,
matrix uMat,
poly &  l,
poly &  u,
poly &  lTimesU 
)

LU-decomposition of a given (m x n)-matrix with performing only those divisions that yield zero remainders.

Given an (m x n) matrix A, the method computes its LDU-decomposition, that is, it computes matrices P, L, D, and U such that

  • P * A = L * D^(-1) * U,
  • P is an (m x m) permutation matrix, i.e., its row/columns form the standard basis of K^m,
  • L is an (m x m) matrix in lower triangular form of full rank,
  • D is an (m x m) diagonal matrix with full rank, and
  • U is an (m x n) matrix in upper row echelon form.
    From these conditions, it follows immediately that also A = P * L * D^(-1) * U, since P is self-inverse.

In contrast to luDecomp, this method only performs those divisions which yield zero remainders. Hence, when e.g. computing over a rational function field and starting with polynomial entries only (or over Q and starting with integer entries only), then any newly computed matrix entry will again be polynomial (or integer).

The method will modify all argument matrices except aMat, so that they will finally contain the matrices P, L, D, and U as specified above. Moreover, in order to further speed up subsequent calls of luSolveViaLDUDecomp, the two denominators and their product will also be returned.

Parameters
[in]aMatthe initial matrix A,
[out]pMatthe row permutation matrix P
[out]lMatthe lower triangular matrix L
[out]dMatthe diagonal matrix D
[out]uMatthe upper row echelon matrix U
[out]lthe product of pivots of L
[out]uthe product of pivots of U
[out]lTimesUthe product of l and u

Definition at line 1343 of file linearAlgebra.cc.

1345 {
1346  int rr = aMat->rows();
1347  int cc = aMat->cols();
1348  /* we use an int array to store all row permutations;
1349  note that we only make use of the entries [1..rr] */
1350  int* permut = new int[rr + 1];
1351  for (int i = 1; i <= rr; i++) permut[i] = i;
1352  /* fill lMat and dMat with the (rr x rr) unit matrix */
1353  unitMatrix(rr, lMat);
1354  unitMatrix(rr, dMat);
1355  uMat = mpNew(rr, cc);
1356  /* copy aMat into uMat: */
1357  for (int r = 1; r <= rr; r++)
1358  for (int c = 1; c <= cc; c++)
1359  MATELEM(uMat, r, c) = pCopy(MATELEM(aMat, r, c));
1360  u = pOne(); l = pOne();
1361 
1362  int col = 1; int row = 1;
1363  while ((col <= cc) & (row < rr))
1364  {
1365  int pivotR; int pivotC; bool pivotValid = false;
1366  while (col <= cc)
1367  {
1368  pivotValid = pivot(uMat, row, rr, col, col, &pivotR, &pivotC);
1369  if (pivotValid) break;
1370  col++;
1371  }
1372  if (pivotValid)
1373  {
1374  if (pivotR != row)
1375  {
1376  swapRows(row, pivotR, uMat);
1377  poly p = MATELEM(dMat, row, row);
1378  MATELEM(dMat, row, row) = MATELEM(dMat, pivotR, pivotR);
1379  MATELEM(dMat, pivotR, pivotR) = p;
1380  swapColumns(row, pivotR, lMat);
1381  swapRows(row, pivotR, lMat);
1382  int temp = permut[row];
1383  permut[row] = permut[pivotR]; permut[pivotR] = temp;
1384  }
1385  /* in gg, we compute the gcd of all non-zero elements in
1386  uMat[row..rr, col];
1387  the next number is the pivot and thus guaranteed to be different
1388  from zero: */
1389  number gg = nCopy(pGetCoeff(MATELEM(uMat, row, col))); number t;
1390  for (int r = row + 1; r <= rr; r++)
1391  {
1392  if (MATELEM(uMat, r, col) != NULL)
1393  {
1394  t = gg;
1395  gg = n_Gcd(t, pGetCoeff(MATELEM(uMat, r, col)),currRing->cf);
1396  nDelete(&t);
1397  }
1398  }
1399  t = nDiv(pGetCoeff(MATELEM(uMat, row, col)), gg);
1400  nNormalize(t); /* this division works without remainder */
1401  if (!nIsOne(t))
1402  {
1403  for (int r = row; r <= rr; r++)
1404  MATELEM(dMat, r, r)=__p_Mult_nn(MATELEM(dMat, r, r), t,currRing);
1405  MATELEM(lMat, row, row)=__p_Mult_nn(MATELEM(lMat, row, row), t,currRing);
1406  }
1407  l = pMult(l, pCopy(MATELEM(lMat, row, row)));
1408  u = pMult(u, pCopy(MATELEM(uMat, row, col)));
1409 
1410  for (int r = row + 1; r <= rr; r++)
1411  {
1412  if (MATELEM(uMat, r, col) != NULL)
1413  {
1414  number g = n_Gcd(pGetCoeff(MATELEM(uMat, row, col)),
1415  pGetCoeff(MATELEM(uMat, r, col)),
1416  currRing->cf);
1417  number f1 = nDiv(pGetCoeff(MATELEM(uMat, r, col)), g);
1418  nNormalize(f1); /* this division works without remainder */
1419  number f2 = nDiv(pGetCoeff(MATELEM(uMat, row, col)), g);
1420  nNormalize(f2); /* this division works without remainder */
1421  pDelete(&MATELEM(uMat, r, col)); MATELEM(uMat, r, col) = NULL;
1422  for (int c = col + 1; c <= cc; c++)
1423  {
1424  poly p = MATELEM(uMat, r, c);
1425  p=__p_Mult_nn(p, f2,currRing);
1426  poly q = pCopy(MATELEM(uMat, row, c));
1427  q=__p_Mult_nn(q, f1,currRing); q = pNeg(q);
1428  MATELEM(uMat, r, c) = pAdd(p, q);
1429  }
1430  number tt = nDiv(g, gg);
1431  nNormalize(tt); /* this division works without remainder */
1432  MATELEM(lMat, r, r)=__p_Mult_nn(MATELEM(lMat, r, r), tt, currRing);
1433  nDelete(&tt);
1434  MATELEM(lMat, r, row) = pCopy(MATELEM(lMat, r, r));
1435  MATELEM(lMat, r, row)=__p_Mult_nn(MATELEM(lMat, r, row), f1,currRing);
1436  nDelete(&f1); nDelete(&f2); nDelete(&g);
1437  }
1438  else
1439  MATELEM(lMat, r, r)=__p_Mult_nn(MATELEM(lMat, r, r), t, currRing);
1440  }
1441  nDelete(&t); nDelete(&gg);
1442  col++; row++;
1443  }
1444  }
1445  /* one factor in the product u might be missing: */
1446  if (row == rr)
1447  {
1448  while ((col <= cc) && (MATELEM(uMat, row, col) == NULL)) col++;
1449  if (col <= cc) u = pMult(u, pCopy(MATELEM(uMat, row, col)));
1450  }
1451 
1452  /* building the permutation matrix from 'permut' and computing l */
1453  pMat = mpNew(rr, rr);
1454  for (int r = 1; r <= rr; r++)
1455  MATELEM(pMat, r, permut[r]) = pOne();
1456  delete[] permut;
1457 
1458  lTimesU = ppMult_qq(l, u);
1459 }
int l
Definition: cfEzgcd.cc:100
int & rows()
Definition: matpol.h:23
int & cols()
Definition: matpol.h:24
static FORCE_INLINE number n_Gcd(number a, number b, const coeffs r)
in Z: return the gcd of 'a' and 'b' in Z/nZ, Z/2^kZ: computed as in the case Z in Z/pZ,...
Definition: coeffs.h:664
bool pivot(const matrix aMat, const int r1, const int r2, const int c1, const int c2, int *bestR, int *bestC, const ring R)
This code computes a score for each non-zero matrix entry in aMat[r1..r2, c1..c2].
#define nIsOne(n)
Definition: numbers.h:25
#define nNormalize(n)
Definition: numbers.h:30
#define __p_Mult_nn(p, n, r)
Definition: p_polys.h:943
#define pMult(p, q)
Definition: polys.h:207

◆ lowerLeftTriangleInverse()

bool lowerLeftTriangleInverse ( const matrix  lMat,
matrix iMat,
bool  diagonalIsOne 
)

Computes the inverse of a given (n x n)-matrix in lower left triangular form.

This method expects an (n x n)-matrix, that is, it must have as many rows as columns. Moreover, lMat[i,j] = 0, at least for all j > i, that ism lMat is in lower left triangular form.
If the argument diagonalIsOne is true, then we know additionally, that lMat[i, i] = 1, for all i. In this case lMat is invertible. Contrariwise, if diagonalIsOne is false, we do not know anything about the diagonal entries. (Note that they may still all be 1.)
In general, there are two cases:
1) The matrix is not invertible. Then the method returns false, and &iMat remains unchanged.
2) The matrix is invertible. Then the method returns true, and the content of iMat is the inverse of lMat.

Returns
true iff lMat is invertible, false otherwise
Parameters
[in]lMat(n x n)-matrix in lower left triangular form
[out]iMatinverse of lMat if invertible
[in]diagonalIsOneif true, then all diagonal entries of lMat are 1

Definition at line 300 of file linearAlgebra.cc.

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 }
int k
Definition: cfEzgcd.cc:99
#define nInvers(a)
Definition: numbers.h:33
#define pNSet(n)
Definition: polys.h:313
#define pNormalize(p)
Definition: polys.h:317

◆ luDecomp()

void luDecomp ( const matrix  aMat,
matrix pMat,
matrix lMat,
matrix uMat,
const ring  r = currRing 
)

LU-decomposition of a given (m x n)-matrix.

Given an (m x n) matrix A, the method computes its LU-decomposition, that is, it computes matrices P, L, and U such that

  • P * A = L * U,
  • P is an (m x m) permutation matrix, i.e., its row/columns form the standard basis of K^m,
  • L is an (m x m) matrix in lower triangular form with all diagonal entries equal to 1,
  • U is an (m x n) matrix in upper row echelon form.
    From these conditions, it follows immediately that also A = P * L * U, since P is self-inverse.

The method will modify all argument matrices except aMat, so that they will finally contain the matrices P, L, and U as specified above.

Parameters
[in]aMatthe initial matrix A,
[out]pMatthe row permutation matrix P
[out]lMatthe lower triangular matrix L
[out]uMatthe upper row echelon matrix U
[in]Rcurrent ring

Definition at line 103 of file linearAlgebra.cc.

105 {
106  int rr = aMat->rows();
107  int cc = aMat->cols();
108  pMat = mpNew(rr, rr);
109  uMat = mp_Copy(aMat,R); /* copy aMat into uMat: */
110 
111  /* we use an int array to store all row permutations;
112  note that we only make use of the entries [1..rr] */
113  int* permut = new int[rr + 1];
114  for (int i = 1; i <= rr; i++) permut[i] = i;
115 
116  /* fill lMat with the (rr x rr) unit matrix */
117  unitMatrix(rr, lMat,R);
118 
119  int bestR; int bestC; int intSwap; poly pSwap; int cOffset = 0;
120  for (int r = 1; r < rr; r++)
121  {
122  if (r > cc) break;
123  while ((r + cOffset <= cc) &&
124  (!pivot(uMat, r, rr, r + cOffset, r + cOffset, &bestR, &bestC, R)))
125  cOffset++;
126  if (r + cOffset <= cc)
127  {
128  /* swap rows with indices r and bestR in permut */
129  intSwap = permut[r];
130  permut[r] = permut[bestR];
131  permut[bestR] = intSwap;
132 
133  /* swap rows with indices r and bestR in uMat;
134  it is sufficient to do this for columns >= r + cOffset*/
135  for (int c = r + cOffset; c <= cc; c++)
136  {
137  pSwap = MATELEM(uMat, r, c);
138  MATELEM(uMat, r, c) = MATELEM(uMat, bestR, c);
139  MATELEM(uMat, bestR, c) = pSwap;
140  }
141 
142  /* swap rows with indices r and bestR in lMat;
143  we must do this only for columns < r */
144  for (int c = 1; c < r; c++)
145  {
146  pSwap = MATELEM(lMat, r, c);
147  MATELEM(lMat, r, c) = MATELEM(lMat, bestR, c);
148  MATELEM(lMat, bestR, c) = pSwap;
149  }
150 
151  /* perform next Gauss elimination step, i.e., below the
152  row with index r;
153  we need to adjust lMat and uMat;
154  we are certain that the matrix entry at [r, r + cOffset]
155  is non-zero: */
156  number pivotElement = pGetCoeff(MATELEM(uMat, r, r + cOffset));
157  poly p;
158  for (int rGauss = r + 1; rGauss <= rr; rGauss++)
159  {
160  p = MATELEM(uMat, rGauss, r + cOffset);
161  if (p != NULL)
162  {
163  number n = n_Div(pGetCoeff(p), pivotElement, R->cf);
164  n_Normalize(n,R->cf);
165 
166  /* filling lMat;
167  old entry was zero, so no need to call pDelete(.) */
168  MATELEM(lMat, rGauss, r) = p_NSet(n_Copy(n,R->cf),R);
169 
170  /* adjusting uMat: */
171  MATELEM(uMat, rGauss, r + cOffset) = NULL; p_Delete(&p,R);
172  n = n_InpNeg(n,R->cf);
173  for (int cGauss = r + cOffset + 1; cGauss <= cc; cGauss++)
174  {
175  MATELEM(uMat, rGauss, cGauss)
176  = p_Add_q(MATELEM(uMat, rGauss, cGauss),
177  pp_Mult_nn(MATELEM(uMat, r, cGauss), n, R), R);
178  p_Normalize(MATELEM(uMat, rGauss, cGauss),R);
179  }
180 
181  n_Delete(&n,R->cf); /* clean up n */
182  }
183  }
184  }
185  }
186 
187  /* building the permutation matrix from 'permut' */
188  for (int r = 1; r <= rr; r++)
189  MATELEM(pMat, r, permut[r]) = p_One(R);
190  delete[] permut;
191 
192  return;
193 }
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition: coeffs.h:451
static FORCE_INLINE number n_InpNeg(number n, const coeffs r)
in-place negation of n MUST BE USED: n = n_InpNeg(n) (no copy is returned)
Definition: coeffs.h:557
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
Definition: coeffs.h:615
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition: coeffs.h:455
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition: coeffs.h:578
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition: matpol.cc:64
void p_Normalize(poly p, const ring r)
Definition: p_polys.cc:3847
poly p_One(const ring r)
Definition: p_polys.cc:1309
poly p_NSet(number n, const ring r)
returns the poly representing the number n, destroys n
Definition: p_polys.cc:1465
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:908
static poly pp_Mult_nn(poly p, number n, const ring r)
Definition: p_polys.h:964
static void p_Delete(poly *p, const ring r)
Definition: p_polys.h:873

◆ luInverse()

bool luInverse ( const matrix  aMat,
matrix iMat,
const ring  R 
)

This code first computes the LU-decomposition of aMat, and then calls the method for inverting a matrix which is given by its LU-decomposition.

Computes the inverse of a given (n x n)-matrix.

Parameters
[in]aMatmatrix to be inverted
[out]iMatinverse of aMat if invertible
[in]Rcurrent ring

Definition at line 200 of file linearAlgebra.cc.

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 }
bool luInverseFromLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, matrix &iMat, const ring R)
This code computes the inverse by inverting lMat and uMat, and then performing two matrix multiplicat...
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix

◆ luInverseFromLUDecomp()

bool luInverseFromLUDecomp ( const matrix  pMat,
const matrix  lMat,
const matrix  uMat,
matrix iMat,
const ring  R 
)

This code computes the inverse by inverting lMat and uMat, and then performing two matrix multiplications.

Computes the inverse of an (n x n)-matrix which is given by its LU- decomposition.

Parameters
[in]pMatpermutation matrix of an LU- decomposition
[in]lMatlower left matrix of an LU- decomposition
[in]uMatupper right matrix of an LU- decomposition
[out]iMatinverse of A if invertible

Definition at line 352 of file linearAlgebra.cc.

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 }
bool upperRightTriangleInverse(const matrix uMat, matrix &iMat, bool diagonalIsOne, const ring R)
Computes the inverse of a given (n x n)-matrix in upper right triangular form.
bool lowerLeftTriangleInverse(const matrix lMat, matrix &iMat, bool diagonalIsOne)
Computes the inverse of a given (n x n)-matrix in lower left triangular form.

◆ luRank()

int luRank ( const matrix  aMat,
const bool  isRowEchelon,
const ring  r = currRing 
)

Computes the rank of a given (m x n)-matrix.

The matrix may already be given in row echelon form, e.g., when the user has before called luDecomp and passes the upper triangular matrix U to luRank. In this case, the second argument can be set to true in order to pass this piece of information to the method. Otherwise, luDecomp will be called first to compute the matrix U. The rank is then read off the matrix U.

Returns
the rank as an int
See also
luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat)
Parameters
[in]aMatinput matrix
[in]isRowEchelonif true then aMat is assumed to be already in row echelon form

Definition at line 230 of file linearAlgebra.cc.

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 }
int rankFromRowEchelonForm(const matrix aMat)

◆ luSolveViaLDUDecomp()

bool luSolveViaLDUDecomp ( const matrix  pMat,
const matrix  lMat,
const matrix  dMat,
const matrix  uMat,
const poly  l,
const poly  u,
const poly  lTimesU,
const matrix  bVec,
matrix xVec,
matrix H 
)

Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LDU-decomposition.

The method expects the LDU-decomposition of A, that is, pMat * A = lMat * dMat^(-1) * uMat, where the argument matrices have the appropriate proteries; see method 'lduDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &dMat, matrix &uMat, poly &l, poly &u, poly &lTimesU)'.
Instead of trying to invert A and return x = A^(-1)*b, this method 1) computes b' = l * pMat * b, 2) solves the simple system L * y = b', 3) computes y' = u * dMat * y, 4) solves the simple system U * y'' = y', 5) computes x = 1/(lTimesU) * y''. Note that steps 1), 2) and 3) will always work, as L is in any case invertible. Moreover, y and thus y' are uniquely determined. Step 4) will only work if and only if y' is in the column span of U. In that case, there may be infinitely many solutions. In contrast to luSolveViaLUDecomp, this algorithm guarantees that all divisions which have to be performed in steps 2) and 4) will work without remainder. This is due to properties of the given LDU- decomposition. Only the final step 5) performs a division of a vector by a member of the ground field. (Suppose the ground field is Q or some rational function field. Then, if A contains only integers or polynomials, respectively, then steps 1) - 4) will keep all data integer or polynomial, respectively. This may speed up computations considerably.) For the outcome, there are three cases:
1) There is no solution. Then the method returns false, and &xVec as well as &H remain unchanged.
2) There is a unique solution. Then the method returns true and H is the 1x1 matrix with zero entry.
3) There are infinitely many solutions. Then the method returns true and some solution of the given original linear system. Furthermore, the columns of H span the solution space of the homogeneous linear system. The dimension of the solution space is then the number of columns of H.

Returns
true if there is at least one solution, false otherwise
Parameters
[in]pMatpermutation matrix of an LDU- decomposition
[in]lMatlower left matrix of an LDU- decomposition
[in]dMatdiagonal matrix of an LDU- decomposition
[in]uMatupper right matrix of an LDU- decomposition
[in]lpivot product l of an LDU decomp.
[in]upivot product u of an LDU decomp.
[in]lTimesUproduct of l and u
[in]bVecright-hand side of the linear system to be solved
[out]xVecsolution of A*x = b
[out]Hmatrix with columns spanning homogeneous solution space

Definition at line 1461 of file linearAlgebra.cc.

1465 {
1466  int m = uMat->rows(); int n = uMat->cols();
1467  matrix cVec = mpNew(m, 1); /* for storing l * pMat * bVec */
1468  matrix yVec = mpNew(m, 1); /* for storing the unique solution of
1469  lMat * yVec = cVec */
1470 
1471  /* compute cVec = l * pMat * bVec but without actual matrix mult. */
1472  for (int r = 1; r <= m; r++)
1473  {
1474  for (int c = 1; c <= m; c++)
1475  {
1476  if (MATELEM(pMat, r, c) != NULL)
1477  { MATELEM(cVec, r, 1) = ppMult_qq(l, MATELEM(bVec, c, 1)); break; }
1478  }
1479  }
1480 
1481  /* solve lMat * yVec = cVec; this will always work as lMat is invertible;
1482  moreover, all divisions are guaranteed to be without remainder */
1483  poly p; poly q;
1484  for (int r = 1; r <= m; r++)
1485  {
1486  p = pNeg(pCopy(MATELEM(cVec, r, 1)));
1487  for (int c = 1; c < r; c++)
1488  p = pAdd(p, ppMult_qq(MATELEM(lMat, r, c), MATELEM(yVec, c, 1) ));
1489  /* The following division is without remainder. */
1490  q = pNSet(nInvers(pGetCoeff(MATELEM(lMat, r, r))));
1491  MATELEM(yVec, r, 1) = pMult(pNeg(p), q);
1492  pNormalize(MATELEM(yVec, r, 1));
1493  }
1494 
1495  /* compute u * dMat * yVec and put result into yVec */
1496  for (int r = 1; r <= m; r++)
1497  {
1498  p = ppMult_qq(u, MATELEM(dMat, r, r));
1499  MATELEM(yVec, r, 1) = pMult(p, MATELEM(yVec, r, 1));
1500  }
1501 
1502  /* determine whether uMat * xVec = yVec is solvable */
1503  bool isSolvable = true;
1504  bool isZeroRow; int nonZeroRowIndex;
1505  for (int r = m; r >= 1; r--)
1506  {
1507  isZeroRow = true;
1508  for (int c = 1; c <= n; c++)
1509  if (MATELEM(uMat, r, c) != NULL) { isZeroRow = false; break; }
1510  if (isZeroRow)
1511  {
1512  if (MATELEM(yVec, r, 1) != NULL) { isSolvable = false; break; }
1513  }
1514  else { nonZeroRowIndex = r; break; }
1515  }
1516 
1517  if (isSolvable)
1518  {
1519  xVec = mpNew(n, 1);
1520  matrix N = mpNew(n, n); int dim = 0;
1521  /* solve uMat * xVec = yVec and determine a basis of the solution
1522  space of the homogeneous system uMat * xVec = 0;
1523  We do not know in advance what the dimension (dim) of the latter
1524  solution space will be. Thus, we start with the possibly too wide
1525  matrix N and later copy the relevant columns of N into H. */
1526  int nonZeroC; int lastNonZeroC = n + 1;
1527  for (int r = nonZeroRowIndex; r >= 1; r--)
1528  {
1529  for (nonZeroC = 1; nonZeroC <= n; nonZeroC++)
1530  if (MATELEM(uMat, r, nonZeroC) != NULL) break;
1531  for (int w = lastNonZeroC - 1; w >= nonZeroC + 1; w--)
1532  {
1533  /* this loop will only be done when the given linear system has
1534  more than one, i.e., infinitely many solutions */
1535  dim++;
1536  /* now we fill two entries of the dim-th column of N */
1537  MATELEM(N, w, dim) = pNeg(pCopy(MATELEM(uMat, r, nonZeroC)));
1538  MATELEM(N, nonZeroC, dim) = pCopy(MATELEM(uMat, r, w));
1539  }
1540  for (int d = 1; d <= dim; d++)
1541  {
1542  /* here we fill the entry of N at [r, d], for each valid vector
1543  that we already have in N;
1544  again, this will only be done when the given linear system has
1545  more than one, i.e., infinitely many solutions */
1546  p = NULL;
1547  for (int c = nonZeroC + 1; c <= n; c++)
1548  if (MATELEM(N, c, d) != NULL)
1549  p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(N, c, d)));
1550  /* The following division may be with remainder but only takes place
1551  when dim > 0. */
1552  q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
1553  MATELEM(N, nonZeroC, d) = pMult(pNeg(p), q);
1554  pNormalize(MATELEM(N, nonZeroC, d));
1555  }
1556  p = pNeg(pCopy(MATELEM(yVec, r, 1)));
1557  for (int c = nonZeroC + 1; c <= n; c++)
1558  if (MATELEM(xVec, c, 1) != NULL)
1559  p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(xVec, c, 1)));
1560  /* The following division is without remainder. */
1561  q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
1562  MATELEM(xVec, nonZeroC, 1) = pMult(pNeg(p), q);
1563  pNormalize(MATELEM(xVec, nonZeroC, 1));
1564  lastNonZeroC = nonZeroC;
1565  }
1566 
1567  /* divide xVec by l*u = lTimesU and put result in xVec */
1568  number z = nInvers(pGetCoeff(lTimesU));
1569  for (int c = 1; c <= n; c++)
1570  {
1571  MATELEM(xVec, c, 1)=__p_Mult_nn(MATELEM(xVec, c, 1), z,currRing);
1572  pNormalize(MATELEM(xVec, c, 1));
1573  }
1574  nDelete(&z);
1575 
1576  if (dim == 0)
1577  {
1578  /* that means the given linear system has exactly one solution;
1579  we return as H the 1x1 matrix with entry zero */
1580  H = mpNew(1, 1);
1581  }
1582  else
1583  {
1584  /* copy the first 'dim' columns of N into H */
1585  H = mpNew(n, dim);
1586  for (int r = 1; r <= n; r++)
1587  for (int c = 1; c <= dim; c++)
1588  MATELEM(H, r, c) = pCopy(MATELEM(N, r, c));
1589  }
1590  /* clean up N */
1591  idDelete((ideal*)&N);
1592  }
1593 
1594  idDelete((ideal*)&cVec);
1595  idDelete((ideal*)&yVec);
1596 
1597  return isSolvable;
1598 }
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
const CanonicalForm & w
Definition: facAbsFact.cc:51
CanonicalForm H
Definition: facAbsFact.cc:60
int dim(ideal I, ring r)

◆ luSolveViaLUDecomp()

bool luSolveViaLUDecomp ( const matrix  pMat,
const matrix  lMat,
const matrix  uMat,
const matrix  bVec,
matrix xVec,
matrix H 
)

Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LU-decomposition.

The method expects the LU-decomposition of A, that is, pMat * A = lMat * uMat, where the argument matrices have the appropriate proteries; see method 'luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat)'.
Instead of trying to invert A and return x = A^(-1)*b, this method 1) computes b' = pMat * b, 2) solves the simple system L * y = b', and then 3) solves the simple system U * x = y. Note that steps 1) and 2) will always work, as L is in any case invertible. Moreover, y is uniquely determined. Step 3) will only work if and only if y is in the column span of U. In that case, there may be infinitely many solutions. Thus, there are three cases:
1) There is no solution. Then the method returns false, and &xVec as well as &H remain unchanged.
2) There is a unique solution. Then the method returns true and H is the 1x1 matrix with zero entry.
3) There are infinitely many solutions. Then the method returns true and some solution of the given original linear system. Furthermore, the columns of H span the solution space of the homogeneous linear system. The dimension of the solution space is then the number of columns of H.

Returns
true if there is at least one solution, false otherwise
Parameters
[in]pMatpermutation matrix of an LU- decomposition
[in]lMatlower left matrix of an LU- decomposition
[in]uMatupper right matrix of an LU- decomposition
[in]bVecright-hand side of the linear system to be solved
[out]xVecsolution of A*x = b
[out]Hmatrix with columns spanning homogeneous solution space

Definition at line 377 of file linearAlgebra.cc.

380 {
381  int m = uMat->rows(); int n = uMat->cols();
382  matrix cVec = mpNew(m, 1); /* for storing pMat * bVec */
383  matrix yVec = mpNew(m, 1); /* for storing the unique solution of
384  lMat * yVec = cVec */
385 
386  /* compute cVec = pMat * bVec but without actual multiplications */
387  for (int r = 1; r <= m; r++)
388  {
389  for (int c = 1; c <= m; c++)
390  {
391  if (MATELEM(pMat, r, c) != NULL)
392  { MATELEM(cVec, r, 1) = pCopy(MATELEM(bVec, c, 1)); break; }
393  }
394  }
395 
396  /* solve lMat * yVec = cVec; this will always work as lMat is invertible;
397  moreover, no divisions are needed, since lMat[i, i] = 1, for all i */
398  for (int r = 1; r <= m; r++)
399  {
400  poly p = pNeg(pCopy(MATELEM(cVec, r, 1)));
401  for (int c = 1; c < r; c++)
402  p = pAdd(p, ppMult_qq(MATELEM(lMat, r, c), MATELEM(yVec, c, 1) ));
403  MATELEM(yVec, r, 1) = pNeg(p);
404  pNormalize(MATELEM(yVec, r, 1));
405  }
406 
407  /* determine whether uMat * xVec = yVec is solvable */
408  bool isSolvable = true;
409  bool isZeroRow;
410  int nonZeroRowIndex = 0 ; // handle case that the matrix is zero
411  for (int r = m; r >= 1; r--)
412  {
413  isZeroRow = true;
414  for (int c = 1; c <= n; c++)
415  if (MATELEM(uMat, r, c) != NULL) { isZeroRow = false; break; }
416  if (isZeroRow)
417  {
418  if (MATELEM(yVec, r, 1) != NULL) { isSolvable = false; break; }
419  }
420  else { nonZeroRowIndex = r; break; }
421  }
422 
423  if (isSolvable)
424  {
425  xVec = mpNew(n, 1);
426  matrix N = mpNew(n, n); int dim = 0;
427  poly p; poly q;
428  /* solve uMat * xVec = yVec and determine a basis of the solution
429  space of the homogeneous system uMat * xVec = 0;
430  We do not know in advance what the dimension (dim) of the latter
431  solution space will be. Thus, we start with the possibly too wide
432  matrix N and later copy the relevant columns of N into H. */
433  int nonZeroC = 0 ;
434  int lastNonZeroC = n + 1;
435 
436  for (int r = nonZeroRowIndex; r >= 1; r--)
437  {
438  for (nonZeroC = 1; nonZeroC <= n; nonZeroC++)
439  if (MATELEM(uMat, r, nonZeroC) != NULL) break;
440 
441  for (int w = lastNonZeroC - 1; w >= nonZeroC + 1; w--)
442  {
443  /* this loop will only be done when the given linear system has
444  more than one, i.e., infinitely many solutions */
445  dim++;
446  /* now we fill two entries of the dim-th column of N */
447  MATELEM(N, w, dim) = pNeg(pCopy(MATELEM(uMat, r, nonZeroC)));
448  MATELEM(N, nonZeroC, dim) = pCopy(MATELEM(uMat, r, w));
449  }
450  for (int d = 1; d <= dim; d++)
451  {
452  /* here we fill the entry of N at [r, d], for each valid vector
453  that we already have in N;
454  again, this will only be done when the given linear system has
455  more than one, i.e., infinitely many solutions */
456  p = NULL;
457  for (int c = nonZeroC + 1; c <= n; c++)
458  if (MATELEM(N, c, d) != NULL)
459  p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(N, c, d)));
460  q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
461  MATELEM(N, nonZeroC, d) = pMult(pNeg(p), q);
462  pNormalize(MATELEM(N, nonZeroC, d));
463  }
464  p = pNeg(pCopy(MATELEM(yVec, r, 1)));
465  for (int c = nonZeroC + 1; c <= n; c++)
466  if (MATELEM(xVec, c, 1) != NULL)
467  p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(xVec, c, 1)));
468  q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
469  MATELEM(xVec, nonZeroC, 1) = pMult(pNeg(p), q);
470  pNormalize(MATELEM(xVec, nonZeroC, 1));
471  lastNonZeroC = nonZeroC;
472  }
473  for (int w = lastNonZeroC - 1; w >= 1; w--)
474  {
475  // remaining variables are free
476  dim++;
477  MATELEM(N, w, dim) = pOne();
478  }
479 
480  if (dim == 0)
481  {
482  /* that means the given linear system has exactly one solution;
483  we return as H the 1x1 matrix with entry zero */
484  H = mpNew(1, 1);
485  }
486  else
487  {
488  /* copy the first 'dim' columns of N into H */
489  H = mpNew(n, dim);
490  for (int r = 1; r <= n; r++)
491  for (int c = 1; c <= dim; c++)
492  MATELEM(H, r, c) = pCopy(MATELEM(N, r, c));
493  }
494  /* clean up N */
495  idDelete((ideal*)&N);
496  }
497  idDelete((ideal*)&cVec);
498  idDelete((ideal*)&yVec);
499 
500  return isSolvable;
501 }

◆ matrixBlock()

void matrixBlock ( const matrix  aMat,
const matrix  bMat,
matrix block 
)

Creates a new matrix which contains the first argument in the top left corner, and the second in the bottom right.

All other entries of the resulting matrix which will be created in the third argument, are zero.

Parameters
[in]aMatthe top left input matrix
[in]bMatthe bottom right input matrix
[out]blockthe new block matrix

Definition at line 805 of file linearAlgebra.cc.

806 {
807  int rowsA = MATROWS(aMat);
808  int rowsB = MATROWS(bMat);
809  int n = rowsA + rowsB;
810  block = mpNew(n, n);
811  for (int i = 1; i <= rowsA; i++)
812  for (int j = 1; j <= rowsA; j++)
813  MATELEM(block, i, j) = pCopy(MATELEM(aMat, i, j));
814  for (int i = 1; i <= rowsB; i++)
815  for (int j = 1; j <= rowsB; j++)
816  MATELEM(block, i + rowsA, j + rowsA) = pCopy(MATELEM(bMat, i, j));
817 }
int j
Definition: facHensel.cc:110
#define block
Definition: scanner.cc:666

◆ mpTrafo()

void mpTrafo ( matrix H,
int  it,
const number  tolerance,
const ring  R 
)

Performs one transformation step on the given matrix H as part of the gouverning QR double shift algorithm.

The method will change the given matrix H side-effect-wise. The resulting matrix H' will be in Hessenberg form. The iteration index is needed, since for the 11th and 21st iteration, the transformation step is different from the usual step, to avoid convergence problems of the gouverning QR double shift process (who is also the only caller of this method).

Parameters
H[in/out] the matrix to be transformed
[in]ititeration index
[in]toleranceaccuracy for square roots

Definition at line 979 of file linearAlgebra.cc.

985 {
986  int n = MATROWS(H);
987  number trace; number det; number tmp1; number tmp2; number tmp3;
988 
989  if ((it != 11) && (it != 21)) /* the standard case */
990  {
991  /* in this case 'trace' will really be the trace of the lowermost
992  (2x2) block of hMat */
993  trace = nInit(0);
994  det = nInit(0);
995  if (MATELEM(H, n - 1, n - 1) != NULL)
996  {
997  tmp1 = nAdd(trace, pGetCoeff(MATELEM(H, n - 1, n - 1)));
998  nDelete(&trace);
999  trace = tmp1;
1000  }
1001  if (MATELEM(H, n, n) != NULL)
1002  {
1003  tmp1 = nAdd(trace, pGetCoeff(MATELEM(H, n, n)));
1004  nDelete(&trace);
1005  trace = tmp1;
1006  }
1007  /* likewise 'det' will really be the determinante of the lowermost
1008  (2x2) block of hMat */
1009  if ((MATELEM(H, n - 1, n - 1 ) != NULL) && (MATELEM(H, n, n) != NULL))
1010  {
1011  tmp1 = nMult(pGetCoeff(MATELEM(H, n - 1, n - 1)),
1012  pGetCoeff(MATELEM(H, n, n)));
1013  tmp2 = nAdd(tmp1, det); nDelete(&tmp1); nDelete(&det);
1014  det = tmp2;
1015  }
1016  if ((MATELEM(H, n - 1, n) != NULL) && (MATELEM(H, n, n - 1) != NULL))
1017  {
1018  tmp1 = nMult(pGetCoeff(MATELEM(H, n - 1, n)),
1019  pGetCoeff(MATELEM(H, n, n - 1)));
1020  tmp2 = nSub(det, tmp1); nDelete(&tmp1); nDelete(&det);
1021  det = tmp2;
1022  }
1023  }
1024  else
1025  {
1026  /* for it = 11 or it = 21, we use special formulae to avoid convergence
1027  problems of the gouverning QR double shift algorithm (who is the only
1028  caller of this method) */
1029  /* trace = 3/2 * (|hMat[n, n-1]| + |hMat[n-1, n-2]|) */
1030  tmp1 = nInit(0);
1031  if (MATELEM(H, n, n - 1) != NULL)
1032  { nDelete(&tmp1); tmp1 = nCopy(pGetCoeff(MATELEM(H, n, n - 1))); }
1033  if (!nGreaterZero(tmp1)) tmp1 = nInpNeg(tmp1);
1034  tmp2 = nInit(0);
1035  if (MATELEM(H, n - 1, n - 2) != NULL)
1036  { nDelete(&tmp2); tmp2 = nCopy(pGetCoeff(MATELEM(H, n - 1, n - 2))); }
1037  if (!nGreaterZero(tmp2)) tmp2 = nInpNeg(tmp2);
1038  tmp3 = nAdd(tmp1, tmp2); nDelete(&tmp1); nDelete(&tmp2);
1039  tmp1 = nInit(3); tmp2 = nInit(2);
1040  trace = nDiv(tmp1, tmp2); nDelete(&tmp1); nDelete(&tmp2);
1041  tmp1 = nMult(tmp3, trace); nDelete(&trace);
1042  trace = tmp1;
1043  /* det = (|hMat[n, n-1]| + |hMat[n-1, n-2]|)^2 */
1044  det = nMult(tmp3, tmp3); nDelete(&tmp3);
1045  }
1046  matrix c = mpNew(n, 1);
1047  trace = nInpNeg(trace);
1048  MATELEM(c,1,1) = pAdd(pAdd(pAdd(ppMult_qq(MATELEM(H,1,1), MATELEM(H,1,1)),
1049  ppMult_qq(MATELEM(H,1,2), MATELEM(H,2,1))),
1050  __pp_Mult_nn(MATELEM(H,1,1), trace, currRing)),
1051  __p_Mult_nn(pOne(), det,currRing));
1052  MATELEM(c,2,1) = pAdd(pMult(pCopy(MATELEM(H,2,1)),
1053  pAdd(pCopy(MATELEM(H,1,1)),
1054  pCopy(MATELEM(H,2,2)))),
1055  __pp_Mult_nn(MATELEM(H,2,1), trace,currRing));
1056  MATELEM(c,3,1) = ppMult_qq(MATELEM(H,2,1), MATELEM(H,3,2));
1057  nDelete(&trace); nDelete(&det);
1058 
1059  /* for applying hessenbergStep, we need to make sure that c[1, 1] is
1060  not zero */
1061  if ((MATELEM(c,1,1) != NULL) &&
1062  ((MATELEM(c,2,1) != NULL) || (MATELEM(c,3,1) != NULL)))
1063  {
1064  matrix uVec; matrix hMat;
1065  tmp1 = hessenbergStep(c, uVec, hMat, tolerance); nDelete(&tmp1);
1066  /* now replace H by hMat * H * hMat: */
1067  matrix wMat = mp_Mult(hMat, H,R); idDelete((ideal*)&H);
1068  matrix H1 = mp_Mult(wMat, hMat,R);
1069  idDelete((ideal*)&wMat); idDelete((ideal*)&hMat);
1070  /* now need to re-establish Hessenberg form of H1 and put it in H */
1071  hessenberg(H1, wMat, H, tolerance,R);
1072  idDelete((ideal*)&wMat); idDelete((ideal*)&H1);
1073  }
1074  else if ((MATELEM(c,1,1) == NULL) && (MATELEM(c,2,1) != NULL))
1075  {
1076  swapRows(1, 2, H);
1077  swapColumns(1, 2, H);
1078  }
1079  else if ((MATELEM(c,1,1) == NULL) && (MATELEM(c,3,1) != NULL))
1080  {
1081  swapRows(1, 3, H);
1082  swapColumns(1, 3, H);
1083  }
1084  else
1085  { /* c is the zero vector or a multiple of e_1;
1086  no hessenbergStep needed */ }
1087 }
void hessenberg(const matrix aMat, matrix &pMat, matrix &hessenbergMat, const number tolerance, const ring R)
Computes the Hessenberg form of a given square matrix.
#define __pp_Mult_nn(p, n, r)
Definition: p_polys.h:974

◆ pivot()

bool pivot ( const matrix  aMat,
const int  r1,
const int  r2,
const int  c1,
const int  c2,
int *  bestR,
int *  bestC,
const ring  R 
)

This code computes a score for each non-zero matrix entry in aMat[r1..r2, c1..c2].

Finds the best pivot element in some part of a given matrix.

If all entries are zero, false is returned, otherwise true. In the latter case, the minimum of all scores is sought, and the row and column indices of the corresponding matrix entry are stored in bestR and bestC.

Parameters
[in]aMatany matrix with number entries
[in]r1lower row index
[in]r2upper row index
[in]c1lower column index
[in]c2upper column index
[out]bestRaddress of row index of best pivot element
[out]bestCaddress of column index of best pivot element
[in]Rcurrent ring

Definition at line 68 of file linearAlgebra.cc.

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 }
int pivotScore(number n, const ring r)
The returned score is based on the implementation of 'nSize' for numbers (, see numbers....

◆ pivotScore()

int pivotScore ( number  n,
const ring  r 
)

The returned score is based on the implementation of 'nSize' for numbers (, see numbers.h): nSize(n) provides a measure for the complexity of n.

Returns a pivot score for any non-zero matrix entry.

Thus, less complex pivot elements will be prefered, and get therefore a smaller pivot score. Consequently, we simply return the value of nSize. An exception to this rule are the ground fields R, long R, and long C: In these, the result of nSize relates to |n|. And since a larger modulus of the pivot element ensures a numerically more stable Gauss elimination, we would like to have a smaller score for larger values of |n|. Therefore, in R, long R, and long C, the negative of nSize will be returned.

Parameters
[in]na non-zero matrix entry
[in]rcurrent ring

Definition at line 50 of file linearAlgebra.cc.

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 }
static FORCE_INLINE int n_Size(number n, const coeffs r)
return a non-negative measure for the complexity of n; return 0 only when n represents zero; (used fo...
Definition: coeffs.h:570
const CanonicalForm int s
Definition: facAbsFact.cc:51
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:519
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:546
static BOOLEAN rField_is_long_R(const ring r)
Definition: ring.h:543

◆ printMatrix()

void printMatrix ( const matrix  m)

Definition at line 522 of file linearAlgebra.cc.

523 {
524  int rr = MATROWS(m); int cc = MATCOLS(m);
525  printf("\n-------------\n");
526  for (int r = 1; r <= rr; r++)
527  {
528  for (int c = 1; c <= cc; c++)
529  printf("%s ", pString(MATELEM(m, r, c)));
530  printf("\n");
531  }
532  printf("-------------\n");
533 }
char * pString(poly p)
Definition: polys.h:306

◆ printNumber()

void printNumber ( const number  z)

Definition at line 506 of file linearAlgebra.cc.

507 {
508  if (nIsZero(z)) printf("number = 0\n");
509  else
510  {
511  poly p = pOne();
512  pSetCoeff(p, nCopy(z));
513  pSetm(p);
514  printf("number = %s\n", pString(p));
515  pDelete(&p);
516  }
517 }

◆ printSolutions()

void printSolutions ( const int  a,
const int  b,
const int  c 
)

Definition at line 575 of file linearAlgebra.cc.

576 {
577  printf("\n------\n");
578  /* build the polynomial a*x^2 + b*x + c: */
579  poly p = NULL; poly q = NULL; poly r = NULL;
580  if (a != 0)
581  { p = pOne(); pSetExp(p, 1, 2); pSetm(p); pSetCoeff(p, nInit(a)); }
582  if (b != 0)
583  { q = pOne(); pSetExp(q, 1, 1); pSetm(q); pSetCoeff(q, nInit(b)); }
584  if (c != 0)
585  { r = pOne(); pSetCoeff(r, nInit(c)); }
586  p = pAdd(p, q); p = pAdd(p, r);
587  printf("poly = %s\n", pString(p));
588  number tol = tenToTheMinus(20);
589  number s1; number s2; int nSol = quadraticSolve(p, s1, s2, tol);
590  nDelete(&tol);
591  printf("solution code = %d\n", nSol);
592  if ((1 <= nSol) && (nSol <= 3))
593  {
594  if (nSol != 3) { printNumber(s1); nDelete(&s1); }
595  else { printNumber(s1); nDelete(&s1); printNumber(s2); nDelete(&s2); }
596  }
597  printf("------\n");
598  pDelete(&p);
599 }
void printNumber(const number z)
int quadraticSolve(const poly p, number &s1, number &s2, const number tolerance)
Returns all solutions of a quadratic polynomial relation with real coefficients.
number tenToTheMinus(const int exponent)
Returns one over the exponent-th power of ten.

◆ qrDS()

bool qrDS ( const int  n,
matrix queue,
int &  queueL,
number *  eigenValues,
int &  eigenValuesL,
const number  tol1,
const number  tol2,
const ring  R 
)

Definition at line 1090 of file linearAlgebra.cc.

1100 {
1101  bool deflationFound = true;
1102  /* we loop until the working queue is empty,
1103  provided we always find deflation */
1104  while (deflationFound && (queueL > 0))
1105  {
1106  /* take out last queue entry */
1107  matrix currentMat = queue[queueL - 1]; queueL--;
1108  int m = MATROWS(currentMat);
1109  if (m == 1)
1110  {
1111  number newEigenvalue;
1112  /* the entry at [1, 1] is the eigenvalue */
1113  if (MATELEM(currentMat, 1, 1) == NULL) newEigenvalue = nInit(0);
1114  else newEigenvalue = nCopy(pGetCoeff(MATELEM(currentMat, 1, 1)));
1115  eigenValues[eigenValuesL++] = newEigenvalue;
1116  }
1117  else if (m == 2)
1118  {
1119  /* there are two eigenvalues which come as zeros of the characteristic
1120  polynomial */
1121  poly p; charPoly(currentMat, p);
1122  number s1; number s2;
1123  int nSol = quadraticSolve(p, s1, s2, tol2); pDelete(&p);
1124  assume(nSol >= 2);
1125  eigenValues[eigenValuesL++] = s1;
1126  /* if nSol = 2, then s1 is a double zero, and s2 is invalid: */
1127  if (nSol == 2) s2 = nCopy(s1);
1128  eigenValues[eigenValuesL++] = s2;
1129  }
1130  else /* m > 2 */
1131  {
1132  /* bring currentMat into Hessenberg form to fasten computations: */
1133  matrix mm1; matrix mm2;
1134  hessenberg(currentMat, mm1, mm2, tol2,R);
1135  idDelete((ideal*)&currentMat); idDelete((ideal*)&mm1);
1136  currentMat = mm2;
1137  int it = 1; bool doLoop = true;
1138  while (doLoop && (it <= 30 * m))
1139  {
1140  /* search for deflation */
1141  number w1; number w2;
1142  number test1; number test2; bool stopCriterion = false; int k;
1143  for (k = 1; k < m; k++)
1144  {
1145  test1 = absValue(MATELEM(currentMat, k + 1, k));
1146  w1 = absValue(MATELEM(currentMat, k, k));
1147  w2 = absValue(MATELEM(currentMat, k + 1, k + 1));
1148  test2 = nMult(tol1, nAdd(w1, w2));
1149  nDelete(&w1); nDelete(&w2);
1150  if (!nGreater(test1, test2)) stopCriterion = true;
1151  nDelete(&test1); nDelete(&test2);
1152  if (stopCriterion) break;
1153  }
1154  if (k < m) /* found deflation at position (k + 1, k) */
1155  {
1156  pDelete(&MATELEM(currentMat, k + 1, k)); /* make this entry zero */
1157  subMatrix(currentMat, 1, k, 1, k, queue[queueL++]);
1158  subMatrix(currentMat, k + 1, m, k + 1, m, queue[queueL++]);
1159  doLoop = false;
1160  }
1161  else /* no deflation found yet */
1162  {
1163  mpTrafo(currentMat, it, tol2,R);
1164  it++; /* try again */
1165  }
1166  }
1167  if (doLoop) /* could not find deflation for currentMat */
1168  {
1169  deflationFound = false;
1170  }
1171  idDelete((ideal*)&currentMat);
1172  }
1173  }
1174  return deflationFound;
1175 }
number absValue(poly p)
void mpTrafo(matrix &H, int it, const number tolerance, const ring R)
Performs one transformation step on the given matrix H as part of the gouverning QR double shift algo...
#define assume(x)
Definition: mod2.h:387
#define nGreater(a, b)
Definition: numbers.h:28

◆ quadraticSolve()

int quadraticSolve ( const poly  p,
number &  s1,
number &  s2,
const number  tolerance 
)

Returns all solutions of a quadratic polynomial relation with real coefficients.

The method assumes that the polynomial is univariate in the first ring variable, and that the current ground field is the complex numbers. The polynomial has degree <= 2. Thus, there may be a) infinitely many zeros, when p == 0; then -1 is returned; b) no zero, when p == constant <> 0; then 0 is returned; c) one zero, when p is linear; then 1 is returned; d) one double zero; then 2 is returned; e) two distinct zeros; then 3 is returned; In the case e), s1 and s2 will contain the two distinct solutions. In cases c) and d) s1 will contain the single/double solution.

Returns
the number of distinct zeros
Parameters
[in]pthe polynomial
[out]s1first zero, if any
[out]s2second zero, if any
[in]toleranceaccuracy

Definition at line 626 of file linearAlgebra.cc.

628 {
629  poly q = pCopy(p);
630  int result;
631 
632  if (q == NULL) result = -1;
633  else
634  {
635  int degree = pGetExp(q, 1);
636  if (degree == 0) result = 0; /* constant polynomial <> 0 */
637  else
638  {
639  number c2 = nInit(0); /* coefficient of var(1)^2 */
640  number c1 = nInit(0); /* coefficient of var(1)^1 */
641  number c0 = nInit(0); /* coefficient of var(1)^0 */
642  if (pGetExp(q, 1) == 2)
643  { nDelete(&c2); c2 = nCopy(pGetCoeff(q)); q = q->next; }
644  if ((q != NULL) && (pGetExp(q, 1) == 1))
645  { nDelete(&c1); c1 = nCopy(pGetCoeff(q)); q = q->next; }
646  if ((q != NULL) && (pGetExp(q, 1) == 0))
647  { nDelete(&c0); c0 = nCopy(pGetCoeff(q)); q = q->next; }
648 
649  if (degree == 1)
650  {
651  c0 = nInpNeg(c0);
652  s1 = nDiv(c0, c1);
653  result = 1;
654  }
655  else
656  {
657  number tmp = nMult(c0, c2);
658  number tmp2 = nAdd(tmp, tmp); nDelete(&tmp);
659  number tmp4 = nAdd(tmp2, tmp2); nDelete(&tmp2);
660  number discr = nSub(nMult(c1, c1), tmp4); nDelete(&tmp4);
661  if (nIsZero(discr))
662  {
663  tmp = nAdd(c2, c2);
664  s1 = nDiv(c1, tmp); nDelete(&tmp);
665  s1 = nInpNeg(s1);
666  result = 2;
667  }
668  else if (nGreaterZero(discr))
669  {
670  realSqrt(discr, tolerance, tmp); /* sqrt of the discriminant */
671  tmp2 = nSub(tmp, c1);
672  tmp4 = nAdd(c2, c2);
673  s1 = nDiv(tmp2, tmp4); nDelete(&tmp2);
674  tmp = nInpNeg(tmp);
675  tmp2 = nSub(tmp, c1); nDelete(&tmp);
676  s2 = nDiv(tmp2, tmp4); nDelete(&tmp2); nDelete(&tmp4);
677  result = 3;
678  }
679  else
680  {
681  discr = nInpNeg(discr);
682  realSqrt(discr, tolerance, tmp); /* sqrt of |discriminant| */
683  tmp2 = nAdd(c2, c2);
684  tmp4 = nDiv(tmp, tmp2); nDelete(&tmp);
685  tmp = nDiv(c1, tmp2); nDelete(&tmp2);
686  tmp = nInpNeg(tmp);
687  s1 = (number)new gmp_complex(((gmp_complex*)tmp)->real(),
688  ((gmp_complex*)tmp4)->real());
689  tmp4 = nInpNeg(tmp4);
690  s2 = (number)new gmp_complex(((gmp_complex*)tmp)->real(),
691  ((gmp_complex*)tmp4)->real());
692  nDelete(&tmp); nDelete(&tmp4);
693  result = 3;
694  }
695  nDelete(&discr);
696  }
697  nDelete(&c0); nDelete(&c1); nDelete(&c2);
698  }
699  }
700  pDelete(&q);
701 
702  return result;
703 }
int degree(const CanonicalForm &f)

◆ rankFromRowEchelonForm()

int rankFromRowEchelonForm ( const matrix  aMat)

Definition at line 217 of file linearAlgebra.cc.

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 }

◆ realSqrt()

bool realSqrt ( const number  n,
const number  tolerance,
number &  root 
)

Computes the square root of a non-negative real number and returns it as a new number.

The method assumes that the current ground field is the complex numbers, and that the argument has imaginary part zero. If the real part is negative, then false is returned, otherwise true. The square root will be computed in the last argument by Heron's iteration formula with the first argument as the starting value. The iteration will stop as soon as two successive values (in the resulting sequence) differ by no more than the given tolerance, which is assumed to be a positive real number.

Returns
the square root of the non-negative argument number
Parameters
[in]nthe input number
[in]toleranceaccuracy of iteration
[out]rootthe root of n

Definition at line 601 of file linearAlgebra.cc.

602 {
603  if (!nGreaterZero(n)) return false;
604  if (nIsZero(n)) return nInit(0);
605 
606  number oneHalf = complexNumber(0.5, 0.0);
607  number nHalf = nMult(n, oneHalf);
608  root = nCopy(n);
609  number nOld = complexNumber(10.0, 0.0);
610  number nDiff = nCopy(nOld);
611 
612  while (nGreater(nDiff, tolerance))
613  {
614  nDelete(&nOld);
615  nOld = root;
616  root = nAdd(nMult(oneHalf, nOld), nDiv(nHalf, nOld));
617  nDelete(&nDiff);
618  nDiff = nSub(nOld, root);
619  if (!nGreaterZero(nDiff)) nDiff = nInpNeg(nDiff);
620  }
621 
622  nDelete(&nOld); nDelete(&nDiff); nDelete(&oneHalf); nDelete(&nHalf);
623  return true;
624 }
number complexNumber(const double r, const double i)
Creates a new complex number from real and imaginary parts given by doubles.

◆ similar()

int similar ( const number *  nn,
const int  nnLength,
const number  n,
const number  tolerance 
)

Tries to find the number n in the array nn[0..nnLength-1].

The method assumes that the ground field is the complex numbers. n and an entry of nn will be regarded equal when the absolute value of their difference is not larger than the given tolerance. In this case, the index of the respective entry of nn is returned, otherwise -1.

Returns
the index of n in nn (up to tolerance) or -1
Parameters
[in]nnarray of numbers to look-up
[in]nnLengthlength of nn
[in]nnumber to loop-up in nn
[in]tolerancetolerance for comparison

Definition at line 1188 of file linearAlgebra.cc.

1194 {
1195  int result = -1;
1196  number tt = nMult(tolerance, tolerance);
1197  number nr = (number)new gmp_complex(((gmp_complex*)n)->real());
1198  number ni = (number)new gmp_complex(((gmp_complex*)n)->imag());
1199  number rr; number ii;
1200  number w1; number w2; number w3; number w4; number w5;
1201  for (int i = 0; i < nnLength; i++)
1202  {
1203  rr = (number)new gmp_complex(((gmp_complex*)nn[i])->real());
1204  ii = (number)new gmp_complex(((gmp_complex*)nn[i])->imag());
1205  w1 = nSub(nr, rr); w2 = nMult(w1, w1);
1206  w3 = nSub(ni, ii); w4 = nMult(w3, w3);
1207  w5 = nAdd(w2, w4);
1208  if (!nGreater(w5, tt)) result = i;
1209  nDelete(&w1); nDelete(&w2); nDelete(&w3); nDelete(&w4);
1210  nDelete(&w5); nDelete(&rr); nDelete(&ii);
1211  if (result != -1) break;
1212  }
1213  nDelete(&tt); nDelete(&nr); nDelete(&ni);
1214  return result;
1215 }

◆ subMatrix()

bool subMatrix ( const matrix  aMat,
const int  rowIndex1,
const int  rowIndex2,
const int  colIndex1,
const int  colIndex2,
matrix subMat 
)

Creates a new matrix which is a submatrix of the first argument, and returns true in case of success.

The method will be successful whenever rowIndex1 <= rowIndex2 and colIndex1 <= colIndex2. In this case and only then true will be returned and the last argument will afterwards contain a copy of the respective submatrix of the first argument. Note that this method may also be used to copy an entire matrix.

Returns
true iff the submatrix could be build
Parameters
[in]aMatthe input matrix
[in]rowIndex1lower row index
[in]rowIndex2higher row index
[in]colIndex1lower column index
[in]colIndex2higher column index
[out]subMatthe new submatrix

Definition at line 733 of file linearAlgebra.cc.

735 {
736  if (rowIndex1 > rowIndex2) return false;
737  if (colIndex1 > colIndex2) return false;
738  int rr = rowIndex2 - rowIndex1 + 1;
739  int cc = colIndex2 - colIndex1 + 1;
740  subMat = mpNew(rr, cc);
741  for (int r = 1; r <= rr; r++)
742  for (int c = 1; c <= cc; c++)
743  MATELEM(subMat, r, c) =
744  pCopy(MATELEM(aMat, rowIndex1 + r - 1, colIndex1 + c - 1));
745  return true;
746 }

◆ swapColumns()

void swapColumns ( int  column1,
int  column2,
matrix aMat 
)

Swaps two columns of a given matrix and thereby modifies it.

The method expects two valid, distinct column indices, i.e. in [1..n], where n is the number of columns in aMat.

Parameters
[in]column1index of first column to swap
[in]column2index of second column to swap
aMat[in/out] matrix subject to swapping

Definition at line 793 of file linearAlgebra.cc.

794 {
795  poly p;
796  int rr = MATROWS(aMat);
797  for (int r = 1; r <= rr; r++)
798  {
799  p = MATELEM(aMat, r, column1);
800  MATELEM(aMat, r, column1) = MATELEM(aMat, r, column2);
801  MATELEM(aMat, r, column2) = p;
802  }
803 }

◆ swapRows()

void swapRows ( int  row1,
int  row2,
matrix aMat 
)

Swaps two rows of a given matrix and thereby modifies it.

The method expects two valid, distinct row indices, i.e. in [1..n], where n is the number of rows in aMat.

Parameters
[in]row1index of first row to swap
[in]row2index of second row to swap
aMat[in/out] matrix subject to swapping

Definition at line 781 of file linearAlgebra.cc.

782 {
783  poly p;
784  int cc = MATCOLS(aMat);
785  for (int c = 1; c <= cc; c++)
786  {
787  p = MATELEM(aMat, row1, c);
788  MATELEM(aMat, row1, c) = MATELEM(aMat, row2, c);
789  MATELEM(aMat, row2, c) = p;
790  }
791 }

◆ tenToTheMinus()

number tenToTheMinus ( const int  exponent)

Returns one over the exponent-th power of ten.

The method assumes that the base ring is the complex numbers.

return one over the exponent-th power of 10

Parameters
[in]exponentthe exponent for 1/10

Definition at line 554 of file linearAlgebra.cc.

557 {
558  number ten = complexNumber(10.0, 0.0);
559  number result = complexNumber(1.0, 0.0);
560  number tmp;
561  /* compute 10^{-exponent} inside result by subsequent divisions by 10 */
562  for (int i = 1; i <= exponent; i++)
563  {
564  tmp = nDiv(result, ten);
565  nDelete(&result);
566  result = tmp;
567  }
568  nDelete(&ten);
569  return result;
570 }
int exponent(const CanonicalForm &f, int q)
int exponent ( const CanonicalForm & f, int q )

◆ unitMatrix()

bool unitMatrix ( const int  n,
matrix unitMat,
const ring  r = currRing 
)

Creates a new matrix which is the (nxn) unit matrix, and returns true in case of success.

The method will be successful whenever n >= 1. In this case and only then true will be returned and the new (nxn) unit matrix will reside inside the second argument.

Returns
true iff the (nxn) unit matrix could be build
Parameters
[in]nsize of the matrix
[out]unitMatthe new (nxn) unit matrix
[in]Rcurrent ring

Definition at line 95 of file linearAlgebra.cc.

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 }

◆ upperRightTriangleInverse()

bool upperRightTriangleInverse ( const matrix  uMat,
matrix iMat,
bool  diagonalIsOne,
const ring  r = currRing 
)

Computes the inverse of a given (n x n)-matrix in upper right triangular form.

This method expects a quadratic matrix, that is, it must have as many rows as columns. Moreover, uMat[i, j] = 0, at least for all i > j, that is, u is in upper right triangular form.
If the argument diagonalIsOne is true, then we know additionally, that uMat[i, i] = 1, for all i. In this case uMat is invertible. Contrariwise, if diagonalIsOne is false, we do not know anything about the diagonal entries. (Note that they may still all be 1.)
In general, there are two cases:
1) The matrix is not invertible. Then the method returns false, and &iMat remains unchanged.
2) The matrix is invertible. Then the method returns true, and the content of iMat is the inverse of uMat.

Returns
true iff uMat is invertible, false otherwise
Parameters
[in]uMat(n x n)-matrix in upper right triangular form
[out]iMatinverse of uMat if invertible
[in]diagonalIsOneif true, then all diagonal entries of uMat are 1
[in]Rcurrent ring

Definition at line 251 of file linearAlgebra.cc.

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 }
static FORCE_INLINE number n_Invers(number a, const coeffs r)
return the multiplicative inverse of 'a'; raise an error if 'a' is not invertible
Definition: coeffs.h:564
#define p_GetCoeff(p, r)
Definition: monomials.h:50
static poly p_Neg(poly p, const ring r)
Definition: p_polys.h:1079
static poly p_Mult_q(poly p, poly q, const ring r)
Definition: p_polys.h:1086
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition: p_polys.h:1123
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition: p_polys.h:818