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