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

Diff for /OpenXM_contrib2/asir2018/engine/f-itv.c between version 1.3 and 1.5

version 1.3, 2019/06/04 07:11:23 version 1.5, 2019/11/12 10:53:22
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2018/engine/f-itv.c,v 1.2 2018/09/28 08:20:28 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2018/engine/f-itv.c,v 1.4 2019/10/17 03:03:12 kondoh Exp $
 */  */
 #if defined(INTERVAL)  #if defined(INTERVAL)
 #include "ca.h"  #include "ca.h"
Line 120  Num tobf(Num,int);
Line 120  Num tobf(Num,int);
   
 #define BFPREC(a) (((BF)(a))->body->_mpfr_prec)  #define BFPREC(a) (((BF)(a))->body->_mpfr_prec)
   
   double mpfr2dblDown(mpfr_t a)
   {
           return mpfr_get_d(a,MPFR_RNDD);
   }
   
   double mpfr2dblUp(mpfr_t a)
   {
           return mpfr_get_d(a,MPFR_RNDU);
   }
   
   
   void toInterval(Num a, int prec, int type, Num *rp)
   {
           if ( ! a ) {
                   *rp = 0;
           } else if (type == EvalIntervalDouble) {
                   if (NID(a)==N_C) {
                           double inf,sup;
                           C z;
                           IntervalDouble re, im;
   
                           if ( ! ((C)a)->r ) {
                                   re = 0;
                           } else {
                                   inf = toRealDown(((C)a)->r);
                                   sup = toRealUp(((C)a)->r);
                                   MKIntervalDouble(inf,sup,re);
               }
                           if ( ! ((C)a)->i ) {
                                   im = 0;
                           } else {
                                   inf = toRealDown(((C)a)->i);
                                   sup = toRealUp(((C)a)->i);
                                   MKIntervalDouble(inf,sup,im);
                           }
                           if ( !re && !im )
                                   z = 0;
                           else
                                   reimtocplx((Num)re,(Num)im,(Num *)&z);
                           *rp = (Num)z;
                   } else {
                           double inf,sup;
                           IntervalDouble c;
   
                           inf = toRealDown(a);
                           sup = toRealUp(a);
   
                           MKIntervalDouble(inf,sup,c);
                           *rp = (Num) c;
                   }
           } else if (type == EvalIntervalBigFloat) {
                   if (NID(a)==N_C) {
                           Num ai,as;
                           Num inf,sup;
                           C z;
                           IntervalBigFloat re, im;
                           int current_roundmode;
   
                           current_roundmode = mpfr_roundmode;
   
                           if ( ! ((C)a)->r )
                                   re = 0;
                           else {
                                   itvtois((Itv)((C)a)->r,&ai,&as);
                           mpfr_roundmode = MPFR_RNDD;
                                   inf = tobf(ai, prec);
                                   mpfr_roundmode = MPFR_RNDU;
                                   sup = tobf(as, prec);
                                   istoitv(inf,sup,(Itv *)&re);
                           }
   
                           if ( ! ((C)a)->i )
                                   im = 0;
                           else {
                                   itvtois((Itv)((C)a)->i,&ai,&as);
                           mpfr_roundmode = MPFR_RNDD;
                                   inf = tobf(ai, prec);
                                   mpfr_roundmode = MPFR_RNDU;
                                   sup = tobf(as, prec);
                                   istoitv(inf,sup,(Itv *)&im);
                           }
   
                   mpfr_roundmode = current_roundmode;
                           reimtocplx((Num)re,(Num)im,(Num *)&z);
                           *rp = (Num)z;
                   } else {
                           Num ai,as;
                           Num inf,sup;
                           IntervalBigFloat c;
                           int current_roundmode;
   
                           itvtois((Itv)a,&ai,&as);
   
                           current_roundmode = mpfr_roundmode;
                   mpfr_roundmode = MPFR_RNDD;
                           inf = tobf(ai, prec);
                           mpfr_roundmode = MPFR_RNDU;
                           sup = tobf(as, prec);
                           istoitv(inf,sup,(Itv *)&c);
                   mpfr_roundmode = current_roundmode;
                           *rp = (Num) c;
                   }
           } else {
                   error("toInterval: not supported types.");
                   *rp = 0;
           }
   }
   
   
 void additvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *rp)  void additvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *rp)
 {  {
   Num  ai,as,bi,bs;    Num  ai,as,bi,bs;
Line 387  void divitvf(IntervalBigFloat a, IntervalBigFloat b, I
Line 495  void divitvf(IntervalBigFloat a, IntervalBigFloat b, I
   
 void pwritvf(Itv a, Num e, Itv *c)  void pwritvf(Itv a, Num e, Itv *c)
 {  {
   int ei;    long ei;
   Itv t;    Itv t;
   
   if ( !e )    if ( !e )
Line 403  void pwritvf(Itv a, Num e, Itv *c)
Line 511  void pwritvf(Itv a, Num e, Itv *c)
     error("pwritvf() : can't calculate a fractional power");      error("pwritvf() : can't calculate a fractional power");
 #endif  #endif
   } else {    } else {
     ei = QTOS((Q)e);      //ei = QTOS((Q)e);
       ei = mpz_get_si(BDY((Q)e));
     if (ei<0) ei = -ei;      if (ei<0) ei = -ei;
     pwritv0f(a,ei,&t);      pwritv0f(a,ei,&t);
     if ( SGN((Q)e) < 0 )  //    if ( SGN((Q)e) < 0 )
       if ( sgnq((Q)e) < 0 )
         divitvf((IntervalBigFloat)ONE,(IntervalBigFloat)t,(IntervalBigFloat *)c); /* bug ?? */          divitvf((IntervalBigFloat)ONE,(IntervalBigFloat)t,(IntervalBigFloat *)c); /* bug ?? */
     else      else
       *c = t;        *c = t;
   }    }
 }  }
   
 void pwritv0f(Itv a, int e, Itv *c)  void pwritv0f(Itv a, long e, Itv *c)
 {  {
   Num inf, sup;    Num inf, sup;
   Num ai,Xmin,Xmax;    Num ai,Xmin,Xmax;
Line 424  void pwritv0f(Itv a, int e, Itv *c)
Line 534  void pwritv0f(Itv a, int e, Itv *c)
   if ( e == 1 )    if ( e == 1 )
     *c = a;      *c = a;
   else {    else {
     STOQ(e,ne);      STOZ(e,ne);
     if ( !(e%2) ) {      if ( !(e%2) ) {
       if ( initvp(0,a) ) {        if ( initvp(0,a) ) {
         Xmin = 0;          Xmin = 0;
Line 533  void miditvf(Itv a, Num *b)
Line 643  void miditvf(Itv a, Num *b)
   else if ( (NID(a) <= N_B) )    else if ( (NID(a) <= N_B) )
     *b = (Num)a;      *b = (Num)a;
   else {    else {
     STOQ(2,TWO);      //STOQ(2,TWO);
     itvtois(a,&ai,&as);      itvtois(a,&ai,&as);
     addbf(ai,as,&t);      addbf(ai,as,&t);
     divbf(t,(Num)TWO,b);      divbf(t,(Num)TWO,b);
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.5
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>