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