00001 00030 #include <itpp/base/gf2mat.h> 00031 #include <itpp/base/specmat.h> 00032 #include <itpp/base/matfunc.h> 00033 #include <itpp/base/converters.h> 00034 #include <iostream> 00035 00036 namespace itpp 00037 { 00038 00039 // ====== IMPLEMENTATION OF THE ALIST CLASS ========== 00040 00041 GF2mat_sparse_alist::GF2mat_sparse_alist(const std::string &fname) 00042 : data_ok(false) 00043 { 00044 read(fname); 00045 } 00046 00047 void GF2mat_sparse_alist::read(const std::string &fname) 00048 { 00049 std::ifstream file; 00050 std::string line; 00051 std::stringstream ss; 00052 00053 file.open(fname.c_str()); 00054 it_assert(file.is_open(), 00055 "GF2mat_sparse_alist::read(): Could not open file \"" 00056 << fname << "\" for reading"); 00057 00058 // parse N and M: 00059 getline(file, line); 00060 ss << line; 00061 ss >> N >> M; 00062 it_assert(!ss.fail(), 00063 "GF2mat_sparse_alist::read(): Wrong alist data (N or M)"); 00064 it_assert((N > 0) && (M > 0), 00065 "GF2mat_sparse_alist::read(): Wrong alist data"); 00066 ss.seekg(0, std::ios::end); 00067 ss.clear(); 00068 00069 // parse max_num_n and max_num_m 00070 getline(file, line); 00071 ss << line; 00072 ss >> max_num_n >> max_num_m; 00073 it_assert(!ss.fail(), 00074 "GF2mat_sparse_alist::read(): Wrong alist data (max_num_{n,m})"); 00075 it_assert((max_num_n >= 0) && (max_num_n <= N) && 00076 (max_num_m >= 0) && (max_num_m <= M), 00077 "GF2mat_sparse_alist::read(): Wrong alist data"); 00078 ss.seekg(0, std::ios::end); 00079 ss.clear(); 00080 00081 // parse weight of each column n 00082 num_nlist.set_size(N); 00083 num_nlist.clear(); 00084 getline(file, line); 00085 ss << line; 00086 for (int i = 0; i < N; i++) { 00087 ss >> num_nlist(i); 00088 it_assert(!ss.fail(), 00089 "GF2mat_sparse_alist::read(): Wrong alist data (num_nlist(" 00090 << i << "))"); 00091 it_assert((num_nlist(i) >= 0) && (num_nlist(i) <= M), 00092 "GF2mat_sparse_alist::read(): Wrong alist data (num_nlist(" 00093 << i << "))"); 00094 } 00095 ss.seekg(0, std::ios::end); 00096 ss.clear(); 00097 00098 // parse weight of each row m 00099 num_mlist.set_size(M); 00100 num_mlist.clear(); 00101 getline(file, line); 00102 ss << line; 00103 for (int i = 0; i < M; i++) { 00104 ss >> num_mlist(i); 00105 it_assert(!ss.fail(), 00106 "GF2mat_sparse_alist::read(): Wrong alist data (num_mlist(" 00107 << i << "))"); 00108 it_assert((num_mlist(i) >= 0) && (num_mlist(i) <= N), 00109 "GF2mat_sparse_alist::read(): Wrong alist data (num_mlist(" 00110 << i << "))"); 00111 } 00112 ss.seekg(0, std::ios::end); 00113 ss.clear(); 00114 00115 // parse coordinates in the n direction with non-zero entries 00116 nlist.set_size(N, max_num_n); 00117 nlist.clear(); 00118 for (int i = 0; i < N; i++) { 00119 getline(file, line); 00120 ss << line; 00121 for (int j = 0; j < num_nlist(i); j++) { 00122 ss >> nlist(i, j); 00123 it_assert(!ss.fail(), 00124 "GF2mat_sparse_alist::read(): Wrong alist data (nlist(" 00125 << i << "," << j << "))"); 00126 it_assert((nlist(i, j) >= 0) && (nlist(i, j) <= M), 00127 "GF2mat_sparse_alist::read(): Wrong alist data (nlist(" 00128 << i << "," << j << "))"); 00129 } 00130 ss.seekg(0, std::ios::end); 00131 ss.clear(); 00132 } 00133 00134 // parse coordinates in the m direction with non-zero entries 00135 mlist.set_size(M, max_num_m); 00136 mlist.clear(); 00137 for (int i = 0; i < M; i++) { 00138 getline(file, line); 00139 ss << line; 00140 for (int j = 0; j < num_mlist(i); j++) { 00141 ss >> mlist(i, j); 00142 it_assert(!ss.fail(), 00143 "GF2mat_sparse_alist::read(): Wrong alist data (mlist(" 00144 << i << "," << j << "))"); 00145 it_assert((mlist(i, j) >= 0) && (mlist(i, j) <= N), 00146 "GF2mat_sparse_alist::read(): Wrong alist data (mlist(" 00147 << i << "," << j << "))"); 00148 } 00149 ss.seekg(0, std::ios::end); 00150 ss.clear(); 00151 } 00152 00153 file.close(); 00154 data_ok = true; 00155 } 00156 00157 void GF2mat_sparse_alist::write(const std::string &fname) const 00158 { 00159 it_assert(data_ok, 00160 "GF2mat_sparse_alist::write(): alist data not ready for writing"); 00161 00162 std::ofstream file(fname.c_str(), std::ofstream::out); 00163 it_assert(file.is_open(), 00164 "GF2mat_sparse_alist::write(): Could not open file \"" 00165 << fname << "\" for writing"); 00166 00167 file << N << " " << M << std::endl; 00168 file << max_num_n << " " << max_num_m << std::endl; 00169 00170 for (int i = 0; i < num_nlist.length() - 1; i++) 00171 file << num_nlist(i) << " "; 00172 file << num_nlist(num_nlist.length() - 1) << std::endl; 00173 00174 for (int i = 0; i < num_mlist.length() - 1; i++) 00175 file << num_mlist(i) << " "; 00176 file << num_mlist(num_mlist.length() - 1) << std::endl; 00177 00178 for (int i = 0; i < N; i++) { 00179 for (int j = 0; j < num_nlist(i) - 1; j++) 00180 file << nlist(i, j) << " "; 00181 file << nlist(i, num_nlist(i) - 1) << std::endl; 00182 } 00183 00184 for (int i = 0; i < M; i++) { 00185 for (int j = 0; j < num_mlist(i) - 1; j++) 00186 file << mlist(i, j) << " "; 00187 file << mlist(i, num_mlist(i) - 1) << std::endl; 00188 } 00189 00190 file.close(); 00191 } 00192 00193 00194 GF2mat_sparse GF2mat_sparse_alist::to_sparse(bool transpose) const 00195 { 00196 GF2mat_sparse sbmat(M, N, max_num_m); 00197 00198 for (int i = 0; i < M; i++) { 00199 for (int j = 0; j < num_mlist(i); j++) { 00200 sbmat.set_new(i, mlist(i, j) - 1, bin(1)); 00201 } 00202 } 00203 sbmat.compact(); 00204 00205 if (transpose) { 00206 return sbmat.transpose(); 00207 } 00208 else { 00209 return sbmat; 00210 } 00211 } 00212 00213 00214 // ---------------------------------------------------------------------- 00215 // WARNING: This method is very slow. Sparse_Mat has to be extended with 00216 // some extra functions to improve the performance of this. 00217 // ---------------------------------------------------------------------- 00218 void GF2mat_sparse_alist::from_sparse(const GF2mat_sparse &sbmat, 00219 bool transpose) 00220 { 00221 if (transpose) { 00222 from_sparse(sbmat.transpose(), false); 00223 } 00224 else { 00225 // check matrix dimension 00226 M = sbmat.rows(); 00227 N = sbmat.cols(); 00228 00229 num_mlist.set_size(M); 00230 num_nlist.set_size(N); 00231 00232 // fill mlist matrix, num_mlist vector and max_num_m 00233 mlist.set_size(0, 0); 00234 int tmp_cols = 0; // initial number of allocated columns 00235 for (int i = 0; i < M; i++) { 00236 ivec temp_row(0); 00237 for (int j = 0; j < N; j++) { 00238 if (sbmat(i, j) == bin(1)) { 00239 temp_row = concat(temp_row, j + 1); 00240 } 00241 } 00242 int trs = temp_row.size(); 00243 if (trs > tmp_cols) { 00244 tmp_cols = trs; 00245 mlist.set_size(M, tmp_cols, true); 00246 } 00247 else if (trs < tmp_cols) { 00248 temp_row.set_size(tmp_cols, true); 00249 } 00250 mlist.set_row(i, temp_row); 00251 num_mlist(i) = trs; 00252 } 00253 max_num_m = max(num_mlist); 00254 00255 // fill nlist matrix, num_nlist vector and max_num_n 00256 nlist.set_size(0, 0); 00257 tmp_cols = 0; // initial number of allocated columns 00258 for (int j = 0; j < N; j++) { 00259 ivec temp_row = sbmat.get_col(j).get_nz_indices() + 1; 00260 int trs = temp_row.size(); 00261 if (trs > tmp_cols) { 00262 tmp_cols = trs; 00263 nlist.set_size(N, tmp_cols, true); 00264 } 00265 else if (trs < tmp_cols) { 00266 temp_row.set_size(tmp_cols, true); 00267 } 00268 nlist.set_row(j, temp_row); 00269 num_nlist(j) = trs; 00270 } 00271 max_num_n = max(num_nlist); 00272 00273 data_ok = true; 00274 } 00275 } 00276 00277 00278 // ---------------------------------------------------------------------- 00279 // Implementation of a dense GF2 matrix class 00280 // ---------------------------------------------------------------------- 00281 00282 GF2mat::GF2mat(int i, int j): nrows(i), ncols(j), 00283 nwords((j >> shift_divisor) + 1) 00284 { 00285 data.set_size(nrows, nwords); 00286 data.clear(); 00287 } 00288 00289 GF2mat::GF2mat(): nrows(1), ncols(1), nwords(1) 00290 { 00291 data.set_size(nrows, nwords); 00292 data.clear(); 00293 } 00294 00295 GF2mat::GF2mat(const bvec &x, bool is_column) 00296 { 00297 if (is_column) { // create column vector 00298 nrows = length(x); 00299 ncols = 1; 00300 nwords = 1; 00301 data.set_size(nrows, nwords); 00302 data.clear(); 00303 for (int i = 0; i < nrows; i++) { 00304 set(i, 0, x(i)); 00305 } 00306 } 00307 else { // create row vector 00308 nrows = 1; 00309 ncols = length(x); 00310 nwords = (ncols >> shift_divisor) + 1; 00311 data.set_size(nrows, nwords); 00312 data.clear(); 00313 for (int i = 0; i < ncols; i++) { 00314 set(0, i, x(i)); 00315 } 00316 } 00317 } 00318 00319 00320 GF2mat::GF2mat(const bmat &X): nrows(X.rows()), ncols(X.cols()) 00321 { 00322 nwords = (ncols >> shift_divisor) + 1; 00323 data.set_size(nrows, nwords); 00324 data.clear(); 00325 for (int i = 0; i < nrows; i++) { 00326 for (int j = 0; j < ncols; j++) { 00327 set(i, j, X(i, j)); 00328 } 00329 } 00330 } 00331 00332 00333 GF2mat::GF2mat(const GF2mat_sparse &X) 00334 { 00335 nrows = X.rows(); 00336 ncols = X.cols(); 00337 nwords = (ncols >> shift_divisor) + 1; 00338 data.set_size(nrows, nwords); 00339 for (int i = 0; i < nrows; i++) { 00340 for (int j = 0; j < nwords; j++) { 00341 data(i, j) = 0; 00342 } 00343 } 00344 00345 for (int j = 0; j < ncols; j++) { 00346 for (int i = 0; i < X.get_col(j).nnz(); i++) { 00347 bin b = X.get_col(j).get_nz_data(i); 00348 set(X.get_col(j).get_nz_index(i), j, b); 00349 } 00350 } 00351 } 00352 00353 GF2mat::GF2mat(const GF2mat_sparse &X, int m1, int n1, int m2, int n2) 00354 { 00355 it_assert(X.rows() > m2, "GF2mat(): indexes out of range"); 00356 it_assert(X.cols() > n2, "GF2mat(): indexes out of range"); 00357 it_assert(m1 >= 0 && n1 >= 0 && m2 >= m1 && n2 >= n1, 00358 "GF2mat::GF2mat(): indexes out of range"); 00359 00360 nrows = m2 - m1 + 1; 00361 ncols = n2 - n1 + 1; 00362 nwords = (ncols >> shift_divisor) + 1; 00363 data.set_size(nrows, nwords); 00364 00365 for (int i = 0; i < nrows; i++) { 00366 for (int j = 0; j < nwords; j++) { 00367 data(i, j) = 0; 00368 } 00369 } 00370 00371 for (int i = 0; i < nrows; i++) { 00372 for (int j = 0; j < ncols; j++) { 00373 bin b = X(i + m1, j + n1); 00374 set(i, j, b); 00375 } 00376 } 00377 } 00378 00379 00380 GF2mat::GF2mat(const GF2mat_sparse &X, const ivec &columns) 00381 { 00382 it_assert(X.cols() > max(columns), 00383 "GF2mat::GF2mat(): index out of range"); 00384 it_assert(min(columns) >= 0, 00385 "GF2mat::GF2mat(): column index must be positive"); 00386 00387 nrows = X.rows(); 00388 ncols = length(columns); 00389 nwords = (ncols >> shift_divisor) + 1; 00390 data.set_size(nrows, nwords); 00391 00392 for (int i = 0; i < nrows; i++) { 00393 for (int j = 0; j < nwords; j++) { 00394 data(i, j) = 0; 00395 } 00396 } 00397 00398 for (int j = 0; j < ncols; j++) { 00399 for (int i = 0; i < X.get_col(columns(j)).nnz(); i++) { 00400 bin b = X.get_col(columns(j)).get_nz_data(i); 00401 set(X.get_col(columns(j)).get_nz_index(i), j, b); 00402 } 00403 } 00404 } 00405 00406 00407 void GF2mat::set_size(int m, int n, bool copy) 00408 { 00409 nrows = m; 00410 ncols = n; 00411 nwords = (ncols >> shift_divisor) + 1; 00412 data.set_size(nrows, nwords, copy); 00413 if (!copy) 00414 data.clear(); 00415 } 00416 00417 00418 GF2mat_sparse GF2mat::sparsify() const 00419 { 00420 GF2mat_sparse Z(nrows, ncols); 00421 for (int i = 0; i < nrows; i++) { 00422 for (int j = 0; j < ncols; j++) { 00423 if (get(i, j) == 1) { 00424 Z.set(i, j, 1); 00425 } 00426 } 00427 } 00428 00429 return Z; 00430 } 00431 00432 bvec GF2mat::bvecify() const 00433 { 00434 it_assert(nrows == 1 || ncols == 1, 00435 "GF2mat::bvecify() matrix must be a vector"); 00436 int n = (nrows == 1 ? ncols : nrows); 00437 bvec result(n); 00438 if (nrows == 1) { 00439 for (int i = 0; i < n; i++) { 00440 result(i) = get(0, i); 00441 } 00442 } 00443 else { 00444 for (int i = 0; i < n; i++) { 00445 result(i) = get(i, 0); 00446 } 00447 } 00448 return result; 00449 } 00450 00451 00452 void GF2mat::set_row(int i, bvec x) 00453 { 00454 it_assert(length(x) == ncols, 00455 "GF2mat::set_row(): dimension mismatch"); 00456 for (int j = 0; j < ncols; j++) { 00457 set(i, j, x(j)); 00458 } 00459 } 00460 00461 void GF2mat::set_col(int j, bvec x) 00462 { 00463 it_assert(length(x) == nrows, 00464 "GF2mat::set_col(): dimension mismatch"); 00465 for (int i = 0; i < nrows; i++) { 00466 set(i, j, x(i)); 00467 } 00468 } 00469 00470 00471 GF2mat gf2dense_eye(int m) 00472 { 00473 GF2mat Z(m, m); 00474 for (int i = 0; i < m; i++) { 00475 Z.set(i, i, 1); 00476 } 00477 return Z; 00478 } 00479 00480 GF2mat GF2mat::get_submatrix(int m1, int n1, int m2, int n2) const 00481 { 00482 it_assert(m1 >= 0 && n1 >= 0 && m2 >= m1 && n2 >= n1 00483 && m2 < nrows && n2 < ncols, 00484 "GF2mat::get_submatrix() index out of range"); 00485 GF2mat result(m2 - m1 + 1, n2 - n1 + 1); 00486 00487 for (int i = m1; i <= m2; i++) { 00488 for (int j = n1; j <= n2; j++) { 00489 result.set(i - m1, j - n1, get(i, j)); 00490 } 00491 } 00492 00493 return result; 00494 } 00495 00496 00497 GF2mat GF2mat::concatenate_vertical(const GF2mat &X) const 00498 { 00499 it_assert(X.ncols == ncols, 00500 "GF2mat::concatenate_vertical(): dimension mismatch"); 00501 it_assert(X.nwords == nwords, 00502 "GF2mat::concatenate_vertical(): dimension mismatch"); 00503 00504 GF2mat result(nrows + X.nrows, ncols); 00505 for (int i = 0; i < nrows; i++) { 00506 for (int j = 0; j < nwords; j++) { 00507 result.data(i, j) = data(i, j); 00508 } 00509 } 00510 00511 for (int i = 0; i < X.nrows; i++) { 00512 for (int j = 0; j < nwords; j++) { 00513 result.data(i + nrows, j) = X.data(i, j); 00514 } 00515 } 00516 00517 return result; 00518 } 00519 00520 GF2mat GF2mat::concatenate_horizontal(const GF2mat &X) const 00521 { 00522 it_assert(X.nrows == nrows, 00523 "GF2mat::concatenate_horizontal(): dimension mismatch"); 00524 00525 GF2mat result(nrows, X.ncols + ncols); 00526 for (int i = 0; i < nrows; i++) { 00527 for (int j = 0; j < ncols; j++) { 00528 result.set(i, j, get(i, j)); 00529 } 00530 } 00531 00532 for (int i = 0; i < nrows; i++) { 00533 for (int j = 0; j < X.ncols; j++) { 00534 result.set(i, j + ncols, X.get(i, j)); 00535 } 00536 } 00537 00538 return result; 00539 } 00540 00541 bvec GF2mat::get_row(int i) const 00542 { 00543 bvec result(ncols); 00544 for (int j = 0; j < ncols; j++) { 00545 result(j) = get(i, j); 00546 } 00547 00548 return result; 00549 } 00550 00551 bvec GF2mat::get_col(int j) const 00552 { 00553 bvec result(nrows); 00554 for (int i = 0; i < nrows; i++) { 00555 result(i) = get(i, j); 00556 } 00557 00558 return result; 00559 } 00560 00561 00562 int GF2mat::T_fact(GF2mat &T, GF2mat &U, ivec &perm) const 00563 { 00564 T = gf2dense_eye(nrows); 00565 U = *this; 00566 00567 perm = zeros_i(ncols); 00568 for (int i = 0; i < ncols; i++) { 00569 perm(i) = i; 00570 } 00571 00572 if (nrows > 250) { // avoid cluttering output ... 00573 it_info_debug("Performing T-factorization of GF(2) matrix... rows: " 00574 << nrows << " cols: " << ncols << " .... " << std::endl); 00575 } 00576 int pdone = 0; 00577 for (int j = 0; j < nrows; j++) { 00578 // Now working on diagonal element j,j 00579 // First try find a row with a 1 in column i 00580 int i1, j1; 00581 for (i1 = j; i1 < nrows; i1++) { 00582 for (j1 = j; j1 < ncols; j1++) { 00583 if (U.get(i1, j1) == 1) { goto found; } 00584 } 00585 } 00586 00587 return j; 00588 00589 found: 00590 U.swap_rows(i1, j); 00591 T.swap_rows(i1, j); 00592 U.swap_cols(j1, j); 00593 00594 int temp = perm(j); 00595 perm(j) = perm(j1); 00596 perm(j1) = temp; 00597 00598 // now subtract row i from remaining rows 00599 for (int i1 = j + 1; i1 < nrows; i1++) { 00600 if (U.get(i1, j) == 1) { 00601 int ptemp = floor_i(100.0 * (i1 + j * nrows) / (nrows * nrows)); 00602 if (nrows > 250 && ptemp > pdone + 10) { 00603 it_info_debug(ptemp << "% done."); 00604 pdone = ptemp; 00605 } 00606 U.add_rows(i1, j); 00607 T.add_rows(i1, j); 00608 } 00609 } 00610 } 00611 return nrows; 00612 } 00613 00614 00615 int GF2mat::T_fact_update_bitflip(GF2mat &T, GF2mat &U, 00616 ivec &perm, int rank, 00617 int r, int c) const 00618 { 00619 // First, update U (before re-triangulization) 00620 int c0 = c; 00621 for (c = 0; c < ncols; c++) { 00622 if (perm(c) == c0) { 00623 goto foundidx; 00624 } 00625 } 00626 it_error("GF2mat::T_fact_update_bitflip() - internal error"); 00627 00628 foundidx: 00629 for (int i = 0; i < nrows; i++) { 00630 if (T.get(i, r) == 1) { 00631 U.addto_element(i, c, 1); 00632 } 00633 } 00634 00635 // first move column c to the end 00636 bvec lastcol = U.get_col(c); 00637 int temp_perm = perm(c); 00638 for (int j = c; j < ncols - 1; j++) { 00639 U.set_col(j, U.get_col(j + 1)); 00640 perm(j) = perm(j + 1); 00641 } 00642 U.set_col(ncols - 1, lastcol); 00643 perm(ncols - 1) = temp_perm; 00644 00645 // then, if the matrix is tall, also move row c to the end 00646 if (nrows >= ncols) { 00647 bvec lastrow_U = U.get_row(c); 00648 bvec lastrow_T = T.get_row(c); 00649 for (int i = c; i < nrows - 1; i++) { 00650 U.set_row(i, U.get_row(i + 1)); 00651 T.set_row(i, T.get_row(i + 1)); 00652 } 00653 U.set_row(nrows - 1, lastrow_U); 00654 T.set_row(nrows - 1, lastrow_T); 00655 00656 // Do Gaussian elimination on the last row 00657 for (int j = c; j < ncols; j++) { 00658 if (U.get(nrows - 1, j) == 1) { 00659 U.add_rows(nrows - 1, j); 00660 T.add_rows(nrows - 1, j); 00661 } 00662 } 00663 } 00664 00665 // Now, continue T-factorization from the point (rank-1,rank-1) 00666 for (int j = rank - 1; j < nrows; j++) { 00667 int i1, j1; 00668 for (i1 = j; i1 < nrows; i1++) { 00669 for (j1 = j; j1 < ncols; j1++) { 00670 if (U.get(i1, j1) == 1) { 00671 goto found; 00672 } 00673 } 00674 } 00675 00676 return j; 00677 00678 found: 00679 U.swap_rows(i1, j); 00680 T.swap_rows(i1, j); 00681 U.swap_cols(j1, j); 00682 00683 int temp = perm(j); 00684 perm(j) = perm(j1); 00685 perm(j1) = temp; 00686 00687 for (int i1 = j + 1; i1 < nrows; i1++) { 00688 if (U.get(i1, j) == 1) { 00689 U.add_rows(i1, j); 00690 T.add_rows(i1, j); 00691 } 00692 } 00693 } 00694 00695 return nrows; 00696 } 00697 00698 bool GF2mat::T_fact_update_addcol(GF2mat &T, GF2mat &U, 00699 ivec &perm, bvec newcol) const 00700 { 00701 int i0 = T.rows(); 00702 int j0 = U.cols(); 00703 it_assert(j0 > 0, "GF2mat::T_fact_update_addcol(): dimension mismatch"); 00704 it_assert(i0 == T.cols(), 00705 "GF2mat::T_fact_update_addcol(): dimension mismatch"); 00706 it_assert(i0 == U.rows(), 00707 "GF2mat::T_fact_update_addcol(): dimension mismatch"); 00708 it_assert(length(perm) == j0, 00709 "GF2mat::T_fact_update_addcol(): dimension mismatch"); 00710 it_assert(U.get(j0 - 1, j0 - 1) == 1, 00711 "GF2mat::T_fact_update_addcol(): dimension mismatch"); 00712 // The following test is VERY TIME-CONSUMING 00713 it_assert_debug(U.row_rank() == j0, 00714 "GF2mat::T_fact_update_addcol(): factorization has incorrect rank"); 00715 00716 bvec z = T * newcol; 00717 GF2mat Utemp = U.concatenate_horizontal(GF2mat(z, 1)); 00718 00719 // start working on position (j0,j0) 00720 int i; 00721 for (i = j0; i < i0; i++) { 00722 if (Utemp.get(i, j0) == 1) { 00723 goto found; 00724 } 00725 } 00726 return (false); // adding the new column would not improve the rank 00727 00728 found: 00729 perm.set_length(j0 + 1, true); 00730 perm(j0) = j0; 00731 00732 Utemp.swap_rows(i, j0); 00733 T.swap_rows(i, j0); 00734 00735 for (int i1 = j0 + 1; i1 < i0; i1++) { 00736 if (Utemp.get(i1, j0) == 1) { 00737 Utemp.add_rows(i1, j0); 00738 T.add_rows(i1, j0); 00739 } 00740 } 00741 00742 U = Utemp; 00743 return (true); // the new column was successfully added 00744 } 00745 00746 00747 00748 00749 GF2mat GF2mat::inverse() const 00750 { 00751 it_assert(nrows == ncols, "GF2mat::inverse(): Matrix must be square"); 00752 00753 // first compute the T-factorization 00754 GF2mat T, U; 00755 ivec perm; 00756 int rank = T_fact(T, U, perm); 00757 it_assert(rank == ncols, "GF2mat::inverse(): Matrix is not full rank"); 00758 00759 // backward substitution 00760 for (int i = ncols - 2; i >= 0; i--) { 00761 for (int j = ncols - 1; j > i; j--) { 00762 if (U.get(i, j) == 1) { 00763 U.add_rows(i, j); 00764 T.add_rows(i, j); 00765 } 00766 } 00767 } 00768 T.permute_rows(perm, 1); 00769 return T; 00770 } 00771 00772 int GF2mat::row_rank() const 00773 { 00774 GF2mat T, U; 00775 ivec perm; 00776 int rank = T_fact(T, U, perm); 00777 return rank; 00778 } 00779 00780 bool GF2mat::is_zero() const 00781 { 00782 for (int i = 0; i < nrows; i++) { 00783 for (int j = 0; j < nwords; j++) { 00784 if (data(i, j) != 0) { 00785 return false; 00786 } 00787 } 00788 } 00789 return true; 00790 } 00791 00792 bool GF2mat::operator==(const GF2mat &X) const 00793 { 00794 if (X.nrows != nrows) { return false; } 00795 if (X.ncols != ncols) { return false; } 00796 it_assert(X.nwords == nwords, "GF2mat::operator==() dimension mismatch"); 00797 00798 for (int i = 0; i < nrows; i++) { 00799 for (int j = 0; j < nwords; j++) { 00800 // if (X.get(i,j)!=get(i,j)) { 00801 if (X.data(i, j) != data(i, j)) { 00802 return false; 00803 } 00804 } 00805 } 00806 return true; 00807 } 00808 00809 00810 void GF2mat::add_rows(int i, int j) 00811 { 00812 it_assert(i >= 0 && i < nrows, "GF2mat::add_rows(): out of range"); 00813 it_assert(j >= 0 && j < nrows, "GF2mat::add_rows(): out of range"); 00814 for (int k = 0; k < nwords; k++) { 00815 data(i, k) ^= data(j, k); 00816 } 00817 } 00818 00819 void GF2mat::swap_rows(int i, int j) 00820 { 00821 it_assert(i >= 0 && i < nrows, "GF2mat::swap_rows(): index out of range"); 00822 it_assert(j >= 0 && j < nrows, "GF2mat::swap_rows(): index out of range"); 00823 for (int k = 0; k < nwords; k++) { 00824 unsigned char temp = data(i, k); 00825 data(i, k) = data(j, k); 00826 data(j, k) = temp; 00827 } 00828 } 00829 00830 void GF2mat::swap_cols(int i, int j) 00831 { 00832 it_assert(i >= 0 && i < ncols, "GF2mat::swap_cols(): index out of range"); 00833 it_assert(j >= 0 && j < ncols, "GF2mat::swap_cols(): index out of range"); 00834 bvec temp = get_col(i); 00835 set_col(i, get_col(j)); 00836 set_col(j, temp); 00837 } 00838 00839 00840 void GF2mat::operator=(const GF2mat &X) 00841 { 00842 nrows = X.nrows; 00843 ncols = X.ncols; 00844 nwords = X.nwords; 00845 data = X.data; 00846 } 00847 00848 GF2mat operator*(const GF2mat &X, const GF2mat &Y) 00849 { 00850 it_assert(X.ncols == Y.nrows, "GF2mat::operator*(): dimension mismatch"); 00851 it_assert(X.nwords > 0, "Gfmat::operator*(): dimension mismatch"); 00852 it_assert(Y.nwords > 0, "Gfmat::operator*(): dimension mismatch"); 00853 00854 /* 00855 // this can be done more efficiently? 00856 GF2mat result(X.nrows,Y.ncols); 00857 for (int i=0; i<X.nrows; i++) { 00858 for (int j=0; j<Y.ncols; j++) { 00859 bin b=0; 00860 for (int k=0; k<X.ncols; k++) { 00861 bin x = X.get(i,k); 00862 bin y = Y.get(k,j); 00863 b ^= (x&y); 00864 } 00865 result.set(i,j,b); 00866 } 00867 } 00868 return result; */ 00869 00870 // is this better? 00871 return mult_trans(X, Y.transpose()); 00872 } 00873 00874 bvec operator*(const GF2mat &X, const bvec &y) 00875 { 00876 it_assert(length(y) == X.ncols, "GF2mat::operator*(): dimension mismatch"); 00877 it_assert(X.nwords > 0, "Gfmat::operator*(): dimension mismatch"); 00878 00879 /* 00880 // this can be done more efficiently? 00881 bvec result = zeros_b(X.nrows); 00882 for (int i=0; i<X.nrows; i++) { 00883 // do the nwords-1 data columns first 00884 for (int j=0; j<X.nwords-1; j++) { 00885 int ind = j<<shift_divisor; 00886 unsigned char r=X.data(i,j); 00887 while (r) { 00888 result(i) ^= (r & y(ind)); 00889 r >>= 1; 00890 ind++; 00891 } 00892 } 00893 // do the last column separately 00894 for (int j=(X.nwords-1)<<shift_divisor; j<X.ncols; j++) { 00895 result(i) ^= (X.get(i,j) & y(j)); 00896 } 00897 } 00898 return result; */ 00899 00900 // is this better? 00901 return (mult_trans(X, GF2mat(y, 0))).bvecify(); 00902 } 00903 00904 GF2mat mult_trans(const GF2mat &X, const GF2mat &Y) 00905 { 00906 it_assert(X.ncols == Y.ncols, "GF2mat::mult_trans(): dimension mismatch"); 00907 it_assert(X.nwords > 0, "GF2mat::mult_trans(): dimension mismatch"); 00908 it_assert(Y.nwords > 0, "GF2mat::mult_trans(): dimension mismatch"); 00909 it_assert(X.nwords == Y.nwords, "GF2mat::mult_trans(): dimension mismatch"); 00910 00911 GF2mat result(X.nrows, Y.nrows); 00912 00913 for (int i = 0; i < X.nrows; i++) { 00914 for (int j = 0; j < Y.nrows; j++) { 00915 bin b = 0; 00916 int kindx = i; 00917 int kindy = j; 00918 for (int k = 0; k < X.nwords; k++) { 00919 unsigned char r = X.data(kindx) & Y.data(kindy); 00920 /* The following can be speeded up by using a small lookup 00921 table for the number of ones and shift only a few times (or 00922 not at all if table is large) */ 00923 while (r) { 00924 b ^= r & 1; 00925 r >>= 1; 00926 }; 00927 kindx += X.nrows; 00928 kindy += Y.nrows; 00929 } 00930 result.set(i, j, b); 00931 } 00932 } 00933 return result; 00934 } 00935 00936 GF2mat GF2mat::transpose() const 00937 { 00938 // CAN BE SPEEDED UP 00939 GF2mat result(ncols, nrows); 00940 00941 for (int i = 0; i < nrows; i++) { 00942 for (int j = 0; j < ncols; j++) { 00943 result.set(j, i, get(i, j)); 00944 } 00945 } 00946 return result; 00947 } 00948 00949 GF2mat operator+(const GF2mat &X, const GF2mat &Y) 00950 { 00951 it_assert(X.nrows == Y.nrows, "GF2mat::operator+(): dimension mismatch"); 00952 it_assert(X.ncols == Y.ncols, "GF2mat::operator+(): dimension mismatch"); 00953 it_assert(X.nwords == Y.nwords, "GF2mat::operator+(): dimension mismatch"); 00954 GF2mat result(X.nrows, X.ncols); 00955 00956 for (int i = 0; i < X.nrows; i++) { 00957 for (int j = 0; j < X.nwords; j++) { 00958 result.data(i, j) = X.data(i, j) ^ Y.data(i, j); 00959 } 00960 } 00961 00962 return result; 00963 } 00964 00965 void GF2mat::permute_cols(ivec &perm, bool I) 00966 { 00967 it_assert(length(perm) == ncols, 00968 "GF2mat::permute_cols(): dimensions do not match"); 00969 00970 GF2mat temp = (*this); 00971 for (int j = 0; j < ncols; j++) { 00972 if (I == 0) { 00973 set_col(j, temp.get_col(perm(j))); 00974 } 00975 else { 00976 set_col(perm(j), temp.get_col(j)); 00977 } 00978 } 00979 } 00980 00981 void GF2mat::permute_rows(ivec &perm, bool I) 00982 { 00983 it_assert(length(perm) == nrows, 00984 "GF2mat::permute_rows(): dimensions do not match"); 00985 00986 GF2mat temp = (*this); 00987 for (int i = 0; i < nrows; i++) { 00988 if (I == 0) { 00989 for (int j = 0; j < nwords; j++) { 00990 data(i, j) = temp.data(perm(i), j); 00991 } 00992 } 00993 else { 00994 for (int j = 0; j < nwords; j++) { 00995 data(perm(i), j) = temp.data(i, j); 00996 } 00997 } 00998 } 00999 } 01000 01001 01002 std::ostream &operator<<(std::ostream &os, const GF2mat &X) 01003 { 01004 int i, j; 01005 os << "---- GF(2) matrix of dimension " << X.nrows << "*" << X.ncols 01006 << " -- Density: " << X.density() << " ----" << std::endl; 01007 01008 for (i = 0; i < X.nrows; i++) { 01009 os << " "; 01010 for (j = 0; j < X.ncols; j++) { 01011 os << X.get(i, j) << " "; 01012 } 01013 os << std::endl; 01014 } 01015 01016 return os; 01017 } 01018 01019 double GF2mat::density() const 01020 { 01021 int no_of_ones = 0; 01022 01023 for (int i = 0; i < nrows; i++) { 01024 for (int j = 0; j < ncols; j++) { 01025 no_of_ones += (get(i, j) == 1 ? 1 : 0); 01026 } 01027 } 01028 01029 return ((double) no_of_ones) / (nrows*ncols); 01030 } 01031 01032 01033 it_file &operator<<(it_file &f, const GF2mat &X) 01034 { 01035 // 3 64-bit unsigned words for: nrows, ncols and nwords + rest for char data 01036 uint64_t bytecount = 3 * sizeof(uint64_t) 01037 + X.nrows * X.nwords * sizeof(char); 01038 f.write_data_header("GF2mat", bytecount); 01039 01040 f.low_level_write(static_cast<uint64_t>(X.nrows)); 01041 f.low_level_write(static_cast<uint64_t>(X.ncols)); 01042 f.low_level_write(static_cast<uint64_t>(X.nwords)); 01043 for (int i = 0; i < X.nrows; i++) { 01044 for (int j = 0; j < X.nwords; j++) { 01045 f.low_level_write(static_cast<char>(X.data(i, j))); 01046 } 01047 } 01048 return f; 01049 } 01050 01051 it_ifile &operator>>(it_ifile &f, GF2mat &X) 01052 { 01053 it_file::data_header h; 01054 01055 f.read_data_header(h); 01056 if (h.type == "GF2mat") { 01057 uint64_t tmp; 01058 f.low_level_read(tmp); 01059 X.nrows = static_cast<int>(tmp); 01060 f.low_level_read(tmp); 01061 X.ncols = static_cast<int>(tmp); 01062 f.low_level_read(tmp); 01063 X.nwords = static_cast<int>(tmp); 01064 X.data.set_size(X.nrows, X.nwords); 01065 for (int i = 0; i < X.nrows; i++) { 01066 for (int j = 0; j < X.nwords; j++) { 01067 char r; 01068 f.low_level_read(r); 01069 X.data(i, j) = static_cast<unsigned char>(r); 01070 } 01071 } 01072 } 01073 else { 01074 it_error("it_ifile &operator>>() - internal error"); 01075 } 01076 01077 return f; 01078 } 01079 01080 } // namespace itpp 01081
Generated on Tue Nov 23 08:47:55 2010 for IT++ by Doxygen 1.6.1