#ifndef MAT_H #define MAT_H #include #include #define LOOP(i,n) for(int i = 0; i < (n); i++) // library for matrix computations in variable rings (based on GMP types) #ifdef QEXT #include "qext.h" #define NUMBER qext_number #define INIT qext_init #define CLEAR qext_clear #define SET qext_set #define SET_INT qext_set_int #define SET_Q qext_set_q #define SET_ZERO(x) qext_set_int((x),0) #define SET_ONE(x) qext_set_int((x),1) #define ADD qext_add #define SUB qext_sub #define NEG qext_neg #define MUL qext_mul #define DIV qext_div #define INV qext_inv #define CMP qext_cmp #define PRINT qext_print #define SNPRINT qext_snprint #define TYPE struct qext_type* #define GETTYPE(x) ((x)->type) #else #define NUMBER mpq_t #define INIT(x,t) mpq_init(x) #define CLEAR mpq_clear #define SET mpq_set #define SET_INT(x,y) mpq_set_si(x,y,1) #define SET_Q SET #define SET_ZERO(x) SET_INT(x,0) #define SET_ONE(x) SET_INT(x,1) #define ADD mpq_add #define SUB mpq_sub #define NEG mpq_neg #define MUL mpq_mul #define DIV mpq_div #define INV mpq_inv #define PRINT(x) gmp_printf("%Qd", x) #define SNPRINT(out, size, x) gmp_snprintf(out, size, "%Qd", x) #define TYPE int #define GETTYPE(x) 0 #endif #define M(m,i,j) ((m)->x[(i)+(m)->n*(j)]) struct _mat{ int n; NUMBER *x; }; typedef struct _mat mat[1]; typedef struct _mat_workspace { mat tmp_mat; NUMBER tmp_num; NUMBER tmp_num2; } mat_workspace; mat_workspace *mat_workspace_init(int n, TYPE t); void mat_workspace_clear(mat_workspace *ws); void mat_init(mat m, int n, TYPE t); void mat_get(NUMBER out, mat m, int i, int j); void mat_set(mat m, int i, int j, NUMBER x); NUMBER *mat_ref(mat m, int i, int j); void mat_zero(mat m); void mat_identity(mat m); void mat_copy(mat to, mat from); void mat_clear(mat m); int mat_same(mat m1, mat m2); void mat_multiply(mat_workspace *ws, mat out, mat in1, mat in2); void mat_det(mat_workspace *ws, NUMBER out, mat in); void mat_pseudoinverse(mat_workspace *ws, mat out, mat in); void mat_trace(NUMBER out, mat in); void mat_print(mat in); #endif