[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.3 and 1.10

version 1.3, 2003/02/14 22:29:08 version 1.10, 2019/11/12 10:52:04
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2000/engine/f-itv.c,v 1.2 2002/01/08 04:14:37 kondoh 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"
 extern long prec;  long get_pariprec();
 #endif  #endif
   
 void ToBf(Num a, BF *rp)  void ToBf(Num a, BF *rp)
 {  {
         GEN     pa, pb, pc;    GEN  pa, pb, pc;
         BF      bn, bd;    BF  bn, bd;
         BF      c;    BF  c;
         long    ltop, lbot;    long  ltop, lbot;
   
         if ( ! a ) {    if ( ! a ) {
                 *rp = 0;      *rp = 0;
         }    }
         else switch ( NID(a) ) {    else switch ( NID(a) ) {
                 case N_Q:      case N_Q:
                         ltop = avma;        ltop = avma;
                         ritopa_i(NM((Q)a), SGN((Q)a), &pa);        ritopa_i(NM((Q)a), SGN((Q)a), &pa);
                         pb = cgetr(prec);        pb = cgetr(get_pariprec());
                         mpaff(pa, pb);        mpaff(pa, pb);
                         if ( INT((Q)a) ) {        if ( INT((Q)a) ) {
                                 lbot = avma;          lbot = avma;
                                 pb = gerepile(ltop, lbot, pb);          pb = gerepile(ltop, lbot, pb);
                                 patori(pb, &c);          patori(pb, &c);
                                 cgiv(pb);          cgiv(pb);
                                 *rp = c;          *rp = c;
                         } else {        } else {
                                 patori(pb, &bn);          patori(pb, &bn);
                                 ritopa_i(DN((Q)a), 1, &pa);          ritopa_i(DN((Q)a), 1, &pa);
                                 pb = cgetr(prec);          pb = cgetr(get_pariprec());
                                 mpaff(pa, pb);          mpaff(pa, pb);
                                 lbot = avma;          lbot = avma;
                                 pb = gerepile(ltop, lbot, pb);          pb = gerepile(ltop, lbot, pb);
                                 patori(pb, &bd);          patori(pb, &bd);
                                 cgiv(pb);          cgiv(pb);
                                 divbf((Num)bn,(Num)bd, (Num *)&c);          divbf((Num)bn,(Num)bd, (Num *)&c);
                                 *rp = c;          *rp = c;
                         }        }
                         break;        break;
                 case N_R:      case N_R:
                         pb = dbltor(BDY((Real)a));        pb = dbltor(BDY((Real)a));
                         patori(pb, &c);        patori(pb, &c);
                         cgiv(pb);        cgiv(pb);
                         *rp = c;        *rp = c;
                         break;        break;
                 case N_B:      case N_B:
                         *rp = (BF)a;        *rp = (BF)a;
                         break;        break;
                 default:      default:
                         error("ToBf : not implemented");        error("ToBf : not implemented");
                         break;        break;
         }    }
 }  }
   
 void double2bf(double d, BF *rp)  void double2bf(double d, BF *rp)
 {  {
         GEN     p;    GEN  p;
   
         p = dbltor(d);    p = dbltor(d);
         patori(p, rp);    patori(p, rp);
         cgiv(p);    cgiv(p);
 }  }
   
 double  bf2double(BF a)  double  bf2double(BF a)
 {  {
         GEN     p;    GEN  p;
   
         ritopa(a, &p);    ritopa(a, &p);
         return rtodbl(p);    return rtodbl(p);
 }  }
   
 void getulp(BF a, BF *au)  void getulp(BF a, BF *au)
 {  {
         GEN     b, c;    GEN  b, c;
         long    lb, i;    long  lb, i;
   
         ritopa(a, &b);    ritopa(a, &b);
         lb = lg(b);    lb = lg(b);
         c = cgetr(lb);    c = cgetr(lb);
         c[1] = b[1];    c[1] = b[1];
         for (i=2;i<lb-1;i++) c[i] = 0;    for (i=2;i<lb-1;i++) c[i] = 0;
         c[i] = 1;    c[i] = 1;
         setsigne(c, 1);    setsigne(c, 1);
         patori(c,au);    patori(c,au);
         cgiv(c);    cgiv(c);
         cgiv(b);    cgiv(b);
 }  }
   
 void addulp(IntervalBigFloat a, IntervalBigFloat *c)  void addulp(IntervalBigFloat a, IntervalBigFloat *c)
 {  {
         Num     ai, as, aiu, asu, inf, sup;    Num  ai, as, aiu, asu, inf, sup;
   
         itvtois((Itv)a,&ai,&as);    itvtois((Itv)a,&ai,&as);
         if ( ai ) {    if ( ai ) {
                 getulp((BF)ai, (BF *)&aiu);      getulp((BF)ai, (BF *)&aiu);
                 subbf(ai,aiu,&inf);      subbf(ai,aiu,&inf);
         } else  inf = ai;    } else  inf = ai;
         if ( as ) {    if ( as ) {
                 getulp((BF)as, (BF *)&asu);      getulp((BF)as, (BF *)&asu);
                 addbf(as,asu,&sup);      addbf(as,asu,&sup);
         } 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);
         Num     inf,sup;  }
         GEN     pa, pb, z;  
         long    ltop, lbot;  
   
         if ( !a )  double mpfr2dblUp(mpfr_t a)
                 *c = b;  {
         else if ( !b )          return mpfr_get_d(a,MPFR_RNDU);
                 *c = a;  }
         else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )  
                 addnum(0,(Num)a,(Num)b,(Num *)c);  
         else {  void toInterval(Num a, int prec, int type, Num *rp)
                 itvtois((Itv)a,&inf,&sup);  {
                 ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);          if ( ! a ) {
                 itvtois((Itv)b,&inf,&sup);                  *rp = 0;
                 ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);          } 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;
     //GEN  pa, pb, z;
     //long  ltop, lbot;
     int current_roundmode;
   
     if ( !a )
       *rp = b;
     else if ( !b )
       *rp = a;
     else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
       addnum(0,(Num)a,(Num)b,(Num *)rp);
     else {
       itvtois((Itv)a,&ai,&as);
       itvtois((Itv)b,&bi,&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
       return;
     }
   }
   
 #if 1  void subitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *rp)
                 addnum(0,ai,bi,&inf);  {
                 addnum(0,as,bs,&sup);    Num  ai,as,bi,bs;
                 istoitv(inf,sup,(Itv *)&tmp);    IntervalBigFloat c;
                 addulp((IntervalBigFloat)tmp, c);  //,mas, mbs;
                 return;    Num  inf,sup;
 #else  //,tmp;
                 ltop = avma;    //GEN  pa, pb, z;
                 ritopa(ai,&pa);    //long  ltop, lbot;
                 ritopa(bi,&pb);    int current_roundmode;
                 lbot = avma;  
                 z = gerepile(ltop,lbot,PariAddDown(pa,pb));  
                 patori(z,&inf); cgiv(z);  
   
         /* MUST check if ai, as, bi, bs are bigfloat. */    if ( !a )
       chsgnnum((Num)b,(Num *)rp);
     else if ( !b )
       *rp = a;
     else if ( (NID(a) < N_IntervalBigFloat) && (NID(b) < N_IntervalBigFloat ) )
       subnum(0,(Num)a,(Num)b,(Num *)rp);
     else {
       itvtois((Itv)a,&ai,&as);
       itvtois((Itv)b,&bi,&bs);
   
                 /*  as + bs = ( - ( (-as) + (-bs) ) ) */      current_roundmode = mpfr_roundmode;
                 chsgnbf(as,&mas);  
                 chsgnbf(bs,&mbs);  
                 ltop = avma;  
                 ritopa(mas,&pa);  
                 ritopa(mbs,&pb);  
                 lbot = avma;  
                 z = gerepile(ltop,lbot,PariAddDown(pa,pb));  
                 patori(z,&tmp); cgiv(z);  
   
                 chsgnbf(tmp,&sup);      mpfr_roundmode = MPFR_RNDD;
                 istoitv(inf,sup,c);      ai = tobf(ai, DEFAULTPREC);
 #endif      bi = tobf(bi, DEFAULTPREC);
         }  
 }  
   
 void subitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)      mpfr_roundmode = MPFR_RNDU;
 {      as = tobf(as, DEFAULTPREC);
         Num     ai,as,bi,bs,mas, mbs;      bs = tobf(bs, DEFAULTPREC);
         Num     inf,sup,tmp;      //subnum(0,as,bi,&sup);
         GEN     pa, pb, z;      subbf(as,bi,&sup);
         long    ltop, lbot;  
   
         if ( !a )      mpfr_roundmode = MPFR_RNDD;
                 chsgnnum((Num)b,(Num *)c);      //subnum(0,ai,bs,&inf);
         else if ( !b )      subbf(ai,bs,&inf);
                 *c = a;  
         else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )  
                 subnum(0,(Num)a,(Num)b,(Num *)c);  
         else {  
                 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);  
 #if 1  
                 subnum(0,ai,bs,&inf);  
                 subnum(0,as,bi,&sup);  
                 istoitv(inf,sup,(Itv *)&tmp);  
                 addulp((IntervalBigFloat)tmp, c);  
 #else  
   
 /* MUST check if ai, as, bi, bs are bigfloat. */      istoitv(inf,sup,(Itv *)&c);
                 /* ai - bs = ai + (-bs)  */      mpfr_roundmode = current_roundmode;
                 chsgnbf(bs,&mbs);      //MKIntervalBigFloat((BF)inf, (BF)sup, c);
                 ltop = avma;      *rp = c;
                 ritopa(ai,&pa);  
                 ritopa(mbs,&pb);  
                 lbot = avma;  
                 z = gerepile(ltop,lbot,PariAddDown(pa,pb));  
                 patori(z,&inf); cgiv(z);  
   
                 /* as - bi = ( - ( bi + (-as) ) ) */      //addulp((IntervalBigFloat)tmp, c);
                 chsgnbf(as,&mas);    }
                 ltop = avma;  
                 ritopa(mas,&pa);  
                 ritopa(bi,&pb);  
                 lbot = avma;  
                 z = gerepile(ltop,lbot,PariAddDown(pa,pb));  
                 patori(z,&tmp); cgiv(z);  
   
                 chsgnbf(tmp,&sup);  
                 istoitv(inf,sup,c);  
 #endif  
         }  
 }  }
   
 void mulitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)  void mulitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
 {  {
         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);  
   
         /* MUST check if ai, as, bi, bs are bigfloat. */      current_roundmode = mpfr_roundmode;
                 chsgnnum(as,&a2);  
                 chsgnnum(bs,&b2);  
   
       mpfr_roundmode = MPFR_RNDU;
       as = tobf(as, DEFAULTPREC);
       bs = tobf(bs, DEFAULTPREC);
   
                 if ( ba = (compnum(0,0,a2) > 0) ) {      mpfr_roundmode = MPFR_RNDD;
                         a1 = ai;      ai = tobf(ai, DEFAULTPREC);
                 } else {      bi = tobf(bi, DEFAULTPREC);
                         a1 = a2;  
                         a2 = ai;  //    itvtois((Itv)a,&inf,&sup);
                 }  //    ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
                 if ( bb = (compnum(0,0,b2) > 0) ) {  //    itvtois((Itv)b,&inf,&sup);
                         b1 = bi;  //    ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
                 } else {  
                         b1 = b2;    /* MUST check if ai, as, bi, bs are bigfloat. */
                         b2 = bi;      chsgnnum(as,&a2);
                 }      chsgnnum(bs,&b2);
                 mulnum(0,a2,b2,&t);  
                 subnum(0,0,t,&c2);  
                 if ( compnum(0,0,b1) > 0 ) {      if ( ba = (compnum(0,0,a2) > 0) ) {
                         mulnum(0,a2,b1,&t);        a1 = ai;
                         subnum(0,0,t,&c1);      } else {
                         if ( compnum(0,0,a1) > 0 ) {        a1 = a2;
                                 mulnum(0,a1,b2,&t);        a2 = ai;
                                 subnum(0,0,t,&p);      }
                                 if ( compnum(0,c1,p) > 0 ) c1 = p;      if ( bb = (compnum(0,0,b2) > 0) ) {
                                 mulnum(0,a1,b1,&t);        b1 = bi;
                                 subnum(0,0,t,&p);      } else {
                                 if ( compnum(0,c2,p) > 0 ) c2 = p;        b1 = b2;
                         }        b2 = bi;
                 } else {      }
                         if ( compnum(0,0,a1) > 0 ) {      mulnum(0,a2,b2,&t);
                                 mulnum(0,a1,b2,&t);      subnum(0,0,t,&c2);
                                 subnum(0,0,t,&c1);      if ( compnum(0,0,b1) > 0 ) {
                         } else {        mulnum(0,a2,b1,&t);
                                 mulnum(0,a1,b1,&c1);        subnum(0,0,t,&c1);
                         }        if ( compnum(0,0,a1) > 0 ) {
                 }          mulnum(0,a1,b2,&t);
                 if ( ba == bb ) {          subnum(0,0,t,&p);
                         subnum(0,0,c2,&t);          if ( compnum(0,c1,p) > 0 ) c1 = p;
                         istoitv(c1,t,(Itv *)&tmp);          mulnum(0,a1,b1,&t);
                 } else {          subnum(0,0,t,&p);
                         subnum(0,0,c1,&t);          if ( compnum(0,c2,p) > 0 ) c2 = p;
                         istoitv(c2,t,(Itv *)&tmp);        }
                 }      } else {
                 addulp((IntervalBigFloat)tmp, c);        if ( compnum(0,0,a1) > 0 ) {
         }          mulnum(0,a1,b2,&t);
           subnum(0,0,t,&c1);
         } else {
           mulnum(0,a1,b1,&c1);
         }
       }
       if ( ba == bb ) {
         subnum(0,0,c2,&t);
         istoitv(c1,t,(Itv *)c);
       } else {
         subnum(0,0,c1,&t);
         istoitv(c2,t,(Itv *)c);
       }
       //addulp((IntervalBigFloat)tmp, c);
     }
       mpfr_roundmode = current_roundmode;
 }  }
   
 int     initvf(Num a, Itv b)  int     initvf(Num a, Itv b)
 {  {
         Num inf, sup;    Num inf, sup;
   
         itvtois(b, &inf, &sup);    itvtois(b, &inf, &sup);
         if ( compnum(0,inf,a) > 0 ) return 0;    if ( compnum(0,inf,a) > 0 ) return 0;
         else if ( compnum(0,a,sup) > 0 ) return 0;    else if ( compnum(0,a,sup) > 0 ) return 0;
         else return 1;    else return 1;
 }  }
   
 int     itvinitvf(Itv a, Itv b)  int     itvinitvf(Itv a, Itv b)
 {  {
         Num ai,as,bi,bs;    Num ai,as,bi,bs;
   
         itvtois(a, &ai, &as);    itvtois(a, &ai, &as);
         itvtois(b, &bi, &bs);    itvtois(b, &bi, &bs);
         if ( compnum(0,ai,bi) >= 0  && compnum(0,bs,as) >= 0 ) return 1;    if ( compnum(0,ai,bi) >= 0  && compnum(0,bs,as) >= 0 ) return 1;
         else return 0;    else return 0;
 }  }
   
 void divitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)  void divitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
 {  {
         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");
         else if ( !a )    else if ( !a )
                 *c = 0;      *c = 0;
         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);
                 if ( ba = (compnum(0,0,a2) > 0) ) {      if ( ba = (compnum(0,0,a2) > 0) ) {
                         a1 = ai;        a1 = ai;
                 } else {      } else {
                         a1 = a2;        a1 = a2;
                         a2 = ai;        a2 = ai;
                 }      }
                 if ( bb = (compnum(0,0,b2) > 0) ) {      if ( bb = (compnum(0,0,b2) > 0) ) {
                         b1 = bi;        b1 = bi;
                 } else {      } else {
                         b1 = b2;        b1 = b2;
                         b2 = bi;        b2 = bi;
                 }      }
                 if ( compnum(0,0,b1) >= 0 ) {      if ( compnum(0,0,b1) >= 0 ) {
                         error("divitv : division by interval including 0.");        error("divitv : division by interval including 0.");
                 } else {      } else {
                         divnum(0,a2,b1,&c2);        divnum(0,a2,b1,&c2);
                         if ( compnum(0,0,a1) > 0 ) {        if ( compnum(0,0,a1) > 0 ) {
                                 divnum(0,a1,b1,&c1);          divnum(0,a1,b1,&c1);
                         } else {        } else {
                                 divnum(0,a1,b2,&t);          divnum(0,a1,b2,&t);
                                 subnum(0,0,t,&c1);          subnum(0,0,t,&c1);
                         }        }
                 }      }
                 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)
 {  {
         int ei;    int ei;
         Itv t;    Itv t;
   
         if ( !e )    if ( !e )
                 *c = (Itv)ONE;      *c = (Itv)ONE;
         else if ( !a )    else if ( !a )
                 *c = 0;      *c = 0;
         else if ( NID(a) <= N_B )    else if ( NID(a) <= N_B )
                 pwrnum(0,(Num)a,e,(Num *)c);      pwrnum(0,(Num)a,e,(Num *)c);
         else if ( !INT(e) ) {    else if ( !INT(e) ) {
 #if defined(PARI) && 0  #if defined(PARI) && 0
                 GEN pa,pe,z;      gpui_ri((Obj)a,(Obj)c,(Obj *)c);
                 int ltop,lbot;  
   
                 ltop = avma; ritopa(a,&pa); ritopa(e,&pe); lbot = avma;  
                 z = gerepile(ltop,lbot,gpui(pa,pe,prec));  
                 patori(z,c); cgiv(z);  
 #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);
                 pwritv0f(a,ei,&t);      if (ei<0) ei = -ei;
                 if ( SGN((Q)e) < 0 )      pwritv0f(a,ei,&t);
                         divbf((Num)ONE,(Num)t,(Num *)c);      if ( SGN((Q)e) < 0 )
                 else          divitvf((IntervalBigFloat)ONE,(IntervalBigFloat)t,(IntervalBigFloat *)c); /* bug ?? */
                         *c = t;      else
         }        *c = t;
     }
 }  }
   
 void pwritv0f(Itv a, int e, Itv *c)  void pwritv0f(Itv a, int e, Itv *c)
 {  {
         Num inf, sup;    Num inf, sup;
         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 {
                 STOQ(e,ne);      STOQ(e,ne);
                 if ( !(e%2) ) {      if ( !(e%2) ) {
                         if ( initvp(0,a) ) {        if ( initvp(0,a) ) {
                                 Xmin = 0;          Xmin = 0;
                                 chsgnnum(INF(a),&ai);          chsgnnum(INF(a),&ai);
                                 if ( compnum(0,ai,SUP(a)) > 0 ) {          if ( compnum(0,ai,SUP(a)) > 0 ) {
                                         Xmax = ai;            Xmax = ai;
                                 } else {          } else {
                                         Xmax = SUP(a);            Xmax = SUP(a);
                                 }          }
                         } else {        } else {
                                 if ( compnum(0,INF(a),0) > 0 ) {          if ( compnum(0,INF(a),0) > 0 ) {
                                         Xmin = INF(a);            Xmin = INF(a);
                                         Xmax = SUP(a);            Xmax = SUP(a);
                                 } else {          } else {
                                         Xmin = SUP(a);            Xmin = SUP(a);
                                         Xmax = INF(a);            Xmax = INF(a);
                                 }          }
                         }        }
                 } else {      } else {
                         Xmin = INF(a);        Xmin = INF(a);
                         Xmax = SUP(a);        Xmax = SUP(a);
                 }      }
                 if ( ! Xmin )   inf = 0;  
                 else            pwrbf(Xmin,(Num)ne,&inf);      current_roundmode = mpfr_roundmode;
                 if ( ! Xmax )   sup = 0;      if ( ! Xmin )  inf = 0;
                 else            pwrbf(Xmax,(Num)ne,&sup);      else {
                 istoitv(inf,sup,(Itv *)&tmp);        mpfr_roundmode = MPFR_RNDD;
                 addulp(tmp, (IntervalBigFloat *)c);        pwrbf(Xmin,(Num)ne,&inf);
         }      }
       if ( ! Xmax )   sup = 0;
       else {
         mpfr_roundmode = MPFR_RNDU;
         pwrbf(Xmax,(Num)ne,&sup);
       }
       istoitv(inf,sup,(Itv *)&tmp);
       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;
     }
 }  }
   
 int cmpitvf(Itv a, Itv b)  int cmpitvf(Itv a, Itv b)
 {  {
         Num     ai,as,bi,bs;    Num  ai,as,bi,bs;
         int     s,t;    int  s,t;
   
         if ( !a ) {    if ( !a ) {
                 if ( !b || (NID(b)<=N_B) )      if ( !b || (NID(b)<=N_B) ) {
                         return compnum(0,(Num)a,(Num)b);        return compnum(0,(Num)a,(Num)b);
                 else      } else {
                         return -1;        itvtois(b,&bi,&bs);
         } else if ( !b ) {        if ( compnum(0,(Num)a,bs) > 0 ) return 1;
                 if ( !a || (NID(a)<=N_B) )        else if ( compnum(0,bi,(Num)a) > 0 ) return -1;
                         return compnum(0,(Num)a,(Num)b);        else  return 0;
                 else      }
                         return initvp((Num)b,a);    } else if ( !b ) {
         } else {      if ( !a || (NID(a)<=N_B) ) {
                 itvtois(a,&ai,&as);        return compnum(0,(Num)a,(Num)b);
                 itvtois(b,&bi,&bs);      } else {
                 s = compnum(0,ai,bs) ;        itvtois(a,&ai,&as);
                 t = compnum(0,bi,as) ;        if ( compnum(0,ai,(Num)b) > 0 ) return 1;
                 if ( s > 0 ) return 1;        else if ( compnum(0,(Num)b,as) > 0 ) return -1;
                 else if ( t > 0 ) return -1;        else  return 0;
                 else  return 0;      }
         }    } else {
       itvtois(a,&ai,&as);
       itvtois(b,&bi,&bs);
       s = compnum(0,ai,bs) ;
       t = compnum(0,bi,as) ;
       if ( s > 0 ) return 1;
       else if ( t > 0 ) return -1;
       else  return 0;
     }
 }  }
   
 void miditvf(Itv a, Num *b)  void miditvf(Itv a, Num *b)
 {  {
         Num     ai,as;    Num  ai,as;
         Num     t;    Num  t;
         Q       TWO;  
   
         if ( ! a ) *b = 0;    if ( ! a ) *b = 0;
         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);
         }    }
 }  }
   
 void cupitvf(Itv a, Itv b, Itv *c)  void cupitvf(Itv a, Itv b, Itv *c)
 {  {
         Num     ai,as,bi,bs;    Num  ai,as,bi,bs;
         Num     inf,sup;    Num  inf,sup;
   
         itvtois(a,&ai,&as);    itvtois(a,&ai,&as);
         itvtois(b,&bi,&bs);    itvtois(b,&bi,&bs);
         if ( compnum(0,ai,bi) > 0 )     inf = bi;    if ( compnum(0,ai,bi) > 0 )  inf = bi;
         else                            inf = ai;    else        inf = ai;
         if ( compnum(0,as,bs) > 0 )     sup = as;    if ( compnum(0,as,bs) > 0 )  sup = as;
         else                            sup = bs;    else        sup = bs;
         istoitv(inf,sup,c);    istoitv(inf,sup,c);
 }  }
   
 void capitvf(Itv a, Itv b, Itv *c)  void capitvf(Itv a, Itv b, Itv *c)
 {  {
         Num     ai,as,bi,bs;    Num  ai,as,bi,bs;
         Num     inf,sup;    Num  inf,sup;
   
         itvtois(a,&ai,&as);    itvtois(a,&ai,&as);
         itvtois(b,&bi,&bs);    itvtois(b,&bi,&bs);
         if ( compnum(0,ai,bi) > 0 )     inf = ai;    if ( compnum(0,ai,bi) > 0 )  inf = ai;
         else                            inf = bi;    else        inf = bi;
         if ( compnum(0,as,bs) > 0 )     sup = bs;    if ( compnum(0,as,bs) > 0 )  sup = bs;
         else                            sup = as;    else        sup = as;
         if ( compnum(0,inf,sup) > 0 )   *c = 0;    if ( compnum(0,inf,sup) > 0 )  *c = 0;
         else                            istoitv(inf,sup,c);    else        istoitv(inf,sup,c);
 }  }
   
   
 void widthitvf(Itv a, Num *b)  void widthitvf(Itv a, Num *b)
 {  {
         Num     ai,as;    Num  ai,as;
   
         if ( ! a ) *b = 0;    if ( ! a ) *b = 0;
         else if ( (NID(a) <= N_B) )    else if ( (NID(a) <= N_B) )
                 *b = (Num)a;      *b = (Num)a;
         else {    else {
                 itvtois(a,&ai,&as);      itvtois(a,&ai,&as);
                 subnum(0,as,ai,b);      subnum(0,as,ai,b);
         }    }
 }  }
   
 void absitvf(Itv a, Num *b)  void absitvf(Itv a, Num *b)
 {  {
         Num     ai,as,bi,bs;    Num  ai,as,bi,bs;
   
         if ( ! a ) *b = 0;    if ( ! a ) *b = 0;
         else if ( (NID(a) <= N_B) ) {    else if ( (NID(a) <= N_B) ) {
                 if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);      if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
                 else *b = (Num)a;      else *b = (Num)a;
         } else {    } else {
                 itvtois(a,&ai,&as);      itvtois(a,&ai,&as);
                 if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);      if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
                 else bi = ai;      else bi = ai;
                 if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);      if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
                 else bs = as;      else bs = as;
                 if ( compnum(0,bi,bs) > 0 )     *b = bi;      if ( compnum(0,bi,bs) > 0 )  *b = bi;
                 else                            *b = bs;      else        *b = bs;
         }    }
 }  }
   
 void distanceitvf(Itv a, Itv b, Num *c)  void distanceitvf(Itv a, Itv b, Num *c)
 {  {
         Num     ai,as,bi,bs;    Num  ai,as,bi,bs;
         Num     di,ds;    Num  di,ds;
         Itv     d;    Itv  d;
   
         itvtois(a,&ai,&as);    itvtois(a,&ai,&as);
         itvtois(b,&bi,&bs);    itvtois(b,&bi,&bs);
         subnum(0,ai,bi,&di);    subnum(0,ai,bi,&di);
         subnum(0,as,bs,&ds);    subnum(0,as,bs,&ds);
         istoitv(di,ds,&d);    istoitv(di,ds,&d);
         absitvf(d,c);    absitvf(d,c);
 }  }
   
 #endif  #endif
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.10
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>