=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v retrieving revision 1.54 retrieving revision 1.75 diff -u -p -r1.54 -r1.75 --- OpenXM_contrib2/asir2000/builtin/array.c 2006/06/17 10:12:06 1.54 +++ OpenXM_contrib2/asir2000/builtin/array.c 2017/09/17 02:34:02 1.75 @@ -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.53 2006/06/12 11:52:10 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.74 2017/09/15 01:52:51 noro Exp $ */ #include "ca.h" #include "base.h" @@ -54,7 +54,9 @@ #include #include +#if !defined(_MSC_VER) #include +#endif #define F4_INTRAT_PERIOD 8 @@ -66,13 +68,15 @@ extern int DP_Print; /* XXX */ -void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(); +void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(), Ptriangleq(); void Pinvmat(); void Pnewbytearray(),Pmemoryplot_to_coord(); void Pgeneric_gauss_elim(); void Pgeneric_gauss_elim_mod(); +void Pindep_rows_mod(); + void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat(); void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(), Pltov(); void Pgeninv_sf_swap(); @@ -95,6 +99,9 @@ void Pmat(); void Pmatc(); void Pnd_det(); void Plu_mat(); +void Pmat_col(); +void Plusolve_prep(); +void Plusolve_main(); struct ftab array_tab[] = { {"lu_mat",Plu_mat,1}, @@ -103,6 +110,7 @@ struct ftab array_tab[] = { {"mat_to_gfmmat",Pmat_to_gfmmat,2}, {"generic_gauss_elim",Pgeneric_gauss_elim,1}, {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2}, + {"indep_rows_mod",Pindep_rows_mod,2}, {"newvect",Pnewvect,-2}, {"vect",Pvect,-99999999}, {"vector",Pnewvect,-2}, @@ -140,9 +148,211 @@ struct ftab array_tab[] = { {"nbpoly_up2",Pnbpoly_up2,2}, {"mat_swap_row_destructive",Pmat_swap_row_destructive,3}, {"mat_swap_col_destructive",Pmat_swap_col_destructive,3}, + {"mat_col",Pmat_col,2}, + {"lusolve_prep",Plusolve_prep,1}, + {"lusolve_main",Plusolve_main,1}, + {"triangleq",Ptriangleq,1}, {0,0,0}, }; +typedef struct _ent { int j; unsigned int e; } ent; + +ent *get_row(FILE *,int *l); +void put_row(FILE *out,int l,ent *a); +void lu_elim(int *l,ent **a,int k,int i,int mul,int mod); +void lu_append(int *,ent **,int *,int,int,int); +void solve_l(int *,ent **,int,int *,int); +void solve_u(int *,ent **,int,int *,int); + + +static int *ul,*ll; +static ent **u,**l; +static int modulus; + +void Plusolve_prep(NODE arg,Q *rp) +{ + char *fname; + FILE *in; + int len,i,rank; + int *rhs; + + fname = BDY((STRING)ARG0(arg)); + in = fopen(fname,"r"); + modulus = getw(in); + len = getw(in); + ul = (int *)MALLOC_ATOMIC(len*sizeof(int)); + u = (ent **)MALLOC(len*sizeof(ent *)); + ll = (int *)MALLOC_ATOMIC(len*sizeof(int)); + l = (ent **)MALLOC(len*sizeof(ent *)); + for ( i = 0; i < len; i++ ) { + u[i] = get_row(in,&ul[i]); + } + for ( i = 0; i < len; i++ ) { + l[i] = get_row(in,&ll[i]); + } + fclose(in); + *rp = ONE; +} + +void Plusolve_main(NODE arg,VECT *rp) +{ + Q *d,*p; + VECT v,r; + int len,i; + int *rhs; + + v = (VECT)ARG0(arg); len = v->len; + d = (Q *)BDY(v); + rhs = (int *)MALLOC_ATOMIC(len*sizeof(int)); + for ( i = 0; i < len; i++ ) rhs[i] = QTOS(d[i]); + solve_l(ll,l,len,rhs,modulus); + solve_u(ul,u,len,rhs,modulus); + NEWVECT(r); r->len = len; + r->body = (pointer *)MALLOC(len*sizeof(pointer)); + p = (Q *)r->body; + for ( i = 0; i < len; i++ ) + STOQ(rhs[i],p[i]); + *rp = r; +} + +ent *get_row(FILE *in,int *l) +{ + int len,i; + ent *a; + + *l = len = getw(in); + a = (ent *)MALLOC_ATOMIC(len*sizeof(ent)); + for ( i = 0; i < len; i++ ) { + a[i].j = getw(in); + a[i].e = getw(in); + } + return a; +} + +void lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int mod) +{ + int i,j,k,s,mul; + unsigned int inv; + int *ll2; + + ll2 = (int *)MALLOC_ATOMIC(n*sizeof(int)); + for ( i = 0; i < n; i++ ) ll2[i] = 0; + for ( i = 0; i < n; i++ ) { + fprintf(stderr,"i=%d\n",i); + inv = invm(u[i][0].e,mod); + for ( k = i+1; k < n; k++ ) + if ( u[k][0].j == n-i ) { + s = u[k][0].e; + DMAR(s,inv,0,mod,mul); + lu_elim(ul,u,k,i,mul,mod); + lu_append(ll,l,ll2,k,i,mul); + } + } +} + +#define INITLEN 10 + +void lu_append(int *l,ent **a,int *l2,int k,int i,int mul) +{ + int len; + ent *p; + + len = l[k]; + if ( !len ) { + a[k] = p = (ent *)MALLOC_ATOMIC(INITLEN*sizeof(ent)); + p[0].j = i; p[0].e = mul; + l[k] = 1; l2[k] = INITLEN; + } else { + if ( l2[k] == l[k] ) { + l2[k] *= 2; + a[k] = REALLOC(a[k],l2[k]*sizeof(ent)); + } + p =a[k]; + p[l[k]].j = i; p[l[k]].e = mul; + l[k]++; + } +} + +/* a[k] = a[k]-mul*a[i] */ + +void lu_elim(int *l,ent **a,int k,int i,int mul,int mod) +{ + ent *ak,*ai,*w; + int lk,li,j,m,p,q,r,s,t,j0; + + ak = a[k]; ai = a[i]; lk = l[k]; li = l[i]; + w = (ent *)alloca((lk+li)*sizeof(ent)); + p = 0; q = 0; j = 0; + mul = mod-mul; + while ( p < lk && q < li ) { + if ( ak[p].j > ai[q].j ) { + w[j] = ak[p]; j++; p++; + } else if ( ak[p].j < ai[q].j ) { + w[j].j = ai[q].j; + t = ai[q].e; + DMAR(t,mul,0,mod,r); + w[j].e = r; + j++; q++; + } else { + t = ai[q].e; s = ak[p].e; + DMAR(t,mul,s,mod,r); + if ( r ) { + w[j].j = ai[q].j; w[j].e = r; j++; + } + p++; q++; + } + } + if ( q == li ) + while ( p < lk ) { + w[j] = ak[p]; j++; p++; + } + else if ( p == lk ) + while ( q < li ) { + w[j].j = ai[q].j; + t = ai[q].e; + DMAR(t,mul,0,mod,r); + w[j].e = r; + j++; q++; + } + if ( j <= lk ) { + for ( m = 0; m < j; m++ ) ak[m] = w[m]; + } else { + a[k] = ak = (ent *)MALLOC_ATOMIC(j*sizeof(ent)); + for ( m = 0; m < j; m++ ) ak[m] = w[m]; + } + l[k] = j; +} + +void solve_l(int *ll,ent **l,int n,int *rhs,int mod) +{ + int j,k,s,len; + ent *p; + + for ( j = 0; j < n; j++ ) { + len = ll[j]; p = l[j]; + for ( k = 0, s = 0; k < len; k++ ) + s = dmar(p[k].e,rhs[p[k].j],s,mod); + rhs[j] -= s; + if ( rhs[j] < 0 ) rhs[j] += mod; + } +} + +void solve_u(int *ul,ent **u,int n,int *rhs,int mod) +{ + int j,k,s,len,inv; + ent *p; + + for ( j = n-1; j >= 0; j-- ) { + len = ul[j]; p = u[j]; + for ( k = 1, s = 0; k < len; k++ ) + s = dmar(p[k].e,rhs[p[k].j],s,mod); + rhs[j] -= s; + if ( rhs[j] < 0 ) rhs[j] += mod; + inv = invm((unsigned int)p[0].e,mod); + rhs[j] = dmar(rhs[j],inv,0,mod); + } +} + int comp_obj(Obj *a,Obj *b) { return arf_comp(CO,*a,*b); @@ -150,6 +360,7 @@ int comp_obj(Obj *a,Obj *b) static FUNC generic_comp_obj_func; static NODE generic_comp_obj_arg; +static NODE generic_comp_obj_option; int generic_comp_obj(Obj *a,Obj *b) { @@ -157,7 +368,7 @@ int generic_comp_obj(Obj *a,Obj *b) BDY(generic_comp_obj_arg)=(pointer)(*a); BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b); - r = (Q)bevalf(generic_comp_obj_func,generic_comp_obj_arg); + r = (Q)bevalf_with_opts(generic_comp_obj_func,generic_comp_obj_arg,generic_comp_obj_option); if ( !r ) return 0; else @@ -204,7 +415,8 @@ void Pqsort(NODE arg,LIST *rp) func = (FUNC)v->priv; } generic_comp_obj_func = func; - MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n); + MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n); + generic_comp_obj_option = current_option; qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj); } if (OID(t) == O_LIST) { @@ -338,8 +550,8 @@ void Psepmat_destructive(NODE arg,LIST *rp) sgn = SGN(ent); divn(nm,mod,&quo,&rem); /* if ( quo != nm && rem != nm ) */ -/* GC_free(nm); */ -/* GC_free(ent); */ +/* GCFREE(nm); */ +/* GCFREE(ent); */ NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]); } MKNODE(n1,mat1,0); MKNODE(n0,mat,n1); @@ -391,11 +603,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); } @@ -403,7 +617,7 @@ void Pnewvect(NODE arg,VECT *rp) } void Pvect(NODE arg,VECT *rp) { - int len,i,r; + int len,i; VECT vect; pointer *vb; NODE tn; @@ -457,15 +671,29 @@ void Pnewbytearray(NODE arg,BYTEARRAY *rp) ac = argc(arg); if ( ac == 1 ) { - /* ARG0(arg) must be a filename */ - asir_assert(ARG0(arg),O_STR,"newbytearray"); - 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); + 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)); @@ -675,6 +903,10 @@ void Pvtol(NODE arg,LIST *rp) pointer *a; int len,i; + if ( OID(ARG0(arg)) == O_LIST ) { + *rp = ARG0(arg); + return; + } asir_assert(ARG0(arg),O_VECT,"vtol"); v = (VECT)ARG0(arg); len = v->len; a = BDY(v); for ( i = len - 1, n = 0; i >= 0; i-- ) { @@ -686,9 +918,18 @@ void Pvtol(NODE arg,LIST *rp) void Pltov(NODE arg,VECT *rp) { NODE n; - VECT v; + VECT v,v0; int len,i; + if ( OID(ARG0(arg)) == O_VECT ) { + v0 = (VECT)ARG0(arg); len = v0->len; + MKVECT(v,len); + for ( i = 0; i < len; i++ ) { + BDY(v)[i] = BDY(v0)[i]; + } + *rp = v; + return; + } asir_assert(ARG0(arg),O_LIST,"ltov"); n = (NODE)BDY((LIST)ARG0(arg)); len = length(n); @@ -901,7 +1142,7 @@ void Pgeneric_gauss_elim(NODE arg,LIST *rp) int *ri,*ci; VECT rind,cind; Q dn,q; - int i,j,k,l,row,col,t,rank; + int i,row,col,t,rank; int is_hensel = 0; char *key; Obj value; @@ -939,6 +1180,46 @@ void Pgeneric_gauss_elim(NODE arg,LIST *rp) MKLIST(*rp,n0); } +void Pindep_rows_mod(NODE arg,VECT *rp) +{ + MAT m,mat; + VECT rind; + Q **tmat; + int **wmat,**row0; + Q *rib; + int *rowstat,*p; + Q q; + int md,i,j,k,l,row,col,t,rank; + + asir_assert(ARG0(arg),O_MAT,"indep_rows_mod"); + asir_assert(ARG1(arg),O_N,"indep_rows_mod"); + 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]; + + rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int)); + for ( i = 0; i < row; i++ ) + for ( j = 0; j < col; j++ ) + if ( q = (Q)tmat[i][j] ) { + t = rem(NM(q),md); + if ( t && SGN(q) < 0 ) + t = (md - t) % md; + wmat[i][j] = t; + } else + wmat[i][j] = 0; + rank = indep_rows_mod(wmat,row,col,md,rowstat); + + MKVECT(rind,rank); + rib = (Q *)rind->body; + for ( j = 0; j < rank; j++ ) { + STOQ(rowstat[j],rib[j]); + } + *rp = rind; +} + /* input : a row x col matrix A A[I] <-> A[I][0]*x_0+A[I][1]*x_1+... @@ -1244,8 +1525,10 @@ RESET: } } +void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm); + /* XXX broken */ -int lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm) +void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm) { Q **a0,**b; Q *aiq; @@ -1347,7 +1630,7 @@ int lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm) } } -int nmat(N **m,int n) +void nmat(N **m,int n) { int i,j; @@ -1542,7 +1825,7 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn } } -int generic_gauss_elim_hensel_dalg(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp) +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; @@ -1566,12 +1849,9 @@ int generic_gauss_elim_hensel_dalg(MAT mat,MAT *nmmat, N wn; Q wq; NumberField nf; - DP *mb; DP m; int col1; - nf = get_numberfield(); - mb = nf->mb; a0 = (Q **)mat->body; row = mat->row; col = mat->col; w = (int **)almat(row,col); @@ -2063,8 +2343,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++ ) @@ -2077,6 +2356,24 @@ void red_by_vect(int m,unsigned int *p,unsigned int *r } } +#if defined(__GNUC__) && SIZEOF_LONG==8 +/* 64bit vector += UNIT vector(normalized) */ + +void red_by_vect64(int m, U64 *p,unsigned int *c,U64 *r,unsigned int hc,int len) +{ + U64 t; + + /* (p[0],c[0]) is normalized */ + *p++ = 0; *c++ = 0; r++; len--; + for ( ; len; len--, r++, p++, c++ ) + if ( *r ) { + t = (*p)+(*r)*hc; + if ( t < *p ) (*c)++; + *p = t; + } +} +#endif + void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len) { *p++ = 0; r++; len--; @@ -2085,6 +2382,22 @@ void red_by_vect_sf(int m,unsigned int *p,unsigned int *p = _addsf(_mulsf(*r,hc),*p); } +extern GZ current_mod_lf; +extern int current_mod_lf_size; + +void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len) +{ + mpz_set_ui(*p++,0); r++; len--; + for ( ; len; len--, r++, p++ ) { + mpz_addmul(*p,*r,hc); +#if 0 + if ( mpz_size(*p) > current_mod_lf_size ) + mpz_mod(*p,*p,BDY(current_mod_lf)); +#endif + } +} + + extern unsigned int **psca; void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind, @@ -2173,6 +2486,103 @@ int generic_gauss_elim_mod(int **mat0,int row,int col, return rank; } +int generic_gauss_elim_mod2(int **mat0,int row,int col,int md,int *colstat,int *rowstat) +{ + int i,j,k,l,inv,a,rank; + unsigned int *t,*pivot,*pk; + unsigned int **mat; + + for ( i = 0; i < row; i++ ) rowstat[i] = i; + mat = (unsigned int **)mat0; + for ( rank = 0, j = 0; j < col; j++ ) { + 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; + } else + colstat[j] = 1; + if ( i != rank ) { + t = mat[i]; mat[i] = mat[rank]; mat[rank] = t; + k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k; + } + pivot = mat[rank]; + inv = invm(pivot[j],md); + for ( k = j, pk = pivot+k; k < col; k++, pk++ ) + if ( *pk ) { + if ( *pk >= (unsigned int)md ) + *pk %= md; + DMAR(*pk,inv,0,md,*pk) + } + for ( i = rank+1; i < row; i++ ) { + t = mat[i]; + if ( a = t[j] ) + red_by_vect(md,t+j,pivot+j,md-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]; + 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] >= (unsigned int)md ) + t[k] %= md; + l++; + } + return rank; +} + +int indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat) +{ + int i,j,k,l,inv,a,rank; + unsigned int *t,*pivot,*pk; + unsigned int **mat; + + for ( i = 0; i < row; i++ ) rowstat[i] = i; + mat = (unsigned int **)mat0; + for ( rank = 0, j = 0; j < col; j++ ) { + for ( i = rank; i < row; i++ ) + mat[i][j] %= md; + for ( i = rank; i < row; i++ ) + if ( mat[i][j] ) + break; + if ( i == row ) continue; + if ( i != rank ) { + t = mat[i]; mat[i] = mat[rank]; mat[rank] = t; + k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k; + } + pivot = mat[rank]; + inv = invm(pivot[j],md); + for ( k = j, pk = pivot+k; k < col; k++, pk++ ) + if ( *pk ) { + if ( *pk >= (unsigned int)md ) + *pk %= md; + DMAR(*pk,inv,0,md,*pk) + } + for ( i = rank+1; i < row; i++ ) { + t = mat[i]; + if ( a = t[j] ) + red_by_vect(md,t+j,pivot+j,md-a,col-j); + } + rank++; + } + 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; @@ -2699,9 +3109,7 @@ void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp) TOGFMMAT(row,col,wmat,*rp); } -void Pgeninvm_swap(arg,rp) -NODE arg; -LIST *rp; +void Pgeninvm_swap(NODE arg,LIST *rp) { MAT m; pointer **mat; @@ -2747,12 +3155,8 @@ LIST *rp; } } -gauss_elim_geninv_mod_swap(mat,row,col,md,invmatp,indexp) -unsigned int **mat; -int row,col; -unsigned int md; -unsigned int ***invmatp; -int **indexp; +int gauss_elim_geninv_mod_swap(unsigned int **mat,int row,int col,unsigned int md, + unsigned int ***invmatp,int **indexp) { int i,j,k,inv,a,n,m; unsigned int *t,*pivot,*s; @@ -2945,10 +3349,10 @@ void inner_product_int(Q *a,Q *b,int n,Q *r) t = wma; wma = sum; sum = t; } } - GC_free(wm); - GC_free(wma); + GCFREE(wm); + GCFREE(wma); if ( !sgn ) { - GC_free(sum); + GCFREE(sum); *r = 0; } else NTOQ(sum,sgn,*r); @@ -3003,10 +3407,10 @@ void inner_product_mat_int_mod(Q **a,int **b,int n,int t = wma; wma = sum; sum = t; } } - GC_free(wm); - GC_free(wma); + GCFREE(wm); + GCFREE(wma); if ( !sgn ) { - GC_free(sum); + GCFREE(sum); *r = 0; } else NTOQ(sum,sgn,*r); @@ -3438,4 +3842,66 @@ void Pnd_det(NODE arg,P *rp) nd_det(0,ARG0(arg),rp); else nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp); +} + +void Pmat_col(NODE arg,VECT *rp) +{ + int i,j,n; + MAT mat; + VECT vect; + + asir_assert(ARG0(arg),O_MAT,"mat_col"); + asir_assert(ARG1(arg),O_N,"mat_col"); + mat = (MAT)ARG0(arg); + j = QTOS((Q)ARG1(arg)); + if ( j < 0 || j >= mat->col) { + error("mat_col : Out of range"); + } + n = mat->row; + MKVECT(vect,n); + for(i=0; iv; + for ( k = i+1; k < n; k++ ) + if ( p[k] ) { + if ( OID(p[k]) == O_N ) return 0; + if ( p[k]->v == v ) p[k] = 0; + } + i++; + } + for ( r = 0, i = 0; i < n; i++ ) { + if ( p[i] ) { + MKNODE(r1,p[i],r); r = r1; + } + } + return r; +} + +void Ptriangleq(NODE arg,LIST *rp) +{ + NODE ret; + + asir_assert(ARG0(arg),O_LIST,"sparseleq"); + ret = triangleq(BDY((LIST)ARG0(arg))); + MKLIST(*rp,ret); }