=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v retrieving revision 1.13 retrieving revision 1.23 diff -u -p -r1.13 -r1.23 --- OpenXM_contrib2/asir2000/builtin/array.c 2001/06/07 05:14:48 1.13 +++ OpenXM_contrib2/asir2000/builtin/array.c 2001/10/01 01:58:01 1.23 @@ -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.12 2001/06/07 04:54:38 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.22 2001/09/17 08:37:30 noro Exp $ */ #include "ca.h" #include "base.h" @@ -67,12 +67,14 @@ 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 **); +void reduce_sp_by_red_mod_compress (int *,CDP *,int *,int,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 Pinvmat(); void Pnewbytearray(); void Pgeneric_gauss_elim_mod(); @@ -90,6 +92,7 @@ void Px962_irredpoly_up2(); void Pirredpoly_up2(); void Pnbpoly_up2(); void Pqsort(); +void Pexponent_vector(); struct ftab array_tab[] = { {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4}, @@ -97,7 +100,10 @@ struct ftab array_tab[] = { {"mat_to_gfmmat",Pmat_to_gfmmat,2}, {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2}, {"newvect",Pnewvect,-2}, + {"vector",Pnewvect,-2}, + {"exponent_vector",Pexponent_vector,-99999999}, {"newmat",Pnewmat,-3}, + {"matrix",Pnewmat,-3}, {"newbytearray",Pnewbytearray,-2}, {"sepmat_destructive",Psepmat_destructive,2}, {"sepvect",Psepvect,2}, @@ -105,6 +111,7 @@ struct ftab array_tab[] = { {"vtol",Pvtol,1}, {"size",Psize,1}, {"det",Pdet,-2}, + {"invmat",Pinvmat,-2}, {"leqm",Pleqm,2}, {"leqm1",Pleqm1,2}, {"geninvm",Pgeninvm,2}, @@ -376,6 +383,13 @@ VECT *rp; *rp = vect; } +void Pexponent_vector(arg,rp) +NODE arg; +DP *rp; +{ + nodetod(arg,rp); +} + void Pnewbytearray(arg,rp) NODE arg; BYTEARRAY *rp; @@ -629,6 +643,41 @@ P *rp; } } +void Pinvmat(arg,rp) +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+... @@ -1292,7 +1341,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++; } } @@ -1379,25 +1428,114 @@ 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 +*/ + +int red_by_compress(m,p,r,ri,hc,len) +int m; +unsigned int *p; +register unsigned int *r; +register unsigned int *ri; +unsigned int hc; +register 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 */ + +int red_by_vect(m,p,r,hc,len) +int m; +unsigned int *p,*r; +unsigned int hc; +int len; +{ + register unsigned int up,lo; + unsigned int 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; + } +} + +extern unsigned int **psca; + +void reduce_sp_by_red_mod_compress (sp,redmat,ind,nred,col,md) +int *sp; +CDP *redmat; +int *ind; +int nred,col; +int md; +{ + int i,j,k,len; + unsigned int *tj; + CDP ri; + unsigned int hc,up,lo,up1,lo1,c; + 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] >= 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 generic_gauss_elim_mod(mat0,row,col,md,colstat) +int **mat0; int row,col,md; int *colstat; { int i,j,k,l,inv,a,rank,zzz; - int *t,*pivot,*pk,*tk; + unsigned int *t,*pivot,*pk,*tk; + 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,35 +1548,14 @@ int *colstat; inv = invm(pivot[j],md); for ( k = j, pk = pivot+k; k < col; k++, pk++ ) if ( *pk ) { + if ( *pk >= 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++; } @@ -1447,33 +1564,19 @@ 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++; - } - } + 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] >= md ) + t[k] %= md; + l++; } return rank; }