=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v retrieving revision 1.24 retrieving revision 1.27 diff -u -p -r1.24 -r1.27 --- OpenXM_contrib2/asir2000/builtin/array.c 2001/10/09 01:36:05 1.24 +++ OpenXM_contrib2/asir2000/builtin/array.c 2003/01/06 01:16:37 1.27 @@ -45,7 +45,7 @@ * 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.26 2002/02/06 00:55:03 noro Exp $ */ #include "ca.h" #include "base.h" @@ -64,10 +64,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 Pgeninv_sf_swap(); void sepvect(); void Pmulmat_gf2n(); void Pbconvmat_gf2n(); @@ -80,11 +82,14 @@ void Pirredpoly_up2(); void Pnbpoly_up2(); void Pqsort(); void Pexponent_vector(); +void Pmat_swap_row_destructive(); +void Pmat_swap_col_destructive(); 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}, {"vector",Pnewvect,-2}, @@ -103,6 +108,7 @@ struct ftab array_tab[] = { {"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 +119,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}, }; @@ -635,6 +643,45 @@ void Pinvmat(NODE arg,LIST *rp) B[I] <-> 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) + C : a vector of length col-rank(A) + B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+... +*/ + void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp) { NODE n0; @@ -949,7 +996,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++ ) @@ -1894,7 +1941,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 +1989,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 +2044,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 +2335,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