/**************************************************************************\

MODULE: mat_zz_pE

SUMMARY:

Defines the class mat_zz_pE.

\**************************************************************************/


#include <NTL/matrix.h>
#include <NTL/vec_vec_lzz_pE.h>

NTL_matrix_decl(zz_pE,vec_zz_pE,vec_vec_zz_pE,mat_zz_pE)
NTL_io_matrix_decl(zz_pE,vec_zz_pE,vec_vec_zz_pE,mat_zz_pE)
NTL_eq_matrix_decl(zz_pE,vec_zz_pE,vec_vec_zz_pE,mat_zz_pE)

void add(mat_zz_pE& X, const mat_zz_pE& A, const mat_zz_pE& B); 
// X = A + B

void sub(mat_zz_pE& X, const mat_zz_pE& A, const mat_zz_pE& B); 
// X = A - B

void negate(mat_zz_pE& X, const mat_zz_pE& A); 
// X = - A

void mul(mat_zz_pE& X, const mat_zz_pE& A, const mat_zz_pE& B); 
// X = A * B

void mul(vec_zz_pE& x, const mat_zz_pE& A, const vec_zz_pE& b); 
// x = A * b

void mul(vec_zz_pE& x, const vec_zz_pE& a, const mat_zz_pE& B); 
// x = a * B

void mul(mat_zz_pE& X, const mat_zz_pE& A, const zz_pE& b);
void mul(mat_zz_pE& X, const mat_zz_pE& A, const zz_p& b);
void mul(mat_zz_pE& X, const mat_zz_pE& A, long b);
// X = A * b

void mul(mat_zz_pE& X, const zz_pE& a, const mat_zz_pE& B);
void mul(mat_zz_pE& X, const zz_p& a, const mat_zz_pE& B);
void mul(mat_zz_pE& X, long a, const mat_zz_pE& B);
// X = a * B


void determinant(zz_pE& d, const mat_zz_pE& A);
zz_pE determinant(const mat_zz_pE& a); 
// d = determinant(A)


void transpose(mat_zz_pE& X, const mat_zz_pE& A);
mat_zz_pE transpose(const mat_zz_pE& A);
// X = transpose of A

void solve(zz_pE& d, vec_zz_pE& X,
           const mat_zz_pE& A, const vec_zz_pE& 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.

void inv(zz_pE& d, mat_zz_pE& X, const mat_zz_pE& A);
// A is an n x n matrix.  Computes d = determinant(A).  If d != 0,
// computes X = A^{-1}.

void sqr(mat_zz_pE& X, const mat_zz_pE& A);
mat_zz_pE sqr(const mat_zz_pE& A);
// X = A*A   

void inv(mat_zz_pE& X, const mat_zz_pE& A);
mat_zz_pE inv(const mat_zz_pE& A);
// X = A^{-1}; error is raised if A is  singular

void power(mat_zz_pE& X, const mat_zz_pE& A, const ZZ& e);
mat_zz_pE power(const mat_zz_pE& A, const ZZ& e);

void power(mat_zz_pE& X, const mat_zz_pE& A, long e);
mat_zz_pE power(const mat_zz_pE& A, long e);
// X = A^e; e may be negative (in which case A must be nonsingular).

void ident(mat_zz_pE& X, long n);
mat_zz_pE ident_mat_zz_pE(long n);
// X = n x n identity matrix

long IsIdent(const mat_zz_pE& A, long n);
// test if A is the n x n identity matrix

void diag(mat_zz_pE& X, long n, const zz_pE& d);
mat_zz_pE diag(long n, const zz_pE& d);
// X = n x n diagonal matrix with d on diagonal

long IsDiag(const mat_zz_pE& A, long n, const zz_pE& d);
// test if X is an  n x n diagonal matrix with d on diagonal




long gauss(mat_zz_pE& M);
long gauss(mat_zz_pE& 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_pE& X, const mat_zz_pE& A);
// The rows of X are computed as basis of A's row space.  X is is row
// echelon form

void kernel(mat_zz_pE& X, const mat_zz_pE& A);
// Computes a basis for the kernel of the map x -> x*A. where x is a
// row vector.



// miscellaneous:

void clear(mat_zz_pE& a);
// x = 0 (dimension unchanged)

long IsZero(const mat_zz_pE& a);
// test if a is the zero matrix (any dimension)


// operator notation:

mat_zz_pE operator+(const mat_zz_pE& a, const mat_zz_pE& b);
mat_zz_pE operator-(const mat_zz_pE& a, const mat_zz_pE& b);
mat_zz_pE operator*(const mat_zz_pE& a, const mat_zz_pE& b);

mat_zz_pE operator-(const mat_zz_pE& a);


// matrix/scalar multiplication:

mat_zz_pE operator*(const mat_zz_pE& a, const zz_pE& b);
mat_zz_pE operator*(const mat_zz_pE& a, const zz_p& b);
mat_zz_pE operator*(const mat_zz_pE& a, long b);

mat_zz_pE operator*(const zz_pE& a, const mat_zz_pE& b);
mat_zz_pE operator*(const zz_p& a, const mat_zz_pE& b);
mat_zz_pE operator*(long a, const mat_zz_pE& b);

// matrix/vector multiplication:

vec_zz_pE operator*(const mat_zz_pE& a, const vec_zz_pE& b);

vec_zz_pE operator*(const vec_zz_pE& a, const mat_zz_pE& b);


// assignment operator notation:

mat_zz_pE& operator+=(mat_zz_pE& x, const mat_zz_pE& a);
mat_zz_pE& operator-=(mat_zz_pE& x, const mat_zz_pE& a);
mat_zz_pE& operator*=(mat_zz_pE& x, const mat_zz_pE& a);

mat_zz_pE& operator*=(mat_zz_pE& x, const zz_pE& a);
mat_zz_pE& operator*=(mat_zz_pE& x, const zz_p& a);
mat_zz_pE& operator*=(mat_zz_pE& x, long a);

vec_zz_pE& operator*=(vec_zz_pE& x, const mat_zz_pE& a);