=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v retrieving revision 1.29 retrieving revision 1.38 diff -u -p -r1.29 -r1.38 --- OpenXM_contrib2/asir2000/builtin/array.c 2003/06/09 16:18:09 1.29 +++ OpenXM_contrib2/asir2000/builtin/array.c 2004/09/21 05:23:13 1.38 @@ -45,13 +45,15 @@ * 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.28 2003/05/29 16:44:59 saito Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.37 2004/09/15 01:43:32 noro Exp $ */ #include "ca.h" #include "base.h" #include "parse.h" #include "inline.h" +#define F4_INTRAT_PERIOD 8 + #if 0 #undef DMAR #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d); @@ -68,7 +70,7 @@ 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 Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(), Pltov(); void Pgeninv_sf_swap(); void sepvect(); void Pmulmat_gf2n(); @@ -87,6 +89,7 @@ void Pmat_swap_col_destructive(); void Pvect(); void Pmat(); void Pmatc(); +void Pnd_det(); struct ftab array_tab[] = { {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4}, @@ -108,8 +111,10 @@ struct ftab array_tab[] = { {"sepvect",Psepvect,2}, {"qsort",Pqsort,-2}, {"vtol",Pvtol,1}, + {"ltov",Pltov,1}, {"size",Psize,1}, {"det",Pdet,-2}, + {"nd_det",Pnd_det,-2}, {"invmat",Pinvmat,-2}, {"leqm",Pleqm,2}, {"leqm1",Pleqm1,2}, @@ -156,12 +161,28 @@ int generic_comp_obj(Obj *a,Obj *b) void Pqsort(NODE arg,VECT *rp) { VECT vect; - NODE n; + NODE n,n1; P p; V v; + FUNC func; + int len,i; + pointer *a; + Obj t; - asir_assert(ARG0(arg),O_VECT,"qsort"); - vect = (VECT)ARG0(arg); + t = ARG0(arg); + if (OID(t) == O_LIST) { + n = (NODE)BDY((LIST)t); + len = length(n); + MKVECT(vect,len); + for ( i = 0; i < len; i++, n = NEXT(n) ) { + BDY(vect)[i] = BDY(n); + } + + }else if (OID(t) != O_VECT) { + error("qsort : invalid argument"); + }else { + vect = (VECT)t; + } if ( argc(arg) == 1 ) qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj); else { @@ -169,13 +190,25 @@ void Pqsort(NODE arg,VECT *rp) if ( !p || OID(p)!=2 ) error("qsort : invalid argument"); v = VR(p); - if ( (int)v->attr != V_SR ) - error("qsort : no such function"); - generic_comp_obj_func = (FUNC)v->priv; + gen_searchf(NAME(v),&func); + if ( !func ) { + if ( (int)v->attr != V_SR ) + error("qsort : no such function"); + func = (FUNC)v->priv; + } + generic_comp_obj_func = func; MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n); qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj); } - *rp = vect; + if (OID(t) == O_LIST) { + a = BDY(vect); + for ( i = len - 1, n = 0; i >= 0; i-- ) { + MKNODE(n1,a[i],n); n = n1; + } + MKLIST((LIST)*rp,n); + }else { + *rp = vect; + } } void PNBmul_gf2n(NODE arg,GF2N *rp) @@ -495,11 +528,21 @@ void Pmat(NODE arg, MAT *rp) } for (row = 0, tn = arg; tn; tn = NEXT(tn), row++); + if ( row == 1 ) { + if ( OID(ARG0(arg)) == O_MAT ) { + *rp=ARG0(arg); + return; + } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) { + error("mat : invalid argument"); + } + } if ( OID(ARG0(arg)) == O_VECT ) { v = ARG0(arg); col = v->len; } else if ( OID(ARG0(arg)) == O_LIST ) { for (col = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), col++); + } else { + error("mat : invalid argument"); } MKMAT(m,row,col); @@ -536,11 +579,21 @@ void Pmatc(NODE arg, MAT *rp) } for (col = 0, tn = arg; tn; tn = NEXT(tn), col++); + if ( col == 1 ) { + if ( OID(ARG0(arg)) == O_MAT ) { + *rp=ARG0(arg); + return; + } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) { + error("matc : invalid argument"); + } + } if ( OID(ARG0(arg)) == O_VECT ) { v = ARG0(arg); row = v->len; } else if ( OID(ARG0(arg)) == O_LIST ) { for (row = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), row++); + } else { + error("matc : invalid argument"); } MKMAT(m,row,col); @@ -576,6 +629,21 @@ void Pvtol(NODE arg,LIST *rp) MKLIST(*rp,n); } +void Pltov(NODE arg,VECT *rp) +{ + NODE n; + VECT v; + int len,i; + + asir_assert(ARG0(arg),O_LIST,"ltov"); + n = (NODE)BDY((LIST)ARG0(arg)); + len = length(n); + MKVECT(v,len); + for ( i = 0; i < len; i++, n = NEXT(n) ) + BDY(v)[i] = BDY(n); + *rp = v; +} + void Premainder(NODE arg,Obj *rp) { Obj a; @@ -933,6 +1001,7 @@ int gauss_elim_mod(int **mat,int row,int col,int md) } struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb; +struct oEGT eg_conv; int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp) { @@ -1051,7 +1120,7 @@ RESET: add_eg(&eg_chrem_split,&tmp0,&tmp1); get_eg(&tmp0); - if ( ind % 16 ) + if ( ind % F4_INTRAT_PERIOD ) ret = 0; else ret = intmtoratm(crmat,m1,*nm,dn); @@ -1190,7 +1259,7 @@ 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 % 16) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) { + if ( !(count % F4_INTRAT_PERIOD) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) { for ( j = k = l = 0; j < col; j++ ) if ( cinfo[j] ) rind[k++] = j; @@ -1555,6 +1624,14 @@ void red_by_vect(int m,unsigned int *p,unsigned int *r } } +void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len) +{ + *p++ = 0; r++; len--; + for ( ; len; len--, r++, p++ ) + if ( *r ) + *p = _addsf(_mulsf(*r,hc),*p); +} + extern unsigned int **psca; void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind, @@ -1643,6 +1720,50 @@ int generic_gauss_elim_mod(int **mat0,int row,int col, return rank; } +int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat) +{ + int i,j,k,l,inv,a,rank; + unsigned int *t,*pivot,*pk; + unsigned int **mat; + + mat = (unsigned int **)mat0; + for ( rank = 0, j = 0; j < col; j++ ) { + for ( i = rank; i < row; i++ ) + if ( mat[i][j] ) + break; + if ( i == row ) { + colstat[j] = 0; + continue; + } else + colstat[j] = 1; + if ( i != rank ) { + t = mat[i]; mat[i] = mat[rank]; mat[rank] = t; + } + pivot = mat[rank]; + inv = _invsf(pivot[j]); + for ( k = j, pk = pivot+k; k < col; k++, pk++ ) + if ( *pk ) + *pk = _mulsf(*pk,inv); + for ( i = rank+1; i < row; i++ ) { + t = mat[i]; + if ( a = t[j] ) + red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j); + } + rank++; + } + for ( j = col-1, l = rank-1; j >= 0; j-- ) + if ( colstat[j] ) { + pivot = mat[l]; + for ( i = 0; i < l; i++ ) { + t = mat[i]; + if ( a = t[j] ) + red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j); + } + l--; + } + return rank; +} + /* LU decomposition; a[i][i] = 1/U[i][i] */ int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm) @@ -2796,4 +2917,12 @@ void printimat(int **mat,int row,int col) } printf("\n"); } +} + +void Pnd_det(NODE arg,P *rp) +{ + if ( argc(arg) == 1 ) + nd_det(0,ARG0(arg),rp); + else + nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp); }