[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.4

version 1.2, 2018/09/28 08:20:28 version 1.4, 2019/10/17 03:03:12
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.3 2019/06/04 07:11:23 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 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;
Line 253  double  ToRealUp(Num a)
Line 293  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:
Line 280  void  Num2double(Num a, double *inf, double *sup)
Line 321  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));
Line 480  void  divitvd(Num a, Num b, IntervalDouble *c)
Line 521  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 539  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;

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.4

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>