=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v retrieving revision 1.13 retrieving revision 1.63 diff -u -p -r1.13 -r1.63 --- OpenXM_contrib2/asir2000/builtin/array.c 2001/06/07 05:14:48 1.13 +++ OpenXM_contrib2/asir2000/builtin/array.c 2013/09/09 07:29:25 1.63 @@ -45,13 +45,21 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.12 2001/06/07 04:54:38 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.62 2013/06/14 05:55:24 ohara Exp $ */ #include "ca.h" #include "base.h" #include "parse.h" #include "inline.h" +#include +#include +#if !defined(_MSC_VER) +#include +#endif + +#define F4_INTRAT_PERIOD 8 + #if 0 #undef DMAR #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d); @@ -59,26 +67,17 @@ extern int DP_Print; /* XXX */ -void inner_product_mat_int_mod(Q **,int **,int,int,int,Q *); -void solve_by_lu_mod(int **,int,int,int **,int); -void solve_by_lu_gfmmat(GFMMAT,unsigned int,unsigned int *,unsigned int *); -int lu_gfmmat(GFMMAT,unsigned int,int *); -void mat_to_gfmmat(MAT,unsigned int,GFMMAT *); -int generic_gauss_elim_mod(int **,int,int,int,int *); -int generic_gauss_elim(MAT ,MAT *,Q *,int **,int **); - -int gauss_elim_mod(int **,int,int,int); -int gauss_elim_mod1(int **,int,int,int); -int gauss_elim_geninv_mod(unsigned int **,int,int,int); -int gauss_elim_geninv_mod_swap(unsigned int **,int,int,unsigned int,unsigned int ***,int **); void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(); -void Pnewbytearray(); +void Pinvmat(); +void Pnewbytearray(),Pmemoryplot_to_coord(); +void Pgeneric_gauss_elim(); void Pgeneric_gauss_elim_mod(); void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat(); -void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(); +void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(), Pltov(); +void Pgeninv_sf_swap(); void sepvect(); void Pmulmat_gf2n(); void Pbconvmat_gf2n(); @@ -90,25 +89,50 @@ void Px962_irredpoly_up2(); void Pirredpoly_up2(); void Pnbpoly_up2(); void Pqsort(); +void Pexponent_vector(); +void Pmat_swap_row_destructive(); +void Pmat_swap_col_destructive(); +void Pvect(); +void Pmat(); +void Pmatc(); +void Pnd_det(); +void Plu_mat(); +void Pmat_col(); +void Plusolve_prep(); +void Plusolve_main(); struct ftab array_tab[] = { + {"lu_mat",Plu_mat,1}, {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4}, {"lu_gfmmat",Plu_gfmmat,2}, {"mat_to_gfmmat",Pmat_to_gfmmat,2}, + {"generic_gauss_elim",Pgeneric_gauss_elim,1}, {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2}, {"newvect",Pnewvect,-2}, + {"vect",Pvect,-99999999}, + {"vector",Pnewvect,-2}, + {"exponent_vector",Pexponent_vector,-99999999}, {"newmat",Pnewmat,-3}, + {"matrix",Pnewmat,-3}, + {"mat",Pmat,-99999999}, + {"matr",Pmat,-99999999}, + {"matc",Pmatc,-99999999}, {"newbytearray",Pnewbytearray,-2}, + {"memoryplot_to_coord",Pmemoryplot_to_coord,1}, {"sepmat_destructive",Psepmat_destructive,2}, {"sepvect",Psepvect,2}, {"qsort",Pqsort,-2}, {"vtol",Pvtol,1}, + {"ltov",Pltov,1}, {"size",Psize,1}, {"det",Pdet,-2}, + {"nd_det",Pnd_det,-2}, + {"invmat",Pinvmat,-2}, {"leqm",Pleqm,2}, {"leqm1",Pleqm1,2}, {"geninvm",Pgeninvm,2}, {"geninvm_swap",Pgeninvm_swap,2}, + {"geninv_sf_swap",Pgeninv_sf_swap,1}, {"remainder",Premainder,2}, {"sremainder",Psremainder,2}, {"mulmat_gf2n",Pmulmat_gf2n,1}, @@ -119,26 +143,224 @@ struct ftab array_tab[] = { {"x962_irredpoly_up2",Px962_irredpoly_up2,2}, {"irredpoly_up2",Pirredpoly_up2,2}, {"nbpoly_up2",Pnbpoly_up2,2}, + {"mat_swap_row_destructive",Pmat_swap_row_destructive,3}, + {"mat_swap_col_destructive",Pmat_swap_col_destructive,3}, + {"mat_col",Pmat_col,2}, + {"lusolve_prep",Plusolve_prep,1}, + {"lusolve_main",Plusolve_main,1}, {0,0,0}, }; -int comp_obj(a,b) -Obj *a,*b; +typedef struct _ent { int j; unsigned int e; } ent; + +ent *get_row(FILE *,int *l); +void put_row(FILE *out,int l,ent *a); +int lu_elim(int *l,ent **a,int k,int i,int mul,int mod); + +static int *ul,*ll; +static ent **u,**l; +static int modulus; + +void Plusolve_prep(NODE arg,Q *rp) { + char *fname; + FILE *in; + int len,i,rank; + int *rhs; + + fname = BDY((STRING)ARG0(arg)); + in = fopen(fname,"r"); + modulus = getw(in); + len = getw(in); + ul = (int *)MALLOC_ATOMIC(len*sizeof(int)); + u = (ent **)MALLOC(len*sizeof(ent *)); + ll = (int *)MALLOC_ATOMIC(len*sizeof(int)); + l = (ent **)MALLOC(len*sizeof(ent *)); + for ( i = 0; i < len; i++ ) { + u[i] = get_row(in,&ul[i]); + } + for ( i = 0; i < len; i++ ) { + l[i] = get_row(in,&ll[i]); + } + fclose(in); + *rp = ONE; +} + +void Plusolve_main(NODE arg,VECT *rp) +{ + Q *d,*p; + VECT v,r; + int len,i; + int *rhs; + + v = (VECT)ARG0(arg); len = v->len; + d = (Q *)BDY(v); + rhs = (int *)MALLOC_ATOMIC(len*sizeof(int)); + for ( i = 0; i < len; i++ ) rhs[i] = QTOS(d[i]); + solve_l(ll,l,len,rhs,modulus); + solve_u(ul,u,len,rhs,modulus); + NEWVECT(r); r->len = len; + r->body = (pointer *)MALLOC(len*sizeof(pointer)); + p = (Q *)r->body; + for ( i = 0; i < len; i++ ) + STOQ(rhs[i],p[i]); + *rp = r; +} + +ent *get_row(FILE *in,int *l) +{ + int len,i; + ent *a; + + *l = len = getw(in); + a = (ent *)MALLOC_ATOMIC(len*sizeof(ent)); + for ( i = 0; i < len; i++ ) { + a[i].j = getw(in); + a[i].e = getw(in); + } + return a; +} + +int lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int mod) +{ + int i,j,k,s,mul; + unsigned int inv; + int *ll2; + + ll2 = (int *)MALLOC_ATOMIC(n*sizeof(int)); + for ( i = 0; i < n; i++ ) ll2[i] = 0; + for ( i = 0; i < n; i++ ) { + fprintf(stderr,"i=%d\n",i); + inv = invm(u[i][0].e,mod); + for ( k = i+1; k < n; k++ ) + if ( u[k][0].j == n-i ) { + s = u[k][0].e; + DMAR(s,inv,0,mod,mul); + lu_elim(ul,u,k,i,mul,mod); + lu_append(ll,l,ll2,k,i,mul); + } + } +} + +#define INITLEN 10 + +lu_append(int *l,ent **a,int *l2,int k,int i,int mul) +{ + int len; + ent *p; + + len = l[k]; + if ( !len ) { + a[k] = p = (ent *)MALLOC_ATOMIC(INITLEN*sizeof(ent)); + p[0].j = i; p[0].e = mul; + l[k] = 1; l2[k] = INITLEN; + } else { + if ( l2[k] == l[k] ) { + l2[k] *= 2; + a[k] = REALLOC(a[k],l2[k]*sizeof(ent)); + } + p =a[k]; + p[l[k]].j = i; p[l[k]].e = mul; + l[k]++; + } +} + +/* a[k] = a[k]-mul*a[i] */ + +int lu_elim(int *l,ent **a,int k,int i,int mul,int mod) +{ + ent *ak,*ai,*w; + int lk,li,j,m,p,q,r,s,t,j0; + + ak = a[k]; ai = a[i]; lk = l[k]; li = l[i]; + w = (ent *)alloca((lk+li)*sizeof(ent)); + p = 0; q = 0; j = 0; + mul = mod-mul; + while ( p < lk && q < li ) { + if ( ak[p].j > ai[q].j ) { + w[j] = ak[p]; j++; p++; + } else if ( ak[p].j < ai[q].j ) { + w[j].j = ai[q].j; + t = ai[q].e; + DMAR(t,mul,0,mod,r); + w[j].e = r; + j++; q++; + } else { + t = ai[q].e; s = ak[p].e; + DMAR(t,mul,s,mod,r); + if ( r ) { + w[j].j = ai[q].j; w[j].e = r; j++; + } + p++; q++; + } + } + if ( q == li ) + while ( p < lk ) { + w[j] = ak[p]; j++; p++; + } + else if ( p == lk ) + while ( q < li ) { + w[j].j = ai[q].j; + t = ai[q].e; + DMAR(t,mul,0,mod,r); + w[j].e = r; + j++; q++; + } + if ( j <= lk ) { + for ( m = 0; m < j; m++ ) ak[m] = w[m]; + } else { + a[k] = ak = (ent *)MALLOC_ATOMIC(j*sizeof(ent)); + for ( m = 0; m < j; m++ ) ak[m] = w[m]; + } + l[k] = j; +} + +int solve_l(int *ll,ent **l,int n,int *rhs,int mod) +{ + int j,k,s,len; + ent *p; + + for ( j = 0; j < n; j++ ) { + len = ll[j]; p = l[j]; + for ( k = 0, s = 0; k < len; k++ ) + s = dmar(p[k].e,rhs[p[k].j],s,mod); + rhs[j] -= s; + if ( rhs[j] < 0 ) rhs[j] += mod; + } +} + +int solve_u(int *ul,ent **u,int n,int *rhs,int mod) +{ + int j,k,s,len,inv; + ent *p; + + for ( j = n-1; j >= 0; j-- ) { + len = ul[j]; p = u[j]; + for ( k = 1, s = 0; k < len; k++ ) + s = dmar(p[k].e,rhs[p[k].j],s,mod); + rhs[j] -= s; + if ( rhs[j] < 0 ) rhs[j] += mod; + inv = invm((unsigned int)p[0].e,mod); + rhs[j] = dmar(rhs[j],inv,0,mod); + } +} + +int comp_obj(Obj *a,Obj *b) +{ return arf_comp(CO,*a,*b); } static FUNC generic_comp_obj_func; static NODE generic_comp_obj_arg; +static NODE generic_comp_obj_option; -int generic_comp_obj(a,b) -Obj *a,*b; +int generic_comp_obj(Obj *a,Obj *b) { Q r; BDY(generic_comp_obj_arg)=(pointer)(*a); BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b); - r = (Q)bevalf(generic_comp_obj_func,generic_comp_obj_arg); + r = (Q)bevalf_with_opts(generic_comp_obj_func,generic_comp_obj_arg,generic_comp_obj_option); if ( !r ) return 0; else @@ -146,19 +368,31 @@ Obj *a,*b; } -void Pqsort(arg,rp) -NODE arg; -VECT *rp; +void Pqsort(NODE arg,LIST *rp) { VECT vect; - char buf[BUFSIZ]; - char *fname; - NODE n; + NODE n,n1; P p; V v; + FUNC func; + int len,i; + pointer *a; + Obj t; - asir_assert(ARG0(arg),O_VECT,"qsort"); - vect = (VECT)ARG0(arg); + t = ARG0(arg); + if (OID(t) == O_LIST) { + n = (NODE)BDY((LIST)t); + len = length(n); + MKVECT(vect,len); + for ( i = 0; i < len; i++, n = NEXT(n) ) { + BDY(vect)[i] = BDY(n); + } + + }else if (OID(t) != O_VECT) { + error("qsort : invalid argument"); + }else { + vect = (VECT)t; + } if ( argc(arg) == 1 ) qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj); else { @@ -166,18 +400,29 @@ VECT *rp; if ( !p || OID(p)!=2 ) error("qsort : invalid argument"); v = VR(p); - if ( (int)v->attr != V_SR ) - error("qsort : no such function"); - generic_comp_obj_func = (FUNC)v->priv; - MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n); + gen_searchf(NAME(v),&func); + if ( !func ) { + if ( (int)v->attr != V_SR ) + error("qsort : no such function"); + func = (FUNC)v->priv; + } + generic_comp_obj_func = func; + MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n); + generic_comp_obj_option = current_option; qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj); } - *rp = vect; + if (OID(t) == O_LIST) { + a = BDY(vect); + for ( i = len - 1, n = 0; i >= 0; i-- ) { + MKNODE(n1,a[i],n); n = n1; + } + MKLIST(*rp,n); + }else { + *rp = (LIST)vect; + } } -void PNBmul_gf2n(arg,rp) -NODE arg; -GF2N *rp; +void PNBmul_gf2n(NODE arg,GF2N *rp) { GF2N a,b; GF2MAT mat; @@ -214,9 +459,7 @@ GF2N *rp; } } -void Pmul_vect_mat_gf2n(arg,rp) -NODE arg; -GF2N *rp; +void Pmul_vect_mat_gf2n(NODE arg,GF2N *rp) { GF2N a; GF2MAT mat; @@ -247,9 +490,7 @@ GF2N *rp; } } -void Pbconvmat_gf2n(arg,rp) -NODE arg; -LIST *rp; +void Pbconvmat_gf2n(NODE arg,LIST *rp) { P p0,p1; int to; @@ -269,9 +510,7 @@ LIST *rp; MKLIST(*rp,n0); } -void Pmulmat_gf2n(arg,rp) -NODE arg; -GF2MAT *rp; +void Pmulmat_gf2n(NODE arg,GF2MAT *rp) { GF2MAT m; @@ -280,9 +519,7 @@ GF2MAT *rp; *rp = m; } -void Psepmat_destructive(arg,rp) -NODE arg; -LIST *rp; +void Psepmat_destructive(NODE arg,LIST *rp) { MAT mat,mat1; int i,j,row,col; @@ -305,25 +542,20 @@ LIST *rp; sgn = SGN(ent); divn(nm,mod,&quo,&rem); /* if ( quo != nm && rem != nm ) */ -/* GC_free(nm); */ -/* GC_free(ent); */ +/* GCFREE(nm); */ +/* GCFREE(ent); */ NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]); } MKNODE(n1,mat1,0); MKNODE(n0,mat,n1); MKLIST(*rp,n0); } -void Psepvect(arg,rp) -NODE arg; -VECT *rp; +void Psepvect(NODE arg,VECT *rp) { sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp); } -void sepvect(v,d,rp) -VECT v; -int d; -VECT *rp; +void sepvect(VECT v,int d,VECT *rp) { int i,j,k,n,q,q1,r; pointer *pv,*pw,*pu; @@ -347,9 +579,7 @@ VECT *rp; } } -void Pnewvect(arg,rp) -NODE arg; -VECT *rp; +void Pnewvect(NODE arg,VECT *rp) { int len,i,r; VECT vect; @@ -365,34 +595,103 @@ VECT *rp; if ( argc(arg) == 2 ) { list = (LIST)ARG1(arg); asir_assert(list,O_LIST,"newvect"); +#if 0 for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ); if ( r > len ) { *rp = vect; return; } +#endif for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) ) vb[i] = (pointer)BDY(tn); } *rp = vect; } -void Pnewbytearray(arg,rp) -NODE arg; -BYTEARRAY *rp; +void Pvect(NODE arg,VECT *rp) { + int len,i; + VECT vect; + pointer *vb; + NODE tn; + + if ( !arg ) { + *rp =0; + return; + } + + for (len = 0, tn = arg; tn; tn = NEXT(tn), len++); + if ( len == 1 ) { + if ( ARG0(arg) != 0 ) { + switch ( OID(ARG0(arg)) ) { + case O_VECT: + *rp = ARG0(arg); + return; + case O_LIST: + for ( len = 0, tn = ARG0(arg); tn; tn = NEXT(tn), len++ ); + MKVECT(vect,len-1); + for ( i = 0, tn = BDY((LIST)ARG0(arg)), vb =BDY(vect); + tn; i++, tn = NEXT(tn) ) + vb[i] = (pointer)BDY(tn); + *rp=vect; + return; + } + } + } + MKVECT(vect,len); + for ( i = 0, tn = arg, vb = BDY(vect); tn; i++, tn = NEXT(tn) ) + vb[i] = (pointer)BDY(tn); + *rp = vect; +} + +void Pexponent_vector(NODE arg,DP *rp) { + nodetod(arg,rp); +} + +void Pnewbytearray(NODE arg,BYTEARRAY *rp) +{ int len,i,r; BYTEARRAY array; unsigned char *vb; char *str; LIST list; NODE tn; + int ac; + struct stat sbuf; + char *fname; + FILE *fp; - asir_assert(ARG0(arg),O_N,"newbytearray"); - len = QTOS((Q)ARG0(arg)); - if ( len < 0 ) - error("newbytearray : invalid size"); - MKBYTEARRAY(array,len); - if ( argc(arg) == 2 ) { + ac = argc(arg); + if ( ac == 1 ) { + if ( !OID((Obj)ARG0(arg)) ) error("newbytearray : invalid argument"); + switch ( OID((Obj)ARG0(arg)) ) { + case O_STR: + fname = BDY((STRING)ARG0(arg)); + fp = fopen(fname,"rb"); + if ( !fp ) error("newbytearray : fopen failed"); + if ( stat(fname,&sbuf) < 0 ) + error("newbytearray : stat failed"); + len = sbuf.st_size; + MKBYTEARRAY(array,len); + fread(BDY(array),len,sizeof(char),fp); + break; + case O_N: + if ( !RATN(ARG0(arg)) ) + error("newbytearray : invalid argument"); + len = QTOS((Q)ARG0(arg)); + if ( len < 0 ) + error("newbytearray : invalid size"); + MKBYTEARRAY(array,len); + break; + default: + error("newbytearray : invalid argument"); + } + } else if ( ac == 2 ) { + asir_assert(ARG0(arg),O_N,"newbytearray"); + len = QTOS((Q)ARG0(arg)); + if ( len < 0 ) + error("newbytearray : invalid size"); + MKBYTEARRAY(array,len); if ( !ARG1(arg) ) error("newbytearray : invalid initialization"); switch ( OID((Obj)ARG1(arg)) ) { @@ -416,14 +715,44 @@ BYTEARRAY *rp; if ( !ARG1(arg) ) error("newbytearray : invalid initialization"); } - } + } else + error("newbytearray : invalid argument"); *rp = array; } -void Pnewmat(arg,rp) -NODE arg; -MAT *rp; +#define MEMORY_GETPOINT(a,len,x,y) (((a)[(len)*(y)+((x)>>3)])&(1<<((x)&7))) + +void Pmemoryplot_to_coord(NODE arg,LIST *rp) { + int len,blen,y,i,j; + unsigned char *a; + NODE r0,r,n; + LIST l; + BYTEARRAY ba; + Q iq,jq; + + asir_assert(ARG0(arg),O_LIST,"memoryplot_to_coord"); + arg = BDY((LIST)ARG0(arg)); + len = QTOS((Q)ARG0(arg)); + blen = (len+7)/8; + y = QTOS((Q)ARG1(arg)); + ba = (BYTEARRAY)ARG2(arg); a = ba->body; + r0 = 0; + for ( j = 0; j < y; j++ ) + for ( i = 0; i < len; i++ ) + if ( MEMORY_GETPOINT(a,blen,i,j) ) { + NEXTNODE(r0,r); + STOQ(i,iq); STOQ(j,jq); + n = mknode(2,iq,jq); + MKLIST(l,n); + BDY(r) = l; + } + if ( r0 ) NEXT(r) = 0; + MKLIST(*rp,r0); +} + +void Pnewmat(NODE arg,MAT *rp) +{ int row,col; int i,j,r,c; NODE tn,sn; @@ -457,10 +786,110 @@ MAT *rp; *rp = m; } -void Pvtol(arg,rp) -NODE arg; -LIST *rp; +void Pmat(NODE arg, MAT *rp) { + int row,col; + int i; + MAT m; + pointer **mb; + pointer *ent; + NODE tn, sn; + VECT v; + + if ( !arg ) { + *rp =0; + return; + } + + for (row = 0, tn = arg; tn; tn = NEXT(tn), row++); + if ( row == 1 ) { + if ( OID(ARG0(arg)) == O_MAT ) { + *rp=ARG0(arg); + return; + } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) { + error("mat : invalid argument"); + } + } + if ( OID(ARG0(arg)) == O_VECT ) { + v = ARG0(arg); + col = v->len; + } else if ( OID(ARG0(arg)) == O_LIST ) { + for (col = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), col++); + } else { + error("mat : invalid argument"); + } + + MKMAT(m,row,col); + for (row = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), row++) { + if ( BDY(tn) == 0 ) { + error("mat : invalid argument"); + } else if ( OID(BDY(tn)) == O_VECT ) { + v = tn->body; + ent = BDY(v); + for (i = 0; i < v->len; i++ ) mb[row][i] = (Obj)ent[i]; + } else if ( OID(BDY(tn)) == O_LIST ) { + for (col = 0, sn = BDY((LIST)BDY(tn)); sn; col++, sn = NEXT(sn) ) + mb[row][col] = (pointer)BDY(sn); + } else { + error("mat : invalid argument"); + } + } + *rp = m; +} + +void Pmatc(NODE arg, MAT *rp) +{ + int row,col; + int i; + MAT m; + pointer **mb; + pointer *ent; + NODE tn, sn; + VECT v; + + if ( !arg ) { + *rp =0; + return; + } + + for (col = 0, tn = arg; tn; tn = NEXT(tn), col++); + if ( col == 1 ) { + if ( OID(ARG0(arg)) == O_MAT ) { + *rp=ARG0(arg); + return; + } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) { + error("matc : invalid argument"); + } + } + if ( OID(ARG0(arg)) == O_VECT ) { + v = ARG0(arg); + row = v->len; + } else if ( OID(ARG0(arg)) == O_LIST ) { + for (row = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), row++); + } else { + error("matc : invalid argument"); + } + + MKMAT(m,row,col); + for (col = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), col++) { + if ( BDY(tn) == 0 ) { + error("matc : invalid argument"); + } else if ( OID(BDY(tn)) == O_VECT ) { + v = tn->body; + ent = BDY(v); + for (i = 0; i < v->len; i++ ) mb[i][col] = (Obj)ent[i]; + } else if ( OID(BDY(tn)) == O_LIST ) { + for (row = 0, sn = BDY((LIST)BDY(tn)); sn; row++, sn = NEXT(sn) ) + mb[row][col] = (pointer)BDY(sn); + } else { + error("matc : invalid argument"); + } + } + *rp = m; +} + +void Pvtol(NODE arg,LIST *rp) +{ NODE n,n1; VECT v; pointer *a; @@ -474,10 +903,23 @@ LIST *rp; MKLIST(*rp,n); } -void Premainder(arg,rp) -NODE arg; -Obj *rp; +void Pltov(NODE arg,VECT *rp) { + NODE n; + VECT v; + int len,i; + + asir_assert(ARG0(arg),O_LIST,"ltov"); + n = (NODE)BDY((LIST)ARG0(arg)); + len = length(n); + MKVECT(v,len); + for ( i = 0; i < len; i++, n = NEXT(n) ) + BDY(v)[i] = BDY(n); + *rp = v; +} + +void Premainder(NODE arg,Obj *rp) +{ Obj a; VECT v,w; MAT m,l; @@ -524,9 +966,7 @@ Obj *rp; } } -void Psremainder(arg,rp) -NODE arg; -Obj *rp; +void Psremainder(NODE arg,Obj *rp) { Obj a; VECT v,w; @@ -575,9 +1015,7 @@ Obj *rp; } } -void Psize(arg,rp) -NODE arg; -LIST *rp; +void Psize(NODE arg,LIST *rp) { int n,m; @@ -596,6 +1034,10 @@ LIST *rp; n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col; STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s); break; + case O_IMAT: + n = ((IMAT)ARG0(arg))->row; m = ((IMAT)ARG0(arg))->col; + STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s); + break; default: error("size : invalid argument"); break; } @@ -603,9 +1045,7 @@ LIST *rp; MKLIST(*rp,t); } -void Pdet(arg,rp) -NODE arg; -P *rp; +void Pdet(NODE arg,P *rp) { MAT m; int n,i,j,mod; @@ -629,36 +1069,130 @@ P *rp; } } +void Pinvmat(NODE arg,LIST *rp) +{ + MAT m,r; + int n,i,j,mod; + P dn; + P **mat,**imat,**w; + NODE nd; + + m = (MAT)ARG0(arg); + asir_assert(m,O_MAT,"invmat"); + if ( m->row != m->col ) + error("invmat : non-square matrix"); + else if ( argc(arg) == 1 ) { + n = m->row; + invmatp(CO,(P **)BDY(m),n,&imat,&dn); + NEWMAT(r); r->row = n; r->col = n; r->body = (pointer **)imat; + nd = mknode(2,r,dn); + MKLIST(*rp,nd); + } else { + n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m); + w = (P **)almat_pointer(n,n); + for ( i = 0; i < n; i++ ) + for ( j = 0; j < n; j++ ) + ptomp(mod,mat[i][j],&w[i][j]); +#if 0 + detmp(CO,mod,w,n,&d); + mptop(d,rp); +#else + error("not implemented yet"); +#endif + } +} + /* input : a row x col matrix A A[I] <-> A[I][0]*x_0+A[I][1]*x_1+... + output : [B,D,R,C] + B : a rank(A) x col-rank(A) matrix + D : the denominator + R : a vector of length rank(A) + C : a vector of length col-rank(A) + B[I] <-> D*x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+... +*/ + +void Pgeneric_gauss_elim(NODE arg,LIST *rp) +{ + NODE n0,opt,p; + MAT m,nm; + int *ri,*ci; + VECT rind,cind; + Q dn,q; + int i,row,col,t,rank; + int is_hensel = 0; + char *key; + Obj value; + + if ( current_option ) { + for ( opt = current_option; opt; opt = NEXT(opt) ) { + p = BDY((LIST)BDY(opt)); + key = BDY((STRING)BDY(p)); + value = (Obj)BDY(NEXT(p)); + if ( !strcmp(key,"hensel") && value ) { + is_hensel = value ? 1 : 0; + break; + } + } + } + asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim"); + m = (MAT)ARG0(arg); + row = m->row; col = m->col; + if ( is_hensel ) + rank = generic_gauss_elim_hensel(m,&nm,&dn,&ri,&ci); + else + rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci); + t = col-rank; + MKVECT(rind,rank); + MKVECT(cind,t); + for ( i = 0; i < rank; i++ ) { + STOQ(ri[i],q); + BDY(rind)[i] = (pointer)q; + } + for ( i = 0; i < t; i++ ) { + STOQ(ci[i],q); + BDY(cind)[i] = (pointer)q; + } + n0 = mknode(4,nm,dn,rind,cind); + MKLIST(*rp,n0); +} + +/* + input : a row x col matrix A + A[I] <-> A[I][0]*x_0+A[I][1]*x_1+... + output : [B,R,C] B : a rank(A) x col-rank(A) matrix R : a vector of length rank(A) C : a vector of length col-rank(A) + RN : a vector of length rank(A) indicating useful rows + B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+... */ -void Pgeneric_gauss_elim_mod(arg,rp) -NODE arg; -LIST *rp; +void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp) { NODE n0; MAT m,mat; - VECT rind,cind; + VECT rind,cind,rnum; Q **tmat; - int **wmat; - Q *rib,*cib; - int *colstat; + int **wmat,**row0; + Q *rib,*cib,*rnb; + int *colstat,*p; Q q; - int md,i,j,k,l,row,col,t,n,rank; + int md,i,j,k,l,row,col,t,rank; asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod"); asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod"); m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg)); row = m->row; col = m->col; tmat = (Q **)m->body; wmat = (int **)almat(row,col); + + row0 = (int **)ALLOCA(row*sizeof(int *)); + for ( i = 0; i < row; i++ ) row0[i] = wmat[i]; + colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); for ( i = 0; i < row; i++ ) for ( j = 0; j < col; j++ ) @@ -671,6 +1205,13 @@ LIST *rp; wmat[i][j] = 0; rank = generic_gauss_elim_mod(wmat,row,col,md,colstat); + MKVECT(rnum,rank); + rnb = (Q *)rnum->body; + for ( i = 0; i < rank; i++ ) + for ( j = 0, p = wmat[i]; j < row; j++ ) + if ( p == row0[j] ) + STOQ(j,rnb[i]); + MKMAT(mat,rank,col-rank); tmat = (Q **)mat->body; for ( i = 0; i < rank; i++ ) @@ -688,13 +1229,11 @@ LIST *rp; } else { STOQ(j,cib[l]); l++; } - n0 = mknode(3,mat,rind,cind); + n0 = mknode(4,mat,rind,cind,rnum); MKLIST(*rp,n0); } -void Pleqm(arg,rp) -NODE arg; -VECT *rp; +void Pleqm(NODE arg,VECT *rp) { MAT m; VECT vect; @@ -733,9 +1272,7 @@ VECT *rp; } } -int gauss_elim_mod(mat,row,col,md) -int **mat; -int row,col,md; +int gauss_elim_mod(int **mat,int row,int col,int md) { int i,j,k,inv,a,n; int *t,*pivot; @@ -773,12 +1310,9 @@ int row,col,md; } struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb; +struct oEGT eg_conv; -int generic_gauss_elim(mat,nm,dn,rindp,cindp) -MAT mat; -MAT *nm; -Q *dn; -int **rindp,**cindp; +int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp) { int **wmat; Q **bmat; @@ -895,7 +1429,7 @@ RESET: add_eg(&eg_chrem_split,&tmp0,&tmp1); get_eg(&tmp0); - if ( ind % 16 ) + if ( ind % F4_INTRAT_PERIOD ) ret = 0; else ret = intmtoratm(crmat,m1,*nm,dn); @@ -930,12 +1464,123 @@ RESET: } } -int generic_gauss_elim_hensel(mat,nmmat,dn,rindp,cindp) -MAT mat; -MAT *nmmat; -Q *dn; -int **rindp,**cindp; +/* XXX broken */ +int lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm) { + Q **a0,**b; + Q *aiq; + N **a; + N *ai; + Q q,q1,dn2,a1,q0,bik; + MAT m; + unsigned int md; + int n,ind,i,j,rank,t,inv,t1,ret,min,k; + int **w; + int *wi,*rinfo0,*rinfo; + N m1,m2,m3,u,s; + + a0 = (Q **)mat->body; + n = mat->row; + if ( n != mat->col ) + error("lu_dec_cr : non-square matrix"); + w = (int **)almat(n,n); + MKMAT(m,n,n); + a = (N **)m->body; + UTON(1,m1); + rinfo0 = 0; + ind = 0; + while ( 1 ) { + md = get_lprime(ind); + /* mat mod md */ + for ( i = 0; i < n; i++ ) + for ( j = 0, aiq = a0[i], wi = w[i]; j < n; j++ ) + if ( q = aiq[j] ) { + t = rem(NM(q),md); + if ( t && SGN(q) < 0 ) + t = (md - t) % md; + wi[j] = t; + } else + wi[j] = 0; + + if ( !lu_mod((unsigned int **)w,n,md,&rinfo) ) continue; + printf("."); fflush(stdout); + if ( !rinfo0 ) + *perm = rinfo0 = rinfo; + else { + for ( i = 0; i < n; i++ ) + if ( rinfo[i] != rinfo0[i] ) break; + if ( i < n ) continue; + } + if ( UNIN(m1) ) { + for ( i = 0; i < n; i++ ) + for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) { + UTON(wi[j],u); ai[j] = u; + } + UTON(md,m1); + } else { + inv = invm(rem(m1,md),md); + UTON(md,m2); muln(m1,m2,&m3); + for ( i = 0; i < n; i++ ) + for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) + if ( ai[i] ) { + /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */ + t = rem(ai[j],md); + if ( wi[j] >= t ) + t = wi[j]-t; + else + t = md-(t-wi[j]); + DMAR(t,inv,0,md,t1) + UTON(t1,u); + muln(m1,u,&s); + addn(ai[j],s,&u); ai[j] = u; + } else if ( wi[j] ) { + /* f3 = m1*(m1 mod m2)^(-1)*f2 */ + DMAR(wi[j],inv,0,md,t) + UTON(t,u); + muln(m1,u,&s); ai[j] = s; + } + m1 = m3; + } + if ( (++ind%8) == 0 ) { + ret = intmtoratm(m,m1,lu,dn); + if ( ret ) { + b = (Q **)lu->body; + mulq(*dn,*dn,&dn2); + for ( i = 0; i < n; i++ ) { + for ( j = 0; j < n; j++ ) { + q = 0; + min = MIN(i,j); + for ( k = 0; k <= min; k++ ) { + bik = k==i ? *dn : b[i][k]; + mulq(bik,b[k][j],&q0); + addq(q,q0,&q1); q = q1; + } + mulq(a0[rinfo0[i]][j],dn2,&q1); + if ( cmpq(q,q1) ) break; + } + if ( j < n ) break; + } + if ( i == n ) + return; + } + } + } +} + +int nmat(N **m,int n) +{ + int i,j; + + for ( i = 0; i < n; i++ ) { + for ( j = 0; j < n; j++ ) { + printn(m[i][j]); printf(" "); + } + printf("\n"); + } +} + +int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp) +{ MAT bmat,xmat; Q **a0,**a,**b,**x,**nm; Q *ai,*bi,*xi; @@ -950,7 +1595,13 @@ int **rindp,**cindp; int *cinfo,*rinfo; int *rind,*cind; int count; - struct oEGT eg_mul,eg_inv,tmp0,tmp1; + int ret; + struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1; + int period; + int *wx,*ptr; + int wxsize,nsize; + N wn; + Q wq; a0 = (Q **)mat->body; row = mat->row; col = mat->col; @@ -968,7 +1619,13 @@ int **rindp,**cindp; } else wi[j] = 0; - rank = find_lhs_and_lu_mod(w,row,col,md,&rinfo,&cinfo); + if ( DP_Print > 3 ) { + fprintf(asir_out,"LU decomposition.."); fflush(asir_out); + } + rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo); + if ( DP_Print > 3 ) { + fprintf(asir_out,"done.\n"); fflush(asir_out); + } a = (Q **)almat_pointer(rank,rank); /* lhs mat */ MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */ for ( j = li = ri = 0; j < col; j++ ) @@ -998,9 +1655,17 @@ int **rindp,**cindp; *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int)); init_eg(&eg_mul); init_eg(&eg_inv); - for ( q = ONE, count = 0; ; count++ ) { - fprintf(stderr,"."); + init_eg(&eg_check); init_eg(&eg_intrat); + period = F4_INTRAT_PERIOD; + nsize = period; + wxsize = rank*ri*nsize; + wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int)); + for ( i = 0; i < wxsize; i++ ) wx[i] = 0; + for ( q = ONE, count = 0; ; ) { + if ( DP_Print > 3 ) + fprintf(stderr,"o"); /* wc = -b mod md */ + get_eg(&tmp0); for ( i = 0; i < rank; i++ ) for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) if ( u = (Q)bi[j] ) { @@ -1010,19 +1675,213 @@ int **rindp,**cindp; wi[j] = t; } else wi[j] = 0; - /* wc = A^(-1)wc; wc is normalized */ - get_eg(&tmp0); - solve_by_lu_mod(w,rank,md,wc,ri); + /* wc = A^(-1)wc; wc is not normalized */ + solve_by_lu_mod(w,rank,md,wc,ri,0); + /* wx += q*wc */ + ptr = wx; + for ( i = 0; i < rank; i++ ) + for ( j = 0, wi = wc[i]; j < ri; j++ ) { + if ( wi[j] ) + muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr); + ptr += nsize; + } + count++; get_eg(&tmp1); add_eg(&eg_inv,&tmp0,&tmp1); - /* x = x-q*wc */ + get_eg(&tmp0); for ( i = 0; i < rank; i++ ) - for ( j = 0, xi = x[i], wi = wc[i]; j < ri; j++ ) { - STOQ(wi[j],u); mulq(q,u,&s); - subq(xi[j],s,&u); xi[j] = u; + for ( j = 0; j < ri; j++ ) { + inner_product_mat_int_mod(a,wc,rank,i,j,&u); + addq(b[i][j],u,&s); + if ( s ) { + t = divin(NM(s),md,&tn); + if ( t ) + error("generic_gauss_elim_hensel:incosistent"); + NTOQ(tn,SGN(s),b[i][j]); + } else + b[i][j] = 0; } + get_eg(&tmp1); + add_eg(&eg_mul,&tmp0,&tmp1); + /* q = q*md */ + mulq(q,mdq,&u); q = u; + if ( count == period ) { + get_eg(&tmp0); + ptr = wx; + for ( i = 0; i < rank; i++ ) + for ( j = 0, xi = x[i]; j < ri; + j++, ptr += nsize ) { + for ( k = nsize-1; k >= 0 && !ptr[k]; k-- ); + if ( k >= 0 ) { + wn = NALLOC(k+1); + PL(wn) = k+1; + for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l]; + NTOQ(wn,1,wq); + subq(xi[j],wq,&u); xi[j] = u; + } + } + ret = intmtoratm_q(xmat,NM(q),*nmmat,dn); + get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1); + if ( ret ) { + rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); + cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); + for ( j = k = l = 0; j < col; j++ ) + if ( cinfo[j] ) + rind[k++] = j; + else + cind[l++] = j; + get_eg(&tmp0); + ret = gensolve_check(mat,*nmmat,*dn,rind,cind); + get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1); + if ( ret ) { + if ( DP_Print > 3 ) { + fprintf(stderr,"\n"); + print_eg("INV",&eg_inv); + print_eg("MUL",&eg_mul); + print_eg("INTRAT",&eg_intrat); + print_eg("CHECK",&eg_check); + fflush(asir_out); + } + *rindp = rind; + *cindp = cind; + for ( j = k = 0; j < col; j++ ) + if ( !cinfo[j] ) + cind[k++] = j; + return rank; + } + } else { + period = period*3/2; + count = 0; + nsize += period; + wxsize += rank*ri*nsize; + wx = (int *)REALLOC(wx,wxsize*sizeof(int)); + for ( i = 0; i < wxsize; i++ ) wx[i] = 0; + } + } + } + } +} + +int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT *nmmat,Q *dn,int **rindp,int **cindp) +{ + MAT bmat,xmat; + Q **a0,**a,**b,**x,**nm; + Q *ai,*bi,*xi; + int row,col; + int **w; + int *wi; + int **wc; + Q mdq,q,s,u; + N tn; + int ind,md,i,j,k,l,li,ri,rank; + unsigned int t; + int *cinfo,*rinfo; + int *rind,*cind; + int count; + int ret; + struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1; + int period; + int *wx,*ptr; + int wxsize,nsize; + N wn; + Q wq; + NumberField nf; + DP m; + int col1; + + a0 = (Q **)mat->body; + row = mat->row; col = mat->col; + w = (int **)almat(row,col); + for ( ind = 0; ; ind++ ) { + md = get_lprime(ind); + STOQ(md,mdq); + for ( i = 0; i < row; i++ ) + for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ ) + if ( q = (Q)ai[j] ) { + t = rem(NM(q),md); + if ( t && SGN(q) < 0 ) + t = (md - t) % md; + wi[j] = t; + } else + wi[j] = 0; + + if ( DP_Print ) { + fprintf(asir_out,"LU decomposition.."); fflush(asir_out); + } + rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo); + if ( DP_Print ) { + fprintf(asir_out,"done.\n"); fflush(asir_out); + } + for ( i = 0; i < col-1; i++ ) { + if ( !cinfo[i] ) { + m = mb[i]; + for ( j = i+1; j < col-1; j++ ) + if ( dp_redble(mb[j],m) ) + cinfo[j] = -1; + } + } + a = (Q **)almat_pointer(rank,rank); /* lhs mat */ + MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */ + for ( j = li = ri = 0; j < col; j++ ) + if ( cinfo[j] > 0 ) { + /* the column is in lhs */ + for ( i = 0; i < rank; i++ ) { + w[i][li] = w[i][j]; + a[i][li] = a0[rinfo[i]][j]; + } + li++; + } else if ( !cinfo[j] ) { + /* the column is in rhs */ + for ( i = 0; i < rank; i++ ) + b[i][ri] = a0[rinfo[i]][j]; + ri++; + } + + /* solve Ax+B=0; A: rank x rank, B: rank x ri */ + MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body; + MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body; + /* use the right part of w as work area */ + wc = (int **)almat(rank,ri); + for ( i = 0; i < rank; i++ ) + wc[i] = w[i]+rank; + *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); + *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); + init_eg(&eg_mul); init_eg(&eg_inv); + init_eg(&eg_check); init_eg(&eg_intrat); + period = F4_INTRAT_PERIOD; + nsize = period; + wxsize = rank*ri*nsize; + wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int)); + for ( i = 0; i < wxsize; i++ ) wx[i] = 0; + for ( q = ONE, count = 0; ; ) { + if ( DP_Print ) + fprintf(stderr,"o"); + /* wc = -b mod md */ get_eg(&tmp0); for ( i = 0; i < rank; i++ ) + for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) + if ( u = (Q)bi[j] ) { + t = rem(NM(u),md); + if ( t && SGN(u) > 0 ) + t = (md - t) % md; + wi[j] = t; + } else + wi[j] = 0; + /* wc = A^(-1)wc; wc is not normalized */ + solve_by_lu_mod(w,rank,md,wc,ri,0); + /* wx += q*wc */ + ptr = wx; + for ( i = 0; i < rank; i++ ) + for ( j = 0, wi = wc[i]; j < ri; j++ ) { + if ( wi[j] ) + muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr); + ptr += nsize; + } + count++; + get_eg(&tmp1); + add_eg(&eg_inv,&tmp0,&tmp1); + get_eg(&tmp0); + for ( i = 0; i < rank; i++ ) for ( j = 0; j < ri; j++ ) { inner_product_mat_int_mod(a,wc,rank,i,j,&u); addq(b[i][j],u,&s); @@ -1038,18 +1897,50 @@ int **rindp,**cindp; add_eg(&eg_mul,&tmp0,&tmp1); /* q = q*md */ mulq(q,mdq,&u); q = u; - if ( !(count % 16) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) { - for ( j = k = l = 0; j < col; j++ ) - if ( cinfo[j] ) - rind[k++] = j; - else - cind[l++] = j; - if ( gensolve_check(mat,*nmmat,*dn,rind,cind) ) { - fprintf(stderr,"\n"); - print_eg("INV",&eg_inv); - print_eg("MUL",&eg_mul); - fflush(asir_out); - return rank; + if ( count == period ) { + get_eg(&tmp0); + ptr = wx; + for ( i = 0; i < rank; i++ ) + for ( j = 0, xi = x[i]; j < ri; + j++, ptr += nsize ) { + for ( k = nsize-1; k >= 0 && !ptr[k]; k-- ); + if ( k >= 0 ) { + wn = NALLOC(k+1); + PL(wn) = k+1; + for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l]; + NTOQ(wn,1,wq); + subq(xi[j],wq,&u); xi[j] = u; + } + } + ret = intmtoratm_q(xmat,NM(q),*nmmat,dn); + get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1); + if ( ret ) { + for ( j = k = l = 0; j < col; j++ ) + if ( cinfo[j] > 0 ) + rind[k++] = j; + else if ( !cinfo[j] ) + cind[l++] = j; + get_eg(&tmp0); + ret = gensolve_check(mat,*nmmat,*dn,rind,cind); + get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1); + if ( ret ) { + if ( DP_Print > 3 ) { + fprintf(stderr,"\n"); + print_eg("INV",&eg_inv); + print_eg("MUL",&eg_mul); + print_eg("INTRAT",&eg_intrat); + print_eg("CHECK",&eg_check); + fflush(asir_out); + } + return rank; + } + } else { + period = period*3/2; + count = 0; + nsize += period; + wxsize += rank*ri*nsize; + wx = (int *)REALLOC(wx,wxsize*sizeof(int)); + for ( i = 0; i < wxsize; i++ ) wx[i] = 0; } } } @@ -1058,13 +1949,10 @@ int **rindp,**cindp; int f4_nocheck; -int gensolve_check(mat,nm,dn,rind,cind) -MAT mat,nm; -Q dn; -int *rind,*cind; +int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind) { int row,col,rank,clen,i,j,k,l; - Q s,t,u; + Q s,t; Q *w; Q *mati,*nmk; @@ -1109,13 +1997,10 @@ int *rind,*cind; /* assuming 0 < c < m */ -int inttorat(c,m,b,sgnp,nmp,dnp) -N c,m,b; -int *sgnp; -N *nmp,*dnp; +int inttorat(N c,N m,N b,int *sgnp,N *nmp,N *dnp) { - Q qq,t,u1,v1,r1,nm; - N q,r,u2,v2,r2; + Q qq,t,u1,v1,r1; + N q,u2,v2,r2; u1 = 0; v1 = ONE; u2 = m; v2 = c; while ( cmpn(v2,b) >= 0 ) { @@ -1134,14 +2019,10 @@ N *nmp,*dnp; /* mat->body = N ** */ -int intmtoratm(mat,md,nm,dn) -MAT mat; -N md; -MAT nm; -Q *dn; +int intmtoratm(MAT mat,N md,MAT nm,Q *dn) { N t,s,b; - Q bound,dn0,dn1,nm1,q,tq; + Q dn0,dn1,nm1,q; int i,j,k,l,row,col; Q **rmat; N **tmat; @@ -1191,14 +2072,10 @@ Q *dn; /* mat->body = Q ** */ -int intmtoratm_q(mat,md,nm,dn) -MAT mat; -N md; -MAT nm; -Q *dn; +int intmtoratm_q(MAT mat,N md,MAT nm,Q *dn) { N t,s,b; - Q bound,dn0,dn1,nm1,q,tq; + Q dn0,dn1,nm1,q; int i,j,k,l,row,col; Q **rmat; Q **tmat; @@ -1250,10 +2127,7 @@ Q *dn; #define ONE_STEP1 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; -void reduce_reducers_mod(mat,row,col,md) -int **mat; -int row,col; -int md; +void reduce_reducers_mod(int **mat,int row,int col,int md) { int i,j,k,l,hc,zzz; int *t,*s,*tj,*ind; @@ -1292,7 +2166,7 @@ int md; ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 } - for ( ; k >= 0; k-- ) { + for ( ; k > 0; k-- ) { if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; } } @@ -1308,10 +2182,7 @@ int md; 2. reduce spolys by the reduced reducers */ -void pre_reduce_mod(mat,row,col,nred,md) -int **mat; -int row,col,nred; -int md; +void pre_reduce_mod(int **mat,int row,int col,int nred,int md) { int i,j,k,l,hc,inv; int *t,*s,*tk,*ind; @@ -1361,14 +2232,10 @@ int md; mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order */ -void reduce_sp_by_red_mod(sp,redmat,ind,nred,col,md) -int *sp,**redmat; -int *ind; -int nred,col; -int md; +void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int col,int md) { int i,j,k,hc,zzz; - int *t,*s,*tj; + int *s,*tj; /* reduce the spolys by redmat */ for ( i = nred-1; i >= 0; i-- ) { @@ -1379,25 +2246,104 @@ int md; hc = md-hc; s = redmat[i]+j; tj = sp+j; - for ( k = col-j; k >= 0; k-- ) { + for ( k = col-j; k > 0; k-- ) { if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; } } } } +/* + mat[i] : compressed reducers (i=0,...,nred-1) + mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order +*/ + +void red_by_compress(int m,unsigned int *p,unsigned int *r, + unsigned int *ri,unsigned int hc,int len) +{ + unsigned int up,lo; + unsigned int dmy; + unsigned int *pj; + + p[*ri] = 0; r++; ri++; + for ( len--; len; len--, r++, ri++ ) { + pj = p+ *ri; + DMA(*r,hc,*pj,up,lo); + if ( up ) { + DSAB(m,up,lo,dmy,*pj); + } else + *pj = lo; + } +} + +/* p -= hc*r */ + +void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len) +{ + unsigned int up,lo,dmy; + + *p++ = 0; r++; len--; + for ( ; len; len--, r++, p++ ) + if ( *r ) { + DMA(*r,hc,*p,up,lo); + if ( up ) { + DSAB(m,up,lo,dmy,*p); + } else + *p = lo; + } +} + +void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len) +{ + *p++ = 0; r++; len--; + for ( ; len; len--, r++, p++ ) + if ( *r ) + *p = _addsf(_mulsf(*r,hc),*p); +} + +extern unsigned int **psca; + +void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind, + int nred,int col,int md) +{ + int i,len; + CDP ri; + unsigned int hc; + unsigned int *usp; + + usp = (unsigned int *)sp; + /* reduce the spolys by redmat */ + for ( i = nred-1; i >= 0; i-- ) { + /* reduce sp by redmat[i] */ + usp[ind[i]] %= md; + if ( hc = usp[ind[i]] ) { + /* sp = sp-hc*redmat[i] */ + hc = md-hc; + ri = redmat[i]; + len = ri->len; + red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len); + } + } + for ( i = 0; i < col; i++ ) + if ( usp[i] >= (unsigned int)md ) + usp[i] %= md; +} + #define ONE_STEP2 if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++; -int generic_gauss_elim_mod(mat,row,col,md,colstat) -int **mat; -int row,col,md; -int *colstat; +int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat) { - int i,j,k,l,inv,a,rank,zzz; - int *t,*pivot,*pk,*tk; + int i,j,k,l,inv,a,rank; + unsigned int *t,*pivot,*pk; + unsigned int **mat; + mat = (unsigned int **)mat0; for ( rank = 0, j = 0; j < col; j++ ) { - for ( i = rank; i < row && !mat[i][j]; i++ ); + for ( i = rank; i < row; i++ ) + mat[i][j] %= md; + for ( i = rank; i < row; i++ ) + if ( mat[i][j] ) + break; if ( i == row ) { colstat[j] = 0; continue; @@ -1410,36 +2356,68 @@ int *colstat; inv = invm(pivot[j],md); for ( k = j, pk = pivot+k; k < col; k++, pk++ ) if ( *pk ) { + if ( *pk >= (unsigned int)md ) + *pk %= md; DMAR(*pk,inv,0,md,*pk) } for ( i = rank+1; i < row; i++ ) { t = mat[i]; - if ( a = t[j] ) { - a = md - a; pk = pivot+j; tk = t+j; - k = col-j; - for ( ; k >= 64; k -= 64 ) { - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - } - for ( ; k >= 0; k -- ) { - if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++; - } + if ( a = t[j] ) + red_by_vect(md,t+j,pivot+j,md-a,col-j); + } + rank++; + } + for ( j = col-1, l = rank-1; j >= 0; j-- ) + if ( colstat[j] ) { + pivot = mat[l]; + for ( i = 0; i < l; i++ ) { + t = mat[i]; + t[j] %= md; + if ( a = t[j] ) + red_by_vect(md,t+j,pivot+j,md-a,col-j); } + l--; } + for ( j = 0, l = 0; l < rank; j++ ) + if ( colstat[j] ) { + t = mat[l]; + for ( k = j; k < col; k++ ) + if ( t[k] >= (unsigned int)md ) + t[k] %= md; + l++; + } + return rank; +} + +int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat) +{ + int i,j,k,l,inv,a,rank; + unsigned int *t,*pivot,*pk; + unsigned int **mat; + + mat = (unsigned int **)mat0; + for ( rank = 0, j = 0; j < col; j++ ) { + for ( i = rank; i < row; i++ ) + if ( mat[i][j] ) + break; + if ( i == row ) { + colstat[j] = 0; + continue; + } else + colstat[j] = 1; + if ( i != rank ) { + t = mat[i]; mat[i] = mat[rank]; mat[rank] = t; + } + pivot = mat[rank]; + inv = _invsf(pivot[j]); + for ( k = j, pk = pivot+k; k < col; k++, pk++ ) + if ( *pk ) + *pk = _mulsf(*pk,inv); + for ( i = rank+1; i < row; i++ ) { + t = mat[i]; + if ( a = t[j] ) + red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j); + } rank++; } for ( j = col-1, l = rank-1; j >= 0; j-- ) @@ -1447,31 +2425,8 @@ int *colstat; pivot = mat[l]; for ( i = 0; i < l; i++ ) { t = mat[i]; - if ( a = t[j] ) { - a = md-a; pk = pivot+j; tk = t+j; - k = col-j; - for ( ; k >= 64; k -= 64 ) { - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2 - } - for ( ; k >= 0; k -- ) { - if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++; - } - } + if ( a = t[j] ) + red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j); } l--; } @@ -1480,13 +2435,10 @@ int *colstat; /* LU decomposition; a[i][i] = 1/U[i][i] */ -int lu_gfmmat(mat,md,perm) -GFMMAT mat; -unsigned int md; -int *perm; +int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm) { int row,col; - int i,j,k,l; + int i,j,k; unsigned int *t,*pivot; unsigned int **a; unsigned int inv,m; @@ -1538,12 +2490,10 @@ int *perm; cinfo[j]=1 <=> j-th column is contained in the LU decomp. */ -int find_lhs_and_lu_mod(a,row,col,md,rinfo,cinfo) -unsigned int **a; -unsigned int md; -int **rinfo,**cinfo; +int find_lhs_and_lu_mod(unsigned int **a,int row,int col, + unsigned int md,int **rinfo,int **cinfo) { - int i,j,k,l,d; + int i,j,k,d; int *rp,*cp; unsigned int *t,*pivot; unsigned int inv,m; @@ -1582,6 +2532,40 @@ int **rinfo,**cinfo; return d; } +int lu_mod(unsigned int **a,int n,unsigned int md,int **rinfo) +{ + int i,j,k; + int *rp; + unsigned int *t,*pivot; + unsigned int inv,m; + + *rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int)); + for ( i = 0; i < n; i++ ) rp[i] = i; + for ( k = 0; k < n; k++ ) { + for ( i = k; i < n && !a[i][k]; i++ ); + if ( i == n ) return 0; + if ( i != k ) { + j = rp[i]; rp[i] = rp[k]; rp[k] = j; + t = a[i]; a[i] = a[k]; a[k] = t; + } + pivot = a[k]; + inv = invm(pivot[k],md); + for ( i = k+1; i < n; i++ ) { + t = a[i]; + if ( m = t[k] ) { + DMAR(inv,m,0,md,t[k]) + for ( j = k+1, m = md - t[k]; j < n; j++ ) + if ( pivot[j] ) { + unsigned int tj; + DMAR(m,pivot[j],t[j],md,tj) + t[j] = tj; + } + } + } + } + return 1; +} + /* Input a : n x n matrix; a result of LU-decomposition @@ -1591,12 +2575,7 @@ int **rinfo,**cinfo; b = a^(-1)b */ -void solve_by_lu_mod(a,n,md,b,l) -int **a; -int n; -int md; -int **b; -int l; +void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize) { unsigned int *y,*c; int i,j,k; @@ -1629,14 +2608,16 @@ int l; DMAR(t,a[i][i],0,md,c[i]) } /* copy c to b[.][k] with normalization */ - for ( i = 0; i < n; i++ ) - b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]); + if ( normalize ) + for ( i = 0; i < n; i++ ) + b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]); + else + for ( i = 0; i < n; i++ ) + b[i][k] = c[i]; } } -void Pleqm1(arg,rp) -NODE arg; -VECT *rp; +void Pleqm1(NODE arg,VECT *rp) { MAT m; VECT vect; @@ -1675,9 +2656,7 @@ VECT *rp; } } -gauss_elim_mod1(mat,row,col,md) -int **mat; -int row,col,md; +int gauss_elim_mod1(int **mat,int row,int col,int md) { int i,j,k,inv,a,n; int *t,*pivot; @@ -1714,9 +2693,7 @@ int row,col,md; return -1; } -void Pgeninvm(arg,rp) -NODE arg; -LIST *rp; +void Pgeninvm(NODE arg,LIST *rp) { MAT m; pointer **mat; @@ -1750,17 +2727,15 @@ LIST *rp; MKMAT(mat1,col,row); MKMAT(mat2,row-col,row); for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ ) for ( j = 0; j < row; j++ ) - STOQ(wmat[i][j+col],tmat[i][j]); + UTOQ(wmat[i][j+col],tmat[i][j]); for ( tmat = (Q **)mat2->body; i < row; i++ ) for ( j = 0; j < row; j++ ) - STOQ(wmat[i][j+col],tmat[i-col][j]); + UTOQ(wmat[i][j+col],tmat[i-col][j]); MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1); } } -int gauss_elim_geninv_mod(mat,row,col,md) -unsigned int **mat; -int row,col,md; +int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md) { int i,j,k,inv,a,n,m; unsigned int *t,*pivot; @@ -1796,9 +2771,7 @@ int row,col,md; return 0; } -void Psolve_by_lu_gfmmat(arg,rp) -NODE arg; -VECT *rp; +void Psolve_by_lu_gfmmat(NODE arg,VECT *rp) { GFMMAT lu; Q *perm,*rhs,*v; @@ -1819,15 +2792,12 @@ VECT *rp; solve_by_lu_gfmmat(lu,md,b,sol); MKVECT(r,n); for ( i = 0, v = (Q *)r->body; i < n; i++ ) - STOQ(sol[i],v[i]); + UTOQ(sol[i],v[i]); *rp = r; } -void solve_by_lu_gfmmat(lu,md,b,x) -GFMMAT lu; -unsigned int md; -unsigned int *b; -unsigned int *x; +void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md, + unsigned int *b,unsigned int *x) { int n; unsigned int **a; @@ -1859,10 +2829,30 @@ unsigned int *x; } } -void Plu_gfmmat(arg,rp) -NODE arg; -LIST *rp; +void Plu_mat(NODE arg,LIST *rp) { + MAT m,lu; + Q dn; + Q *v; + int n,i; + int *iperm; + VECT perm; + NODE n0; + + asir_assert(ARG0(arg),O_MAT,"lu_mat"); + m = (MAT)ARG0(arg); + n = m->row; + MKMAT(lu,n,n); + lu_dec_cr(m,lu,&dn,&iperm); + MKVECT(perm,n); + for ( i = 0, v = (Q *)perm->body; i < n; i++ ) + STOQ(iperm[i],v[i]); + n0 = mknode(3,lu,dn,perm); + MKLIST(*rp,n0); +} + +void Plu_gfmmat(NODE arg,LIST *rp) +{ MAT m; GFMMAT mm; unsigned int md; @@ -1872,8 +2862,8 @@ LIST *rp; VECT perm; NODE n0; - asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat"); - asir_assert(ARG1(arg),O_N,"mat_to_gfmmat"); + asir_assert(ARG0(arg),O_MAT,"lu_gfmmat"); + asir_assert(ARG1(arg),O_N,"lu_gfmmat"); m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg)); mat_to_gfmmat(m,md,&mm); row = m->row; @@ -1891,9 +2881,7 @@ LIST *rp; MKLIST(*rp,n0); } -void Pmat_to_gfmmat(arg,rp) -NODE arg; -GFMMAT *rp; +void Pmat_to_gfmmat(NODE arg,GFMMAT *rp) { MAT m; unsigned int md; @@ -1904,10 +2892,7 @@ GFMMAT *rp; mat_to_gfmmat(m,md,rp); } -void mat_to_gfmmat(m,md,rp) -MAT m; -unsigned int md; -GFMMAT *rp; +void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp) { unsigned int **wmat; unsigned int t; @@ -2033,14 +3018,108 @@ int **indexp; return 0; } +void Pgeninv_sf_swap(NODE arg,LIST *rp) +{ + MAT m; + GFS **mat,**tmat; + Q *tvect; + GFS q; + int **wmat,**invmat; + int *index; + unsigned int t; + int i,j,row,col,status; + MAT mat1; + VECT vect1; + NODE node1,node2; + + asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap"); + m = (MAT)ARG0(arg); + row = m->row; col = m->col; mat = (GFS **)m->body; + wmat = (int **)almat(row,col+row); + for ( i = 0; i < row; i++ ) { + bzero((char *)wmat[i],(col+row)*sizeof(int)); + for ( j = 0; j < col; j++ ) + if ( q = (GFS)mat[i][j] ) + wmat[i][j] = FTOIF(CONT(q)); + wmat[i][col+i] = _onesf(); + } + status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index); + if ( status > 0 ) + *rp = 0; + else { + MKMAT(mat1,col,col); + for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ ) + for ( j = 0; j < col; j++ ) + if ( t = invmat[i][j] ) { + MKGFS(IFTOF(t),tmat[i][j]); + } + MKVECT(vect1,row); + for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ ) + STOQ(index[i],tvect[i]); + MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1); + } +} + +int gauss_elim_geninv_sf_swap(int **mat,int row,int col, + int ***invmatp,int **indexp) +{ + int i,j,k,inv,a,n,m,u; + int *t,*pivot,*s; + int *index; + int **invmat; + + n = col; m = row+col; + *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int)); + for ( i = 0; i < row; i++ ) + index[i] = i; + for ( j = 0; j < n; j++ ) { + for ( i = j; i < row && !mat[i][j]; i++ ); + if ( i == row ) { + *indexp = 0; *invmatp = 0; return 1; + } + if ( i != j ) { + t = mat[i]; mat[i] = mat[j]; mat[j] = t; + k = index[i]; index[i] = index[j]; index[j] = k; + } + pivot = mat[j]; + inv = _invsf(pivot[j]); + for ( k = j; k < m; k++ ) + if ( pivot[k] ) + pivot[k] = _mulsf(pivot[k],inv); + for ( i = j+1; i < row; i++ ) { + t = mat[i]; + if ( a = t[j] ) + for ( k = j, a = _chsgnsf(a); k < m; k++ ) + if ( pivot[k] ) { + u = _mulsf(pivot[k],a); + t[k] = _addsf(u,t[k]); + } + } + } + for ( j = n-1; j >= 0; j-- ) { + pivot = mat[j]; + for ( i = j-1; i >= 0; i-- ) { + t = mat[i]; + if ( a = t[j] ) + for ( k = j, a = _chsgnsf(a); k < m; k++ ) + if ( pivot[k] ) { + u = _mulsf(pivot[k],a); + t[k] = _addsf(u,t[k]); + } + } + } + *invmatp = invmat = (int **)almat(col,col); + for ( i = 0; i < col; i++ ) + for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ ) + s[j] = t[col+index[j]]; + return 0; +} + void _addn(N,N,N); int _subn(N,N,N); void _muln(N,N,N); -void inner_product_int(a,b,n,r) -Q *a,*b; -int n; -Q *r; +void inner_product_int(Q *a,Q *b,int n,Q *r) { int la,lb,i; int sgn,sgn1; @@ -2082,10 +3161,10 @@ Q *r; t = wma; wma = sum; sum = t; } } - GC_free(wm); - GC_free(wma); + GCFREE(wm); + GCFREE(wma); if ( !sgn ) { - GC_free(sum); + GCFREE(sum); *r = 0; } else NTOQ(sum,sgn,*r); @@ -2093,11 +3172,7 @@ Q *r; /* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */ -void inner_product_mat_int_mod(a,b,n,k,l,r) -Q **a; -int **b; -int n,k,l; -Q *r; +void inner_product_mat_int_mod(Q **a,int **b,int n,int k,int l,Q *r) { int la,lb,i; int sgn,sgn1; @@ -2144,18 +3219,16 @@ Q *r; t = wma; wma = sum; sum = t; } } - GC_free(wm); - GC_free(wma); + GCFREE(wm); + GCFREE(wma); if ( !sgn ) { - GC_free(sum); + GCFREE(sum); *r = 0; } else NTOQ(sum,sgn,*r); } -void Pmul_mat_vect_int(arg,rp) -NODE arg; -VECT *rp; +void Pmul_mat_vect_int(NODE arg,VECT *rp) { MAT mat; VECT vect,r; @@ -2166,14 +3239,13 @@ VECT *rp; row = mat->row; col = mat->col; MKVECT(r,row); - for ( i = 0; i < row; i++ ) - inner_product_int(mat->body[i],vect->body,col,&r->body[i]); + for ( i = 0; i < row; i++ ) { + inner_product_int((Q *)mat->body[i],(Q *)vect->body,col,(Q *)&r->body[i]); + } *rp = r; } -void Pnbpoly_up2(arg,rp) -NODE arg; -GF2N *rp; +void Pnbpoly_up2(NODE arg,GF2N *rp) { int m,type,ret; UP2 r; @@ -2187,11 +3259,9 @@ GF2N *rp; *rp = 0; } -void Px962_irredpoly_up2(arg,rp) -NODE arg; -GF2N *rp; +void Px962_irredpoly_up2(NODE arg,GF2N *rp) { - int m,type,ret,w; + int m,ret,w; GF2N prev; UP2 r; @@ -2207,18 +3277,16 @@ GF2N *rp; bzero((char *)r->b,w*sizeof(unsigned int)); } } - ret = _generate_irreducible_polynomial(r,m,type); + ret = _generate_irreducible_polynomial(r,m); if ( ret == 0 ) MKGF2N(r,*rp); else *rp = 0; } -void Pirredpoly_up2(arg,rp) -NODE arg; -GF2N *rp; +void Pirredpoly_up2(NODE arg,GF2N *rp) { - int m,type,ret,w; + int m,ret,w; GF2N prev; UP2 r; @@ -2234,13 +3302,55 @@ GF2N *rp; bzero((char *)r->b,w*sizeof(unsigned int)); } } - ret = _generate_good_irreducible_polynomial(r,m,type); + ret = _generate_good_irreducible_polynomial(r,m); if ( ret == 0 ) MKGF2N(r,*rp); else *rp = 0; } +void Pmat_swap_row_destructive(NODE arg, MAT *m) +{ + int i1,i2; + pointer *t; + MAT mat; + + asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive"); + asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive"); + asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive"); + mat = (MAT)ARG0(arg); + i1 = QTOS((Q)ARG1(arg)); + i2 = QTOS((Q)ARG2(arg)); + if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row ) + error("mat_swap_row_destructive : Out of range"); + t = mat->body[i1]; + mat->body[i1] = mat->body[i2]; + mat->body[i2] = t; + *m = mat; +} + +void Pmat_swap_col_destructive(NODE arg, MAT *m) +{ + int j1,j2,i,n; + pointer *mi; + pointer t; + MAT mat; + + asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive"); + asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive"); + asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive"); + mat = (MAT)ARG0(arg); + j1 = QTOS((Q)ARG1(arg)); + j2 = QTOS((Q)ARG2(arg)); + if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col ) + error("mat_swap_col_destructive : Out of range"); + n = mat->row; + for ( i = 0; i < n; i++ ) { + mi = mat->body[i]; + t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t; + } + *m = mat; +} /* * f = type 'type' normal polynomial of degree m if exists * IEEE P1363 A.7.2 @@ -2514,9 +3624,7 @@ PENTA: return 1; } -printqmat(mat,row,col) -Q **mat; -int row,col; +void printqmat(Q **mat,int row,int col) { int i,j; @@ -2528,9 +3636,7 @@ int row,col; } } -printimat(mat,row,col) -int **mat; -int row,col; +void printimat(int **mat,int row,int col) { int i,j; @@ -2540,4 +3646,33 @@ int row,col; } printf("\n"); } +} + +void Pnd_det(NODE arg,P *rp) +{ + if ( argc(arg) == 1 ) + nd_det(0,ARG0(arg),rp); + else + nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp); +} + +void Pmat_col(NODE arg,VECT *rp) +{ + int i,j,n; + MAT mat; + VECT vect; + + asir_assert(ARG0(arg),O_MAT,"mat_col"); + asir_assert(ARG1(arg),O_N,"mat_col"); + mat = (MAT)ARG0(arg); + j = QTOS((Q)ARG1(arg)); + if ( j < 0 || j >= mat->col) { + error("mat_col : Out of range"); + } + n = mat->row; + MKVECT(vect,n); + for(i=0; i