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