?? mat_gf2.c
字號:
#include <NTL/mat_GF2.h>#include <NTL/vec_long.h>#include <NTL/new.h>NTL_START_IMPLmat_GF2::mat_GF2(const mat_GF2& a) { _mat_GF2__numcols = 0; SetDims(a.NumRows(), a.NumCols()); _mat_GF2__rep = a._mat_GF2__rep; } mat_GF2& mat_GF2::operator=(const mat_GF2& a) { SetDims(a.NumRows(), a.NumCols()); _mat_GF2__rep = a._mat_GF2__rep; return *this;} mat_GF2::mat_GF2(INIT_SIZE_TYPE, long n, long m) { _mat_GF2__numcols = 0; SetDims(n, m); } void mat_GF2::kill() { _mat_GF2__numcols = 0; _mat_GF2__rep.kill(); } void mat_GF2::SetDims(long n, long m) { if (n < 0 || m < 0) Error("SetDims: bad args"); if (m != _mat_GF2__numcols) { _mat_GF2__rep.kill(); _mat_GF2__numcols = m; } long oldmax = _mat_GF2__rep.MaxLength(); long i; _mat_GF2__rep.SetLength(n); for (i = oldmax; i < n; i++) _mat_GF2__rep[i].FixLength(m); } void conv(mat_GF2& x, const vec_vec_GF2& a) { long n = a.length(); if (n == 0) { x.SetDims(0, 0); return; } long m = a[0].length(); long i; for (i = 1; i < n; i++) if (a[i].length() != m) Error("nonrectangular matrix"); x.SetDims(n, m); for (i = 0; i < n; i++) x[i] = a[i]; } void swap(mat_GF2& X, mat_GF2& Y) { swap(X._mat_GF2__numcols, Y._mat_GF2__numcols); swap(X._mat_GF2__rep, Y._mat_GF2__rep); } long operator==(const mat_GF2& a, const mat_GF2& b) { if (a.NumCols() != b.NumCols()) return 0; if (a.NumRows() != b.NumRows()) return 0; long n = a.NumRows(); long i; for (i = 0; i < n; i++) if (a[i] != b[i]) return 0; return 1; } long operator!=(const mat_GF2& a, const mat_GF2& b) { return !(a == b); } istream& operator>>(istream& s, mat_GF2& x) { vec_vec_GF2 buf; s >> buf; conv(x, buf); return s; } ostream& operator<<(ostream& s, const mat_GF2& a) { long n = a.NumRows(); long i; s << "["; for (i = 0; i < n; i++) { s << a[i]; s << "\n"; } s << "]"; return s; } void add(mat_GF2& X, const mat_GF2& A, const mat_GF2& B) { long n = A.NumRows(); long m = A.NumCols(); if (B.NumRows() != n || B.NumCols() != m) Error("matrix add: dimension mismatch"); X.SetDims(n, m); long mw = (m + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG; long i; for (i = 0; i < n; i++) { _ntl_ulong *xp = X[i].rep.elts(); const _ntl_ulong *ap = A[i].rep.elts(); const _ntl_ulong *bp = B[i].rep.elts(); long j; for (j = 0; j < mw; j++) xp[j] = ap[j] ^ bp[j]; }} staticvoid mul_aux(vec_GF2& x, const mat_GF2& A, const vec_GF2& b) { long n = A.NumRows(); long l = A.NumCols(); if (l != b.length()) Error("matrix mul: dimension mismatch"); x.SetLength(n); long i; for (i = 0; i < n; i++) { x.put(i, A[i] * b); } } void mul(vec_GF2& x, const mat_GF2& A, const vec_GF2& b) { if (&b == &x || A.position(b) != -1) { vec_GF2 tmp; mul_aux(tmp, A, b); x = tmp; } else mul_aux(x, A, b);} staticvoid mul_aux(vec_GF2& x, const vec_GF2& a, const mat_GF2& B) { long n = B.NumRows(); long l = B.NumCols(); if (n != a.length()) Error("matrix mul: dimension mismatch"); x.SetLength(l); clear(x); const _ntl_ulong *ap = a.rep.elts(); _ntl_ulong a_mask = 1; _ntl_ulong *xp = x.rep.elts(); long lw = (l + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG; long i; for (i = 0; i < n; i++) { if (*ap & a_mask) { const _ntl_ulong *bp = B[i].rep.elts(); long j; for (j = 0; j < lw; j++) xp[j] ^= bp[j]; } a_mask <<= 1; if (!a_mask) { a_mask = 1; ap++; } } } void mul(vec_GF2& x, const vec_GF2& a, const mat_GF2& B){ if (&a == &x || B.position(a) != -1) { vec_GF2 tmp; mul_aux(tmp, a, B); x = tmp; } else mul_aux(x, a, B);} void mul_aux(mat_GF2& X, const mat_GF2& A, const mat_GF2& B) { long n = A.NumRows(); long l = A.NumCols(); long m = B.NumCols(); if (l != B.NumRows()) Error("matrix mul: dimension mismatch"); X.SetDims(n, m); long i; for (i = 1; i <= n; i++) { mul_aux(X(i), A(i), B); } } void mul(mat_GF2& X, const mat_GF2& A, const mat_GF2& B) { if (&X == &A || &X == &B) { mat_GF2 tmp; mul_aux(tmp, A, B); X = tmp; } else mul_aux(X, A, B); } void ident(mat_GF2& X, long n) { X.SetDims(n, n); clear(X); long i; for (i = 0; i < n; i++) X.put(i, i, to_GF2(1));} void det(GF2& d, const mat_GF2& M_in){ long k, n; long i, j; long pos; n = M_in.NumRows(); if (M_in.NumCols() != n) Error("determinant: nonsquare matrix"); if (n == 0) { set(d); return; } mat_GF2 M; M = M_in; long wn = (n + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG; for (k = 0; k < n; k++) { long wk = k/NTL_BITS_PER_LONG; long bk = k - wk*NTL_BITS_PER_LONG; _ntl_ulong k_mask = 1UL << bk; pos = -1; for (i = k; i < n; i++) { if (M[i].rep.elts()[wk] & k_mask) { pos = i; break; } } if (pos != -1) { if (k != pos) { swap(M[pos], M[k]); } _ntl_ulong *y = M[k].rep.elts(); for (i = k+1; i < n; i++) { // M[i] = M[i] + M[k]*M[i,k] if (M[i].rep.elts()[wk] & k_mask) { _ntl_ulong *x = M[i].rep.elts(); for (j = wk; j < wn; j++) x[j] ^= y[j]; } } } else { clear(d); return; } } set(d); return;}staticlong IsUnitVector(const vec_GF2& a, long i){ long wi = i/NTL_BITS_PER_LONG; long bi = i - wi*NTL_BITS_PER_LONG; const _ntl_ulong *p = a.rep.elts(); long wdlen = a.rep.length(); long j; for (j = 0; j < wi; j++) if (p[j] != 0) return 0; if (p[wi] != (1UL << bi)) return 0; for (j = wi+1; j < wdlen; j++) if (p[j] != 0) return 0; return 1;}long IsIdent(const mat_GF2& A, long n){ if (A.NumRows() != n || A.NumCols() != n) return 0; if (n == 0) return 1; long i; for (i = 0; i < n; i++) if (!IsUnitVector(A[i], i)) return 0; return 1;}void AddToCol(mat_GF2& x, long j, const vec_GF2& a)// add a to column j of x// ALIAS RESTRICTION: a should not alias any row of x{ long n = x.NumRows(); long m = x.NumCols(); if (a.length() != n || j < 0 || j >= m) Error("AddToCol: bad args"); long wj = j/NTL_BITS_PER_LONG; long bj = j - wj*NTL_BITS_PER_LONG; _ntl_ulong j_mask = 1UL << bj; const _ntl_ulong *ap = a.rep.elts(); _ntl_ulong a_mask = 1; long i; for (i = 0; i < n; i++) { if (*ap & a_mask) x[i].rep.elts()[wj] ^= j_mask; a_mask <<= 1; if (!a_mask) { a_mask = 1; ap++; } }}void transpose_aux(mat_GF2& X, const mat_GF2& A){ long n = A.NumRows(); long m = A.NumCols(); X.SetDims(m, n); clear(X); long i; for (i = 0; i < n; i++) AddToCol(X, i, A[i]);} void transpose(mat_GF2& X, const mat_GF2& A){ if (&X == &A) { mat_GF2 tmp; transpose_aux(tmp, A); X = tmp; } else transpose_aux(X, A);} void solve(GF2& d, vec_GF2& X, const mat_GF2& A, const vec_GF2& b){ long n = A.NumRows(); if (A.NumCols() != n) Error("solve: nonsquare matrix"); if (b.length() != n) Error("solve: dimension mismatch"); if (n == 0) { X.SetLength(0); set(d); return; } long i, j, k, pos; mat_GF2 M; M.SetDims(n, n+1); for (i = 0; i < n; i++) { AddToCol(M, i, A[i]); } AddToCol(M, n, b); long wn = ((n+1) + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG; for (k = 0; k < n; k++) { long wk = k/NTL_BITS_PER_LONG; long bk = k - wk*NTL_BITS_PER_LONG; _ntl_ulong k_mask = 1UL << bk; pos = -1; for (i = k; i < n; i++) { if (M[i].rep.elts()[wk] & k_mask) { pos = i; break; } } if (pos != -1) { if (k != pos) { swap(M[pos], M[k]); } _ntl_ulong *y = M[k].rep.elts(); for (i = k+1; i < n; i++) { // M[i] = M[i] + M[k]*M[i,k] if (M[i].rep.elts()[wk] & k_mask) { _ntl_ulong *x = M[i].rep.elts(); for (j = wk; j < wn; j++) x[j] ^= y[j]; } } } else { clear(d); return; } } vec_GF2 XX; XX.SetLength(n+1); XX.put(n, 1); for (i = n-1; i >= 0; i--) { XX.put(i, XX*M[i]); } XX.SetLength(n); X = XX; set(d); return;}void inv(GF2& d, mat_GF2& X, const mat_GF2& A){ long n = A.NumRows(); if (A.NumCols() != n) Error("solve: nonsquare matrix"); if (n == 0) { X.SetDims(0, 0); set(d); } long i, j, k, pos; mat_GF2 M; M.SetDims(n, 2*n); vec_GF2 aa; aa.SetLength(2*n); for (i = 0; i < n; i++) { aa = A[i]; aa.SetLength(2*n); aa.put(n+i, 1); M[i] = aa; } long wn = ((2*n) + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG; for (k = 0; k < n; k++) { long wk = k/NTL_BITS_PER_LONG; long bk = k - wk*NTL_BITS_PER_LONG; _ntl_ulong k_mask = 1UL << bk; pos = -1; for (i = k; i < n; i++) { if (M[i].rep.elts()[wk] & k_mask) { pos = i; break; } } if (pos != -1) { if (k != pos) { swap(M[pos], M[k]); } _ntl_ulong *y = M[k].rep.elts(); for (i = k+1; i < n; i++) { // M[i] = M[i] + M[k]*M[i,k] if (M[i].rep.elts()[wk] & k_mask) { _ntl_ulong *x = M[i].rep.elts(); for (j = wk; j < wn; j++) x[j] ^= y[j]; } } } else { clear(d); return; } } vec_GF2 XX; XX.SetLength(2*n); X.SetDims(n, n); clear(X); for (j = 0; j < n; j++) { XX.SetLength(n+j+1); clear(XX); XX.put(n+j, to_GF2(1)); for (i = n-1; i >= 0; i--) { XX.put(i, XX*M[i]); } XX.SetLength(n); AddToCol(X, j, XX); } set(d); return;}long gauss(mat_GF2& M, long w){ long k, l; long i, j; long pos; long n = M.NumRows(); long m = M.NumCols(); if (w < 0 || w > m) Error("gauss: bad args"); long wm = (m + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG; l = 0; for (k = 0; k < w && l < n; k++) { long wk = k/NTL_BITS_PER_LONG; long bk = k - wk*NTL_BITS_PER_LONG; _ntl_ulong k_mask = 1UL << bk; pos = -1; for (i = l; i < n; i++) { if (M[i].rep.elts()[wk] & k_mask) { pos = i; break; } } if (pos != -1) { if (l != pos) swap(M[pos], M[l]); _ntl_ulong *y = M[l].rep.elts(); for (i = l+1; i < n; i++) { // M[i] = M[i] + M[l]*M[i,k] if (M[i].rep.elts()[wk] & k_mask) { _ntl_ulong *x = M[i].rep.elts(); for (j = wk; j < wm; j++) x[j] ^= y[j]; } } l++; } } return l;}long gauss(mat_GF2& M){ return gauss(M, M.NumCols());}void image(mat_GF2& X, const mat_GF2& A){ mat_GF2 M; M = A; long r = gauss(M); M.SetDims(r, M.NumCols()); X = M;}void kernel(mat_GF2& X, const mat_GF2& A){ long m = A.NumRows(); long n = A.NumCols(); mat_GF2 M; long r; transpose(M, A); r = gauss(M); X.SetDims(m-r, m); clear(X); long i, j, k; vec_long D; D.SetLength(m); for (j = 0; j < m; j++) D[j] = -1; j = -1; for (i = 0; i < r; i++) { do { j++; } while (M.get(i, j) == 0); D[j] = i; } for (k = 0; k < m-r; k++) { vec_GF2& v = X[k]; long pos = 0; for (j = m-1; j >= 0; j--) { if (D[j] == -1) { if (pos == k) { v[j] = 1; // v.put(j, to_GF2(1)); } pos++; } else { v[j] = v*M[D[j]]; // v.put(j, v*M[D[j]]); } } }} void mul(mat_GF2& X, const mat_GF2& A, GF2 b){ X = A; if (b == 0) clear(X);}void diag(mat_GF2& X, long n, GF2 d) { if (d == 1) ident(X, n); else { X.SetDims(n, n); clear(X); }} long IsDiag(const mat_GF2& A, long n, GF2 d){ if (A.NumRows() != n || A.NumCols() != n) return 0; if (d == 1) return IsIdent(A, n); else return IsZero(A);}long IsZero(const mat_GF2& a){ long n = a.NumRows(); long i; for (i = 0; i < n; i++) if (!IsZero(a[i])) return 0; return 1;}void clear(mat_GF2& x){ long n = x.NumRows(); long i; for (i = 0; i < n; i++) clear(x[i]);}mat_GF2 operator+(const mat_GF2& a, const mat_GF2& b){ mat_GF2 res; add(res, a, b); NTL_OPT_RETURN(mat_GF2, res);}mat_GF2 operator*(const mat_GF2& a, const mat_GF2& b){ mat_GF2 res; mul_aux(res, a, b); NTL_OPT_RETURN(mat_GF2, res);}mat_GF2 operator-(const mat_GF2& a, const mat_GF2& b){ mat_GF2 res; add(res, a, b); NTL_OPT_RETURN(mat_GF2, res);}vec_GF2 operator*(const mat_GF2& a, const vec_GF2& b){ vec_GF2 res; mul_aux(res, a, b); NTL_OPT_RETURN(vec_GF2, res);}vec_GF2 operator*(const vec_GF2& a, const mat_GF2& b){ vec_GF2 res; mul_aux(res, a, b); NTL_OPT_RETURN(vec_GF2, res);}void inv(mat_GF2& X, const mat_GF2& A){ GF2 d; inv(d, X, A); if (d == 0) Error("inv: non-invertible matrix");}void power(mat_GF2& X, const mat_GF2& A, const ZZ& e){ if (A.NumRows() != A.NumCols()) Error("power: non-square matrix"); if (e == 0) { ident(X, A.NumRows()); return; } mat_GF2 T1, T2; long i, k; k = NumBits(e); T1 = A; for (i = k-2; i >= 0; i--) { sqr(T2, T1); if (bit(e, i)) mul(T1, T2, A); else T1 = T2; } if (e < 0) inv(X, T1); else X = T1;}NTL_END_IMPL?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -