OpenXM_contrib2/asir2000/engine/nd.c - diff

Return to nd.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / engine

Diff for /OpenXM_contrib2/asir2000/engine/nd.c between version 1.53 and 1.56

version 1.53, 2003/08/29 07:37:30

version 1.56, 2003/09/04 08:35:09

Line 1

/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.52 2003/08/27 02:21:16 noro Exp $ */

/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.55 2003/09/03 07:33:35 noro Exp $ */

#include "ca.h"

#include "inline.h"

Line 120 static int nm_adv;

static int nmv_adv;

static int nd_dcomp;

extern int Top,Reverse,dp_nelim;

extern int Top,Reverse,dp_nelim,do_weyl;

/* fundamental macros */

#define TD(d) (d[0])

Line 185 if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT

/* macro for increasing pointer to NMV */

#define NMV_ADV(m) (m = (NMV)(((char *)m)+nmv_adv))

#define NMV_PREV(m) (m = (NMV)(((char *)m)-nmv_adv))

/* external functions */

void GC_gcollect();

Line 233 INLINE int ndl_block_compare(unsigned int *d1,unsigned

Line 234 INLINE int ndl_block_compare(unsigned int *d1,unsigned

INLINE int ndl_equal(unsigned int *d1,unsigned int *d2);

INLINE void ndl_copy(unsigned int *d1,unsigned int *d2);

INLINE void ndl_add(unsigned int *d1,unsigned int *d2,unsigned int *d);

INLINE void ndl_addto(unsigned int *d1,unsigned int *d2);

INLINE void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d);

INLINE int ndl_hash_value(unsigned int *d);

Line 288 ND nd_add_q(ND p1,ND p2);

Line 290 ND nd_add_q(ND p1,ND p2);

INLINE int nd_length(ND p);

/* NDV functions */

ND weyl_ndv_mul_nm(int mod,NM m0,NDV p);

void weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *tab,int tlen);

void ndv_mul_c(int mod,NDV p,int mul);

void ndv_mul_c_q(NDV p,Q mul);

void ndv_realloc(NDV p,int obpe,int oadv,EPOS oepos);

ND ndv_mul_nm(int mod,NDV p,NM m0);

ND ndv_mul_nm(int mod,NM m0,NDV p);

void ndv_dehomogenize(NDV p,struct order_spec *spec);

void ndv_removecont(int mod,NDV p);

void ndv_print(NDV p);

Line 661 INLINE void ndl_add(unsigned int *d1,unsigned int *d2,

Line 665 INLINE void ndl_add(unsigned int *d1,unsigned int *d2,

#endif

}

/* d1 += d2 */

INLINE void ndl_addto(unsigned int *d1,unsigned int *d2)

{

int i;

#if 1

switch ( nd_wpd ) {

case 2:

TD(d1) += TD(d2);

d1[1] += d2[1];

break;

case 3:

TD(d1) += TD(d2);

d1[1] += d2[1];

d1[2] += d2[2];

break;

default:

for ( i = 0; i < nd_wpd; i++ ) d1[i] += d2[i];

break;

}

#else

for ( i = 0; i < nd_wpd; i++ ) d1[i] += d2[i];

#endif

}

INLINE void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d)

{

int i;

Line 1077 int nd_nf(int mod,ND g,NDV *ps,int full,ND *rp)

Line 1106 int nd_nf(int mod,ND g,NDV *ps,int full,ND *rp)

chsgnq(cg,&CQ(mul));

nd_mul_c_q(d,cred); nd_mul_c_q(g,cred);

}

g = nd_add(mod,g,ndv_mul_nm(mod,p,mul));

g = nd_add(mod,g,ndv_mul_nm(mod,mul,p));

sugar = MAX(sugar,SG(p)+TD(DL(mul)));

