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.30 and 1.44

version 1.30, 2003/06/10 16:54:13

version 1.44, 2005/01/12 10:38:07

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.29 2003/06/09 16:18:09 saito Exp $

* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.43 2004/12/18 16:50:10 saito 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 68 void Pgeneric_gauss_elim();

Line 70 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 108 struct ftab array_tab[] = {

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

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

Line 161 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");

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

}else {

*rp = vect;

}

void PNBmul_gf2n(NODE arg,GF2N *rp)

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

Line 629 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 712 void Psize(NODE arg,LIST *rp)

Line 760 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 953 int gauss_elim_mod(int **mat,int row,int col,int md)

Line 1005 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 1071 RESET:

Line 1124 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 1122 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn

Line 1175 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 1170 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn

Line 1229 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 > 3 )

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 1182 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn

Line 1249 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 1210 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn

Line 1279 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 1575 void red_by_vect(int m,unsigned int *p,unsigned int *r

Line 1676 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 1663 int generic_gauss_elim_mod(int **mat0,int row,int col,

Line 1772 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 1771 int find_lhs_and_lu_mod(unsigned int **a,int row,int c

Line 1924 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 1804 void solve_by_lu_mod(int **a,int n,int md,int **b,int

Line 1957 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 2816 void printimat(int **mat,int row,int col)

Line 2973 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>