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