68bool pivot(
const matrix aMat,
const int r1,
const int r2,
const int c1,
69 const int c2,
int* bestR,
int* bestC,
const ring
R)
71 int bestScore;
int score;
bool foundBestScore =
false; poly matEntry;
73 for (
int c = c1; c <= c2; c++)
75 for (
int r = r1; r <= r2; r++)
81 if ((!foundBestScore) || (score < bestScore))
87 foundBestScore =
true;
92 return foundBestScore;
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);
106 int rr = aMat->
rows();
107 int cc = aMat->
cols();
108 pMat =
mpNew(rr, rr);
113 int* permut =
new int[rr + 1];
114 for (
int i = 1;
i <= rr;
i++) permut[
i] =
i;
119 int bestR;
int bestC;
int intSwap; poly pSwap;
int cOffset = 0;
120 for (
int r = 1; r < rr; r++)
123 while ((r + cOffset <= cc) &&
124 (!
pivot(uMat, r, rr, r + cOffset, r + cOffset, &bestR, &bestC,
R)))
126 if (r + cOffset <= cc)
130 permut[r] = permut[bestR];
131 permut[bestR] = intSwap;
135 for (
int c = r + cOffset; c <= cc; c++)
139 MATELEM(uMat, bestR, c) = pSwap;
144 for (
int c = 1; c < r; c++)
148 MATELEM(lMat, bestR, c) = pSwap;
158 for (
int rGauss = r + 1; rGauss <= rr; rGauss++)
160 p =
MATELEM(uMat, rGauss, r + cOffset);
173 for (
int cGauss = r + cOffset + 1; cGauss <= cc; cGauss++)
188 for (
int r = 1; r <= rr; r++)
220 int rr = aMat->
rows();
int cc = aMat->
cols();
221 int r = 1;
int c = 1;
222 while ((r <= rr) && (c <= cc))
225 else { rank++; r++; }
252 bool diagonalIsOne,
const ring
R)
254 int d = uMat->
rows();
258 bool invertible = diagonalIsOne;
262 for (
int r = 1; r <= d; r++)
275 for (
int r = d; r >= 1; r--)
281 for (
int c = r + 1; c <= d; c++)
284 for (
int k = r + 1;
k <= c;
k++)
303 int d = lMat->
rows(); poly
p; poly q;
306 bool invertible = diagonalIsOne;
310 for (
int r = 1; r <= d; r++)
323 for (
int c = d; c >= 1; c--)
329 for (
int r = c + 1; r <= d; r++)
332 for (
int k = c;
k <= r - 1;
k++)
381 int m = uMat->
rows();
int n = uMat->
cols();
387 for (
int r = 1; r <=
m; r++)
389 for (
int c = 1; c <=
m; c++)
398 for (
int r = 1; r <=
m; r++)
401 for (
int c = 1; c < r; c++)
408 bool isSolvable =
true;
410 int nonZeroRowIndex = 0 ;
411 for (
int r =
m; r >= 1; r--)
414 for (
int c = 1; c <= n; c++)
415 if (
MATELEM(uMat, r, c) !=
NULL) { isZeroRow =
false;
break; }
418 if (
MATELEM(yVec, r, 1) !=
NULL) { isSolvable =
false;
break; }
420 else { nonZeroRowIndex = r;
break; }
434 int lastNonZeroC = n + 1;
436 for (
int r = nonZeroRowIndex; r >= 1; r--)
438 for (nonZeroC = 1; nonZeroC <= n; nonZeroC++)
441 for (
int w = lastNonZeroC - 1;
w >= nonZeroC + 1;
w--)
450 for (
int d = 1; d <=
dim; d++)
457 for (
int c = nonZeroC + 1; c <= n; c++)
465 for (
int c = nonZeroC + 1; c <= n; c++)
471 lastNonZeroC = nonZeroC;
473 for (
int w = lastNonZeroC - 1;
w >= 1;
w--)
490 for (
int r = 1; r <= n; r++)
491 for (
int c = 1; c <=
dim; c++)
508 if (
nIsZero(z)) printf(
"number = 0\n");
514 printf(
"number = %s\n",
pString(
p));
525 printf(
"\n-------------\n");
526 for (
int r = 1; r <= rr; r++)
528 for (
int c = 1; c <= cc; c++)
532 printf(
"-------------\n");
577 printf(
"\n------\n");
591 printf(
"solution code = %d\n", nSol);
592 if ((1 <= nSol) && (nSol <= 3))
601bool realSqrt(
const number n,
const number tolerance, number &root)
607 number nHalf =
nMult(n, oneHalf);
610 number nDiff =
nCopy(nOld);
618 nDiff =
nSub(nOld, root);
627 const number tolerance)
639 number c2 =
nInit(0);
640 number c1 =
nInit(0);
641 number c0 =
nInit(0);
657 number tmp =
nMult(c0, c2);
710 for (
int r = 1; r <= rr; r++)
734 const int colIndex1,
const int colIndex2,
matrix &subMat)
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++)
744 pCopy(
MATELEM(aMat, rowIndex1 + r - 1, colIndex1 + c - 1));
750 if (
MATROWS(aMat) != 2)
return false;
751 if (
MATCOLS(aMat) != 2)
return false;
752 number
b =
nInit(0); number t;
785 for (
int c = 1; c <= cc; c++)
797 for (
int r = 1; r <= rr; r++)
809 int n = rowsA + rowsB;
811 for (
int i = 1;
i <= rowsA;
i++)
812 for (
int j = 1;
j <= rowsA;
j++)
814 for (
int i = 1;
i <= rowsB;
i++)
815 for (
int j = 1;
j <= rowsB;
j++)
841 const number tolerance
846 number vNorm;
realSqrt(vNormSquared, tolerance, vNorm);
849 bool v1Sign =
true;
if (
nGreaterZero(v1)) v1Sign =
false;
851 number v1Abs =
nCopy(v1);
if (v1Sign) v1Abs =
nInpNeg(v1Abs);
852 number t1 =
nDiv(v1Abs, vNorm);
853 number one =
nInit(1);
855 number denominator;
realSqrt(t2, tolerance, denominator);
nDelete(&t2);
857 t1 =
nDiv(v1Abs, vNorm);
862 for (
int r = 2; r <= rr; r++)
880 pMat =
mpNew(rr, rr);
882 for (
int r = 1; r <= rr; r++)
883 for (
int c = 1; c <= rr; c++)
900 t1 = vNormSquared;
if (v1Sign) t1 =
nInpNeg(t1);
901 t2 =
nMult(v1, vNorm);
910 const number tolerance,
const ring
R)
914 subMatrix(aMat, 1, n, 1, n, hessenbergMat);
915 for (
int c = 1; c <= n; c++)
918 int r1 = 0;
int r2 = 0;
919 for (
int r = c + 1; r <= n; r++)
923 else if (r2 == 0) { r2 = r;
break; }
956 pTmp =
mp_Mult(pTmpFull, hessenbergMat,
R);
958 hessenbergMat =
mp_Mult(pTmp, pTmpFull,
R);
962 for (
int r = c + 2; r <= n; r++)
982 const number tolerance,
987 number trace; number det; number
tmp1; number
tmp2; number tmp3;
989 if ((it != 11) && (it != 21))
1094 number* eigenValues,
1101 bool deflationFound =
true;
1104 while (deflationFound && (queueL > 0))
1107 matrix currentMat = queue[queueL - 1]; queueL--;
1111 number newEigenvalue;
1115 eigenValues[eigenValuesL++] = newEigenvalue;
1122 number s1; number s2;
1125 eigenValues[eigenValuesL++] = s1;
1127 if (nSol == 2) s2 =
nCopy(s1);
1128 eigenValues[eigenValuesL++] = s2;
1137 int it = 1;
bool doLoop =
true;
1138 while (doLoop && (it <= 30 *
m))
1141 number w1; number w2;
1142 number test1; number test2;
bool stopCriterion =
false;
int k;
1143 for (
k = 1;
k <
m;
k++)
1150 if (!
nGreater(test1, test2)) stopCriterion =
true;
1152 if (stopCriterion)
break;
1157 subMatrix(currentMat, 1,
k, 1,
k, queue[queueL++]);
1169 deflationFound =
false;
1174 return deflationFound;
1192 const number tolerance
1196 number tt =
nMult(tolerance, tolerance);
1199 number rr; number ii;
1200 number w1; number w2; number w3; number w4; number w5;
1201 for (
int i = 0;
i < nnLength;
i++)
1220 const poly f0,
const poly g0,
const int d, poly &
f, poly &
g)
1229 poly
p = f0; poly matEntry; number c;
1246 for (
int row = 2; row <= n + 1; row++)
1247 for (
int col = 2; col <=
m; col++)
1249 if (col > row)
break;
1252 for (
int row = n + 2; row <= n +
m; row++)
1253 for (
int col = row - n; col <=
m; col++)
1255 for (
int row = 2; row <=
m + 1; row++)
1256 for (
int col =
m + 2; col <=
m + n; col++)
1258 if (col -
m > row)
break;
1261 for (
int row =
m + 2; row <= n +
m; row++)
1262 for (
int col = row; col <=
m + n; col++)
1276 for (
int xExp = 1; xExp <= d; xExp++)
1285 int bIsZeroVector =
true;
1293 if (matEntry !=
NULL) bIsZeroVector =
false;
1304 poly fNew =
NULL; poly gNew =
NULL;
1305 for (
int row = 1; row <=
m; row++)
1313 gNew =
pAdd(gNew,
p);
1316 for (
int row =
m + 1; row <=
m + n; row++)
1324 fNew =
pAdd(fNew,
p);
1344 matrix &uMat, poly &
l, poly &u, poly &lTimesU)
1346 int rr = aMat->
rows();
1347 int cc = aMat->
cols();
1350 int* permut =
new int[rr + 1];
1351 for (
int i = 1;
i <= rr;
i++) permut[
i] =
i;
1355 uMat =
mpNew(rr, cc);
1357 for (
int r = 1; r <= rr; r++)
1358 for (
int c = 1; c <= cc; c++)
1362 int col = 1;
int row = 1;
1363 while ((col <= cc) & (row < rr))
1365 int pivotR;
int pivotC;
bool pivotValid =
false;
1368 pivotValid =
pivot(uMat, row, rr, col, col, &pivotR, &pivotC);
1369 if (pivotValid)
break;
1382 int temp = permut[row];
1383 permut[row] = permut[pivotR]; permut[pivotR] = temp;
1390 for (
int r = row + 1; r <= rr; r++)
1403 for (
int r = row; r <= rr; r++)
1410 for (
int r = row + 1; r <= rr; r++)
1422 for (
int c = col + 1; c <= cc; c++)
1430 number tt =
nDiv(
g, gg);
1448 while ((col <= cc) && (
MATELEM(uMat, row, col) ==
NULL)) col++;
1453 pMat =
mpNew(rr, rr);
1454 for (
int r = 1; r <= rr; r++)
1463 const poly
l,
const poly u,
const poly lTimesU,
1466 int m = uMat->
rows();
int n = uMat->
cols();
1472 for (
int r = 1; r <=
m; r++)
1474 for (
int c = 1; c <=
m; c++)
1484 for (
int r = 1; r <=
m; r++)
1487 for (
int c = 1; c < r; c++)
1496 for (
int r = 1; r <=
m; r++)
1503 bool isSolvable =
true;
1504 bool isZeroRow;
int nonZeroRowIndex;
1505 for (
int r =
m; r >= 1; r--)
1508 for (
int c = 1; c <= n; c++)
1509 if (
MATELEM(uMat, r, c) !=
NULL) { isZeroRow =
false;
break; }
1512 if (
MATELEM(yVec, r, 1) !=
NULL) { isSolvable =
false;
break; }
1514 else { nonZeroRowIndex = r;
break; }
1526 int nonZeroC;
int lastNonZeroC = n + 1;
1527 for (
int r = nonZeroRowIndex; r >= 1; r--)
1529 for (nonZeroC = 1; nonZeroC <= n; nonZeroC++)
1531 for (
int w = lastNonZeroC - 1;
w >= nonZeroC + 1;
w--)
1540 for (
int d = 1; d <=
dim; d++)
1547 for (
int c = nonZeroC + 1; c <= n; c++)
1557 for (
int c = nonZeroC + 1; c <= n; c++)
1564 lastNonZeroC = nonZeroC;
1569 for (
int c = 1; c <= n; c++)
1586 for (
int r = 1; r <= n; r++)
1587 for (
int c = 1; c <=
dim; c++)
const CanonicalForm CFMap CFMap & N
gmp_complex numbers based on
Coefficient rings, fields and other domains suitable for Singular polynomials.
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
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,...
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
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)
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...
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...
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
const CanonicalForm int s
const Variable & v
< [in] a sqrfree bivariate poly
#define idDelete(H)
delete an ideal
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-decomposit...
void hessenberg(const matrix aMat, matrix &pMat, matrix &hessenbergMat, const number tolerance, const ring R)
Computes the Hessenberg form of a given square matrix.
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...
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.
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...
number complexNumber(const double r, const double i)
Creates a new complex number from real and imaginary parts given by doubles.
void swapColumns(int column1, int column2, matrix &aMat)
Swaps two columns of a given matrix and thereby modifies it.
void printNumber(const number z)
bool qrDS(const int, matrix *queue, int &queueL, number *eigenValues, int &eigenValuesL, const number tol1, const number tol2, const ring R)
int luRank(const matrix aMat, const bool isRowEchelon, const ring R)
Computes the rank of a given (m x n)-matrix.
int rankFromRowEchelonForm(const matrix aMat)
int quadraticSolve(const poly p, number &s1, number &s2, const number tolerance)
Returns all solutions of a quadratic polynomial relation with real coefficients.
void printMatrix(const matrix m)
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 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 ...
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,...
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 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.
int pivotScore(number n, const ring r)
The returned score is based on the implementation of 'nSize' for numbers (, see numbers....
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].
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].
bool charPoly(const matrix aMat, poly &charPoly)
Computes the characteristic polynomial from a quadratic (2x2) matrix and returns true in case of succ...
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 remai...
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 matr...
void swapRows(int row1, int row2, matrix &aMat)
Swaps two rows of a given matrix and thereby modifies it.
number tenToTheMinus(const int exponent)
Returns one over the exponent-th power of ten.
bool lowerLeftTriangleInverse(const matrix lMat, matrix &iMat, bool diagonalIsOne)
Computes the inverse of a given (n x n)-matrix in lower left triangular form.
number euclideanNormSquared(const matrix aMat)
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.
void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat, const ring R)
LU-decomposition of a given (m x n)-matrix.
void printSolutions(const int a, const int b, const int c)
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
matrix mp_Mult(matrix a, matrix b, const ring R)
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
#define MATELEM(mat, i, j)
1-based access to matrix
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
void p_Normalize(poly p, const ring r)
poly p_NSet(number n, const ring r)
returns the poly representing the number n, destroys n
long p_Deg(poly a, const ring r)
static poly p_Neg(poly p, const ring r)
static poly p_Add_q(poly p, poly q, const ring r)
static poly p_Mult_q(poly p, poly q, const ring r)
#define __pp_Mult_nn(p, n, r)
static poly pp_Mult_nn(poly p, number n, const ring r)
static void p_Delete(poly *p, const ring r)
static poly pp_Mult_qq(poly p, poly q, const ring r)
static poly p_Copy(poly p, const ring r)
returns a copy of p
#define __p_Mult_nn(p, n, r)
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
#define pLmDelete(p)
assume p != NULL, deletes Lm(p)->coef and Lm(p)
#define pGetExp(p, i)
Exponent.
#define pCopy(p)
return a copy of the poly
static BOOLEAN rField_is_R(const ring r)
static BOOLEAN rField_is_long_C(const ring r)
static BOOLEAN rField_is_long_R(const ring r)
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix