=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v retrieving revision 1.3 retrieving revision 1.19 diff -u -p -r1.3 -r1.19 --- OpenXM_contrib2/asir2000/builtin/array.c 2000/04/20 02:20:15 1.3 +++ OpenXM_contrib2/asir2000/builtin/array.c 2001/09/17 02:47:07 1.19 @@ -1,14 +1,63 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.2 2000/03/14 05:25:43 noro Exp $ */ +/* + * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED + * All rights reserved. + * + * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited, + * non-exclusive and royalty-free license to use, copy, modify and + * redistribute, solely for non-commercial and non-profit purposes, the + * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and + * conditions of this Agreement. For the avoidance of doubt, you acquire + * only a limited right to use the SOFTWARE hereunder, and FLL or any + * third party developer retains all rights, including but not limited to + * copyrights, in and to the SOFTWARE. + * + * (1) FLL does not grant you a license in any way for commercial + * purposes. You may use the SOFTWARE only for non-commercial and + * non-profit purposes only, such as academic, research and internal + * business use. + * (2) The SOFTWARE is protected by the Copyright Law of Japan and + * international copyright treaties. If you make copies of the SOFTWARE, + * with or without modification, as permitted hereunder, you shall affix + * to all such copies of the SOFTWARE the above copyright notice. + * (3) An explicit reference to this SOFTWARE and its copyright owner + * shall be made on your publication or presentation in any form of the + * results obtained by use of the SOFTWARE. + * (4) In the event that you modify the SOFTWARE, you shall notify FLL by + * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification + * for such modification or the source code of the modified part of the + * SOFTWARE. + * + * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL + * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND + * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS + * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES' + * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY + * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY. + * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT, + * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY + * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL + * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES + * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES + * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY + * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF + * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART + * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY + * 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.18 2001/09/17 01:18:34 noro Exp $ +*/ #include "ca.h" #include "base.h" #include "parse.h" #include "inline.h" -/* + +#if 0 #undef DMAR #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d); -*/ +#endif -extern int Print; /* XXX */ +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); @@ -18,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 Pnewbytearray(); void Pgeneric_gauss_elim_mod(); @@ -40,6 +91,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}, @@ -47,7 +99,11 @@ 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}, {"qsort",Pqsort,-2}, @@ -308,7 +364,7 @@ VECT *rp; asir_assert(ARG0(arg),O_N,"newvect"); len = QTOS((Q)ARG0(arg)); - if ( len <= 0 ) + if ( len < 0 ) error("newvect : invalid size"); MKVECT(vect,len); if ( argc(arg) == 2 ) { @@ -325,6 +381,57 @@ VECT *rp; *rp = vect; } +void Pexponent_vector(arg,rp) +NODE arg; +DP *rp; +{ + nodetod(arg,rp); +} + +void Pnewbytearray(arg,rp) +NODE arg; +BYTEARRAY *rp; +{ + int len,i,r; + BYTEARRAY array; + unsigned char *vb; + char *str; + LIST list; + NODE tn; + + 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 ) { + if ( !ARG1(arg) ) + error("newbytearray : invalid initialization"); + switch ( OID((Obj)ARG1(arg)) ) { + case O_LIST: + list = (LIST)ARG1(arg); + asir_assert(list,O_LIST,"newbytearray"); + for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ); + if ( r <= len ) { + for ( i = 0, tn = BDY(list), vb = BDY(array); tn; + i++, tn = NEXT(tn) ) + vb[i] = (unsigned char)QTOS((Q)BDY(tn)); + } + break; + case O_STR: + str = BDY((STRING)ARG1(arg)); + r = strlen(str); + if ( r <= len ) + bcopy(str,BDY(array),r); + break; + default: + if ( !ARG1(arg) ) + error("newbytearray : invalid initialization"); + } + } + *rp = array; +} + void Pnewmat(arg,rp) NODE arg; MAT *rp; @@ -339,7 +446,7 @@ MAT *rp; asir_assert(ARG0(arg),O_N,"newmat"); asir_assert(ARG1(arg),O_N,"newmat"); row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg)); - if ( row <= 0 || col <= 0 ) + if ( row < 0 || col < 0 ) error("newmat : invalid size"); MKMAT(m,row,col); if ( argc(arg) == 3 ) { @@ -663,8 +770,10 @@ int row,col,md; t = mat[i]; if ( i != j && (a = t[j]) ) for ( k = j, a = md - a; k <= n; k++ ) { + unsigned int tk; /* t[k] = dmar(pivot[k],a,t[k],md); */ - DMAR(pivot[k],a,t[k],md,t[k]) + DMAR(pivot[k],a,t[k],md,tk) + t[k] = tk; } } } @@ -675,7 +784,7 @@ int row,col,md; return -1; } -struct oEGT eg_mod,eg_elim,eg_chrem,eg_gschk,eg_intrat,eg_symb; +struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb; int generic_gauss_elim(mat,nm,dn,rindp,cindp) MAT mat; @@ -708,10 +817,10 @@ int **rindp,**cindp; colstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int)); for ( ind = 0; ; ind++ ) { - if ( Print ) { + if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); } - md = lprime[ind]; + md = get_lprime(ind); get_eg(&tmp0); for ( i = 0; i < row; i++ ) for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ ) @@ -745,13 +854,13 @@ RESET: } } else { if ( rank < rank0 ) { - if ( Print ) { + if ( DP_Print ) { fprintf(asir_out,"lower rank matrix; continuing...\n"); fflush(asir_out); } continue; } else if ( rank > rank0 ) { - if ( Print ) { + if ( DP_Print ) { fprintf(asir_out,"higher rank matrix; resetting...\n"); fflush(asir_out); } @@ -759,7 +868,7 @@ RESET: } else { for ( j = 0; (jrow; col = mat->col; w = (int **)almat(row,col); for ( ind = 0; ; ind++ ) { - md = lprime[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++ ) @@ -938,7 +1050,7 @@ int **rindp,**cindp; add_eg(&eg_mul,&tmp0,&tmp1); /* q = q*md */ mulq(q,mdq,&u); q = u; - if ( !(count % 2) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) { + if ( !(count % 16) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) { for ( j = k = l = 0; j < col; j++ ) if ( cinfo[j] ) rind[k++] = j; @@ -1148,16 +1260,278 @@ Q *dn; return 1; } -int generic_gauss_elim_mod(mat,row,col,md,colstat) +#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; +{ + int i,j,k,l,hc,zzz; + int *t,*s,*tj,*ind; + + /* reduce the reducers */ + ind = (int *)ALLOCA(row*sizeof(int)); + for ( i = 0; i < row; i++ ) { + t = mat[i]; + for ( j = 0; j < col && !t[j]; j++ ); + /* register the position of the head term */ + ind[i] = j; + for ( l = i-1; l >= 0; l-- ) { + /* reduce mat[i] by mat[l] */ + if ( hc = t[ind[l]] ) { + /* mat[i] = mat[i]-hc*mat[l] */ + j = ind[l]; + s = mat[l]+j; + tj = t+j; + hc = md-hc; + k = col-j; + for ( ; k >= 64; k -= 64 ) { + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1 + } + for ( ; k > 0; k-- ) { + if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; + } + } + } + } +} + +/* + mat[i] : reducers (i=0,...,nred-1) + spolys (i=nred,...,row-1) + mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order + 1. reduce the reducers + 2. reduce spolys by the reduced reducers +*/ + +void pre_reduce_mod(mat,row,col,nred,md) +int **mat; +int row,col,nred; +int md; +{ + int i,j,k,l,hc,inv; + int *t,*s,*tk,*ind; + +#if 1 + /* reduce the reducers */ + ind = (int *)ALLOCA(row*sizeof(int)); + for ( i = 0; i < nred; i++ ) { + /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */ + t = mat[i]; + for ( j = 0; j < col && !t[j]; j++ ); + /* register the position of the head term */ + ind[i] = j; + inv = invm(t[j],md); + for ( k = j; k < col; k++ ) + if ( t[k] ) + DMAR(t[k],inv,0,md,t[k]) + for ( l = i-1; l >= 0; l-- ) { + /* reduce mat[i] by mat[l] */ + if ( hc = t[ind[l]] ) { + /* mat[i] = mat[i]-hc*mat[l] */ + for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k; + k < col; k++, tk++, s++ ) + if ( *s ) + DMAR(*s,hc,*tk,md,*tk) + } + } + } + /* reduce the spolys */ + for ( i = nred; i < row; i++ ) { + t = mat[i]; + for ( l = nred-1; l >= 0; l-- ) { + /* reduce mat[i] by mat[l] */ + if ( hc = t[ind[l]] ) { + /* mat[i] = mat[i]-hc*mat[l] */ + for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k; + k < col; k++, tk++, s++ ) + if ( *s ) + DMAR(*s,hc,*tk,md,*tk) + } + } + } +#endif +} +/* + mat[i] : reducers (i=0,...,nred-1) + 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; +{ + int i,j,k,hc,zzz; + int *t,*s,*tj; + + /* reduce the spolys by redmat */ + for ( i = nred-1; i >= 0; i-- ) { + /* reduce sp by redmat[i] */ + if ( hc = sp[ind[i]] ) { + /* sp = sp-hc*redmat[i] */ + j = ind[i]; + hc = md-hc; + s = redmat[i]+j; + tj = sp+j; + for ( k = col-j; k > 0; k-- ) { + if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++; + } + } + } +} + +/* + rlist : reducers list + ht(BDY(rlist)) < ht(BDY(NEXT(rlist)) < ... w.r.t. the term order +*/ + +void reduce_reducers_mod_compress(rlist,nred,at,col,md,redmatp,indredp) +NODE rlist; +int nred; +DL *at; +int col,md; +CDP **redmatp; +int **indredp; +{ + CDP *redmat; + CDP t; + int *indred,*w; + int i,k; + NODE r; + + *redmatp = redmat = (CDP *)CALLOC(nred,sizeof(CDP)); + *indredp = indred = (int *)CALLOC(nred,sizeof(int)); + w = (int *)CALLOC(col,sizeof(int)); + + _dpmod_to_vect_compress(BDY(rlist),at,&redmat[0]); + indred[0] = redmat[0]->body[0].index; + + for ( i = 1, r = NEXT(rlist); i < nred; i++, r = NEXT(r) ) { + bzero(w,col*sizeof(int)); + _dpmod_to_vect(BDY(r),at,w); + reduce_sp_by_red_mod_compress(w,redmat,indred,i,col,md); + compress_vect(w,col,&redmat[i]); + indred[i] = redmat[i]->body[0].index; + } +} + +/* + 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,hc,len) +int m; +unsigned int *p; +register struct oCM *r; +unsigned int hc; +register int len; +{ + unsigned int up,lo; + unsigned int dmy; + unsigned int *pj; + + p[r->index] = 0; r++; + for ( len--; len; len--, r++ ) { + pj = p+r->index; + DMA0(r->c,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 ) { + DMA0(*r,hc,*p,up,lo); + if ( up ) { + DSAB(m,up,lo,dmy,*p); + } else + *p = lo; + } +} + +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; + struct oCM *rib; + + 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,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(mat0,row,col,md,colstat) +int **mat0; int row,col,md; int *colstat; { - int i,j,k,l,inv,a,rank; - int *t,*pivot; + int i,j,k,l,inv,a,rank,zzz; + 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; @@ -1168,17 +1542,16 @@ int *colstat; } pivot = mat[rank]; inv = invm(pivot[j],md); - for ( k = j; k < col; k++ ) - if ( pivot[k] ) { - DMAR(pivot[k],inv,0,md,pivot[k]) + 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] ) - for ( k = j, a = md - a; k < col; k++ ) - if ( pivot[k] ) { - DMAR(pivot[k],a,t[k],md,t[k]) - } + red_by_vect(md,t+j,pivot+j,md-a,col-j); } rank++; } @@ -1187,14 +1560,20 @@ int *colstat; pivot = mat[l]; for ( i = 0; i < l; i++ ) { t = mat[i]; + t[j] %= md; if ( a = t[j] ) - for ( k = j, a = md-a; k < col; k++ ) - if ( pivot[k] ) { - DMAR(pivot[k],a,t[k],md,t[k]) - } + 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; } @@ -1233,7 +1612,10 @@ int *perm; DMAR(inv,m,0,md,t[k]) for ( j = k+1, m = md - t[k]; j < col; j++ ) if ( pivot[j] ) { - DMAR(m,pivot[j],t[j],md,t[j]) + unsigned int tj; + + DMAR(m,pivot[j],t[j],md,tj) + t[j] = tj; } } } @@ -1288,7 +1670,9 @@ int **rinfo,**cinfo; DMAR(inv,m,0,md,t[k]) for ( j = k+1, m = md - t[k]; j < col; j++ ) if ( pivot[j] ) { - DMAR(m,pivot[j],t[j],md,t[j]) + unsigned int tj; + DMAR(m,pivot[j],t[j],md,tj) + t[j] = tj; } } } @@ -2237,7 +2621,7 @@ int row,col; for ( i = 0; i < row; i++ ) { for ( j = 0; j < col; j++ ) { - printnum(mat[i][j]); printf(" "); + printnum((Num)mat[i][j]); printf(" "); } printf("\n"); }