[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.2 and 1.3

version 1.2, 2018/09/28 08:20:28 version 1.3, 2019/06/04 07:11:23
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2018/engine/f-itv.c,v 1.1 2018/09/19 05:45:07 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2018/engine/f-itv.c,v 1.2 2018/09/28 08:20:28 noro Exp $
 */  */
 #if defined(INTERVAL)  #if defined(INTERVAL)
 #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"
 #include "itv-pari.h"  #include "itv-pari.h"
Line 109  void addulp(IntervalBigFloat a, IntervalBigFloat *c)
Line 110  void addulp(IntervalBigFloat a, IntervalBigFloat *c)
   } else  sup = as;    } else  sup = as;
   istoitv(inf,sup, (Itv *)c);    istoitv(inf,sup, (Itv *)c);
 }  }
   #endif
   
 void additvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)  // in builtin/bfaux.c
   extern int mpfr_roundmode;
   
   // in engine/bf.c
   Num tobf(Num,int);
   
   #define BFPREC(a) (((BF)(a))->body->_mpfr_prec)
   
   
   void additvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *rp)
 {  {
   Num  ai,as,bi,bs,mas,mbs,tmp;    Num  ai,as,bi,bs;
     IntervalBigFloat c;
   //,mas,mbs;
   //,tmp;
   Num  inf,sup;    Num  inf,sup;
   GEN  pa, pb, z;    //GEN  pa, pb, z;
   long  ltop, lbot;    //long  ltop, lbot;
     int current_roundmode;
   
   if ( !a )    if ( !a )
     *c = b;      *rp = b;
   else if ( !b )    else if ( !b )
     *c = a;      *rp = a;
   else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )    else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
     addnum(0,(Num)a,(Num)b,(Num *)c);      addnum(0,(Num)a,(Num)b,(Num *)rp);
   else {    else {
     itvtois((Itv)a,&inf,&sup);      itvtois((Itv)a,&ai,&as);
     ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);      itvtois((Itv)b,&bi,&bs);
     itvtois((Itv)b,&inf,&sup);  
     ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);      current_roundmode = mpfr_roundmode;
   
       mpfr_roundmode = MPFR_RNDD;
       ai = tobf(ai, DEFAULTPREC);
       bi = tobf(bi, DEFAULTPREC);
       //addnum(0,ai,bi,&inf);
       addbf(ai,bi,&inf);
   
       mpfr_roundmode = MPFR_RNDU;
       as = tobf(as, DEFAULTPREC);
       bs = tobf(bs, DEFAULTPREC);
       //addnum(0,as,bs,&sup);
       addbf(as,bs,&sup);
   
       istoitv(inf,sup,(Itv *)&c);
       mpfr_roundmode = current_roundmode;
       //MKIntervalBigFloat((BF)inf, (BF)sup, c);
       *rp = c;
 #if 0  #if 0
 printexpr(CO, ai);  printexpr(CO, ai);
 printexpr(CO, as);  printexpr(CO, as);
 printexpr(CO, bi);  printexpr(CO, bi);
 printexpr(CO, bs);  printexpr(CO, bs);
 #endif  #endif
   
     addnum(0,ai,bi,&inf);  
     addnum(0,as,bs,&sup);  
     istoitv(inf,sup,(Itv *)&tmp);  
     addulp((IntervalBigFloat)tmp, c);  
     return;      return;
   }    }
 }  }
   
 void subitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)  void subitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *rp)
 {  {
   Num  ai,as,bi,bs,mas, mbs;    Num  ai,as,bi,bs;
   Num  inf,sup,tmp;    IntervalBigFloat c;
   GEN  pa, pb, z;  //,mas, mbs;
   long  ltop, lbot;    Num  inf,sup;
   //,tmp;
     //GEN  pa, pb, z;
     //long  ltop, lbot;
     int current_roundmode;
   
   if ( !a )    if ( !a )
     chsgnnum((Num)b,(Num *)c);      chsgnnum((Num)b,(Num *)rp);
   else if ( !b )    else if ( !b )
     *c = a;      *rp = a;
   else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )    else if ( (NID(a) < N_IntervalBigFloat) && (NID(b) < N_IntervalBigFloat ) )
     subnum(0,(Num)a,(Num)b,(Num *)c);      subnum(0,(Num)a,(Num)b,(Num *)rp);
   else {    else {
     itvtois((Itv)a,&inf,&sup);      itvtois((Itv)a,&ai,&as);
     ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);      itvtois((Itv)b,&bi,&bs);
     itvtois((Itv)b,&inf,&sup);  
     ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);      current_roundmode = mpfr_roundmode;
     subnum(0,ai,bs,&inf);  
     subnum(0,as,bi,&sup);      mpfr_roundmode = MPFR_RNDD;
     istoitv(inf,sup,(Itv *)&tmp);      ai = tobf(ai, DEFAULTPREC);
     addulp((IntervalBigFloat)tmp, c);      bi = tobf(bi, DEFAULTPREC);
   
       mpfr_roundmode = MPFR_RNDU;
       as = tobf(as, DEFAULTPREC);
       bs = tobf(bs, DEFAULTPREC);
       //subnum(0,as,bi,&sup);
       subbf(as,bi,&sup);
   
       mpfr_roundmode = MPFR_RNDD;
       //subnum(0,ai,bs,&inf);
       subbf(ai,bs,&inf);
   
       istoitv(inf,sup,(Itv *)&c);
       mpfr_roundmode = current_roundmode;
       //MKIntervalBigFloat((BF)inf, (BF)sup, c);
       *rp = c;
   
       //addulp((IntervalBigFloat)tmp, c);
   }    }
 }  }
   
