[BACK]Return to d-itv.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2018 / engine

Diff for /OpenXM_contrib2/asir2018/engine/d-itv.c between version 1.2 and 1.5

version 1.2, 2018/09/28 08:20:28 version 1.5, 2019/11/12 10:53:22
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2018/engine/d-itv.c,v 1.1 2018/09/19 05:45:07 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2018/engine/d-itv.c,v 1.4 2019/10/17 03:03:12 kondoh Exp $
 */  */
 #if defined(INTERVAL)  #if defined(INTERVAL)
 #include <float.h>  #include <float.h>
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
   #if 0
 #if defined(PARI)  #if defined(PARI)
 #include "genpari.h"  #include "genpari.h"
 #endif  #endif
   #endif
   
   double  addulpd(double);
   double  subulpd(double);
   extern int mpfr_roundmode;
   Num tobf(Num,int);
   
   
 #if defined(ITVDEBUG)  #if defined(ITVDEBUG)
 void printbinint(int d)  void printbinint(int d)
 {  {
Line 31  void printbinint(int d)
Line 39  void printbinint(int d)
 }  }
 #endif  #endif
   
   #if 0
 double NatToRealUp(N a, int *expo)  double NatToRealUp(N a, int *expo)
 {  {
   int *p;    int *p;
Line 181  static double  Q2doubleUp(Q a)
Line 190  static double  Q2doubleUp(Q a)
   
 static double  PARI2doubleDown(BF a)  static double  PARI2doubleDown(BF a)
 {  {
         GEN     p;          //GEN     p;
           Num     p;
   double  d;    double  d;
   
         ritopa(a, &p);          ritopa(a, &p);
Line 194  static double  PARI2doubleUp(BF a)
Line 204  static double  PARI2doubleUp(BF a)
 {  {
   return PARI2doubleDown(a);    return PARI2doubleDown(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  subulpd(double d)
 {  {
   double inf;    double inf;
Line 216  double  addulpd(double d)
Line 255  double  addulpd(double d)
   return (-sup);    return (-sup);
 }  }
   
 double  ToRealDown(Num a)  double  toRealDown(Num a)
 {  {
   double inf;    double inf;
   
Line 227  double  ToRealDown(Num a)
Line 266  double  ToRealDown(Num a)
     case N_R:      case N_R:
       inf = subulpd(BDY((Real)a)); break;        inf = subulpd(BDY((Real)a)); break;
     case N_B:      case N_B:
       inf = PARI2doubleDown((BF)a); break;                  inf = mpfr_get_d(BDY((BF)a),MPFR_RNDD);
                   break;
     case N_IP:      case N_IP:
       inf = ToRealDown(INF((Itv)a));        inf = toRealDown(INF((Itv)a));
       break;        break;
     case N_IntervalDouble:      case N_IntervalDouble:
       inf = INF((IntervalDouble)a); break;        inf = INF((IntervalDouble)a); break;
           case N_IntervalBigFloat:
                   inf = mpfr_get_d(BDY((BF)INF((IntervalBigFloat)a)),MPFR_RNDD);
                   break;
     case N_A:      case N_A:
     default:      default:
       inf = 0.0;        inf = 0.0;
       error("ToRealDown: not supported operands.");        error("toRealDown: not supported operands.");
       break;        break;
   }    }
   return inf;    return inf;
 }  }
   
 double  ToRealUp(Num a)  double  toRealUp(Num a)
 {  {
   double sup;    double sup;
   
Line 253  double  ToRealUp(Num a)
Line 296  double  ToRealUp(Num a)
     case N_R:      case N_R:
       sup = addulpd(BDY((Real)a)); break;        sup = addulpd(BDY((Real)a)); break;
     case N_B:      case N_B:
       sup = PARI2doubleUp((BF)a); break;                  sup = mpfr_get_d(BDY((BF)a),MPFR_RNDU);
                   break;
     case N_IP:      case N_IP:
       sup = ToRealUp(SUP((Itv)a)); break;        sup = toRealUp(SUP((Itv)a)); break;
     case N_IntervalDouble:      case N_IntervalDouble:
       sup = SUP((IntervalDouble)a); break;        sup = SUP((IntervalDouble)a); break;
           case N_IntervalBigFloat:
                   sup = mpfr_get_d(BDY((BF)INF((IntervalBigFloat)a)),MPFR_RNDU);
                   break;
     case N_A:      case N_A:
     default:      default:
       sup = 0.0;        sup = 0.0;
       error("ToRealUp: not supported operands.");        error("toRealUp: not supported operands.");
       break;        break;
   }    }
   return sup;    return sup;
Line 270  double  ToRealUp(Num a)
Line 317  double  ToRealUp(Num a)
   
 void  Num2double(Num a, double *inf, double *sup)  void  Num2double(Num a, double *inf, double *sup)
 {  {
     *inf = 0.0;
     *sup = 0.0;
     if (a && NUM(a) )
   switch ( NID(a) ) {    switch ( NID(a) ) {
     case N_Q:      case N_Q:
       *inf = Q2doubleDown((Q)a);        *inf = Q2doubleDown((Q)a);
Line 280  void  Num2double(Num a, double *inf, double *sup)
Line 330  void  Num2double(Num a, double *inf, double *sup)
       *sup = BDY((Real)a);        *sup = BDY((Real)a);
       break;        break;
     case N_B:      case N_B:
       *inf = PARI2doubleDown((BF)a);        *inf = mpfr_get_d(BDY((BF)a), MPFR_RNDD);
       *sup = PARI2doubleUp((BF)a);        *sup = mpfr_get_d(BDY((BF)a), MPFR_RNDU);
       break;        break;
     case N_IP:      case N_IP:
       *inf = ToRealDown(INF((Itv)a));        *inf = toRealDown(INF((Itv)a));
       *sup = ToRealUp(SUP((Itv)a));        *sup = toRealUp(SUP((Itv)a));
       break;        break;
     case N_IntervalDouble:      case N_IntervalDouble:
       *inf = INF((IntervalDouble)a);        *inf = INF((IntervalDouble)a);
Line 293  void  Num2double(Num a, double *inf, double *sup)
Line 343  void  Num2double(Num a, double *inf, double *sup)
       break;        break;
     case N_A:      case N_A:
     default:      default:
       *inf = 0.0;  
       *sup = 0.0;  
       error("Num2double: not supported operands.");        error("Num2double: not supported operands.");
       break;        break;
   }    }
Line 480  void  divitvd(Num a, Num b, IntervalDouble *c)
Line 528  void  divitvd(Num a, Num b, IntervalDouble *c)
   
 void  pwritvd(Num a, Num e, IntervalDouble *c)  void  pwritvd(Num a, Num e, IntervalDouble *c)
 {  {
   int  ei;    long  ei;
   IntervalDouble  t;    IntervalDouble  t;
   
   if ( !e ) {    if ( !e ) {
Line 498  void  pwritvd(Num a, Num e, IntervalDouble *c)
Line 546  void  pwritvd(Num a, Num e, IntervalDouble *c)
     error("pwritvd : can't calculate a fractional power");      error("pwritvd : can't calculate a fractional power");
 #endif  #endif
   } else {    } else {
     ei = ZTOS((Q)e);      //ei = QTOS((Q)e);
           ei = mpz_get_si(BDY((Z)e));
       if (ei<0) ei = -ei;
     pwritv0d((IntervalDouble)a,ei,&t);      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);        divnum(0,(Num)ONE,(Num)t,(Num *)c);
     else      else
       *c = t;        *c = t;
   }    }
 }  }
   
 void  pwritv0d(IntervalDouble a, int e, IntervalDouble *c)  void  pwritv0d(IntervalDouble a, long e, IntervalDouble *c)
 {  {
   double inf, sup;    double inf, sup;
   double t, Xmin, Xmax;    double t, Xmin, Xmax;
Line 673  void  absitvd(IntervalDouble a, Num *b)
Line 724  void  absitvd(IntervalDouble a, Num *b)
     MKReal(t,rp);      MKReal(t,rp);
     *b = (Num)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)  void  distanceitvd(IntervalDouble a, IntervalDouble b, Num *c)
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.5
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>