/**************************************************************************\ MODULE: mat_zz_p SUMMARY: Defines the class mat_zz_p. Note that the modulus p need not be a prime, except as indicated below. IMPLEMENTATION NOTES: Starting with NTL version 9.7.0 (and 9.7.1), many of the routines here have been optimized to take better advantage of specific hardware features available on 64-bit Intel CPU's. Currently, the mul, inv, determinant, solve, gauss, kernel, and image routines are fastest for p up to 23-bits long (assuming the CPU supports AVX instructions). After that, performance degrades in three stages: stage 1: up to 28-bits; stage 2: up to 31-bits; stage 3: 32-bits and up. For primes up to 23-bits, AVX floating point instructions are used. After that, ordinary integer arithmetic is used. In a future version, I may exploit AVX2 integer instructions to get better stage 2 performance. And in the more distant future, AVX512 instructions will be used, when they become available. On older Intel machines, or non-Intel machines that have "long long" support, one still gets optimizations corresponding to the three stages above. On 32-bit machines, one still gets three stages, just with smaller crossover points. \**************************************************************************/ #include #include "vec_vec_zz_p.h" typedef Mat mat_zz_p; // backward compatibility void add(mat_zz_p& X, const mat_zz_p& A, const mat_zz_p& B); // X = A + B void sub(mat_zz_p& X, const mat_zz_p& A, const mat_zz_p& B); // X = A - B void mul(mat_zz_p& X, const mat_zz_p& A, const mat_zz_p& B); // X = A * B void mul(vec_zz_p& x, const mat_zz_p& A, const vec_zz_p& b); // x = A * b void mul(vec_zz_p& x, const vec_zz_p& a, const mat_zz_p& B); // x = a * B void mul(mat_zz_p& X, const mat_zz_p& A, zz_p b); void mul(mat_zz_p& X, const mat_zz_p& A, long b); // X = A * b void mul(mat_zz_p& X, zz_p a, const mat_zz_p& B); void mul(mat_zz_p& X, long a, const mat_zz_p& B); // X = a * B void transpose(mat_zz_p& X, const mat_zz_p& A); mat_zz_p transpose(const mat_zz_p& A); // X = transpose of A void determinant(zz_p& d, const mat_zz_p& A); zz_p determinant(const mat_zz_p& a); // d = determinant(A) void solve(zz_p& d, vec_zz_p& x, const mat_zz_p& A, const vec_zz_p& b); // A is an n x n matrix, b is a length n vector. Computes d = determinant(A). // If d != 0, solves x*A = b (so x and b are treated as a row vectors). void solve(zz_p& d, const mat_zz_p& A, vec_zz_p& x, const vec_zz_p& b); // A is an n x n matrix, b is a length n vector. Computes d = determinant(A). // If d != 0, solves A*x = b (so x and b are treated as a column vectors). void inv(zz_p& d, mat_zz_p& X, const mat_zz_p& A); // A is an n x n matrix. Computes d = determinant(A). If d != 0, // computes X = A^{-1}. void inv(mat_zz_p& X, const mat_zz_p& A); mat_zz_p inv(const mat_zz_p& A); // X = A^{-1}; error is raised if A is singular void power(mat_zz_p& X, const mat_zz_p& A, const ZZ& e); mat_zz_p power(const mat_zz_p& A, const ZZ& e); void power(mat_zz_p& X, const mat_zz_p& A, long e); mat_zz_p power(const mat_zz_p& A, long e); // X = A^e; e may be negative (in which case A must be nonsingular). // NOTE: the routines determinant, solve, inv, and power (with negative // exponent) all require that the modulus p is prime: during elimination, if a // non-zero pivot element does not have an inverse, and error is raised. The // following "relaxed" versions of these routines will also work with prime // powers, if the optional parameter relax is true (which is the default). // However, note that in these relaxed routines, if a computed determinant // value is zero, this may not be the true determinant: all that you can assume // is that the true determinant is not invertible mod p. If the parameter // relax==false, then these routines behave identically to their "unrelaxed" // counterparts. void relaxed_determinant(zz_p& d, const mat_zz_p& A, bool relax=true); zz_p relaxed_determinant(const mat_zz_p& a, bool relax=true); void relaxed_solve(zz_p& d, vec_zz_p& x, const mat_zz_p& A, const vec_zz_p& b, bool relax=true); void relaxed_solve(zz_p& d, const mat_zz_p& A, vec_zz_p& x, const vec_zz_p& b, bool relax=true); void relaxed_inv(zz_p& d, mat_zz_p& X, const mat_zz_p& A, bool relax=true); void relaxed_inv(mat_zz_p& X, const mat_zz_p& A, bool relax=true); mat_zz_p relaxed_inv(const mat_zz_p& A, bool relax=true); void relaxed_power(mat_zz_p& X, const mat_zz_p& A, const ZZ& e, bool relax=true); mat_zz_p relaxed_power(const mat_zz_p& A, const ZZ& e, bool relax=true); void relaxed_power(mat_zz_p& X, const mat_zz_p& A, long e, bool relax=true); mat_zz_p relaxed_power(const mat_zz_p& A, long e, bool relax=true); void sqr(mat_zz_p& X, const mat_zz_p& A); mat_zz_p sqr(const mat_zz_p& A); // X = A*A void ident(mat_zz_p& X, long n); mat_zz_p ident_mat_zz_p(long n); // X = n x n identity matrix long IsIdent(const mat_zz_p& A, long n); // test if A is the n x n identity matrix void diag(mat_zz_p& X, long n, zz_p d); mat_zz_p diag(long n, zz_p d); // X = n x n diagonal matrix with d on diagonal long IsDiag(const mat_zz_p& A, long n, zz_p d); // test if X is an n x n diagonal matrix with d on diagonal void random(mat_zz_p& x, long n, long m); // x = random n x m matrix mat_zz_p random_mat_zz_p(long n, long m); long gauss(mat_zz_p& M); long gauss(mat_zz_p& M, long w); // Performs unitary row operations so as to bring M into row echelon // form. If the optional argument w is supplied, stops when first w // columns are in echelon form. The return value is the rank (or the // rank of the first w columns). void image(mat_zz_p& X, const mat_zz_p& A); // The rows of X are computed as basis of A's row space. X is is row // echelon form void kernel(mat_zz_p& X, const mat_zz_p& A); // Computes a basis for the kernel of the map x -> x*A. where x is a // row vector. // NOTE: the gauss, image, and kernel routines all require that // the modulus p is prime. // miscellaneous: void clear(mat_zz_p& a); // x = 0 (dimension unchanged) long IsZero(const mat_zz_p& a); // test if a is the zero matrix (any dimension) // operator notation: mat_zz_p operator+(const mat_zz_p& a, const mat_zz_p& b); mat_zz_p operator-(const mat_zz_p& a, const mat_zz_p& b); mat_zz_p operator*(const mat_zz_p& a, const mat_zz_p& b); mat_zz_p operator-(const mat_zz_p& a); // matrix/scalar multiplication: mat_zz_p operator*(const mat_zz_p& a, zz_p b); mat_zz_p operator*(const mat_zz_p& a, long b); mat_zz_p operator*(zz_p a, const mat_zz_p& b); mat_zz_p operator*(long a, const mat_zz_p& b); // matrix/vector multiplication: vec_zz_p operator*(const mat_zz_p& a, const vec_zz_p& b); vec_zz_p operator*(const vec_zz_p& a, const mat_zz_p& b); // assignment operator notation: mat_zz_p& operator+=(mat_zz_p& x, const mat_zz_p& a); mat_zz_p& operator-=(mat_zz_p& x, const mat_zz_p& a); mat_zz_p& operator*=(mat_zz_p& x, const mat_zz_p& a); mat_zz_p& operator*=(mat_zz_p& x, zz_p a); mat_zz_p& operator*=(mat_zz_p& x, long a); vec_zz_p& operator*=(vec_zz_p& x, const mat_zz_p& a);