Line 173  void mulitvf(IntervalBigFloat a, IntervalBigFloat b, I
Line 221  void mulitvf(IntervalBigFloat a, IntervalBigFloat b, I
   Num  ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t,tmp;    Num  ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t,tmp;
   Num  inf, sup;    Num  inf, sup;
   int  ba,bb;    int  ba,bb;
     int current_roundmode;
   
   if ( !a || !b )    if ( !a || !b )
     *c = 0;      *c = 0;
   else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )    else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
     mulnum(0,(Num)a,(Num)b,(Num *)c);      mulnum(0,(Num)a,(Num)b,(Num *)c);
   else {    else {
     itvtois((Itv)a,&inf,&sup);      itvtois((Itv)a,&ai,&as);
     ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);      itvtois((Itv)b,&bi,&bs);
     itvtois((Itv)b,&inf,&sup);  
     ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);  
   
       current_roundmode = mpfr_roundmode;
   
       mpfr_roundmode = MPFR_RNDU;
       as = tobf(as, DEFAULTPREC);
       bs = tobf(bs, DEFAULTPREC);
   
       mpfr_roundmode = MPFR_RNDD;
       ai = tobf(ai, DEFAULTPREC);
       bi = tobf(bi, DEFAULTPREC);
   
   //    itvtois((Itv)a,&inf,&sup);
   //    ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
   //    itvtois((Itv)b,&inf,&sup);
   //    ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
   
   /* MUST check if ai, as, bi, bs are bigfloat. */    /* MUST check if ai, as, bi, bs are bigfloat. */
     chsgnnum(as,&a2);      chsgnnum(as,&a2);
     chsgnnum(bs,&b2);      chsgnnum(bs,&b2);
Line 224  void mulitvf(IntervalBigFloat a, IntervalBigFloat b, I
Line 286  void mulitvf(IntervalBigFloat a, IntervalBigFloat b, I
     }      }
     if ( ba == bb ) {      if ( ba == bb ) {
       subnum(0,0,c2,&t);        subnum(0,0,c2,&t);
       istoitv(c1,t,(Itv *)&tmp);        istoitv(c1,t,(Itv *)c);
     } else {      } else {
       subnum(0,0,c1,&t);        subnum(0,0,c1,&t);
       istoitv(c2,t,(Itv *)&tmp);        istoitv(c2,t,(Itv *)c);
     }      }
     addulp((IntervalBigFloat)tmp, c);      //addulp((IntervalBigFloat)tmp, c);
   }    }
       mpfr_roundmode = current_roundmode;
 }  }
   
 int     initvf(Num a, Itv b)  int     initvf(Num a, Itv b)
