version 1.111, 2004/09/21 07:19:01 |
version 1.117, 2004/12/03 08:57:30 |
|
|
/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.110 2004/09/21 05:23:14 noro Exp $ */ |
/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.116 2004/12/01 12:36:17 noro Exp $ */ |
|
|
#include "nd.h" |
#include "nd.h" |
|
|
|
|
ND _nd_free_list; |
ND _nd_free_list; |
ND_pairs _ndp_free_list; |
ND_pairs _ndp_free_list; |
|
|
|
static int nd_nalg; |
#if 0 |
#if 0 |
static int ndv_alloc; |
static int ndv_alloc; |
#endif |
#endif |
Line 41 static int nd_found,nd_create,nd_notfirst; |
|
Line 42 static int nd_found,nd_create,nd_notfirst; |
|
static int nmv_adv; |
static int nmv_adv; |
static int nd_demand; |
static int nd_demand; |
|
|
|
UINT *nd_det_compute_bound(NDV **dm,int n,int j); |
|
void nd_det_reconstruct(NDV **dm,int n,int j,NDV d); |
|
ND nd_pseudo_monic(int m,ND p); |
|
|
void nd_free_private_storage() |
void nd_free_private_storage() |
{ |
{ |
_nm_free_list = 0; |
_nm_free_list = 0; |
Line 724 int ndl_disjoint(UINT *d1,UINT *d2) |
|
Line 729 int ndl_disjoint(UINT *d1,UINT *d2) |
|
#endif |
#endif |
} |
} |
|
|
int ndl_check_bound2(int index,UINT *d2) |
int ndl_check_bound(UINT *d1,UINT *d2) |
{ |
{ |
UINT u2; |
UINT u2; |
UINT *d1; |
|
int i,j,ind,k; |
int i,j,ind,k; |
|
|
d1 = nd_bound[index]; |
|
ind = 0; |
ind = 0; |
#if USE_UNROLL |
#if USE_UNROLL |
switch ( nd_bpe ) { |
switch ( nd_bpe ) { |
Line 819 int ndl_check_bound2(int index,UINT *d2) |
|
Line 822 int ndl_check_bound2(int index,UINT *d2) |
|
#endif |
#endif |
} |
} |
|
|
|
int ndl_check_bound2(int index,UINT *d2) |
|
{ |
|
return ndl_check_bound(nd_bound[index],d2); |
|
} |
|
|
INLINE int ndl_hash_value(UINT *d) |
INLINE int ndl_hash_value(UINT *d) |
{ |
{ |
int i; |
int i; |
Line 1574 NODE nd_gb(int m,int ishomo,int checkonly) |
|
Line 1582 NODE nd_gb(int m,int ishomo,int checkonly) |
|
NODE r,g,t; |
NODE r,g,t; |
ND_pairs d; |
ND_pairs d; |
ND_pairs l; |
ND_pairs l; |
ND h,nf,s,head; |
ND h,nf,s,head,nf1; |
NDV nfv; |
NDV nfv; |
Q q,num,den; |
Q q,num,den; |
union oNDC dn; |
union oNDC dn; |
|
|
if ( checkonly ) return 0; |
if ( checkonly ) return 0; |
if ( DP_Print ) { printf("+"); fflush(stdout); } |
if ( DP_Print ) { printf("+"); fflush(stdout); } |
nd_removecont(m,nf); |
nd_removecont(m,nf); |
|
if ( nd_nalg ) { |
|
nf1 = nd_pseudo_monic(m,nf); nd_free(nf); |
|
stat = nd_nf(m,nf1,nd_ps,1,0,&nf); |
|
if ( stat ) { |
|
NEXT(l) = d; d = l; |
|
d = nd_reconstruct(0,d); |
|
goto again; |
|
} |
|
} |
nfv = ndtondv(m,nf); nd_free(nf); |
nfv = ndtondv(m,nf); nd_free(nf); |
nh = ndv_newps(m,nfv,0); |
nh = ndv_newps(m,nfv,0); |
d = update_pairs(d,g,nh); |
d = update_pairs(d,g,nh); |
Line 2616 void nd_mul_c(int mod,ND p,int mul) |
|
Line 2633 void nd_mul_c(int mod,ND p,int mul) |
|
int c,c1; |
int c,c1; |
|
|
if ( !p ) return; |
if ( !p ) return; |
|
if ( mul == 1 ) return; |
if ( mod == -1 ) |
if ( mod == -1 ) |
for ( m = BDY(p); m; m = NEXT(m) ) |
for ( m = BDY(p); m; m = NEXT(m) ) |
CM(m) = _mulsf(CM(m),mul); |
CM(m) = _mulsf(CM(m),mul); |
Line 2631 void nd_mul_c_q(ND p,Q mul) |
|
Line 2649 void nd_mul_c_q(ND p,Q mul) |
|
Q c; |
Q c; |
|
|
if ( !p ) return; |
if ( !p ) return; |
|
if ( UNIQ(mul) ) return; |
for ( m = BDY(p); m; m = NEXT(m) ) { |
for ( m = BDY(p); m; m = NEXT(m) ) { |
mulq(CQ(m),mul,&c); CQ(m) = c; |
mulq(CQ(m),mul,&c); CQ(m) = c; |
} |
} |
Line 3709 void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec |
|
Line 3728 void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec |
|
int stat,nvar,max,e; |
int stat,nvar,max,e; |
union oNDC dn; |
union oNDC dn; |
|
|
|
if ( !f ) { |
|
*rp = 0; |
|
return; |
|
} |
pltovl(v,&vv); |
pltovl(v,&vv); |
for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); |
for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); |
|
|
Line 4887 void nd_det(int mod,MAT f,P *rp) |
|
Line 4910 void nd_det(int mod,MAT f,P *rp) |
|
NDV d,s,mij,mjj; |
NDV d,s,mij,mjj; |
ND u; |
ND u; |
NMV nmv; |
NMV nmv; |
|
UINT *bound; |
PGeoBucket bucket; |
PGeoBucket bucket; |
struct order_spec *ord; |
struct order_spec *ord; |
|
|
Line 4904 void nd_det(int mod,MAT f,P *rp) |
|
Line 4928 void nd_det(int mod,MAT f,P *rp) |
|
e = getdeg(tv->v,(P)m[i][j]); |
e = getdeg(tv->v,(P)m[i][j]); |
max = MAX(e,max); |
max = MAX(e,max); |
} |
} |
nd_setup_parameters(nvar,1024); |
nd_setup_parameters(nvar,max); |
dm = (NDV **)almat_pointer(n,n); |
dm = (NDV **)almat_pointer(n,n); |
for ( i = 0, max = 0; i < n; i++ ) |
for ( i = 0, max = 0; i < n; i++ ) |
for ( j = 0; j < n; j++ ) { |
for ( j = 0; j < n; j++ ) { |
Line 4938 void nd_det(int mod,MAT f,P *rp) |
|
Line 4962 void nd_det(int mod,MAT f,P *rp) |
|
} |
} |
sgn = -sgn; |
sgn = -sgn; |
} |
} |
|
bound = nd_det_compute_bound(dm,n,j); |
|
if ( ndl_check_bound(bound,bound) ) |
|
nd_det_reconstruct(dm,n,j,d); |
|
|
for ( i = j+1, mj = dm[j], mjj = mj[j]; i < n; i++ ) { |
for ( i = j+1, mj = dm[j], mjj = mj[j]; i < n; i++ ) { |
if ( DP_Print ) fprintf(stderr," i=%d\n ",i); |
/* if ( DP_Print ) fprintf(stderr," i=%d\n ",i); */ |
mi = dm[i]; mij = mi[j]; |
mi = dm[i]; mij = mi[j]; |
if ( mod ) |
if ( mod ) |
ndv_mul_c(mod,mij,mod-1); |
ndv_mul_c(mod,mij,mod-1); |
else |
else |
ndv_mul_c_q(mij,mone); |
ndv_mul_c_q(mij,mone); |
for ( k = j+1; k < n; k++ ) { |
for ( k = j+1; k < n; k++ ) { |
if ( DP_Print ) fprintf(stderr,"k=%d ",k); |
/* if ( DP_Print ) fprintf(stderr,"k=%d ",k); */ |
bucket = create_pbucket(); |
bucket = create_pbucket(); |
if ( mi[k] ) { |
if ( mi[k] ) { |
nmv = BDY(mjj); len = LEN(mjj); |
nmv = BDY(mjj); len = LEN(mjj); |
Line 4965 void nd_det(int mod,MAT f,P *rp) |
|
Line 4993 void nd_det(int mod,MAT f,P *rp) |
|
u = nd_quo(mod,bucket,d); |
u = nd_quo(mod,bucket,d); |
mi[k] = ndtondv(mod,u); |
mi[k] = ndtondv(mod,u); |
} |
} |
if ( DP_Print ) fprintf(stderr,"\n",k); |
/* if ( DP_Print ) fprintf(stderr,"\n",k); */ |
} |
} |
d = mjj; |
d = mjj; |
} |
} |
Line 5025 ND ndv_mul_nmv_trunc(int mod,NMV m0,NDV p,UINT *d) |
|
Line 5053 ND ndv_mul_nmv_trunc(int mod,NMV m0,NDV p,UINT *d) |
|
return r; |
return r; |
} |
} |
} |
} |
|
} |
|
|
|
void nd_det_reconstruct(NDV **dm,int n,int j,NDV d) |
|
{ |
|
int i,obpe,oadv,h,k,l; |
|
static NM prev_nm_free_list; |
|
EPOS oepos; |
|
|
|
obpe = nd_bpe; |
|
oadv = nmv_adv; |
|
oepos = nd_epos; |
|
if ( obpe < 2 ) nd_bpe = 2; |
|
else if ( obpe < 3 ) nd_bpe = 3; |
|
else if ( obpe < 4 ) nd_bpe = 4; |
|
else if ( obpe < 5 ) nd_bpe = 5; |
|
else if ( obpe < 6 ) nd_bpe = 6; |
|
else if ( obpe < 8 ) nd_bpe = 8; |
|
else if ( obpe < 10 ) nd_bpe = 10; |
|
else if ( obpe < 16 ) nd_bpe = 16; |
|
else if ( obpe < 32 ) nd_bpe = 32; |
|
else error("nd_det_reconstruct : exponent too large"); |
|
|
|
nd_setup_parameters(nd_nvar,0); |
|
prev_nm_free_list = _nm_free_list; |
|
_nm_free_list = 0; |
|
for ( k = j; k < n; k++ ) |
|
for (l = j; l < n; l++ ) |
|
ndv_realloc(dm[k][l],obpe,oadv,oepos); |
|
ndv_realloc(d,obpe,oadv,oepos); |
|
prev_nm_free_list = 0; |
|
#if 0 |
|
GC_gcollect(); |
|
#endif |
|
} |
|
|
|
UINT *nd_det_compute_bound(NDV **dm,int n,int j) |
|
{ |
|
UINT *d0,*d1,*d,*t,*r; |
|
int k,l; |
|
|
|
d0 = (UINT *)ALLOCA(nd_wpd*sizeof(UINT)); |
|
d1 = (UINT *)ALLOCA(nd_wpd*sizeof(UINT)); |
|
for ( k = 0; k < nd_wpd; k++ ) d0[k] = 0; |
|
for ( k = j; k < n; k++ ) |
|
for ( l = j; l < n; l++ ) |
|
if ( dm[k][l] ) { |
|
d = ndv_compute_bound(dm[k][l]); |
|
ndl_lcm(d,d0,d1); |
|
t = d1; d1 = d0; d0 = t; |
|
} |
|
r = (UINT *)ALLOCA(nd_wpd*sizeof(UINT)); |
|
for ( k = 0; k < nd_wpd; k++ ) r[k] = d0[k]; |
|
return r; |
|
} |
|
|
|
DL nd_separate_d(UINT *d,UINT *trans) |
|
{ |
|
int n,ntrans,td,i,e; |
|
DL a; |
|
|
|
n = nd_nvar; ntrans = n-nd_nalg; |
|
ndl_zero(trans); |
|
td = 0; |
|
for ( i = 0; i < ntrans; i++ ) { |
|
e = GET_EXP(d,i); |
|
PUT_EXP(trans,i,e); |
|
td += MUL_WEIGHT(e,i); |
|
} |
|
TD(trans) = td; |
|
if ( nd_blockmask) ndl_weight_mask(trans); |
|
NEWDL(a,nd_nalg); |
|
td = 0; |
|
for ( ; i < n; i++ ) { |
|
e = GET_EXP(d,i); |
|
a->d[i-ntrans] = e; |
|
td += e; |
|
} |
|
a->td = td; |
|
return a; |
|
} |
|
|
|
ND nd_pseudo_monic(int mod,ND p) |
|
{ |
|
UINT *trans,*t; |
|
DL alg; |
|
MP mp0,mp; |
|
NM m,m0,m1; |
|
DL dl; |
|
DP nm; |
|
NDV ndv; |
|
DAlg lc,inv; |
|
ND s,c; |
|
int n,ntrans,i,e,td; |
|
|
|
n = nd_nvar; ntrans = n-nd_nalg; |
|
NEWNM(m0); |
|
NEWNM(m1); |
|
alg = nd_separate_d(HDL(p),DL(m0)); |
|
mp0 = 0; NEXTMP(mp0,mp); mp->c = (P)HCQ(p); mp->dl = alg; |
|
if ( !mp->dl->td ) |
|
return p; |
|
for ( m = NEXT(BDY(p)); m; m = NEXT(m) ) { |
|
alg = nd_separate_d(DL(m),DL(m1)); |
|
if ( !ndl_equal(DL(m0),DL(m1)) ) |
|
break; |
|
NEXTMP(mp0,mp); mp->c = (P)CQ(m); mp->dl = alg; |
|
} |
|
NEXT(mp) = 0; |
|
MKDP(nd_nalg,mp0,nm); |
|
MKDAlg(nm,ONE,lc); |
|
invdalg(lc,&inv); |
|
ndv = ndtondv(0,p); |
|
for ( s = 0, mp = BDY(inv->nm); mp; mp = NEXT(mp) ) { |
|
CQ(m0) = (Q)mp->c; |
|
dl = mp->dl; |
|
for ( td = 0, i = ntrans; i < n; i++ ) { |
|
e = dl->d[i-ntrans]; |
|
ndl_zero(DL(m0)); |
|
PUT_EXP(DL(m0),i,e); |
|
td += MUL_WEIGHT(e,i); |
|
} |
|
TD(DL(m0)) = td; |
|
if ( nd_blockmask) ndl_weight_mask(trans); |
|
s = nd_add(0,s,ndv_mul_nm(0,m0,ndv)); |
|
} |
|
ndv_free(ndv); |
|
return s; |
} |
} |