===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2018/engine/d-itv.c,v
retrieving revision 1.3
retrieving revision 1.5
diff -u -p -r1.3 -r1.5
--- OpenXM_contrib2/asir2018/engine/d-itv.c	2019/06/04 07:11:23	1.3
+++ OpenXM_contrib2/asir2018/engine/d-itv.c	2019/11/12 10:53:22	1.5
@@ -1,5 +1,5 @@
 /*
- * $OpenXM: OpenXM_contrib2/asir2018/engine/d-itv.c,v 1.2 2018/09/28 08:20:28 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2018/engine/d-itv.c,v 1.4 2019/10/17 03:03:12 kondoh Exp $
 */
 #if defined(INTERVAL)
 #include <float.h>
@@ -11,6 +11,12 @@
 #endif
 #endif
 
+double  addulpd(double);
+double  subulpd(double);
+extern int mpfr_roundmode;
+Num tobf(Num,int);
+
+
 #if defined(ITVDEBUG)
 void printbinint(int d)
 {
@@ -33,6 +39,7 @@ void printbinint(int d)
 }
 #endif
 
+#if 0
 double NatToRealUp(N a, int *expo)
 {
   int *p;
@@ -181,7 +188,6 @@ static double  Q2doubleUp(Q a)
   }
 }
 
-#if 0
 static double  PARI2doubleDown(BF a)
 {
         //GEN     p;
@@ -200,6 +206,34 @@ static double  PARI2doubleUp(BF a)
 }
 #endif
 
+static double  Q2doubleDown(Q a)
+{
+	Num b;
+	double rd;
+	int current_roundmode=0;
+
+	current_roundmode = mpfr_roundmode;
+ 	mpfr_roundmode = MPFR_RNDD;
+	b = tobf((Num)a,60);
+ 	mpfr_roundmode = current_roundmode;
+	rd = mpfr_get_d(BDY((BF)b),MPFR_RNDD);
+	return rd;
+}
+
+static double  Q2doubleUp(Q a)
+{
+	Num b;
+	double rd;
+	int current_roundmode=0;
+
+	current_roundmode = mpfr_roundmode;
+ 	mpfr_roundmode = MPFR_RNDU;
+	b = tobf((Num)a,60);
+ 	mpfr_roundmode = current_roundmode;
+	rd = mpfr_get_d(BDY((BF)b),MPFR_RNDU);
+	return rd;
+}
+
 double  subulpd(double d)
 {
   double inf;
@@ -221,7 +255,7 @@ double  addulpd(double d)
   return (-sup);
 }
 
-double  ToRealDown(Num a)
+double  toRealDown(Num a)
 {
   double inf;
 
@@ -232,25 +266,26 @@ double  ToRealDown(Num a)
     case N_R:
       inf = subulpd(BDY((Real)a)); break;
     case N_B:
-      //inf = PARI2doubleDown((BF)a); break;
-	inf = 0;
-      error("ToRealDown: not supported operands.");
-	break;
+		inf = mpfr_get_d(BDY((BF)a),MPFR_RNDD);
+		break;
     case N_IP:
-      inf = ToRealDown(INF((Itv)a));
+      inf = toRealDown(INF((Itv)a));
       break;
     case N_IntervalDouble:
       inf = INF((IntervalDouble)a); break;
+	case N_IntervalBigFloat:
+		inf = mpfr_get_d(BDY((BF)INF((IntervalBigFloat)a)),MPFR_RNDD);
+		break;
     case N_A:
     default:
       inf = 0.0;
-      error("ToRealDown: not supported operands.");
+      error("toRealDown: not supported operands.");
       break;
   }
   return inf;
 }
 