Line 258  void divitvf(IntervalBigFloat a, IntervalBigFloat b, I
Line 321  void divitvf(IntervalBigFloat a, IntervalBigFloat b, I
   Num  ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t,tmp;    Num  ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t,tmp;
   Num  inf, sup;    Num  inf, sup;
   int  ba,bb;    int  ba,bb;
     int current_roundmode;
   
   if ( !b )    if ( !b )
     error("divitv : division by 0");      error("divitv : division by 0");
Line 266  void divitvf(IntervalBigFloat a, IntervalBigFloat b, I
Line 330  void divitvf(IntervalBigFloat a, IntervalBigFloat b, I
   else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )    else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
     divnum(0,(Num)a,(Num)b,(Num *)c);      divnum(0,(Num)a,(Num)b,(Num *)c);
   else {    else {
     itvtois((Itv)a,&inf,&sup);      itvtois((Itv)a,&ai,&as);
     ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);      itvtois((Itv)b,&bi,&bs);
     itvtois((Itv)b,&inf,&sup);  
     ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);      current_roundmode = mpfr_roundmode;
   
       mpfr_roundmode = MPFR_RNDU;
       as = tobf(as, DEFAULTPREC);
       bs = tobf(bs, DEFAULTPREC);
   
       mpfr_roundmode = MPFR_RNDD;
       ai = tobf(ai, DEFAULTPREC);
       bi = tobf(bi, DEFAULTPREC);
   
   //    itvtois((Itv)a,&inf,&sup);
   //    ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
   //    itvtois((Itv)b,&inf,&sup);
   //    ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
 /* MUST check if ai, as, bi, bs are bigfloat. */  /* MUST check if ai, as, bi, bs are bigfloat. */
     chsgnnum(as,&a2);      chsgnnum(as,&a2);
     chsgnnum(bs,&b2);      chsgnnum(bs,&b2);
Line 298  void divitvf(IntervalBigFloat a, IntervalBigFloat b, I
Line 375  void divitvf(IntervalBigFloat a, IntervalBigFloat b, I
     }      }
     if ( ba == bb ) {      if ( ba == bb ) {
       subnum(0,0,c2,&t);        subnum(0,0,c2,&t);
       istoitv(c1,t,(Itv *)&tmp);        istoitv(c1,t,(Itv *)c);
     } else {      } else {
       subnum(0,0,c1,&t);        subnum(0,0,c1,&t);
       istoitv(c2,t,(Itv *)&tmp);        istoitv(c2,t,(Itv *)c);
     }      }
     addulp((IntervalBigFloat)tmp, c);      //addulp((IntervalBigFloat)tmp, c);
   }    }
       mpfr_roundmode = current_roundmode;
 }  }
   
 void pwritvf(Itv a, Num e, Itv *c)  void pwritvf(Itv a, Num e, Itv *c)
