[2cfffe] | 1 | |
---|
| 2 | /**************************************************************************\ |
---|
| 3 | |
---|
| 4 | MODULE: mat_ZZ |
---|
| 5 | |
---|
| 6 | SUMMARY: |
---|
| 7 | |
---|
| 8 | Defines the class mat_ZZ. |
---|
| 9 | |
---|
| 10 | \**************************************************************************/ |
---|
| 11 | |
---|
| 12 | |
---|
| 13 | #include <NTL/matrix.h> |
---|
| 14 | #include <NTL/vec_vec_ZZ.h> |
---|
| 15 | |
---|
| 16 | NTL_matrix_decl(ZZ,vec_ZZ,vec_vec_ZZ,mat_ZZ) |
---|
| 17 | NTL_io_matrix_decl(ZZ,vec_ZZ,vec_vec_ZZ,mat_ZZ) |
---|
| 18 | NTL_eq_matrix_decl(ZZ,vec_ZZ,vec_vec_ZZ,mat_ZZ) |
---|
| 19 | |
---|
| 20 | void add(mat_ZZ& X, const mat_ZZ& A, const mat_ZZ& B); |
---|
| 21 | // X = A + B |
---|
| 22 | |
---|
| 23 | void sub(mat_ZZ& X, const mat_ZZ& A, const mat_ZZ& B); |
---|
| 24 | // X = A - B |
---|
| 25 | |
---|
| 26 | void negate(mat_ZZ& X, const mat_ZZ& A); |
---|
| 27 | // X = - A |
---|
| 28 | |
---|
| 29 | void mul(mat_ZZ& X, const mat_ZZ& A, const mat_ZZ& B); |
---|
| 30 | // X = A * B |
---|
| 31 | |
---|
| 32 | void mul(vec_ZZ& x, const mat_ZZ& A, const vec_ZZ& b); |
---|
| 33 | // x = A * b |
---|
| 34 | |
---|
| 35 | void mul(vec_ZZ& x, const vec_ZZ& a, const mat_ZZ& B); |
---|
| 36 | // x = a * B |
---|
| 37 | |
---|
| 38 | void mul(mat_ZZ& X, const mat_ZZ& A, const ZZ& b); |
---|
| 39 | void mul(mat_ZZ& X, const mat_ZZ& A, long b); |
---|
| 40 | // X = A * b |
---|
| 41 | |
---|
| 42 | void mul(mat_ZZ& X, const ZZ& a, const mat_ZZ& B); |
---|
| 43 | void mul(mat_ZZ& X, long a, const mat_ZZ& B); |
---|
| 44 | // X = a * B |
---|
| 45 | |
---|
| 46 | |
---|
| 47 | |
---|
| 48 | void determinant(ZZ& d, const mat_ZZ& A, long deterministic=0); |
---|
| 49 | ZZ determinant(const mat_ZZ& a, long deterministic=0); |
---|
| 50 | // d = determinant(A). If !deterministic, a randomized strategy may |
---|
| 51 | // be used that errs with probability at most 2^{-80}. |
---|
| 52 | |
---|
| 53 | |
---|
| 54 | |
---|
| 55 | void solve(ZZ& d, vec_ZZ& x, |
---|
| 56 | const mat_ZZ& A, const vec_ZZ& b, |
---|
| 57 | long deterministic=0) |
---|
| 58 | // computes d = determinant(A) and solves x*A = b*d if d != 0; A must |
---|
| 59 | // be a square matrix and have compatible dimensions with b. If |
---|
| 60 | // !deterministic, the computation of d may use a randomized strategy |
---|
| 61 | // that errs with probability 2^{-80}. |
---|
| 62 | |
---|
| 63 | |
---|
| 64 | |
---|
| 65 | void solve1(ZZ& d, vec_ZZ& x, const mat_ZZ& A, const vec_ZZ& b); |
---|
| 66 | // A must be a square matrix. |
---|
| 67 | // If A is singular, this routine sets d = 0 and returns. |
---|
| 68 | // Otherwise, it computes d, x such that x*A == b*d, |
---|
| 69 | // such that d > 0 and minimal. |
---|
| 70 | // Note that d is a positive divisor of the determinant, |
---|
| 71 | // and is not in general equal to the determinant. |
---|
| 72 | // The routine is deterministic, and uses either a Hensel lifting |
---|
| 73 | // strategy. |
---|
| 74 | |
---|
| 75 | // For backward compatability, there is also a routine called |
---|
| 76 | // HenselSolve1 that simply calls solve1. |
---|
| 77 | |
---|
| 78 | |
---|
| 79 | void inv(ZZ& d, mat_ZZ& X, const mat_ZZ& A, long deterministic=0); |
---|
| 80 | // computes d = determinant(A) and solves X*A = I*d if d != 0; A must |
---|
| 81 | // be a square matrix. If !deterministic, the computation of d may |
---|
| 82 | // use a randomized strategy that errs with probability 2^{-80}. |
---|
| 83 | |
---|
| 84 | |
---|
| 85 | // NOTE: See LLL.txt for routines that compute the kernel and |
---|
| 86 | // image of an integer matrix. |
---|
| 87 | |
---|
| 88 | // NOTE: See HNF.txt for a routine that computes Hermite Normal Forms. |
---|
| 89 | |
---|
| 90 | void sqr(mat_ZZ& X, const mat_ZZ& A); |
---|
| 91 | mat_ZZ sqr(const mat_ZZ& A); |
---|
| 92 | // X = A*A |
---|
| 93 | |
---|
| 94 | void inv(mat_ZZ& X, const mat_ZZ& A); |
---|
| 95 | mat_ZZ inv(const mat_ZZ& A); |
---|
| 96 | // X = A^{-1}; error is raised if |det(A)| != 1. |
---|
| 97 | |
---|
| 98 | void power(mat_ZZ& X, const mat_ZZ& A, const ZZ& e); |
---|
| 99 | mat_ZZ power(const mat_ZZ& A, const ZZ& e); |
---|
| 100 | |
---|
| 101 | void power(mat_ZZ& X, const mat_ZZ& A, long e); |
---|
| 102 | mat_ZZ power(const mat_ZZ& A, long e); |
---|
| 103 | // X = A^e; e may be negative (in which case A must be nonsingular). |
---|
| 104 | |
---|
| 105 | |
---|
| 106 | |
---|
| 107 | void ident(mat_ZZ& X, long n); |
---|
| 108 | mat_ZZ ident_mat_ZZ(long n); |
---|
| 109 | // X = n x n identity matrix |
---|
| 110 | |
---|
| 111 | long IsIdent(const mat_ZZ& A, long n); |
---|
| 112 | // test if A is the n x n identity matrix |
---|
| 113 | |
---|
| 114 | void diag(mat_ZZ& X, long n, const ZZ& d); |
---|
| 115 | mat_ZZ diag(long n, const ZZ& d); |
---|
| 116 | // X = n x n diagonal matrix with d on diagonal |
---|
| 117 | |
---|
| 118 | long IsDiag(const mat_ZZ& A, long n, const ZZ& d); |
---|
| 119 | // test if X is an n x n diagonal matrix with d on diagonal |
---|
| 120 | |
---|
| 121 | |
---|
| 122 | void transpose(mat_ZZ& X, const mat_ZZ& A); |
---|
| 123 | mat_ZZ transpose(const mat_ZZ& A); |
---|
| 124 | // X = transpose of A |
---|
| 125 | |
---|
| 126 | |
---|
| 127 | long CRT(mat_ZZ& a, ZZ& prod, const mat_zz_p& A); |
---|
| 128 | // Incremental Chinese Remaindering: If p is the current zz_p modulus with |
---|
| 129 | // (p, prod) = 1; Computes a' such that a' = a mod prod and a' = A mod p, |
---|
| 130 | // with coefficients in the interval (-p*prod/2, p*prod/2]; |
---|
| 131 | // Sets a := a', prod := p*prod, and returns 1 if a's value changed. |
---|
| 132 | |
---|
| 133 | |
---|
| 134 | |
---|
| 135 | // miscellaneous: |
---|
| 136 | |
---|
| 137 | void clear(mat_ZZ& a); |
---|
| 138 | // x = 0 (dimension unchanged) |
---|
| 139 | |
---|
| 140 | long IsZero(const mat_ZZ& a); |
---|
| 141 | // test if a is the zero matrix (any dimension) |
---|
| 142 | |
---|
| 143 | |
---|
| 144 | // operator notation: |
---|
| 145 | |
---|
| 146 | mat_ZZ operator+(const mat_ZZ& a, const mat_ZZ& b); |
---|
| 147 | mat_ZZ operator-(const mat_ZZ& a, const mat_ZZ& b); |
---|
| 148 | mat_ZZ operator*(const mat_ZZ& a, const mat_ZZ& b); |
---|
| 149 | |
---|
| 150 | mat_ZZ operator-(const mat_ZZ& a); |
---|
| 151 | |
---|
| 152 | |
---|
| 153 | // matrix/scalar multiplication: |
---|
| 154 | |
---|
| 155 | mat_ZZ operator*(const mat_ZZ& a, const ZZ& b); |
---|
| 156 | mat_ZZ operator*(const mat_ZZ& a, long b); |
---|
| 157 | |
---|
| 158 | mat_ZZ operator*(const ZZ& a, const mat_ZZ& b); |
---|
| 159 | mat_ZZ operator*(long a, const mat_ZZ& b); |
---|
| 160 | |
---|
| 161 | // matrix/vector multiplication: |
---|
| 162 | |
---|
| 163 | vec_ZZ operator*(const mat_ZZ& a, const vec_ZZ& b); |
---|
| 164 | |
---|
| 165 | vec_ZZ operator*(const vec_ZZ& a, const mat_ZZ& b); |
---|
| 166 | |
---|
| 167 | |
---|
| 168 | |
---|
| 169 | // assignment operator notation: |
---|
| 170 | |
---|
| 171 | mat_ZZ& operator+=(mat_ZZ& x, const mat_ZZ& a); |
---|
| 172 | mat_ZZ& operator-=(mat_ZZ& x, const mat_ZZ& a); |
---|
| 173 | mat_ZZ& operator*=(mat_ZZ& x, const mat_ZZ& a); |
---|
| 174 | |
---|
| 175 | mat_ZZ& operator*=(mat_ZZ& x, const ZZ& a); |
---|
| 176 | mat_ZZ& operator*=(mat_ZZ& x, long a); |
---|
| 177 | |
---|
| 178 | vec_ZZ& operator*=(vec_ZZ& x, const mat_ZZ& a); |
---|
| 179 | |
---|
| 180 | |
---|