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

Diff for /OpenXM_contrib2/asir2000/engine/f-itv.c between version 1.8 and 1.10

version 1.8, 2018/03/29 01:32:51 version 1.10, 2019/11/12 10:52:04
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2000/engine/f-itv.c,v 1.7 2009/03/27 14:42:29 ohara Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/f-itv.c,v 1.9 2019/06/04 07:11:22 kondoh 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)
   
   double mpfr2dblDown(mpfr_t a)
 {  {
   Num  ai,as,bi,bs,mas,mbs,tmp;          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)
   {
     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 329  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 394  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 429  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 438  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 483  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 508  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 = QTOS((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 527  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;
Line 367  void pwritv0f(Itv a, int e, Itv *c)
Line 555  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;
   }    }
 }  }
   
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.10
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>