Line 322  void pwritvf(Itv a, Num e, Itv *c)
Line 400  void pwritvf(Itv a, Num e, Itv *c)
 #if defined(PARI) && 0  #if defined(PARI) && 0
     gpui_ri((Obj)a,(Obj)c,(Obj *)c);      gpui_ri((Obj)a,(Obj)c,(Obj *)c);
 #else  #else
     error("pwritv : can't calculate a fractional power");      error("pwritvf() : can't calculate a fractional power");
 #endif  #endif
   } else {    } else {
     ei = ZTOS((Q)e);      ei = QTOS((Q)e);
       if (ei<0) ei = -ei;
     pwritv0f(a,ei,&t);      pwritv0f(a,ei,&t);
     if ( SGN((Q)e) < 0 )      if ( SGN((Q)e) < 0 )
       divbf((Num)ONE,(Num)t,(Num *)c);          divitvf((IntervalBigFloat)ONE,(IntervalBigFloat)t,(IntervalBigFloat *)c); /* bug ?? */
     else      else
       *c = t;        *c = t;
   }    }
Line 340  void pwritv0f(Itv a, int e, Itv *c)
Line 419  void pwritv0f(Itv a, int e, Itv *c)
   Num ai,Xmin,Xmax;    Num ai,Xmin,Xmax;
   IntervalBigFloat tmp;    IntervalBigFloat tmp;
   Q  ne;    Q  ne;
     int current_roundmode;
   
   if ( e == 1 )    if ( e == 1 )
     *c = a;      *c = a;
   else {    else {
     STOZ(e,ne);      STOQ(e,ne);
     if ( !(e%2) ) {      if ( !(e%2) ) {
       if ( initvp(0,a) ) {        if ( initvp(0,a) ) {
         Xmin = 0;          Xmin = 0;
Line 367  void pwritv0f(Itv a, int e, Itv *c)
Line 447  void pwritv0f(Itv a, int e, Itv *c)
       Xmin = INF(a);        Xmin = INF(a);
       Xmax = SUP(a);        Xmax = SUP(a);
     }      }
   
       current_roundmode = mpfr_roundmode;
     if ( ! Xmin )  inf = 0;      if ( ! Xmin )  inf = 0;
     else    pwrbf(Xmin,(Num)ne,&inf);      else {
         mpfr_roundmode = MPFR_RNDD;
         pwrbf(Xmin,(Num)ne,&inf);
       }
     if ( ! Xmax )   sup = 0;      if ( ! Xmax )   sup = 0;
     else            pwrbf(Xmax,(Num)ne,&sup);      else {
         mpfr_roundmode = MPFR_RNDU;
         pwrbf(Xmax,(Num)ne,&sup);
       }
     istoitv(inf,sup,(Itv *)&tmp);      istoitv(inf,sup,(Itv *)&tmp);
     addulp(tmp, (IntervalBigFloat *)c);      mpfr_roundmode = current_roundmode;
       *c = (Itv)tmp;
    //   addulp(tmp, (IntervalBigFloat *)c);
   }    }
 }  }
   
 void chsgnitvf(Itv a, Itv *c)  void chsgnitvf(IntervalBigFloat a, IntervalBigFloat *rp)
 {  {
   Num inf,sup;    Num inf,sup;
     IntervalBigFloat c;
     int current_roundmode;
   
   if ( !a )    if ( !a )
     *c = 0;      *rp = 0;
   else if ( NID(a) <= N_B )    else if ( NID(a) < N_IntervalBigFloat )
     chsgnnum((Num)a,(Num *)c);      chsgnnum((Num)a,(Num *)rp);
   else {    else {
     chsgnnum(INF((Itv)a),&inf);      current_roundmode = mpfr_roundmode;
     chsgnnum(SUP((Itv)a),&sup);  
     istoitv(inf,sup,c);      mpfr_roundmode = MPFR_RNDU;
       chsgnnum((Num)INF(a),&sup);
       mpfr_roundmode = MPFR_RNDD;
       chsgnnum((Num)SUP(a),&inf);
       //MKIntervalBigFloat((BF)inf,(BF)sup,c);
       istoitv(inf,sup,(Itv *)&c);
       *rp = c;
       mpfr_roundmode = current_roundmode;
   }    }
 }  }
   
Line 434  void miditvf(Itv a, Num *b)
Line 533  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 {
     STOZ(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.2  
changed lines
  Added in v.1.3
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>