=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v retrieving revision 1.47 retrieving revision 1.58 diff -u -p -r1.47 -r1.58 --- OpenXM_contrib2/asir2000/builtin/array.c 2005/11/27 00:07:05 1.47 +++ OpenXM_contrib2/asir2000/builtin/array.c 2009/03/03 10:04:10 1.58 @@ -45,13 +45,19 @@ * 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.46 2005/02/08 18:06:05 saito Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.57 2009/02/03 00:39:23 noro Exp $ */ #include "ca.h" #include "base.h" #include "parse.h" #include "inline.h" +#include +#include +#if !defined(_MSC_VER) +#include +#endif + #define F4_INTRAT_PERIOD 8 #if 0 @@ -64,7 +70,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(); @@ -90,8 +96,10 @@ void Pvect(); void Pmat(); void Pmatc(); void Pnd_det(); +void Plu_mat(); struct ftab array_tab[] = { + {"lu_mat",Plu_mat,1}, {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4}, {"lu_gfmmat",Plu_gfmmat,2}, {"mat_to_gfmmat",Pmat_to_gfmmat,2}, @@ -107,6 +115,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}, @@ -384,11 +393,13 @@ void Pnewvect(NODE arg,VECT *rp) if ( argc(arg) == 2 ) { list = (LIST)ARG1(arg); asir_assert(list,O_LIST,"newvect"); +#if 0 for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ); if ( r > len ) { *rp = vect; return; } +#endif for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) ) vb[i] = (pointer)BDY(tn); } @@ -443,13 +454,42 @@ void Pnewbytearray(NODE arg,BYTEARRAY *rp) char *str; LIST list; NODE tn; + int ac; + struct stat sbuf; + char *fname; + FILE *fp; - 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 ) { + ac = argc(arg); + if ( ac == 1 ) { + if ( !OID((Obj)ARG0(arg)) ) error("newbytearray : invalid argument"); + switch ( OID((Obj)ARG0(arg)) ) { + case O_STR: + fname = BDY((STRING)ARG0(arg)); + fp = fopen(fname,"rb"); + if ( !fp ) error("newbytearray : fopen failed"); + if ( stat(fname,&sbuf) < 0 ) + error("newbytearray : stat failed"); + len = sbuf.st_size; + MKBYTEARRAY(array,len); + fread(BDY(array),len,sizeof(char),fp); + break; + case O_N: + if ( !RATN(ARG0(arg)) ) + error("newbytearray : invalid argument"); + len = QTOS((Q)ARG0(arg)); + if ( len < 0 ) + error("newbytearray : invalid size"); + MKBYTEARRAY(array,len); + break; + default: + error("newbytearray : invalid argument"); + } + } else if ( ac == 2 ) { + asir_assert(ARG0(arg),O_N,"newbytearray"); + len = QTOS((Q)ARG0(arg)); + if ( len < 0 ) + error("newbytearray : invalid size"); + MKBYTEARRAY(array,len); if ( !ARG1(arg) ) error("newbytearray : invalid initialization"); switch ( OID((Obj)ARG1(arg)) ) { @@ -473,10 +513,42 @@ void Pnewbytearray(NODE arg,BYTEARRAY *rp) if ( !ARG1(arg) ) error("newbytearray : invalid initialization"); } - } + } else + error("newbytearray : invalid argument"); *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; + unsigned 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; @@ -842,17 +914,34 @@ void Pinvmat(NODE arg,LIST *rp) 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); @@ -1173,6 +1262,121 @@ RESET: } } +/* XXX broken */ +int lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm) +{ + Q **a0,**b; + Q *aiq; + N **a; + N *ai; + Q q,q1,dn2,a1,q0,bik; + MAT m; + unsigned int md; + int n,ind,i,j,rank,t,inv,t1,ret,min,k; + int **w; + int *wi,*rinfo0,*rinfo; + N m1,m2,m3,u,s; + + a0 = (Q **)mat->body; + n = mat->row; + if ( n != mat->col ) + error("lu_dec_cr : non-square matrix"); + w = (int **)almat(n,n); + MKMAT(m,n,n); + a = (N **)m->body; + UTON(1,m1); + rinfo0 = 0; + ind = 0; + while ( 1 ) { + md = get_lprime(ind); + /* mat mod md */ + for ( i = 0; i < n; i++ ) + for ( j = 0, aiq = a0[i], wi = w[i]; j < n; j++ ) + if ( q = aiq[j] ) { + t = rem(NM(q),md); + if ( t && SGN(q) < 0 ) + t = (md - t) % md; + wi[j] = t; + } else + wi[j] = 0; + + if ( !lu_mod((unsigned int **)w,n,md,&rinfo) ) continue; + printf("."); fflush(stdout); + if ( !rinfo0 ) + *perm = rinfo0 = rinfo; + else { + for ( i = 0; i < n; i++ ) + if ( rinfo[i] != rinfo0[i] ) break; + if ( i < n ) continue; + } + if ( UNIN(m1) ) { + for ( i = 0; i < n; i++ ) + for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) { + UTON(wi[j],u); ai[j] = u; + } + UTON(md,m1); + } else { + inv = invm(rem(m1,md),md); + UTON(md,m2); muln(m1,m2,&m3); + for ( i = 0; i < n; i++ ) + for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) + if ( ai[i] ) { + /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */ + t = rem(ai[j],md); + if ( wi[j] >= t ) + t = wi[j]-t; + else + t = md-(t-wi[j]); + DMAR(t,inv,0,md,t1) + UTON(t1,u); + muln(m1,u,&s); + addn(ai[j],s,&u); ai[j] = u; + } else if ( wi[j] ) { + /* f3 = m1*(m1 mod m2)^(-1)*f2 */ + DMAR(wi[j],inv,0,md,t) + UTON(t,u); + muln(m1,u,&s); ai[j] = s; + } + m1 = m3; + } + if ( (++ind%8) == 0 ) { + ret = intmtoratm(m,m1,lu,dn); + if ( ret ) { + b = (Q **)lu->body; + mulq(*dn,*dn,&dn2); + for ( i = 0; i < n; i++ ) { + for ( j = 0; j < n; j++ ) { + q = 0; + min = MIN(i,j); + for ( k = 0; k <= min; k++ ) { + bik = k==i ? *dn : b[i][k]; + mulq(bik,b[k][j],&q0); + addq(q,q0,&q1); q = q1; + } + mulq(a0[rinfo0[i]][j],dn2,&q1); + if ( cmpq(q,q1) ) break; + } + if ( j < n ) break; + } + if ( i == n ) + return; + } + } + } +} + +int nmat(N **m,int n) +{ + int i,j; + + for ( i = 0; i < n; i++ ) { + for ( j = 0; j < n; j++ ) { + printn(m[i][j]); printf(" "); + } + printf("\n"); + } +} + int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp) { MAT bmat,xmat; @@ -1213,7 +1417,13 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn } else wi[j] = 0; + if ( DP_Print > 3 ) { + 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 > 3 ) { + 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++ ) @@ -1311,6 +1521,8 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn ret = intmtoratm_q(xmat,NM(q),*nmmat,dn); get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1); if ( ret ) { + rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); + cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); for ( j = k = l = 0; j < col; j++ ) if ( cinfo[j] ) rind[k++] = j; @@ -1328,6 +1540,11 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn print_eg("CHECK",&eg_check); fflush(asir_out); } + *rindp = rind; + *cindp = cind; + for ( j = k = 0; j < col; j++ ) + if ( !cinfo[j] ) + cind[k++] = j; return rank; } } else { @@ -1343,6 +1560,191 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn } } +int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT *nmmat,Q *dn,int **rindp,int **cindp) +{ + MAT bmat,xmat; + Q **a0,**a,**b,**x,**nm; + Q *ai,*bi,*xi; + int row,col; + int **w; + int *wi; + int **wc; + Q mdq,q,s,u; + N tn; + int ind,md,i,j,k,l,li,ri,rank; + unsigned int t; + int *cinfo,*rinfo; + int *rind,*cind; + int count; + 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; + NumberField nf; + DP m; + int col1; + + a0 = (Q **)mat->body; + row = mat->row; col = mat->col; + w = (int **)almat(row,col); + for ( ind = 0; ; 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++ ) + if ( q = (Q)ai[j] ) { + t = rem(NM(q),md); + if ( t && SGN(q) < 0 ) + t = (md - t) % md; + wi[j] = t; + } 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); + } + for ( i = 0; i < col-1; i++ ) { + if ( !cinfo[i] ) { + m = mb[i]; + for ( j = i+1; j < col-1; j++ ) + if ( dp_redble(mb[j],m) ) + cinfo[j] = -1; + } + } + 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++ ) + if ( cinfo[j] > 0 ) { + /* the column is in lhs */ + for ( i = 0; i < rank; i++ ) { + w[i][li] = w[i][j]; + a[i][li] = a0[rinfo[i]][j]; + } + li++; + } else if ( !cinfo[j] ) { + /* the column is in rhs */ + for ( i = 0; i < rank; i++ ) + b[i][ri] = a0[rinfo[i]][j]; + ri++; + } + + /* solve Ax+B=0; A: rank x rank, B: rank x ri */ + MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body; + MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body; + /* use the right part of w as work area */ + wc = (int **)almat(rank,ri); + for ( i = 0; i < rank; i++ ) + wc[i] = w[i]+rank; + *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int)); + *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int)); + init_eg(&eg_mul); init_eg(&eg_inv); + 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 ) + 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] ) { + t = rem(NM(u),md); + if ( t && SGN(u) > 0 ) + t = (md - t) % md; + wi[j] = t; + } else + wi[j] = 0; + /* 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, 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++ ) { + inner_product_mat_int_mod(a,wc,rank,i,j,&u); + addq(b[i][j],u,&s); + if ( s ) { + t = divin(NM(s),md,&tn); + if ( t ) + error("generic_gauss_elim_hensel:incosistent"); + NTOQ(tn,SGN(s),b[i][j]); + } else + b[i][j] = 0; + } + get_eg(&tmp1); + add_eg(&eg_mul,&tmp0,&tmp1); + /* q = q*md */ + mulq(q,mdq,&u); q = u; + 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] > 0 ) + rind[k++] = j; + else if ( !cinfo[j] ) + 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; + } + } + } + } +} + int f4_nocheck; int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind) @@ -1676,8 +2078,7 @@ void red_by_compress(int m,unsigned int *p,unsigned in void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len) { - register unsigned int up,lo; - unsigned int dmy; + unsigned int up,lo,dmy; *p++ = 0; r++; len--; for ( ; len; len--, r++, p++ ) @@ -1929,6 +2330,40 @@ int find_lhs_and_lu_mod(unsigned int **a,int row,int c return d; } +int lu_mod(unsigned int **a,int n,unsigned int md,int **rinfo) +{ + int i,j,k; + int *rp; + unsigned int *t,*pivot; + unsigned int inv,m; + + *rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int)); + for ( i = 0; i < n; i++ ) rp[i] = i; + for ( k = 0; k < n; k++ ) { + for ( i = k; i < n && !a[i][k]; i++ ); + if ( i == n ) return 0; + if ( i != k ) { + j = rp[i]; rp[i] = rp[k]; rp[k] = j; + t = a[i]; a[i] = a[k]; a[k] = t; + } + pivot = a[k]; + inv = invm(pivot[k],md); + for ( i = k+1; i < n; i++ ) { + t = a[i]; + if ( m = t[k] ) { + DMAR(inv,m,0,md,t[k]) + for ( j = k+1, m = md - t[k]; j < n; j++ ) + if ( pivot[j] ) { + unsigned int tj; + DMAR(m,pivot[j],t[j],md,tj) + t[j] = tj; + } + } + } + } + return 1; +} + /* Input a : n x n matrix; a result of LU-decomposition @@ -2192,6 +2627,28 @@ void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md, } } +void Plu_mat(NODE arg,LIST *rp) +{ + MAT m,lu; + Q dn; + Q *v; + int n,i; + int *iperm; + VECT perm; + NODE n0; + + asir_assert(ARG0(arg),O_MAT,"lu_mat"); + m = (MAT)ARG0(arg); + n = m->row; + MKMAT(lu,n,n); + lu_dec_cr(m,lu,&dn,&iperm); + MKVECT(perm,n); + for ( i = 0, v = (Q *)perm->body; i < n; i++ ) + STOQ(iperm[i],v[i]); + n0 = mknode(3,lu,dn,perm); + MKLIST(*rp,n0); +} + void Plu_gfmmat(NODE arg,LIST *rp) { MAT m; @@ -2203,8 +2660,8 @@ void Plu_gfmmat(NODE arg,LIST *rp) VECT perm; NODE n0; - asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat"); - asir_assert(ARG1(arg),O_N,"mat_to_gfmmat"); + asir_assert(ARG0(arg),O_MAT,"lu_gfmmat"); + asir_assert(ARG1(arg),O_N,"lu_gfmmat"); m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg)); mat_to_gfmmat(m,md,&mm); row = m->row;