OpenXM_contrib2/asir2000/builtin/array.c - diff

Return to array.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Diff for /OpenXM_contrib2/asir2000/builtin/array.c between version 1.29 and 1.49

version 1.29, 2003/06/09 16:18:09

version 1.49, 2005/12/21 23:18:15

Line 45

* 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.48 2005/11/27 05:37:53 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);

Line 62 extern int DP_Print; /* XXX */

Line 64 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();

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();

Line 87 void Pmat_swap_col_destructive();

Line 89 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},

Line 104 struct ftab array_tab[] = {

Line 107 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},

{"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},

Line 153 int generic_comp_obj(Obj *a,Obj *b)

Line 159 int generic_comp_obj(Obj *a,Obj *b)

}

void Pqsort(NODE arg,VECT *rp)

void Pqsort(NODE arg,LIST *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");

t = ARG0(arg);

vect = (VECT)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 {

Line 169 void Pqsort(NODE arg,VECT *rp)

Line 191 void Pqsort(NODE arg,VECT *rp)

if ( !p || OID(p)!=2 )

error("qsort : invalid argument");

v = VR(p);

if ( (int)v->attr != V_SR )

gen_searchf(NAME(v),&func);

error("qsort : no such function");

if ( !func ) {

generic_comp_obj_func = (FUNC)v->priv;

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(*rp,n);

}else {

*rp = (LIST)vect;

}

void PNBmul_gf2n(NODE arg,GF2N *rp)

Line 444 void Pnewbytearray(NODE arg,BYTEARRAY *rp)

Line 478 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;

Line 495 void Pmat(NODE arg, MAT *rp)

Line 560 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);

Line 536 void Pmatc(NODE arg, MAT *rp)

Line 611 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);

Line 576 void Pvtol(NODE arg,LIST *rp)

Line 661 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;

Line 692 void Psize(NODE arg,LIST *rp)

Line 792 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;

}

Line 760 void Pinvmat(NODE arg,LIST *rp)

Line 864 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);

Line 803 void Pgeneric_gauss_elim(NODE arg,LIST *rp)

Line 925 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]}+...

Line 810 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)

Line 934 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;

int **wmat,**row0;

Q *rib,*cib;

Q *rib,*cib,*rnb;

int *colstat;

int *colstat,*p;

Q q;

int md,i,j,k,l,row,col,t,rank;

Line 823 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)

Line 947 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++ )

Line 835 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)

Line 963 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++ )

Line 852 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)

Line 987 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);

}

Line 933 int gauss_elim_mod(int **mat,int row,int col,int md)

Line 1068 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)

{

Line 1051 RESET:

Line 1187 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);

Line 1102 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn

Line 1238 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;

Line 1120 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn

Line 1262 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++ )

Line 1150 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn

Line 1298 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);

for ( q = ONE, count = 0; ; count++ ) {

init_eg(&eg_check); init_eg(&eg_intrat);

fprintf(stderr,".");

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] ) {

Line 1162 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn

Line 1318 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 */

/* wc = A^(-1)wc; wc is not normalized */

get_eg(&tmp0);

solve_by_lu_mod(w,rank,md,wc,ri,0);

solve_by_lu_mod(w,rank,md,wc,ri);

/* wx += q*wc */

get_eg(&tmp1);

ptr = wx;

add_eg(&eg_inv,&tmp0,&tmp1);

/* x = x-q*wc */

for ( i = 0; i < rank; i++ )

for ( j = 0, xi = x[i], wi = wc[i]; j < ri; j++ ) {

for ( j = 0, wi = wc[i]; j < ri; j++ ) {

STOQ(wi[j],u); mulq(q,u,&s);

if ( wi[j] )

subq(xi[j],s,&u); xi[j] = u;

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++ ) {

Line 1190 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn

Line 1348 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 == period ) {

for ( j = k = l = 0; j < col; j++ )

get_eg(&tmp0);

if ( cinfo[j] )

ptr = wx;

rind[k++] = j;

for ( i = 0; i < rank; i++ )

else

for ( j = 0, xi = x[i]; j < ri;

cind[l++] = j;

j++, ptr += nsize ) {

if ( gensolve_check(mat,*nmmat,*dn,rind,cind) ) {

for ( k = nsize-1; k >= 0 && !ptr[k]; k-- );

fprintf(stderr,"\n");

if ( k >= 0 ) {

print_eg("INV",&eg_inv);

wn = NALLOC(k+1);

print_eg("MUL",&eg_mul);

PL(wn) = k+1;

fflush(asir_out);

for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l];

return rank;

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;

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;

}

Line 1555 void red_by_vect(int m,unsigned int *p,unsigned int *r

Line 1745 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,

Line 1643 int generic_gauss_elim_mod(int **mat0,int row,int col,

Line 1841 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)

Line 1751 int find_lhs_and_lu_mod(unsigned int **a,int row,int c

Line 1993 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;

Line 1784 void solve_by_lu_mod(int **a,int n,int md,int **b,int

Line 2026 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++ )

if ( normalize )

b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);

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];

}

Line 2796 void printimat(int **mat,int row,int col)

Line 3042 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);

}

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>