version 1.27, 2003/01/06 01:16:37 |
version 1.38, 2004/09/21 05:23:13 |
|
|
* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, |
* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, |
* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. |
* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. |
* |
* |
* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.26 2002/02/06 00:55:03 noro Exp $ |
* $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.37 2004/09/15 01:43:32 noro Exp $ |
*/ |
*/ |
#include "ca.h" |
#include "ca.h" |
#include "base.h" |
#include "base.h" |
#include "parse.h" |
#include "parse.h" |
#include "inline.h" |
#include "inline.h" |
|
|
|
#define F4_INTRAT_PERIOD 8 |
|
|
#if 0 |
#if 0 |
#undef DMAR |
#undef DMAR |
#define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d); |
#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 Pgeneric_gauss_elim_mod(); |
|
|
void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat(); |
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 Pgeninv_sf_swap(); |
void sepvect(); |
void sepvect(); |
void Pmulmat_gf2n(); |
void Pmulmat_gf2n(); |
|
|
void Pexponent_vector(); |
void Pexponent_vector(); |
void Pmat_swap_row_destructive(); |
void Pmat_swap_row_destructive(); |
void Pmat_swap_col_destructive(); |
void Pmat_swap_col_destructive(); |
|
void Pvect(); |
|
void Pmat(); |
|
void Pmatc(); |
|
void Pnd_det(); |
|
|
struct ftab array_tab[] = { |
struct ftab array_tab[] = { |
{"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4}, |
{"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4}, |
Line 92 struct ftab array_tab[] = { |
|
Line 98 struct ftab array_tab[] = { |
|
{"generic_gauss_elim",Pgeneric_gauss_elim,1}, |
{"generic_gauss_elim",Pgeneric_gauss_elim,1}, |
{"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2}, |
{"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2}, |
{"newvect",Pnewvect,-2}, |
{"newvect",Pnewvect,-2}, |
|
{"vect",Pvect,-99999999}, |
{"vector",Pnewvect,-2}, |
{"vector",Pnewvect,-2}, |
{"exponent_vector",Pexponent_vector,-99999999}, |
{"exponent_vector",Pexponent_vector,-99999999}, |
{"newmat",Pnewmat,-3}, |
{"newmat",Pnewmat,-3}, |
{"matrix",Pnewmat,-3}, |
{"matrix",Pnewmat,-3}, |
|
{"mat",Pmat,-99999999}, |
|
{"matr",Pmat,-99999999}, |
|
{"matc",Pmatc,-99999999}, |
{"newbytearray",Pnewbytearray,-2}, |
{"newbytearray",Pnewbytearray,-2}, |
{"sepmat_destructive",Psepmat_destructive,2}, |
{"sepmat_destructive",Psepmat_destructive,2}, |
{"sepvect",Psepvect,2}, |
{"sepvect",Psepvect,2}, |
{"qsort",Pqsort,-2}, |
{"qsort",Pqsort,-2}, |
{"vtol",Pvtol,1}, |
{"vtol",Pvtol,1}, |
|
{"ltov",Pltov,1}, |
{"size",Psize,1}, |
{"size",Psize,1}, |
{"det",Pdet,-2}, |
{"det",Pdet,-2}, |
|
{"nd_det",Pnd_det,-2}, |
{"invmat",Pinvmat,-2}, |
{"invmat",Pinvmat,-2}, |
{"leqm",Pleqm,2}, |
{"leqm",Pleqm,2}, |
{"leqm1",Pleqm1,2}, |
{"leqm1",Pleqm1,2}, |
Line 149 int generic_comp_obj(Obj *a,Obj *b) |
|
Line 161 int generic_comp_obj(Obj *a,Obj *b) |
|
void Pqsort(NODE arg,VECT *rp) |
void Pqsort(NODE arg,VECT *rp) |
{ |
{ |
VECT vect; |
VECT vect; |
NODE n; |
NODE n,n1; |
P p; |
P p; |
V v; |
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 ) |
if ( argc(arg) == 1 ) |
qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj); |
qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj); |
else { |
else { |
Line 162 void Pqsort(NODE arg,VECT *rp) |
|
Line 190 void Pqsort(NODE arg,VECT *rp) |
|
if ( !p || OID(p)!=2 ) |
if ( !p || OID(p)!=2 ) |
error("qsort : invalid argument"); |
error("qsort : invalid argument"); |
v = VR(p); |
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); |
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); |
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) |
void PNBmul_gf2n(NODE arg,GF2N *rp) |
Line 355 void Pnewvect(NODE arg,VECT *rp) |
|
Line 395 void Pnewvect(NODE arg,VECT *rp) |
|
*rp = vect; |
*rp = vect; |
} |
} |
|
|
|
void Pvect(NODE arg,VECT *rp) { |
|
int len,i,r; |
|
VECT vect; |
|
pointer *vb; |
|
NODE tn; |
|
|
|
if ( !arg ) { |
|
*rp =0; |
|
return; |
|
} |
|
|
|
for (len = 0, tn = arg; tn; tn = NEXT(tn), len++); |
|
if ( len == 1 ) { |
|
if ( ARG0(arg) != 0 ) { |
|
switch ( OID(ARG0(arg)) ) { |
|
case O_VECT: |
|
*rp = ARG0(arg); |
|
return; |
|
case O_LIST: |
|
for ( len = 0, tn = ARG0(arg); tn; tn = NEXT(tn), len++ ); |
|
MKVECT(vect,len-1); |
|
for ( i = 0, tn = BDY((LIST)ARG0(arg)), vb =BDY(vect); |
|
tn; i++, tn = NEXT(tn) ) |
|
vb[i] = (pointer)BDY(tn); |
|
*rp=vect; |
|
return; |
|
} |
|
} |
|
} |
|
MKVECT(vect,len); |
|
for ( i = 0, tn = arg, vb = BDY(vect); tn; i++, tn = NEXT(tn) ) |
|
vb[i] = (pointer)BDY(tn); |
|
*rp = vect; |
|
} |
|
|
void Pexponent_vector(NODE arg,DP *rp) |
void Pexponent_vector(NODE arg,DP *rp) |
{ |
{ |
nodetod(arg,rp); |
nodetod(arg,rp); |
Line 437 void Pnewmat(NODE arg,MAT *rp) |
|
Line 512 void Pnewmat(NODE arg,MAT *rp) |
|
*rp = m; |
*rp = m; |
} |
} |
|
|
|
void Pmat(NODE arg, MAT *rp) |
|
{ |
|
int row,col; |
|
int i; |
|
MAT m; |
|
pointer **mb; |
|
pointer *ent; |
|
NODE tn, sn; |
|
VECT v; |
|
|
|
if ( !arg ) { |
|
*rp =0; |
|
return; |
|
} |
|
|
|
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); |
|
for (row = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), row++) { |
|
if ( BDY(tn) == 0 ) { |
|
error("mat : invalid argument"); |
|
} else if ( OID(BDY(tn)) == O_VECT ) { |
|
v = tn->body; |
|
ent = BDY(v); |
|
for (i = 0; i < v->len; i++ ) mb[row][i] = (Obj)ent[i]; |
|
} else if ( OID(BDY(tn)) == O_LIST ) { |
|
for (col = 0, sn = BDY((LIST)BDY(tn)); sn; col++, sn = NEXT(sn) ) |
|
mb[row][col] = (pointer)BDY(sn); |
|
} else { |
|
error("mat : invalid argument"); |
|
} |
|
} |
|
*rp = m; |
|
} |
|
|
|
void Pmatc(NODE arg, MAT *rp) |
|
{ |
|
int row,col; |
|
int i; |
|
MAT m; |
|
pointer **mb; |
|
pointer *ent; |
|
NODE tn, sn; |
|
VECT v; |
|
|
|
if ( !arg ) { |
|
*rp =0; |
|
return; |
|
} |
|
|
|
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); |
|
for (col = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), col++) { |
|
if ( BDY(tn) == 0 ) { |
|
error("matc : invalid argument"); |
|
} else if ( OID(BDY(tn)) == O_VECT ) { |
|
v = tn->body; |
|
ent = BDY(v); |
|
for (i = 0; i < v->len; i++ ) mb[i][col] = (Obj)ent[i]; |
|
} else if ( OID(BDY(tn)) == O_LIST ) { |
|
for (row = 0, sn = BDY((LIST)BDY(tn)); sn; row++, sn = NEXT(sn) ) |
|
mb[row][col] = (pointer)BDY(sn); |
|
} else { |
|
error("matc : invalid argument"); |
|
} |
|
} |
|
*rp = m; |
|
} |
|
|
void Pvtol(NODE arg,LIST *rp) |
void Pvtol(NODE arg,LIST *rp) |
{ |
{ |
NODE n,n1; |
NODE n,n1; |
Line 452 void Pvtol(NODE arg,LIST *rp) |
|
Line 629 void Pvtol(NODE arg,LIST *rp) |
|
MKLIST(*rp,n); |
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) |
void Premainder(NODE arg,Obj *rp) |
{ |
{ |
Obj a; |
Obj a; |
Line 809 int gauss_elim_mod(int **mat,int row,int col,int md) |
|
Line 1001 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_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) |
int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp) |
{ |
{ |
|
|
add_eg(&eg_chrem_split,&tmp0,&tmp1); |
add_eg(&eg_chrem_split,&tmp0,&tmp1); |
|
|
get_eg(&tmp0); |
get_eg(&tmp0); |
if ( ind % 16 ) |
if ( ind % F4_INTRAT_PERIOD ) |
ret = 0; |
ret = 0; |
else |
else |
ret = intmtoratm(crmat,m1,*nm,dn); |
ret = intmtoratm(crmat,m1,*nm,dn); |
Line 1066 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn |
|
Line 1259 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn |
|
add_eg(&eg_mul,&tmp0,&tmp1); |
add_eg(&eg_mul,&tmp0,&tmp1); |
/* q = q*md */ |
/* q = q*md */ |
mulq(q,mdq,&u); q = u; |
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++ ) |
for ( j = k = l = 0; j < col; j++ ) |
if ( cinfo[j] ) |
if ( cinfo[j] ) |
rind[k++] = j; |
rind[k++] = j; |
Line 1431 void red_by_vect(int m,unsigned int *p,unsigned int *r |
|
Line 1624 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; |
extern unsigned int **psca; |
|
|
void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind, |
void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind, |
Line 1519 int generic_gauss_elim_mod(int **mat0,int row,int col, |
|
Line 1720 int generic_gauss_elim_mod(int **mat0,int row,int col, |
|
return 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; |
|
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] */ |
/* LU decomposition; a[i][i] = 1/U[i][i] */ |
|
|
int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm) |
int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm) |
Line 2672 void printimat(int **mat,int row,int col) |
|
Line 2917 void printimat(int **mat,int row,int col) |
|
} |
} |
printf("\n"); |
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); |
} |
} |