version 1.33, 2005/11/16 23:42:53 |
version 1.37, 2005/11/26 01:28:12 |
|
|
* 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/engine/dist.c,v 1.32 2004/06/15 16:14:50 ohara Exp $ |
* $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.36 2005/11/25 07:18:32 noro Exp $ |
*/ |
*/ |
#include "ca.h" |
#include "ca.h" |
|
|
Line 1784 int ni_next(int *a,int n) |
|
Line 1784 int ni_next(int *a,int n) |
|
int i,j,k,kj; |
int i,j,k,kj; |
|
|
/* find the first nonzero a[j] */ |
/* find the first nonzero a[j] */ |
for ( j = 0; a[j] == 0; j++ ); |
for ( j = 0; j < n && a[j] == 0; j++ ); |
/* find the first zero a[k] after a[j] */ |
/* find the first zero a[k] after a[j] */ |
for ( k = j; k < n && a[k] == 1; k++ ); |
for ( k = j; k < n && a[k] == 1; k++ ); |
if ( k == n ) return 0; |
if ( k == n ) return 0; |
Line 1837 NBM mul_nbm(NBM a,NBM b) |
|
Line 1837 NBM mul_nbm(NBM a,NBM b) |
|
return m; |
return m; |
} |
} |
|
|
|
NBP nbmtonbp(NBM m) |
|
{ |
|
NODE n; |
|
NBP u; |
|
|
|
MKNODE(n,m,0); |
|
MKNBP(u,n); |
|
return u; |
|
} |
|
|
|
/* a=c*x*rest -> a0= x*rest, ah=x, ar=rest */ |
|
|
|
Q separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) |
|
{ |
|
int i,d1; |
|
NBM t; |
|
|
|
if ( !a->d ) error("separate_nbm : invalid argument"); |
|
|
|
NEWNBM(t); t->d = a->d; t->b = a->b; t->c = ONE; |
|
*a0 = nbmtonbp(t); |
|
|
|
NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = ONE; |
|
if ( NBM_GET(a->b,0) ) NBM_SET(t->b,0); |
|
else NBM_CLR(t->b,0); |
|
*ah = nbmtonbp(t); |
|
|
|
d1 = a->d-1; |
|
NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = ONE; |
|
for ( i = 0; i < d1; i++ ) { |
|
if ( NBM_GET(a->b,i+1) ) NBM_SET(t->b,i); |
|
else NBM_CLR(t->b,i); |
|
} |
|
*ar = nbmtonbp(t); |
|
return a->c; |
|
} |
|
|
|
NBP make_xky(int k) |
|
{ |
|
int k1,i; |
|
NBM t; |
|
|
|
NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = ONE; |
|
k1 = k-1; |
|
for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i); |
|
NBM_CLR(t->b,i); |
|
return nbmtonbp(t); |
|
} |
|
|
|
/* a=c*x^(k-1)*y*rest -> a0= x^(k-1)*y*rest, ah=x^(k-1)*y, ar=rest */ |
|
|
|
Q separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar) |
|
{ |
|
int i,d1,k,k1; |
|
NBM t; |
|
|
|
if ( !a->d ) |
|
error("separate_nbm : invalid argument"); |
|
for ( i = 0; i < a->d && NBM_GET(a->b,i); i++ ); |
|
if ( i == a->d ) |
|
error("separate_nbm : invalid argument"); |
|
k1 = i; |
|
k = i+1; |
|
|
|
NEWNBM(t); t->d = a->d; t->b = a->b; t->c = ONE; |
|
*a0 = nbmtonbp(t); |
|
|
|
NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = ONE; |
|
for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i); |
|
NBM_CLR(t->b,i); |
|
*ah = nbmtonbp(t); |
|
|
|
d1 = a->d-k; |
|
NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = ONE; |
|
for ( i = 0; i < d1; i++ ) { |
|
if ( NBM_GET(a->b,i+k) ) NBM_SET(t->b,i); |
|
else NBM_CLR(t->b,i); |
|
} |
|
*ar = nbmtonbp(t); |
|
|
|
return a->c; |
|
} |
|
#if 0 |
NBP shuffle_mul_nbm(NBM a,NBM b) |
NBP shuffle_mul_nbm(NBM a,NBM b) |
{ |
{ |
int ad,bd,d,i,ai,bi,bit,s; |
int ad,bd,d,i,ai,bi,bit,s; |
Line 1892 NBP shuffle_mul_nbm(NBM a,NBM b) |
|
Line 1975 NBP shuffle_mul_nbm(NBM a,NBM b) |
|
MKNBP(u,r); |
MKNBP(u,r); |
return u; |
return u; |
} |
} |
|
#else |
|
void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp); |
|
void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp); |
|
|
|
NBP shuffle_mul_nbm(NBM a,NBM b) |
|
{ |
|
NBP u,a0,ah,ar,b0,bh,br,a1,b1,t; |
|
Q ac,bc,c; |
|
|
|
if ( !a->d || !b->d ) |
|
u = nbmtonbp(mul_nbm(a,b)); |
|
else { |
|
ac = separate_nbm(a,&a0,&ah,&ar); |
|
bc = separate_nbm(b,&b0,&bh,&br); |
|
mulq(ac,bc,&c); |
|
shuffle_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1); |
|
shuffle_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1); |
|
addnbp(CO,a1,b1,&t); mulnbp(CO,(NBP)c,t,&u); |
|
} |
|
return u; |
|
} |
|
#endif |
|
|
|
int nbmtoxky(NBM a,int *b) |
|
{ |
|
int d,i,j,k; |
|
int *p; |
|
|
|
d = a->d; p = a->b; |
|
for ( i = j = 0, k = 1; i < d; i++ ) { |
|
if ( !NBM_GET(p,i) ) { |
|
b[j++] = k; |
|
k = 1; |
|
} else k++; |
|
} |
|
return j; |
|
} |
|
|
|
#if 0 |
|
NBP harmonic_mul_nbm(NBM a,NBM b) |
|
{ |
|
int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1; |
|
int i,j,k,ia,ib,s; |
|
int *wa,*wb,*w,*wab,*wa1,*wmb; |
|
Q c,c1; |
|
NBM wm,tm; |
|
NODE r,t1,t,p; |
|
NBP u; |
|
|
|
da = a->d; db = b->d; d = da+db; |
|
wa = (int *)ALLOCA(da*sizeof(int)); |
|
wb = (int *)ALLOCA(db*sizeof(int)); |
|
la = nbmtoxky(a,wa); |
|
lb = nbmtoxky(b,wb); |
|
mulq(a->c,b->c,&c); |
|
/* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */ |
|
/* lmax : total length */ |
|
lmax = la+lb; |
|
lmin = la>lb?la:lb; |
|
w = (int *)ALLOCA(lmax*sizeof(int)); |
|
/* position of a+b */ |
|
wab = (int *)ALLOCA(lmax*sizeof(int)); |
|
/* position of a */ |
|
wa1 = (int *)ALLOCA(lmax*sizeof(int)); |
|
NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b; |
|
for ( l = lmin, r = 0; l <= lmax; l++ ) { |
|
lab = lmax - l; |
|
la1 = la - lab; |
|
lb1 = lb - lab; |
|
lab1 = l-lab; |
|
/* partion l into three parts: a, b, a+b */ |
|
/* initialize wab */ |
|
for ( i = 0; i < lab; i++ ) wab[i] = 1; |
|
for ( ; i < l; i++ ) wab[i] = 0; |
|
do { |
|
/* initialize wa1 */ |
|
for ( i = 0; i < la1; i++ ) wa1[i] = 1; |
|
for ( ; i < lab1; i++ ) wa1[i] = 0; |
|
do { |
|
ia = 0; ib = 0; |
|
for ( i = j = 0; i < l; i++ ) |
|
if ( wab[i] ) w[i] = wa[ia++]+wb[ib++]; |
|
else if ( wa1[j++] ) w[i] = wa[ia++]; |
|
else w[i] = wb[ib++]; |
|
for ( i = j = 0; i < l; i++ ) { |
|
for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j); |
|
NBM_CLR(wmb,j); j++; |
|
} |
|
wm->d = j; wm->c = c; |
|
for ( p = 0, t = r; t; p = t, t = NEXT(t) ) { |
|
tm = (NBM)BDY(t); |
|
s = comp_nbm(tm,wm); |
|
if ( s < 0 ) { |
|
/* insert */ |
|
MKNODE(t1,wm,t); |
|
if ( !p ) r = t1; |
|
else NEXT(p) = t1; |
|
NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b; |
|
break; |
|
} else if ( s == 0 ) { |
|
/* add coefs */ |
|
addq(tm->c,c,&c1); |
|
if ( c1 ) tm->c = c1; |
|
else NEXT(p) = NEXT(t); |
|
break; |
|
} |
|
} |
|
if ( !t ) { |
|
/* append */ |
|
MKNODE(t1,wm,t); |
|
if ( !p ) r = t1; |
|
else NEXT(p) = t1; |
|
NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b; |
|
} |
|
} while ( ni_next(wa1,lab1) ); |
|
} while ( ni_next(wab,l) ); |
|
} |
|
MKNBP(u,r); |
|
return u; |
|
} |
|
#else |
|
NBP harmonic_mul_nbm(NBM a,NBM b) |
|
{ |
|
NBP u,a0,ah,ar,b0,bh,br,a1,b1,t,s,abk,ab1; |
|
Q ac,bc,c; |
|
|
|
if ( !a->d || !b->d ) |
|
u = nbmtonbp(mul_nbm(a,b)); |
|
else { |
|
mulq(a->c,b->c,&c); |
|
ac = separate_xky_nbm(a,&a0,&ah,&ar); |
|
bc = separate_xky_nbm(b,&b0,&bh,&br); |
|
mulq(ac,bc,&c); |
|
harmonic_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1); |
|
harmonic_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1); |
|
abk = make_xky(((NBM)BDY(BDY(ah)))->d+((NBM)BDY(BDY(bh)))->d); |
|
harmonic_mulnbp(CO,ar,br,&t); mulnbp(CO,abk,t,&ab1); |
|
addnbp(CO,a1,b1,&t); addnbp(CO,t,ab1,&s); mulnbp(CO,(NBP)c,s,&u); |
|
} |
|
return u; |
|
|
|
} |
|
#endif |
|
|
void addnbp(VL vl,NBP p1,NBP p2, NBP *rp) |
void addnbp(VL vl,NBP p1,NBP p2, NBP *rp) |
{ |
{ |
NODE b1,b2,br,br0; |
NODE b1,b2,br,br0; |
Line 1922 void addnbp(VL vl,NBP p1,NBP p2, NBP *rp) |
|
Line 2148 void addnbp(VL vl,NBP p1,NBP p2, NBP *rp) |
|
NEXTNODE(br0,br); BDY(br) = BDY(b2); |
NEXTNODE(br0,br); BDY(br) = BDY(b2); |
b2 = NEXT(b2); break; |
b2 = NEXT(b2); break; |
} |
} |
if ( !br0 ) |
} |
if ( b1 ) |
if ( !br0 ) |
br0 = b1; |
if ( b1 ) |
else if ( b2 ) |
br0 = b1; |
br0 = b2; |
|
else { |
|
*rp = 0; |
|
return; |
|
} |
|
else if ( b1 ) |
|
NEXT(br) = b1; |
|
else if ( b2 ) |
else if ( b2 ) |
|
br0 = b2; |
|
else { |
|
*rp = 0; |
|
return; |
|
} |
|
else if ( b1 ) |
|
NEXT(br) = b1; |
|
else if ( b2 ) |
NEXT(br) = b2; |
NEXT(br) = b2; |
else |
else |
NEXT(br) = 0; |
NEXT(br) = 0; |
MKNBP(*rp,br0); |
MKNBP(*rp,br0); |
} |
|
} |
} |
} |
} |
|
|
Line 1959 void chsgnnbp(NBP p,NBP *rp) |
|
Line 2185 void chsgnnbp(NBP p,NBP *rp) |
|
NEXTNODE(r0,r); |
NEXTNODE(r0,r); |
m = (NBM)BDY(b); |
m = (NBM)BDY(b); |
NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnq(m->c,&m1->c); |
NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnq(m->c,&m1->c); |
BDY(r) = m; |
BDY(r) = m1; |
} |
} |
if ( r0 ) NEXT(r) = 0; |
if ( r0 ) NEXT(r) = 0; |
MKNBP(*rp,r0); |
MKNBP(*rp,r0); |
Line 1970 void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp); |
|
Line 2196 void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp); |
|
|
|
void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp) |
void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp) |
{ |
{ |
NODE b; |
NODE b,n; |
NBP r,t,s; |
NBP r,t,s; |
|
NBM m; |
|
|
if ( !p1 || !p2 ) *rp = 0; |
if ( !p1 || !p2 ) { |
else if ( length(BDY(p1)) < length(BDY(p2)) ) { |
*rp = 0; return; |
|
} |
|
if ( OID(p1) != O_NBP ) { |
|
if ( !NUM(p1) || !RATN(p1) ) |
|
error("mulnbp : invalid argument"); |
|
NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p1; |
|
MKNODE(n,m,0); MKNBP(p1,n); |
|
} |
|
if ( OID(p2) != O_NBP ) { |
|
if ( !NUM(p2) || !RATN(p2) ) |
|
error("mulnbp : invalid argument"); |
|
NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p2; |
|
MKNODE(n,m,0); MKNBP(p2,n); |
|
} |
|
if ( length(BDY(p1)) < length(BDY(p2)) ) { |
for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) { |
for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) { |
mulnbmnbp(vl,(NBM)BDY(b),p2,&t); |
mulnbmnbp(vl,(NBM)BDY(b),p2,&t); |
addnbp(vl,r,t,&s); r = s; |
addnbp(vl,r,t,&s); r = s; |
Line 2046 void pwrnbp(VL vl,NBP a,Q q,NBP *c) |
|
Line 2287 void pwrnbp(VL vl,NBP a,Q q,NBP *c) |
|
} |
} |
} |
} |
|
|
void shuffle_mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp); |
|
void shuffle_mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp); |
|
void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp); |
|
|
|
void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp) |
void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp) |
{ |
{ |
NODE b; |
NODE b1,b2,n; |
NBP r,t,s; |
NBP r,t,s; |
|
NBM m; |
|
|
if ( !p1 || !p2 ) *rp = 0; |
if ( !p1 || !p2 ) { |
else if ( length(BDY(p1)) < length(BDY(p2)) ) { |
*rp = 0; return; |
for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) { |
|
shuffle_mulnbmnbp(vl,(NBM)BDY(b),p2,&t); |
|
addnbp(vl,r,t,&s); r = s; |
|
} |
|
*rp = r; |
|
} else { |
|
for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) { |
|
shuffle_mulnbpnbm(vl,p1,(NBM)BDY(b),&t); |
|
addnbp(vl,r,t,&s); r = s; |
|
} |
|
*rp = r; |
|
} |
} |
} |
if ( OID(p1) != O_NBP ) { |
|
if ( !NUM(p1) || !RATN(p1) ) |
void shuffle_mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp) |
error("shuffle_mulnbp : invalid argument"); |
{ |
NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p1; |
NODE b; |
MKNODE(n,m,0); MKNBP(p1,n); |
NBP t,s,r; |
} |
|
if ( OID(p2) != O_NBP ) { |
if ( !p ) *rp = 0; |
if ( !NUM(p2) || !RATN(p2) ) |
else { |
error("shuffle_mulnbp : invalid argument"); |
r = 0; |
NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p2; |
for ( b = BDY(p); b; b = NEXT(b) ) { |
MKNODE(n,m,0); MKNBP(p2,n); |
t = shuffle_mul_nbm(m,(NBM)BDY(b)); |
} |
|
for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) ) |
|
for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) { |
|
t = shuffle_mul_nbm(m,(NBM)BDY(b2)); |
addnbp(vl,r,t,&s); r = s; |
addnbp(vl,r,t,&s); r = s; |
} |
} |
*rp = r; |
*rp = r; |
} |
|
} |
} |
|
|
void shuffle_mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp) |
void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp) |
{ |
{ |
NODE b; |
NODE b1,b2,n; |
NBP t,s,r; |
NBP r,t,s; |
|
NBM m; |
|
|
if ( !p ) *rp = 0; |
if ( !p1 || !p2 ) { |
else { |
*rp = 0; return; |
r = 0; |
} |
for ( b = BDY(p); b; b = NEXT(b) ) { |
if ( OID(p1) != O_NBP ) { |
t = shuffle_mul_nbm((NBM)BDY(b),m); |
if ( !NUM(p1) || !RATN(p1) ) |
|
error("harmonic_mulnbp : invalid argument"); |
|
NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p1; |
|
MKNODE(n,m,0); MKNBP(p1,n); |
|
} |
|
if ( OID(p2) != O_NBP ) { |
|
if ( !NUM(p2) || !RATN(p2) ) |
|
error("harmonic_mulnbp : invalid argument"); |
|
NEWNBM(m); m->d = 0; m->b = 0; m->c = (Q)p2; |
|
MKNODE(n,m,0); MKNBP(p2,n); |
|
} |
|
for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) ) |
|
for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) { |
|
t = harmonic_mul_nbm(m,(NBM)BDY(b2)); |
addnbp(vl,r,t,&s); r = s; |
addnbp(vl,r,t,&s); r = s; |
} |
} |
*rp = r; |
*rp = r; |
} |
|
} |
} |