=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v retrieving revision 1.24 retrieving revision 1.46 diff -u -p -r1.24 -r1.46 --- OpenXM_contrib2/asir2000/builtin/array.c 2001/10/09 01:36:05 1.24 +++ OpenXM_contrib2/asir2000/builtin/array.c 2005/02/08 18:06:05 1.46 @@ -45,13 +45,15 @@ * 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.23 2001/10/01 01:58:01 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.45 2005/01/23 14:03:47 noro Exp $ */ #include "ca.h" #include "base.h" #include "parse.h" #include "inline.h" +#define F4_INTRAT_PERIOD 8 + #if 0 #undef DMAR #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d); @@ -64,10 +66,12 @@ void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet( void Pinvmat(); void Pnewbytearray(); +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(); @@ -80,29 +84,43 @@ 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(); struct ftab array_tab[] = { {"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}, {"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}, @@ -113,6 +131,8 @@ 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}, {0,0,0}, }; @@ -138,15 +158,31 @@ int generic_comp_obj(Obj *a,Obj *b) } -void Pqsort(NODE arg,VECT *rp) +void Pqsort(NODE arg,LIST *rp) { VECT vect; - 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 { @@ -154,13 +190,25 @@ void Pqsort(NODE arg,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; + 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); 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(NODE arg,GF2N *rp) @@ -347,6 +395,41 @@ void Pnewvect(NODE arg,VECT *rp) *rp = vect; } +void Pvect(NODE arg,VECT *rp) { + int len,i,r; + 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); @@ -429,6 +512,108 @@ void Pnewmat(NODE arg,MAT *rp) *rp = m; } +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; @@ -444,6 +629,21 @@ void Pvtol(NODE arg,LIST *rp) MKLIST(*rp,n); } +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; @@ -560,6 +760,10 @@ void Psize(NODE arg,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; } @@ -628,6 +832,46 @@ void Pinvmat(NODE arg,LIST *rp) 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; + MAT m,nm; + int *ri,*ci; + VECT rind,cind; + Q dn,q; + int i,j,k,l,row,col,t,rank; + + asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim"); + m = (MAT)ARG0(arg); + row = m->row; col = m->col; + 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) @@ -762,6 +1006,7 @@ int gauss_elim_mod(int **mat,int row,int col,int 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 mat,MAT *nm,Q *dn,int **rindp,int **cindp) { @@ -880,7 +1125,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); @@ -931,7 +1176,13 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn 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; @@ -949,7 +1200,7 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn } else wi[j] = 0; - rank = find_lhs_and_lu_mod(w,row,col,md,&rinfo,&cinfo); + rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo); 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++ ) @@ -979,9 +1230,17 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn *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] ) { @@ -991,17 +1250,19 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn 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); - get_eg(&tmp1); - add_eg(&eg_inv,&tmp0,&tmp1); - /* x = x-q*wc */ + /* 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, 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, 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++ ) { @@ -1019,18 +1280,50 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn 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] ) + 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); + } + 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; } } } @@ -1384,6 +1677,14 @@ void red_by_vect(int m,unsigned int *p,unsigned int *r } } +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, @@ -1472,6 +1773,50 @@ int generic_gauss_elim_mod(int **mat0,int row,int col, 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-- ) + if ( colstat[j] ) { + pivot = mat[l]; + for ( i = 0; i < l; i++ ) { + t = mat[i]; + if ( a = t[j] ) + red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j); + } + l--; + } + return rank; +} + /* LU decomposition; a[i][i] = 1/U[i][i] */ int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm) @@ -1580,7 +1925,7 @@ int find_lhs_and_lu_mod(unsigned int **a,int row,int c b = a^(-1)b */ -void solve_by_lu_mod(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; @@ -1613,8 +1958,12 @@ void solve_by_lu_mod(int **a,int n,int md,int **b,int 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]; } } @@ -1894,7 +2243,9 @@ void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp) TOGFMMAT(row,col,wmat,*rp); } -void Pgeninvm_swap(NODE arg,LIST *rp) +void Pgeninvm_swap(arg,rp) +NODE arg; +LIST *rp; { MAT m; pointer **mat; @@ -1940,8 +2291,12 @@ void Pgeninvm_swap(NODE arg,LIST *rp) } } -int gauss_elim_geninv_mod_swap(unsigned int **mat,int row,int col, - unsigned int md,unsigned int ***invmatp,int **indexp) +gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp) +unsigned int **mat; +int row,col; +unsigned int md; +unsigned int ***invmatp; +int **indexp; { int i,j,k,inv,a,n,m; unsigned int *t,*pivot,*s; @@ -1991,6 +2346,103 @@ int gauss_elim_geninv_mod_swap(unsigned int **mat,int 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); @@ -2185,6 +2637,48 @@ void Pirredpoly_up2(NODE arg,GF2N *rp) *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 @@ -2480,4 +2974,12 @@ void printimat(int **mat,int row,int 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); }