-double  ToRealUp(Num a)
+double  toRealUp(Num a)
 {
   double sup;
 
@@ -261,18 +296,19 @@ double  ToRealUp(Num a)
     case N_R:
       sup = addulpd(BDY((Real)a)); break;
     case N_B:
-      //sup = PARI2doubleUp((BF)a); break;
-      sup = 0;
-      error("ToRealUp: not supported operands.");
-      break;
+		sup = mpfr_get_d(BDY((BF)a),MPFR_RNDU);
+		break;
     case N_IP:
-      sup = ToRealUp(SUP((Itv)a)); break;
+      sup = toRealUp(SUP((Itv)a)); break;
     case N_IntervalDouble:
       sup = SUP((IntervalDouble)a); break;
+	case N_IntervalBigFloat:
+		sup = mpfr_get_d(BDY((BF)INF((IntervalBigFloat)a)),MPFR_RNDU);
+		break;
     case N_A:
     default:
       sup = 0.0;
-      error("ToRealUp: not supported operands.");
+      error("toRealUp: not supported operands.");
       break;
   }
   return sup;
@@ -281,6 +317,9 @@ double  ToRealUp(Num a)
 
 void  Num2double(Num a, double *inf, double *sup)
 {
+  *inf = 0.0;
+  *sup = 0.0;
+  if (a && NUM(a) )
   switch ( NID(a) ) {
     case N_Q:
       *inf = Q2doubleDown((Q)a);
@@ -291,14 +330,12 @@ void  Num2double(Num a, double *inf, double *sup)
       *sup = BDY((Real)a);
       break;
     case N_B:
-      //*inf = PARI2doubleDown((BF)a);
-      //*sup = PARI2doubleUp((BF)a);
       *inf = mpfr_get_d(BDY((BF)a), MPFR_RNDD);
       *sup = mpfr_get_d(BDY((BF)a), MPFR_RNDU);
       break;
     case N_IP:
-      *inf = ToRealDown(INF((Itv)a));
-      *sup = ToRealUp(SUP((Itv)a));
+      *inf = toRealDown(INF((Itv)a));
+      *sup = toRealUp(SUP((Itv)a));
       break;
     case N_IntervalDouble:
       *inf = INF((IntervalDouble)a);
@@ -306,8 +343,6 @@ void  Num2double(Num a, double *inf, double *sup)
       break;
     case N_A:
     default:
-      *inf = 0.0;
-      *sup = 0.0;
       error("Num2double: not supported operands.");
       break;
   }
@@ -493,7 +528,7 @@ void  divitvd(Num a, Num b, IntervalDouble *c)
 
 void  pwritvd(Num a, Num e, IntervalDouble *c)
 {
-  int  ei;
+  long  ei;
   IntervalDouble  t;
 
   if ( !e ) {
@@ -511,17 +546,19 @@ void  pwritvd(Num a, Num e, IntervalDouble *c)
     error("pwritvd : can't calculate a fractional power");
 #endif
   } else {
-    ei = QTOS((Q)e);
+    //ei = QTOS((Q)e);
+	ei = mpz_get_si(BDY((Z)e));
     if (ei<0) ei = -ei;
     pwritv0d((IntervalDouble)a,ei,&t);
-    if ( SGN((Q)e) < 0 )
+//    if ( SGN((Q)e) < 0 )
+    if ( sgnq((Q)e) < 0 )
       divnum(0,(Num)ONE,(Num)t,(Num *)c);
     else
       *c = t;
   }
 }
 
-void  pwritv0d(IntervalDouble a, int e, IntervalDouble *c)
+void  pwritv0d(IntervalDouble a, long e, IntervalDouble *c)
 {
   double inf, sup;
   double t, Xmin, Xmax;
@@ -687,6 +724,24 @@ void  absitvd(IntervalDouble a, Num *b)
     MKReal(t,rp);
     *b = (Num)rp;
   }
+}
+void  absintvald(IntervalDouble a, IntervalDouble *b)
+{
+  double  ai,as,inf,sup;
+
+  ai = INF(a);
+  as = SUP(a);
+  if ( as < 0 ) {
+    sup = -ai;
+    inf = -as;
+  } else if (ai < 0) {
+    inf = 0.0;
+  sup = MAX(as, -ai);
+  } else {
+    inf = ai;
+    sup = as;
+  }
+  MKIntervalDouble(inf,sup,*b);
 }
 
 void  distanceitvd(IntervalDouble a, IntervalDouble b, Num *c)