===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/parse/puref.c,v
retrieving revision 1.6
retrieving revision 1.14
diff -u -p -r1.6 -r1.14
--- OpenXM_contrib2/asir2000/parse/puref.c	2004/03/11 09:52:56	1.6
+++ OpenXM_contrib2/asir2000/parse/puref.c	2018/03/28 05:27:22	1.14
@@ -45,11 +45,13 @@
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *
- * $OpenXM: OpenXM_contrib2/asir2000/parse/puref.c,v 1.5 2003/05/20 06:15:01 noro Exp $ 
+ * $OpenXM: OpenXM_contrib2/asir2000/parse/puref.c,v 1.13 2018/03/27 06:29:19 noro Exp $ 
 */
 #include "ca.h"
 #include "parse.h"
 
+void instoobj(PFINS ins,Obj *rp);
+
 NODE pflist;
 
 void searchpf(char *name,FUNC *fp)
@@ -284,6 +286,59 @@ void derivr(VL vl,Obj a,V v,Obj *b)
 	}
 }
 
+void simple_derivr(VL vl,Obj a,V v,Obj *b)
+{
+	Obj r,s,t,u,nm,dn;
+
+	if ( !a || NUM(a) ) 
+		*b = 0;
+	else
+		switch ( OID(a) ) {
+			case O_P:
+				pderivr(vl,a,v,b); break;
+			case O_R:
+				nm = (Obj)NM((R)a); dn = (Obj)DN((R)a);
+				/* (nm/dn)' = nm'/dn - dn'/dn*nm/dn */
+				pderivr(vl,nm,v,&s); divr(vl,s,dn,&u); reductr(vl,u,&t);
+				pderivr(vl,dn,v,&s); divr(vl,s,dn,&u); reductr(vl,u,&s); mulr(vl,s,a,&u);
+				subr(vl,t,u,&s); reductr(vl,s,b);
+				break;
+			default:
+				error("simple_derivr : invalid argument");
+	}
+}
+
+int obj_is_dependent(Obj a,V v)
+{
+	if ( !a || OID(a) <= O_N ) return 0;
+	else if ( OID(a) == O_P ) return poly_is_dependent((P)a,v);
+	else if ( OID(a) == O_R ) return poly_is_dependent(NM((R)a),v) 
+		|| poly_is_dependent(DN((R)a),v);
+	else
+		error("obj_is_dependent : not implemented");
+}
+
+int poly_is_dependent(P p,V v)
+{
+	DCP dc;
+
+	if ( !p || OID(p) <= O_N ) return 0;
+	else if ( v == VR(p) ) return 1;
+	else {
+		for ( dc = DC(p); dc; dc = NEXT(dc) )
+			if ( poly_is_dependent(COEF(dc),v) ) return 1;
+		return 0;
+	}
+}
+
+void gen_pwrr(VL vl,Obj a,Obj d,Obj *r)
+{
+	if ( INT(d) )
+		pwrr(vl,a,d,r);
+	else
+		mkpow(vl,a,d,r);
+}
+
 void substr(VL vl,int partial,Obj a,V v,Obj b,Obj *c)
 {
 	Obj nm,dn,t;
@@ -338,7 +393,8 @@ void substpr(VL vl,int partial,Obj p,V v0,Obj p0,Obj *
 		for ( c = 0, dc = DC((P)p); dc; dc = NEXT(dc) ) {
 			substpr(vl,partial,(Obj)COEF(dc),v0,p0,&t);
 			if ( DEG(dc) ) {
-				pwrp(vl,x,DEG(dc),(P *)&s); mulr(vl,s,t,&m); 
+				gen_pwrr(vl,(Obj)x,(Obj)DEG(dc),&s);
+				mulr(vl,s,t,&m); 
 				addr(vl,m,c,&a); c = a;
 			} else {
 				addr(vl,t,c,&a); c = a;
@@ -352,7 +408,8 @@ void substpr(VL vl,int partial,Obj p,V v0,Obj p0,Obj *
 		else
 			c = (Obj)COEF(dc);
 		for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc = NEXT(dc) ) {
-				subq(d,DEG(dc),(Q *)&t); pwrr(vl,p0,t,&s); mulr(vl,s,c,&m); 
+				subq(d,DEG(dc),(Q *)&t);
+				gen_pwrr(vl,p0,t,&s); mulr(vl,s,c,&m); 
 				if ( !partial )
 					substpr(vl,partial,(Obj)COEF(dc),v0,p0,&t); 
 				else
@@ -360,7 +417,8 @@ void substpr(VL vl,int partial,Obj p,V v0,Obj p0,Obj *
 				addr(vl,m,t,&c); 
 		}
 		if ( d ) {
-			pwrr(vl,p0,(Obj)d,&t); mulr(vl,t,c,&m); 
+			gen_pwrr(vl,p0,(Obj)d,&t);
+			mulr(vl,t,c,&m); 
 			c = m;
 		}
 		*pr = c;
@@ -441,7 +499,8 @@ void evalv(VL vl,V v,int prec,Obj *rp)
 void evalins(PFINS ins,int prec,Obj *rp)
 {
 	PF pf;
-	PFAD ad;
+	PFINS tins;
+	PFAD ad,tad;
 	int i;
 	Q q;
 	V v;
@@ -449,14 +508,18 @@ void evalins(PFINS ins,int prec,Obj *rp)
 	NODE n0,n;
 
 	pf = ins->pf; ad = ins->ad;
+	tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
+	tins->pf = pf; tad = tins->ad;
+	for ( i = 0; i < pf->argc; i++ ) {
+		tad[i].d = ad[i].d; evalr(CO,ad[i].arg,prec,&tad[i].arg);
+ 	}
 	for ( i = 0; i < pf->argc; i++ )
-		if ( ad[i].d || (ad[i].arg && !NUM(ad[i].arg)) )
-			break;
+		if ( tad[i].d || (tad[i].arg && !NUM(tad[i].arg)) ) break;
 	if ( (i != pf->argc) || !pf->pari ) {
-		instov(ins,&v); MKV(v,x); *rp = (Obj)x;
+		instoobj(tins,rp);
 	} else {
 		for ( n0 = 0, i = 0; i < pf->argc; i++ ) {
-			NEXTNODE(n0,n); BDY(n) = (pointer)ad[i].arg;
+			NEXTNODE(n0,n); BDY(n) = (pointer)tad[i].arg;
 		}
 		if ( prec ) {
 			NEXTNODE(n0,n); STOQ(prec,q); BDY(n) = (pointer)q;
@@ -471,16 +534,27 @@ void devalr(VL vl,Obj a,Obj *c)
 {
 	Obj nm,dn;
 	double d;
-	Real r;
+	Real r,re,im;
+	C z;
+	int nid;
 
 	if ( !a )
 		*c = 0;
 	else {
 		switch ( OID(a) ) {
 			case O_N:
-				d = ToReal(a);
-				MKReal(d,r);
-				*c = (Obj)r;
+				nid = NID((Num)a);
+				if ( nid == N_C ) {
+				  d = ToReal(((C)a)->r); MKReal(d,re);
+				  d = ToReal(((C)a)->i); MKReal(d,im);
+				  reimtocplx(re,im,&z);
+				  *c = (Obj)z;
+				} else if ( nid == N_Q || nid == N_R || nid == N_B ) {
+				  d = ToReal(a);
+				  MKReal(d,r);
+				  *c = (Obj)r;
+				} else
+				  error("devalr : unsupported");
 				break;
 			case O_P:
 				devalp(vl,(P)a,(P *)c); break;
@@ -561,8 +635,9 @@ void devalv(VL vl,V v,Obj *rp)
 
 void devalins(PFINS ins,Obj *rp)
 {
+	PFINS tins;
 	PF pf;
-	PFAD ad;
+	PFAD ad,tad;
 	int i;
 	Real r;
 	double d;
@@ -570,32 +645,88 @@ void devalins(PFINS ins,Obj *rp)
 	P x;
 
 	pf = ins->pf; ad = ins->ad;
+	tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
+	tins->pf = pf; tad = tins->ad;
+	for ( i = 0; i < pf->argc; i++ ) {
+		tad[i].d = ad[i].d; devalr(CO,ad[i].arg,&tad[i].arg);
+ 	}
 	for ( i = 0; i < pf->argc; i++ )
-		if ( ad[i].d || (ad[i].arg && !NUM(ad[i].arg)) )
-			break;
+		if ( tad[i].d || (tad[i].arg && !NUM(tad[i].arg)) ) break;
 	if ( (i != pf->argc) || !pf->libm ) {
-		instov(ins,&v); MKV(v,x); *rp = (Obj)x;
+		instoobj(tins,rp);
 	} else {
+		for ( i = 0; i < pf->argc; i++ )
+			if ( tad[i].arg && NID((Num)tad[i].arg) == N_C )
+				error("devalins : not supported");
 		switch ( pf->argc ) {
 			case 0:
 				d = (*pf->libm)(); break;
 			case 1:
-				d = (*pf->libm)(ToReal(ad[0].arg)); break;
+				d = (*pf->libm)(ToReal(tad[0].arg)); break;
 			case 2:
-				d = (*pf->libm)(ToReal(ad[0].arg),ToReal(ad[1].arg)); break;
+				d = (*pf->libm)(ToReal(tad[0].arg),ToReal(tad[1].arg)); break;
 			case 3:
-				d = (*pf->libm)(ToReal(ad[0].arg),ToReal(ad[1].arg),
-					ToReal(ad[2].arg)); break;
+				d = (*pf->libm)(ToReal(tad[0].arg),ToReal(tad[1].arg),
+					ToReal(tad[2].arg)); break;
 			case 4:
-				d = (*pf->libm)(ToReal(ad[0].arg),ToReal(ad[1].arg),
-					ToReal(ad[2].arg),ToReal(ad[3].arg)); break;
+				d = (*pf->libm)(ToReal(tad[0].arg),ToReal(tad[1].arg),
+					ToReal(tad[2].arg),ToReal(tad[3].arg)); break;
 			default:
-				error("devalv : not supported");
+				error("devalins : not supported");
 		}
 		MKReal(d,r); *rp = (Obj)r;
 	}
 }
 
+extern int evalef,bigfloat;
+
+void simplify_elemfunc_ins(PFINS ins,Obj *rp)
+{
+  if ( evalef ) {
+    if ( bigfloat ) evalins(ins,0,rp); 
+    else devalins(ins,rp);
+  } else instoobj(ins,rp);
+}
+
+void simplify_factorial_ins(PFINS ins,Obj *rp)
+{
+  PFAD ad;
+  Obj a;
+  Q q;
+
+  ad = ins->ad;
+  a = ad[0].arg;
+  if ( !ad[0].d && INT(a) && ( !a || (PL(NM((Q)a)) == 1 && SGN((Q)a) > 0) ) ) {
+    factorial(QTOS((Q)a),&q);
+    *rp = (Obj)q;
+  } else simplify_elemfunc_ins(ins,rp);
+}
+
+void simplify_abs_ins(PFINS ins,Obj *rp)
+{
+  PFAD ad;
+  Obj a;
+  Q q;
+  double t;
+  Real r;
+  struct oNODE arg0;
+
+  ad = ins->ad;
+  a = ad[0].arg;
+  if ( !ad[0].d && NUM(a) && (!a || RATN(a)) ) {
+    if ( !a || SGN((Q)a) > 0 ) *rp = (Obj)a;
+    else {
+      chsgnq((Q)a,&q); *rp = (Obj)q;
+    }
+  } else if ( !ad[0].d && REAL(a) ) {
+    t = fabs(((Real)a)->body);
+    MKReal(t,r); *rp = (Obj)r;
+  } else if ( !ad[0].d && BIGFLOAT(a) ) {
+    arg0.body = (pointer)a; arg0.next = 0;
+    mp_abs(&arg0,rp);
+  } else simplify_elemfunc_ins(ins,rp);
+}
+
 void simplify_ins(PFINS ins,Obj *rp)
 {
 	V v;
@@ -604,17 +735,20 @@ void simplify_ins(PFINS ins,Obj *rp)
 	if ( ins->pf->simplify )
 		(*ins->pf->simplify)(ins,rp);
 	else {
-		instov(ins,&v); MKV(v,t); *rp = (Obj)t;
+		instoobj(ins,rp);
 	}
 }
 
-void instov(PFINS ins,V *vp)
+void instoobj(PFINS ins,Obj *rp)
 {
-	V v;
+	V v,newv;
+	P t;
 
 	NEWV(v); NAME(v) = 0; 
 	v->attr = (pointer)V_PF; v->priv = (pointer)ins;
-	appendpfins(v,vp);
+	appendpfins(v,&newv);
+	MKV(newv,t);
+	*rp = (Obj)t;
 }
 
 void substfr(VL vl,Obj a,PF u,PF f,Obj *c)
@@ -658,7 +792,8 @@ void substfp(VL vl,Obj p,PF u,PF f,Obj *pr)
 			for ( c = 0; dc; dc = NEXT(dc) ) {
 				substfp(vl,(Obj)COEF(dc),u,f,&t);
 				if ( DEG(dc) ) {
-					pwrp(vl,x,DEG(dc),(P *)&s); mulr(vl,s,t,&m); 
+					gen_pwrr(vl,(Obj)x,(Obj)DEG(dc),&s);
+					mulr(vl,s,t,&m); 
 					addr(vl,m,c,&a); c = a;
 				} else {
 					addr(vl,t,c,&a); c = a;
@@ -668,11 +803,12 @@ void substfp(VL vl,Obj p,PF u,PF f,Obj *pr)
 			substfv(vl,v,u,f,&p0);
 			substfp(vl,(Obj)COEF(dc),u,f,&c); 
 			for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc = NEXT(dc) ) {
-					subq(d,DEG(dc),(Q *)&t); pwrr(vl,p0,t,&s); mulr(vl,s,c,&m); 
+					subq(d,DEG(dc),(Q *)&t); 
+					gen_pwrr(vl,p0,t,&s); mulr(vl,s,c,&m); 
 					substfp(vl,(Obj)COEF(dc),u,f,&t); addr(vl,m,t,&c); 
 			}
 			if ( d ) {
-				pwrr(vl,p0,(Obj)d,&t); mulr(vl,t,c,&m); 
+				gen_pwrr(vl,p0,(Obj)d,&t); mulr(vl,t,c,&m); 
 				c = m;
 			}
 		}
@@ -711,6 +847,6 @@ void substfv(VL vl,V v,PF u,PF f,Obj *c)
 		for ( i = 0; i < f->argc; i++ ) {
 			tad[i].d = ad[i].d; substfr(vl,ad[i].arg,u,f,&tad[i].arg);
 		}
-		instov(tins,&v); MKV(v,t); *c = (Obj)t;
+		instoobj(tins,c);
 	}
 }