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

Diff for /OpenXM_contrib2/asir2018/engine/p-itv.c between version 1.3 and 1.6

version 1.3, 2019/06/04 07:11:23 version 1.6, 2019/11/12 10:53:22
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2018/engine/p-itv.c,v 1.2 2018/09/28 08:20:28 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2018/engine/p-itv.c,v 1.5 2019/11/07 08:50:49 kondoh Exp $
 */  */
 #if defined(INTERVAL)  #if defined(INTERVAL)
 #include "ca.h"  #include "ca.h"
Line 15  Num tobf(Num,int);
Line 15  Num tobf(Num,int);
 //int     initvp(Num , Itv );  //int     initvp(Num , Itv );
   
 extern int mpfr_roundmode;  extern int mpfr_roundmode;
   
 extern int zerorewrite;  extern int zerorewrite;
   
 void itvtois(Itv a, Num *inf, Num *sup)  void itvtois(Itv a, Num *inf, Num *sup)
Line 24  void itvtois(Itv a, Num *inf, Num *sup)
Line 23  void itvtois(Itv a, Num *inf, Num *sup)
     *inf = *sup = (Num)0;      *inf = *sup = (Num)0;
   else if ( NID(a) <= N_B ) {    else if ( NID(a) <= N_B ) {
     *inf = (Num)a; *sup = (Num)a;      *inf = (Num)a; *sup = (Num)a;
     } else if ( ITVD(a) ) {
       double dinf, dsup;
       Real rinf, rsup;
       dinf = INF((IntervalDouble)a); dsup = SUP((IntervalDouble)a);
       MKReal(dinf, rinf); MKReal(dsup, rsup);
       *inf = (Num)rinf; *sup = (Num)rsup;
   } else {    } else {
     *inf = INF(a);      *inf = INF(a);
     *sup = SUP(a);      *sup = SUP(a);
Line 32  void itvtois(Itv a, Num *inf, Num *sup)
Line 37  void itvtois(Itv a, Num *inf, Num *sup)
   
 void istoitv(Num inf, Num sup, Itv *rp)  void istoitv(Num inf, Num sup, Itv *rp)
 {  {
   Num  i, s;    Num  ii, ss;
   Num  ni, ns;    Num  ni, ns;
   Itv c;    Itv c;
   int  type=0;    int  type=0;
Line 50  void istoitv(Num inf, Num sup, Itv *rp)
Line 55  void istoitv(Num inf, Num sup, Itv *rp)
   
   current_roundmode = mpfr_roundmode;    current_roundmode = mpfr_roundmode;
   if ( !ns ) {    if ( !ns ) {
     s = 0;      ss = 0;
   }    }
   else if ( NID( ns )==N_B ) {    else if ( NID( ns )==N_B ) {
     type = 1;      type = 1;
   
     mpfr_roundmode = MPFR_RNDU;      mpfr_roundmode = MPFR_RNDU;
     s = tobf(ns, DEFAULTPREC);      ss = tobf(ns, DEFAULTPREC);
     //ToBf(sup, (BF *)&s);      //ToBf(sup, (BF *)&s);
   } else {    } else {
     s = ns;      ss = ns;
   }    }
   if ( !ni ) {    if ( !ni ) {
     i = 0;      ii = 0;
   }    }
   else if ( NID( ni )==N_B ) {    else if ( NID( ni )==N_B ) {
     type = 1;      type = 1;
   
     mpfr_roundmode = MPFR_RNDD;      mpfr_roundmode = MPFR_RNDD;
     i = tobf(inf, DEFAULTPREC);      ii = tobf(ni, DEFAULTPREC);
     //ToBf(inf, (BF *)&i);      //ToBf(inf, (BF *)&i);
   } else {    } else {
     i = ni;      ii = ni;
   }    }
   mpfr_roundmode = current_roundmode;    mpfr_roundmode = current_roundmode;
   
Line 82  void istoitv(Num inf, Num sup, Itv *rp)
Line 87  void istoitv(Num inf, Num sup, Itv *rp)
   } else {    } else {
     NEWItvP(c);      NEWItvP(c);
   }    }
   INF(c) = i; SUP(c) = s;    INF(c) = ii; SUP(c) = ss;
   
   if (zerorewrite && initvp(0,c) ) {    *rp = c;
     *rp = 0;  
     zerorewriteCount++;    ZEROREWRITE
   } else {  
     *rp = c;  
   }  
 }  }
   
 void additvp(Itv a, Itv b, Itv *c)  void additvp(Itv a, Itv b, Itv *c)
Line 272  void divitvp(Itv a, Itv b, Itv *c)
Line 274  void divitvp(Itv a, Itv b, Itv *c)
   
 void pwritvp(Itv a, Num e, Itv *c)  void pwritvp(Itv a, Num e, Itv *c)
 {  {
   int ei;    long ei;
   Itv t;    Itv t;
   
   if ( !e )    if ( !e )
Line 288  void pwritvp(Itv a, Num e, Itv *c)
Line 290  void pwritvp(Itv a, Num e, Itv *c)
     error("pwritv : can't calculate a fractional power");      error("pwritv : 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));
     pwritv0p(a,ei,&t);      pwritv0p(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 pwritv0p(Itv a, int e, Itv *c)  void pwritv0p(Itv a, long e, Itv *c)
 {  {
   Num inf, sup;    Num inf, sup;
   Num ai,Xmin,Xmax;    Num ai,Xmin,Xmax;
Line 306  void pwritv0p(Itv a, int e, Itv *c)
Line 310  void pwritv0p(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 394  void miditvp(Itv a, Num *b)
Line 398  void miditvp(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);      //STOZ(2,TWO);
     itvtois(a,&ai,&as);      itvtois(a,&ai,&as);
     addnum(0,ai,as,&t);      addnum(0,ai,as,&t);
     divnum(0,t,(Num)TWO,b);      divnum(0,t,(Num)TWO,b);
Line 462  void absitvp(Itv a, Num *b)
Line 466  void absitvp(Itv a, Num *b)
     else        *b = bs;      else        *b = bs;
   }    }
 }  }
   
   void absintvalp(Itv a, Itv *c)
   {
     Num  ai,as,mai;
     Num  inf,sup;
   
     itvtois(a,&ai,&as);
     if ( compnum(0,as,0 ) < 0 ) {
       chsgnnum(ai, &sup);
       chsgnnum(as, &inf);
     } else if ( compnum(0,ai,0) < 0 ) {
           inf = 0;
       chsgnnum(ai, &mai);
       if ( compnum(0,as,mai ) > 0 )       sup = as;
       else                                                        sup = mai;
     } else {
       inf = ai;
       sup = as;
     }
     istoitv(inf,sup,c);
   }
   
   
 void distanceitvp(Itv a, Itv b, Num *c)  void distanceitvp(Itv a, Itv b, Num *c)
 {  {
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.6
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>