if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) {

nd_removecont2(d,g);

Line 1160 int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp

Line 1189 int nd_nf_pbucket(int mod,ND g,NDV *ps,int full,ND *rp

g = bucket->body[hindex];

gmag = (double)p_mag((P)HCQ(g));

}

red = ndv_mul_nm(mod,p,mul);

red = ndv_mul_nm(mod,mul,p);

bucket->body[hindex] = nd_remove_head(g);

red = nd_remove_head(red);

add_pbucket(mod,bucket,red);

Line 1249 int nd_nf_direct(int mod,ND g,BaseSet base,int full,ND

Line 1278 int nd_nf_direct(int mod,ND g,BaseSet base,int full,ND

chsgnq(cg,&CQ(mul));

nd_mul_c_q(d,cred); nd_mul_c_q(g,cred);

}

g = nd_add(mod,g,ndv_mul_nm(mod,p,mul));

g = nd_add(mod,g,ndv_mul_nm(mod,mul,p));

sugar = MAX(sugar,SG(p)+TD(DL(mul)));

if ( !mod && hmag && g && ((double)(p_mag((P)HCQ(g))) > hmag) ) {

nd_removecont2(d,g);

Line 1605 NODE nd_reduceall(int m,NODE f)

Line 1634 NODE nd_reduceall(int m,NODE f)

for ( n = 0, t = f; t; t = NEXT(t), n++ );

ps = (NDV *)ALLOCA(n*sizeof(NDV));

bound = (unsigned int **)ALLOCA(n*sizeof(unsigned int *));

for ( i = 0, t = f; i < n; i++, t = NEXT(t) ) {

for ( i = 0, t = f; i < n; i++, t = NEXT(t) ) ps[i] = (NDV)BDY(t);

ps[i] = (NDV)BDY(t);

bound[i] = ndv_compute_bound(ps[i]);

}

qsort(ps,n,sizeof(NDV),(int (*)(const void *,const void *))ndv_compare);

for ( i = 0; i < n; i++ ) bound[i] = ndv_compute_bound(ps[i]);

base.ps = (NDV *)ALLOCA((n-1)*sizeof(NDV));

base.bound = (unsigned int **)ALLOCA((n-1)*sizeof(unsigned int *));

base.len = n-1;

Line 1661 ND_pairs update_pairs( ND_pairs d, NODE /* of index */

Line 1688 ND_pairs update_pairs( ND_pairs d, NODE /* of index */

d1 = nd_newpairs(g,t);

d1 = crit_M(d1);

d1 = crit_F(d1);

prev = 0; cur = head = d1;

if ( do_weyl )

while ( cur ) {

head = d1;

if ( crit_2( cur->i1,cur->i2 ) ) {

else {

remove = cur;

prev = 0; cur = head = d1;

if ( !prev ) head = cur = NEXT(cur);

while ( cur ) {

else cur = NEXT(prev) = NEXT(cur);

if ( crit_2( cur->i1,cur->i2 ) ) {

FREENDP(remove);

remove = cur;

} else {

if ( !prev ) head = cur = NEXT(cur);

prev = cur; cur = NEXT(cur);

else cur = NEXT(prev) = NEXT(cur);

FREENDP(remove);

} else {

prev = cur; cur = NEXT(cur);

}

if ( !d )

Line 1733 ND_pairs crit_B( ND_pairs d, int s )

Line 1764 ND_pairs crit_B( ND_pairs d, int s )

return head;

}

/* XXX : check is necessary */

ND_pairs crit_M( ND_pairs d1 )

{

ND_pairs e,d2,d3,dd,p;

Line 2133 void dltondl(int n,DL dl,unsigned int *r)

Line 2162 void dltondl(int n,DL dl,unsigned int *r)

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

r[j+1] = ndl_weight_mask(r,j);

} else {

if ( nd_isrlex )

for ( i = 0; i < n; i++ ) PUT_EXP(r,i,d[i]);

for ( i = 0; i < n; i++ ) PUT_EXP(r,n-1-i,d[i]);

else

for ( i = 0; i < n; i++ ) PUT_EXP(r,i,d[i]);

TD(r) = ndl_weight(r);

}

Line 2166 DL ndltodl(int n,unsigned int *ndl)

Line 2192 DL ndltodl(int n,unsigned int *ndl)

s += ord_l;

}

} else {

if ( nd_isrlex )

for ( i = 0; i < n; i++ ) d[i] = GET_EXP(ndl,i);

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

d[i] = GET_EXP(ndl,n-1-i);

else

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

d[i] = GET_EXP(ndl,i);

}

return dl;

}

Line 2244 void ndl_print(unsigned int *dl)

Line 2265 void ndl_print(unsigned int *dl)

printf(s==n-1?"%d":"%d,",GET_EXP(dl,s));

}

} else {

if ( nd_isrlex )

for ( i = 0; i < n; i++ ) printf(i==n-1?"%d":"%d,",GET_EXP(dl,i));

for ( i = 0; i < n; i++ ) printf(i==n-1?"%d":"%d,",GET_EXP(dl,n-1-i));

else

for ( i = 0; i < n; i++ ) printf(i==n-1?"%d":"%d,",GET_EXP(dl,i));

}

printf(">>");

}

Line 2550 void nd_setup_parameters() {

Line 2568 void nd_setup_parameters() {

nd_exporigin = nd_get_exporigin(nd_ord);

nd_wpd = nd_exporigin+elen;

nd_epos = (EPOS)MALLOC_ATOMIC(nd_nvar*sizeof(struct oEPOS));

for ( i = 0; i < nd_nvar; i++ ) {

if ( nd_isrlex ) {

nd_epos[i].i = nd_exporigin + i/nd_epw;

for ( i = 0; i < nd_nvar; i++ ) {

nd_epos[i].s = (nd_epw-(i%nd_epw)-1)*nd_bpe;

nd_epos[i].i = nd_exporigin + (nd_nvar-1-i)/nd_epw;

nd_epos[i].s = (nd_epw-((nd_nvar-1-i)%nd_epw)-1)*nd_bpe;

}

} else {

for ( i = 0; i < nd_nvar; i++ ) {

nd_epos[i].i = nd_exporigin + i/nd_epw;

nd_epos[i].s = (nd_epw-(i%nd_epw)-1)*nd_bpe;

}

if ( nd_bpe < 32 ) {

nd_mask0 = (1<<nd_bpe)-1;

Line 2701 void ndl_reconstruct(int obpe,EPOS oepos,unsigned int

Line 2726 void ndl_reconstruct(int obpe,EPOS oepos,unsigned int

s += ord_l;

}

} else {

if ( nd_isrlex )

for ( i = 0; i < n; i++ ) {

ei = GET_EXP_OLD(d,i);

ei = GET_EXP_OLD(d,n-1-i);

PUT_EXP(r,i,ei);

PUT_EXP(r,n-1-i,ei);

}

else

for ( i = 0; i < n; i++ ) {

ei = GET_EXP_OLD(d,i);

PUT_EXP(r,i,ei);

}

Line 2752 int nd_sp(int mod,int trace,ND_pairs p,ND *rp)

Line 2771 int nd_sp(int mod,int trace,ND_pairs p,ND *rp)

NEWNM(m);

CQ(m) = HCQ(p2);

ndl_sub(lcm,HDL(p1),DL(m));

if ( ndl_check_bound2(p->i1,DL(m)) ) return 0;

if ( ndl_check_bound2(p->i1,DL(m)) )

t1 = ndv_mul_nm(mod,p1,m);

return 0;

t1 = ndv_mul_nm(mod,m,p1);

if ( mod ) CM(m) = mod-HCM(p1);

else chsgnq(HCQ(p1),&CQ(m));

ndl_sub(lcm,HDL(p2),DL(m));

Line 2761 int nd_sp(int mod,int trace,ND_pairs p,ND *rp)

Line 2781 int nd_sp(int mod,int trace,ND_pairs p,ND *rp)

nd_free(t1);

return 0;

}

t2 = ndv_mul_nm(mod,p2,m);

t2 = ndv_mul_nm(mod,m,p2);

*rp = nd_add(mod,t1,t2);

FREENM(m);

return 1;

Line 2792 void ndv_mul_c_q(NDV p,Q mul)

Line 2812 void ndv_mul_c_q(NDV p,Q mul)

}

ND ndv_mul_nm(int mod,NDV p,NM m0)

ND weyl_ndv_mul_nm(int mod,NM m0,NDV p) {

int n2,i,j,l,n,tlen;

unsigned int *d0;

NM *tab,*psum;

ND s,r;

NM t;

NMV m1;

if ( !p ) return 0;

n = NV(p); n2 = n>>1;

d0 = DL(m0);

l = LEN(p);

for ( i = 0, tlen = 1; i < n2; i++ ) tlen *= (GET_EXP(d0,n2+i)+1);

tab = (NM *)ALLOCA(tlen*sizeof(NM));

psum = (NM *)ALLOCA(tlen*sizeof(NM));

for ( i = 0; i < tlen; i++ ) psum[i] = 0;

m1 = (NMV)(((char *)BDY(p))+nmv_adv*(l-1));

for ( i = l-1; i >= 0; i--, NMV_PREV(m1) ) {

/* m0(NM) * m1(NMV) => tab(NM) */

weyl_mul_nm_nmv(n,mod,m0,m1,tab,tlen);

for ( j = 0; j < tlen; j++ ) {

if ( tab[j] ) {

NEXT(tab[j]) = psum[j]; psum[j] = tab[j];

}

for ( i = tlen-1, r = 0; i >= 0; i-- )

if ( psum[i] ) {

for ( j = 0, t = psum[i]; t; t = NEXT(t), j++ );

MKND(n,psum[i],j,s);

r = nd_add(mod,r,s);

}

if ( r ) SG(r) = SG(p)+TD(d0);

return r;

}

/* product of monomials */

/* XXX block order is not handled correctly */

void weyl_mul_nm_nmv(int n,int mod,NM m0,NMV m1,NM *tab,int tlen)

{

int i,n2,j,s,curlen,homo,h,a,b,k,l,u,min;

unsigned int *d0,*d1,*d,*dt,*ctab;

Q *ctab_q;

Q q,q1;

unsigned int c0,c1,c;

NM *p;

NM m,t;

for ( i = 0; i < tlen; i++ ) tab[i] = 0;

if ( !m0 || !m1 ) return;

d0 = DL(m0); d1 = DL(m1); n2 = n>>1;

NEWNM(m); d = DL(m);

if ( mod ) {

c0 = CM(m0); c1 = CM(m1); DMAR(c0,c1,0,mod,c); CM(m) = c;

} else

mulq(CQ(m0),CQ(m1),&CQ(m));

for ( i = 0; i < nd_wpd; i++ ) d[i] = 0;

homo = n&1 ? 1 : 0;

if ( homo ) {

/* offset of h-degree */

h = GET_EXP(d0,n-1)+GET_EXP(d1,n-1);

PUT_EXP(DL(m),n-1,h);

TD(DL(m)) = h;

/* XXX other weights */

}

tab[0] = m;

curlen = 1;

NEWNM(m); d = DL(m);

for ( i = 0; i < n2; i++ ) {

a = GET_EXP(d0,i); b = GET_EXP(d1,n2+i);

k = GET_EXP(d0,n2+i); l = GET_EXP(d1,i);

/* xi^a*(Di^k*xi^l)*Di^b */

a += l; b += k;

s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);

if ( !k || !l ) {

for ( j = 0; j < curlen; j++ )

if ( t = tab[j] ) {

dt = DL(t);

PUT_EXP(dt,i,a); PUT_EXP(dt,n2+i,b); TD(dt) += s;

/* XXX other weights */

}

curlen *= k+1;

continue;

}

min = MIN(k,l);

if ( mod ) {

ctab = (unsigned int *)ALLOCA((min+1)*sizeof(unsigned int));

mkwcm(k,l,mod,ctab);

} else {

ctab_q = (Q *)ALLOCA((min+1)*sizeof(Q));

mkwc(k,l,ctab_q);

}

p = tab+curlen;

for ( j = 1; j <= min; j++ ) {

for ( u = 0; u < nd_wpd; u++ ) d[u] = 0;

PUT_EXP(d,i,a-j); PUT_EXP(d,n2+i,b-j);

h = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i);

if ( homo ) {

TD(d) = s;

PUT_EXP(d,n-1,s-h);

} else TD(d) = h;

/* XXX other weights */

if ( mod ) c = ctab[j];

else q = ctab_q[j];

for ( u = 0; u < curlen; u++, p++ ) {

if ( tab[u] ) {

NEWNM(t);

ndl_add(DL(tab[u]),d,DL(t));

if ( mod ) {

c0 = CM(tab[u]); DMAR(c0,c,0,mod,c1); CM(t) = c1;

} else

mulq(CQ(tab[u]),q,&CQ(t));

*p = t;

}

/* destructive for j = 0 */

for ( u = 0; u < nd_wpd; u++ ) d[u] = 0;

PUT_EXP(d,i,a); PUT_EXP(d,n2+i,b);

h = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);

if ( homo ) {

TD(d) = s;

PUT_EXP(d,n-1,s-h);

} else TD(d) = h;

/* XXX other weights */

if ( mod ) c = ctab[0];

else q = ctab_q[0];

p = tab;

for ( u = 0; u < curlen; u++, p++ ) {

if ( tab[u] ) {

ndl_addto(DL(tab[u]),d);

if ( mod ) {

c0 = CM(tab[u]); DMAR(c0,c,0,mod,c1); CM(tab[u]) = c1;

} else {

mulq(CQ(tab[u]),q,&q1); CQ(tab[u]) = q1;

}

curlen *= k+1;

}

FREENM(m);

}

ND ndv_mul_nm(int mod,NM m0,NDV p)

{

NM mr,mr0;

NMV m;

unsigned int *d,*dt,*dm;

Line 2802 ND ndv_mul_nm(int mod,NDV p,NM m0)

Line 2966 ND ndv_mul_nm(int mod,NDV p,NM m0)

ND r;

if ( !p ) return 0;

else if ( do_weyl )

return weyl_ndv_mul_nm(mod,m0,p);

else {

n = NV(p); m = BDY(p);

d = DL(m0);

Line 2838 void ndv_realloc(NDV p,int obpe,int oadv,EPOS oepos)

Line 3004 void ndv_realloc(NDV p,int obpe,int oadv,EPOS oepos)

int len,i,k;

#define NMV_OPREV(m) (m = (NMV)(((char *)m)-oadv))

#define NMV_PREV(m) (m = (NMV)(((char *)m)-nmv_adv))

if ( p ) {

m = BDY(p); len = LEN(p);

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