=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v retrieving revision 1.40 retrieving revision 1.49 diff -u -p -r1.40 -r1.49 --- OpenXM_contrib2/asir2000/builtin/array.c 2004/12/02 13:53:31 1.40 +++ OpenXM_contrib2/asir2000/builtin/array.c 2005/12/21 23:18:15 1.49 @@ -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.39 2004/12/01 12:55:19 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.48 2005/11/27 05:37:53 noro Exp $ */ #include "ca.h" #include "base.h" @@ -64,7 +64,7 @@ extern int DP_Print; /* XXX */ void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(); void Pinvmat(); -void Pnewbytearray(); +void Pnewbytearray(),Pmemoryplot_to_coord(); void Pgeneric_gauss_elim(); void Pgeneric_gauss_elim_mod(); @@ -107,6 +107,7 @@ struct ftab array_tab[] = { {"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}, @@ -158,7 +159,7 @@ int generic_comp_obj(Obj *a,Obj *b) } -void Pqsort(NODE arg,VECT *rp) +void Pqsort(NODE arg,LIST *rp) { VECT vect; NODE n,n1; @@ -205,9 +206,9 @@ void Pqsort(NODE arg,VECT *rp) for ( i = len - 1, n = 0; i >= 0; i-- ) { MKNODE(n1,a[i],n); n = n1; } - MKLIST((LIST)*rp,n); + MKLIST(*rp,n); }else { - *rp = vect; + *rp = (LIST)vect; } } @@ -477,6 +478,37 @@ void Pnewbytearray(NODE arg,BYTEARRAY *rp) *rp = array; } +#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; + 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; @@ -760,6 +792,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; } @@ -828,26 +864,44 @@ 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,R,C] + 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] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+... + 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; + NODE n0,opt,p; MAT m,nm; int *ri,*ci; VECT rind,cind; Q dn,q; int i,j,k,l,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; - rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci); + 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); @@ -871,6 +925,8 @@ void Pgeneric_gauss_elim(NODE arg,LIST *rp) 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]}+... */ @@ -878,11 +934,11 @@ 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,rank; @@ -891,6 +947,10 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp) 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++ ) @@ -903,6 +963,13 @@ void Pgeneric_gauss_elim_mod(NODE arg,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++ ) @@ -920,7 +987,7 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp) } else { STOQ(j,cib[l]); l++; } - n0 = mknode(3,mat,rind,cind); + n0 = mknode(4,mat,rind,cind,rnum); MKLIST(*rp,n0); } @@ -1171,8 +1238,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; @@ -1190,7 +1262,13 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn } 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); + } 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++ ) @@ -1220,11 +1298,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); + init_eg(&eg_check); init_eg(&eg_intrat); period = F4_INTRAT_PERIOD; - for ( q = ONE, count = 0; ; count++ ) { + 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,"."); + 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] ) { @@ -1234,17 +1318,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++ ) { @@ -1262,24 +1348,52 @@ 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 % period) ) - if ( intmtoratm_q(xmat,NM(q),*nmmat,dn) ) { + 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; - if ( gensolve_check(mat,*nmmat,*dn,rind,cind) ) { - if ( DP_Print ) { + 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 *=2; + } 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; + } + } } } } @@ -1879,7 +1993,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; @@ -1912,8 +2026,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]; } }