===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/dist.c,v
retrieving revision 1.34
retrieving revision 1.37
diff -u -p -r1.34 -r1.37
--- OpenXM_contrib2/asir2000/engine/dist.c	2005/11/24 08:16:03	1.34
+++ OpenXM_contrib2/asir2000/engine/dist.c	2005/11/26 01:28:12	1.37
@@ -45,7 +45,7 @@
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *
- * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.33 2005/11/16 23:42:53 noro Exp $ 
+ * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.36 2005/11/25 07:18:32 noro Exp $ 
 */
 #include "ca.h"
 
@@ -1784,7 +1784,7 @@ int ni_next(int *a,int n)
 	int i,j,k,kj;
 
 	/* 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] */
 	for ( k = j; k < n && a[k] == 1; k++ );
 	if ( k == n ) return 0;
@@ -1837,6 +1837,89 @@ NBM mul_nbm(NBM a,NBM b)
 	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)
 {
 	int ad,bd,d,i,ai,bi,bit,s;
@@ -1892,7 +1975,29 @@ NBP shuffle_mul_nbm(NBM a,NBM b)
 	MKNBP(u,r);
 	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;
@@ -1903,12 +2008,12 @@ int nbmtoxky(NBM a,int *b)
 		if ( !NBM_GET(p,i) ) {
 			b[j++] = k;
 			k = 1;
-			i++;
 		} 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;
@@ -1990,7 +2095,30 @@ NBP harmonic_mul_nbm(NBM a,NBM b)
 	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)
 {
 	NODE b1,b2,br,br0;
@@ -2068,11 +2196,26 @@ void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
 
 void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
 {
-	NODE b;
+	NODE b,n;
 	NBP r,t,s;
+	NBM m;
 
-	if ( !p1 || !p2 ) *rp = 0;
-	else if ( length(BDY(p1)) < length(BDY(p2)) ) {
+	if ( !p1 || !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) ) {
 			mulnbmnbp(vl,(NBM)BDY(b),p2,&t);
 			addnbp(vl,r,t,&s); r = s;
@@ -2146,34 +2289,58 @@ void pwrnbp(VL vl,NBP a,Q q,NBP *c)
 
 void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
 {
-	NODE b1,b2;
+	NODE b1,b2,n;
 	NBP r,t,s;
 	NBM m;
 
-	if ( !p1 || !p2 ) *rp = 0;
-	else {
-		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;
-			}
-		*rp = r;
+	if ( !p1 || !p2 ) {
+		*rp = 0; return;
 	}
+	if ( OID(p1) != O_NBP ) {
+		if ( !NUM(p1) || !RATN(p1) ) 
+			error("shuffle_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("shuffle_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 = shuffle_mul_nbm(m,(NBM)BDY(b2));
+			addnbp(vl,r,t,&s); r = s;
+		}
+	*rp = r;
 }
 
 void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
 {
-	NODE b1,b2;
+	NODE b1,b2,n;
 	NBP r,t,s;
 	NBM m;
 
-	if ( !p1 || !p2 ) *rp = 0;
-	else {
-		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;
-			}
-		*rp = r;
+	if ( !p1 || !p2 ) {
+		*rp = 0; return;
 	}
+	if ( OID(p1) != O_NBP ) {
+		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;
+		}
+	*rp = r;
 }