00001
00002
00003
00004 #include"matrices.h"
00005 #include"gcd.h"
00006 #include"powerest.h"
00007
00008 #ifndef GAUSS_H_
00009 #define GAUSS_H_
00010
00017
00028
00029 #ifdef no_default_arguments
00030 template <class field_item>
00031 void rref(const matrix<field_item>& A,
00032 matrix<field_item>& B, matrix<field_item>& Q,
00033 vector<index>& basis, field_item& det,
00034 int movie, ostr& movie_stream);
00035 #else
00036 template <class field_item>
00037 void rref(const matrix<field_item>& A,
00038 matrix<field_item>& B, matrix<field_item>& Q,
00039 vector<index>& basis, field_item& det,
00040 int movie = 0, ostr& movie_stream = cout);
00041 #endif
00042
00047
00048 template <class field_item>
00049 matrix<field_item> inv(const matrix<field_item>& A);
00050
00055
00056 template <class field_item>
00057 field_item det(const matrix<field_item>& A);
00058
00059
00069
00070 #ifdef no_default_arguments
00071 template <class int_item>
00072 void rref_int(const matrix<int_item>& A,
00073 matrix<int_item>& B, matrix<int_item>& Q,
00074 vector<index>& basis, int_item& det,
00075 int movie, ostr& movie_stream);
00076 #else
00077 template <class int_item>
00078 void rref_int(const matrix<int_item>& A,
00079 matrix<int_item>& B, matrix<int_item>& Q,
00080 vector<index>& basis, int_item& det,
00081 int movie = 0, ostr& movie_stream = cout);
00082 #endif
00083
00084
00085
00089
00090 template <class int_item>
00091 int_item det_int(const matrix<int_item>& A);
00092
00093
00094
00104
00105 #ifdef no_default_arguments
00106 template <class int_item>
00107 void smith(const matrix<int_item>& A,
00108 matrix<int_item>& B, matrix<int_item>& P,
00109 matrix<int_item>& Q, index& rank, int_item& det,
00110 int movie, ostr& movie_stream);
00111 #else
00112 template <class int_item>
00113 void smith(const matrix<int_item>& A,
00114 matrix<int_item>& B, matrix<int_item>& P,
00115 matrix<int_item>& Q, index& rank, int_item& det,
00116 int movie = 0, ostr& movie_stream = cout);
00117 #endif
00118
00119
00120
00121
00122 #ifdef no_default_arguments
00123 template <class field_item>
00124 void rref(const matrix<field_item>& A,
00125 matrix<field_item>& B, matrix<field_item>& Q,
00126 vector<index>& basis, field_item& det,
00127 int movie, ostr& movie_stream)
00128 #else
00129 template <class field_item>
00130 void rref(const matrix<field_item>& A,
00131 matrix<field_item>& B, matrix<field_item>& Q,
00132 vector<index>& basis, field_item& det,
00133 int movie, ostr& movie_stream)
00134 #endif
00135
00136
00137 {
00138 det = 1;
00139 index m = A.get_m();
00140 index n = A.get_n();
00141 Q = eye(field_item(1), m);
00142 B = A;
00143 basis.resize(0);
00144
00145 if (movie)
00146 {
00147 movie_stream << "Reduced row echelon form\n";
00148 print2(movie_stream, B, Q);
00149 }
00150
00151 for (index i = 0, j = 0; i < m && j < n; j++)
00152 {
00153
00154
00155 field_item max = 0;
00156 index i_max = i;
00157
00158
00159 if (movie)
00160 movie_stream
00161 << "Find the biggest (in absolute value) entry in the column "
00162 << j << endl;
00163 for (index ii = i; ii < m; ii++)
00164 {
00165 if (absolute_value(B(ii, j)) > max)
00166 {
00167 max = absolute_value(B(ii, j));
00168 i_max = ii;
00169 }
00170 }
00171 if (max == field_item(0))
00172 {
00173
00174 if (movie)
00175 movie_stream << "Column " << j << " is negligible\n";
00176 det = 0;
00177 continue;
00178 }
00179
00180 if (movie)
00181 movie_stream << "Pivot item: (" << i_max << ", " << j << ")\n";
00182
00183 if (i != i_max)
00184 {
00185 B.swap_rows(i, i_max);
00186 Q.swap_rows(i, i_max);
00187 det = -det;
00188 if (movie)
00189 {
00190 movie_stream << "Swap rows: " << i << ", " << i_max << endl;
00191 print2(movie_stream, B, Q);
00192 }
00193 }
00194
00195 field_item pivot = B(i, j);
00196 B.div_row(i, pivot);
00197 Q.div_row(i, pivot);
00198 det = det * pivot;
00199
00200 for (index ii = 0; ii < m; ii++)
00201 if (ii != i)
00202 {
00203 field_item alpha = B(ii, j);
00204 B.add_row(ii, i, -alpha);
00205 Q.add_row(ii, i, -alpha);
00206 }
00207 if (movie)
00208 {
00209 movie_stream << "Elimination in column " << j << endl;
00210 print2(movie_stream, B, Q);
00211 }
00212
00213 basis.join(j);
00214 i++;
00215 }
00216 }
00217
00218 template <class field_item>
00219 field_item det(const matrix<field_item>& A)
00220 {
00221
00222 if (A.get_m() != A.get_n())
00223 {
00224 matrix_error (st_matrix_must_be_square, "matrix must be square");
00225 return 0;
00226 }
00227
00228 matrix<field_item> B;
00229 matrix<field_item> Q;
00230 vector<index> basis;
00231
00232 field_item det;
00233
00234 rref(A, B, Q, basis, det, 0, cout);
00235
00236 return det;
00237 }
00238
00239
00240
00241
00242 template <class field_item>
00243 matrix<field_item> inv(const matrix<field_item>& A)
00244 {
00245
00246 if (A.get_m() != A.get_n())
00247 {
00248 matrix_error (st_matrix_must_be_square, "matrix must be square");
00249 return 0;
00250 }
00251
00252 matrix<field_item> B;
00253 matrix<field_item> Q;
00254 vector<index> basis;
00255
00256 field_item det;
00257
00258 rref(A, B, Q, basis, det, 0, cout);
00259
00260 if (det == 0)
00261 matrix_error (st_matrix_is_singular, "matrix is singular");
00262
00263 return Q;
00264 }
00265
00266
00267 #ifdef no_default_arguments
00268 template <class int_item>
00269 void rref_int(const matrix<int_item>& A,
00270 matrix<int_item>& B, matrix<int_item>& Q,
00271 vector<index>& basis, int_item& det,
00272 int movie, ostr& movie_stream)
00273 #else
00274 template <class int_item>
00275 void rref_int(const matrix<int_item>& A,
00276 matrix<int_item>& B, matrix<int_item>& Q,
00277 vector<index>& basis, int_item& det,
00278 int movie, ostr& movie_stream)
00279 #endif
00280 {
00281 det = 1;
00282 int_item det_denom = 1;
00283 index m = A.get_m();
00284 index n = A.get_n();
00285 Q = eye(int_item(1), m);
00286 B = A;
00287 basis.resize(0);
00288
00289 if (movie)
00290 {
00291 movie_stream << "Reduced row echelon form for integral matrix\n";
00292 print2(movie_stream, B, Q);
00293 }
00294
00295 for (index i = 0, j = 0; i < m && j < n; j++)
00296 {
00297
00298
00299 int_item max = 0;
00300 index i_max = i;
00301
00302
00303 if (movie)
00304 movie_stream
00305 << "Find non-zero entry in the column "
00306 << j << endl;
00307 for (index ii = i; ii < m; ii++)
00308 {
00309 if (absolute_value(B(ii, j)) > max)
00310 {
00311 max = absolute_value(B(ii, j));
00312 i_max = ii;
00313 }
00314 }
00315 if (max == 0)
00316 {
00317
00318 if (movie)
00319 movie_stream << "Column " << j << " is negligible\n";
00320 continue;
00321 }
00322
00323 if (movie)
00324 movie_stream << "Pivot item: (" << i_max << ", " << j << ")\n";
00325
00326 if (i != i_max)
00327 {
00328 B.swap_rows(i, i_max);
00329 Q.swap_rows(i, i_max);
00330 det = -det;
00331 if (movie)
00332 {
00333 movie_stream << "Swap rows: " << i << ", " << i_max << endl;
00334 print2(movie_stream, B, Q);
00335 }
00336 }
00337
00338 int_item pivot = B(i, j);
00339
00340 if (pivot < 0)
00341 {
00342 B.mult_row(i, -1);
00343 Q.mult_row(i, -1);
00344 pivot = -pivot;
00345 det = -det;
00346 }
00347
00348 for (index ii = 0; ii < m; ii++)
00349 if (ii != i)
00350 {
00351 int_item alpha = B(ii, j);
00352 int_item delta = gcd(pivot, alpha);
00353 int_item b_ii = -alpha / delta;
00354 int_item b_i = pivot / delta;
00355 B.mult_row(ii, b_i);
00356 B.add_row(ii, i, b_ii);
00357 Q.mult_row(ii, b_i);
00358 Q.add_row(ii, i, b_ii);
00359 det_denom = det_denom * b_i;
00360 alpha = gcd(gcd(B.row(ii)), gcd(Q.row(ii)));
00361 B.div_row(ii, alpha);
00362 Q.div_row(ii, alpha);
00363 det = det * alpha;
00364 }
00365 if (movie)
00366 {
00367 movie_stream << "Elimination in column " << j << endl;
00368 print2(movie_stream, B, Q);
00369 }
00370
00371
00372 basis.join(j);
00373 i++;
00374 }
00375
00376 for (index i = 0; i < basis.get_n(); i++)
00377 det = det * B(i, basis[i]);
00378 det = det/det_denom;
00379 }
00380
00381
00382 template <class int_item>
00383 int_item det_int(const matrix<int_item>& A)
00384 {
00385
00386 if (A.get_m() != A.get_n())
00387 {
00388 matrix_error (st_matrix_must_be_square,
00389 "matrix must be square");
00390 return 0;
00391 }
00392
00393 matrix<int_item> B;
00394 matrix<int_item> Q;
00395 vector<index> basis;
00396 int_item d;
00397
00398 rref_int(A, B, Q, basis, d, 0, cout);
00399
00400 return d;
00401 }
00402
00403 #ifdef no_default_arguments
00404 template <class int_item>
00405 void smith(const matrix<int_item>& A,
00406 matrix<int_item>& B, matrix<int_item>& P,
00407 matrix<int_item>& Q, index& rank, int_item& det,
00408 int movie, ostr& movie_stream)
00409 #else
00410 template <class int_item>
00411 void smith(const matrix<int_item>& A,
00412 matrix<int_item>& B, matrix<int_item>& P,
00413 matrix<int_item>& Q, index& rank, int_item& det,
00414 int movie, ostr& movie_stream)
00415 #endif
00416 {
00417 det = 1;
00418 index m = A.get_m();
00419 index n = A.get_n();
00420 Q = eye(int_item(1), n);
00421 P = eye(int_item(1), m);
00422 B = A;
00423
00424 if (movie)
00425 {
00426 movie_stream << "Smith's normal diagonal form\n";
00427 print3(movie_stream, Q, P, B);
00428 }
00429
00430 for (index corner = 0; corner < m && corner < n; corner++)
00431 {
00432 index new_i, new_j, i, j;
00433 int_item pivot_item;
00434
00435 if (!find_pivot_item(B, corner, new_i, new_j))
00436
00437 break;
00438
00439 if (movie)
00440 movie_stream << "The smallest (in absolute value) non-zero entry: ("
00441 << new_i << ", " << new_j << ")\n";
00442
00443 while(1)
00444 {
00445 do{
00446 i = new_i;
00447 j = new_j;
00448 pivot_item = B(i, j);
00449
00450 if (movie)
00451 movie_stream << "Pivoting item: (" << i << ", " << j << ")\n";
00452
00453 for (index jj = 0; jj < n; jj++)
00454 if (jj != j)
00455 {
00456 int_item p = quotient(B(i, jj), pivot_item);
00457 B.add_col(jj, j, -p);
00458 Q.add_col(jj, j, -p);
00459 }
00460 for (index ii = 0; ii < m; ii++)
00461 if (ii != i)
00462 {
00463 int_item p = quotient(B(ii, j), pivot_item);
00464 B.add_row(ii, i, -p);
00465 P.add_row(ii, i, -p);
00466 }
00467
00468 if (movie)
00469 {
00470 movie_stream << "After pivoting step:\n";
00471 print3(movie_stream, Q, P, B);
00472 }
00473
00474 new_j = find_pivot_item_in_row(B, corner, i);
00475 if(new_j == j)
00476 new_i = find_pivot_item_in_col(B, corner, j);
00477 else
00478 new_i = i;
00479 }while (i != new_i || j != new_j);
00480
00481
00482
00483 index k, l;
00484 if (!find_non_divisor_of_item(B, corner, pivot_item, k, l))
00485 break;
00486
00487 B.add_row(i, k, 1);
00488 P.add_row(i, k, 1);
00489
00490 if (movie)
00491 {
00492 movie_stream << "Non-divisor entry: (" << k << ", " << l << ")\n";
00493 movie_stream << "Add the " << k << "-th row to the " << i
00494 << "-th row\n";
00495 print3(movie_stream, Q, P, B);
00496 }
00497
00498 }
00499
00500 if (i != corner)
00501 {
00502 B.swap_rows(i, corner);
00503 P.swap_rows(i, corner);
00504 det = -det;
00505 }
00506 if (j != corner)
00507 {
00508 B.swap_cols(j, corner);
00509 Q.swap_cols(j, corner);
00510 det = -det;
00511 }
00512 if (B(corner, corner) < 0)
00513 {
00514 B.mult_row(corner, -1);
00515 P.mult_row(corner, -1);
00516 det = -det;
00517 }
00518 if (movie)
00519 print3(movie_stream, Q, P, B);
00520
00521 }
00522
00523 rank = 0;
00524
00525 for (index i = 0; i < m && i < n; i++)
00526 if(B(i, i) != 0)
00527 {
00528 det = det * B(i, i);
00529 rank++;
00530 }
00531 else
00532 break;
00533 }
00534
00535 template <class int_item>
00536 int find_pivot_item(matrix<int_item>& B, index corner,
00537 index& min_i, index& min_j)
00538 {
00539 int_item min = 0;
00540 for (index i = corner; i < B.get_m(); i++)
00541 for (index j = corner; j < B.get_n(); j++)
00542 if (B(i, j) != 0 && (min == 0 || absolute_value(B(i, j)) < min))
00543 {
00544 min = absolute_value (B(i, j));
00545 min_i = i;
00546 min_j = j;
00547 }
00548 if (min == 0)
00549 return 0;
00550 return 1;
00551 }
00552
00553
00554 template <class int_item>
00555 index find_pivot_item_in_row(matrix<int_item>& B, index corner,
00556 index i)
00557 {
00558 int_item min = 0;
00559 index min_j = 0;
00560
00561 for (index j = corner; j < B.get_n(); j++)
00562 if (B(i, j) != 0 && (min == 0 || absolute_value(B(i, j)) < min))
00563 {
00564 min = absolute_value (B(i, j));
00565 min_j = j;
00566 }
00567
00568 return min_j;
00569 }
00570
00571 template <class int_item>
00572 index find_pivot_item_in_col(matrix<int_item>& B, index corner,
00573 index j)
00574 {
00575 int_item min = 0;
00576 index min_i = 0;
00577
00578 for (index i = corner; i < B.get_m(); i++)
00579 if (B(i, j) != 0 && (min == 0 || absolute_value(B(i, j)) < min))
00580 {
00581 min = absolute_value (B(i, j));
00582 min_i = i;
00583 }
00584
00585 return min_i;
00586 }
00587
00588 template <class int_item>
00589 int find_non_divisor_of_item(matrix<int_item>& B, index corner,
00590 int_item item, index& k, index& l)
00591 {
00592 for (index i = corner; i < B.get_m(); i++)
00593 for (index j = corner; j < B.get_n(); j++)
00594 if (B(i, j) % item != 0)
00595 {
00596 k = i;
00597 l = j;
00598 return 1;
00599 }
00600 return 0;
00601 }
00602
00603
00604
00605 #endif GAUSS_H_