[BACK]Return to itvnum.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Diff for /OpenXM_contrib2/asir2000/builtin/itvnum.c between version 1.10 and 1.13

version 1.10, 2016/06/29 08:16:11 version 1.13, 2019/11/12 10:52:04
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2000/builtin/itvnum.c,v 1.9 2015/08/14 13:51:54 fujimoto Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/itvnum.c,v 1.12 2019/06/04 07:11:22 kondoh Exp $
  */   */
   
 #include "ca.h"  #include "ca.h"
Line 9 
Line 9 
 #include "../plot/ifplot.h"  #include "../plot/ifplot.h"
 #endif  #endif
   
 #if defined(INTERVAL)  #include <mpfr.h>
   #include <mpfi.h>
   // in engine/bf.c
   Num tobf(Num,int);
   
   #if defined(INTERVAL)
 static void Pitv(NODE, Obj *);  static void Pitv(NODE, Obj *);
 static void Pitvd(NODE, Obj *);  static void Pitvd(NODE, Obj *);
 static void Pitvbf(NODE, Obj *);  static void Pitvbf(NODE, Obj *);
Line 24  static void Pcup(NODE, Obj *);
Line 28  static void Pcup(NODE, Obj *);
 static void Pcap(NODE, Obj *);  static void Pcap(NODE, Obj *);
 static void Pwidth(NODE, Obj *);  static void Pwidth(NODE, Obj *);
 static void Pdistance(NODE, Obj *);  static void Pdistance(NODE, Obj *);
 static void Pitvversion(Q *);  static void Pitvversion(NODE, Q *);
 void miditvp(Itv,Num *);  static void PzeroRewriteMode(NODE, Obj *);
 void absitvp(Itv,Num *);  static void PzeroRewriteCountClear(NODE, Obj *);
 int initvd(Num,IntervalDouble);  static void PzeroRewriteCount(NODE, Obj *);
 int initvp(Num,Itv);  //void miditvp(Itv,Num *);
 int itvinitvp(Itv,Itv);  //void absitvp(Itv,Num *);
   //int initvd(Num,IntervalDouble);
   //int initvp(Num,Itv);
   //int itvinitvp(Itv,Itv);
   static void Pevalitv(NODE, Obj *);
   static void Pevalitvd(NODE, Obj *);
   void Ppi_itvd(NODE, Obj *);
   void Pe_itvd(NODE, Obj *);
   void Psinitv(NODE, Obj *);
   void Psinitvd(NODE, Obj *);
 #endif  #endif
 static void Pprintmode(NODE, Obj *);  static void Pprintmode(NODE, Obj *);
   
Line 37  static void Pprintmode(NODE, Obj *);
Line 50  static void Pprintmode(NODE, Obj *);
 static void ccalc(double **, struct canvas *, int);  static void ccalc(double **, struct canvas *, int);
 static void Pifcheck(NODE, Obj *);  static void Pifcheck(NODE, Obj *);
   
 #if     defined(__osf__) && 0  #if  defined(__osf__) && 0
 int     end;  int  end;
 #endif  #endif
   
 struct ftab interval_tab[] = {  struct ftab interval_tab[] = {
         {"printmode",Pprintmode,1},    {"printmode",Pprintmode,1},
 #if defined(INTERVAL)  #if defined(INTERVAL)
         {"itvd",Pitvd,-2},    {"itvd",Pitvd,-2},
         {"intvald",Pitvd,-2},    {"intvald",Pitvd,-2},
         {"itv",Pitv,-2},    {"itv",Pitv,-2},
         {"intval",Pitv,-2},    {"intval",Pitv,-2},
         {"itvbf",Pitvbf,-2},    {"itvbf",Pitvbf,-2},
         {"intvalbf",Pitvbf,-2},    {"intvalbf",Pitvbf,-2},
         {"inf",Pinf,1},    {"inf",Pinf,1},
         {"sup",Psup,1},    {"sup",Psup,1},
         {"absintval",Pabsitv,1},    {"absintval",Pabsitv,1},
         {"disintval",Pdisjitv,2},    {"disintval",Pdisjitv,2},
         {"inintval",Pinitv,2},    {"inintval",Pinitv,2},
         {"cup",Pcup,2},    {"cup",Pcup,2},
         {"cap",Pcap,2},    {"cap",Pcap,2},
         {"mid",Pmid,1},    {"mid",Pmid,1},
         {"width",Pwidth,1},    {"width",Pwidth,1},
         {"diam",Pwidth,1},    {"diam",Pwidth,1},
         {"distance",Pdistance,2},    {"distance",Pdistance,2},
         {"iversion",Pitvversion,0},    {"iversion",Pitvversion,-1},
     {"intvalversion",Pitvversion,-1},
     {"zerorewritemode",PzeroRewriteMode,-1},
     {"zeroRewriteMode",PzeroRewriteMode,-1},
     {"zeroRewriteCountClear",PzeroRewriteCountClear,-1},
     {"zeroRewriteCount",PzeroRewriteCount,-1},
   /* eval */
     {"evalitv",   Pevalitv,       -2},
     {"evalitvd",  Pevalitvd,      1},
   /* math */
     {"piitvd",    Pitvbf_pi,      -1},
     {"eitvd",             Pitvbf_e,       -1},
   
     {"piitv",             Pitvbf_pi,      -1},
     {"eitv",              Pitvbf_e,       -1},
   #if 0
     {"factorialitv",Pfactorialitv,1},
     {"factorialitvd",Pfactorialitvd,1},
   #endif
   
     {"absitv",    Pitvbf_abs,     -2},
     {"absitvd",   Pitvbf_abs,     -2},
   
     {"logitv",    Pitvbf_log,     -2},
     {"logitvd",   Pitvbf_log,     -2},
     {"expitv",    Pitvbf_exp,     -2},
     {"expitvd",   Pitvbf_exp,     -2},
     {"powitv",    Pitvbf_pow,     -3},
     {"powitvd",   Pitvbf_pow,     -3},
   
     {"sinitv",    Pitvbf_sin,     -2},
     {"sinitvd",   Pitvd_sin,      -2},
   
     {"cositv",    Pitvbf_cos,     -2},
     {"cositvd",   Pitvd_cos,      -2},
     {"tanitv",    Pitvbf_tan,     -2},
     {"tanitvd",   Pitvd_tan,      -2},
     {"asinitv",   Pitvbf_asin,    -2},
     {"asinitvd",  Pitvd_asin,     -2},
     {"acositv",   Pitvbf_acos,    -2},
     {"acositvd",  Pitvd_acos,     -2},
     {"atanitv",   Pitvbf_atan,    -2},
     {"atanitvd",  Pitvd_atan,     -2},
     {"sinhitv",   Pitvbf_sinh,    -2},
     {"sinhitvd",  Pitvd_sinh,     -2},
     {"coshitv",   Pitvbf_cosh,    -2},
     {"coshitvd",  Pitvd_cosh,     -2},
     {"tanhitv",   Pitvbf_tanh,    -2},
     {"tanhitvd",  Pitvd_tanh,     -2},
     {"asinhitv",  Pitvbf_asinh,   -2},
     {"asinhitvd", Pitvd_asinh,    -2},
     {"acoshitv",  Pitvbf_acosh,   -2},
     {"acoshitvd", Pitvd_acosh,    -2},
     {"atanhitv",  Pitvbf_atanh,   -2},
     {"atanhitvd", Pitvd_atanh,    -2},
   
 /* plot time check */  /* plot time check */
         {"ifcheck",Pifcheck,-7},    {"ifcheck",Pifcheck,-7},
 #endif  #endif
         {0,0,0},    {0,0,0},
 };  };
   
   extern int mpfr_roundmode;
   
 #if defined(INTERVAL)  #if defined(INTERVAL)
   
 /* plot time check */  /* plot time check */
 static void  static void
 Pifcheck(NODE arg, Obj *rp)  Pifcheck(NODE arg, Obj *rp)
 {  {
         Q m2,p2,s_id;    Q m2,p2,s_id;
         NODE defrange;    NODE defrange;
         LIST xrange,yrange,range[2],list,geom;    LIST xrange,yrange,range[2],list,geom;
         VL vl,vl0;    VL vl,vl0;
         V v[2],av[2];    V v[2],av[2];
         int ri,i,j,sign;    int ri,i,j,sign;
         P poly;    P poly;
         P var;    P var;
         NODE n,n0;    NODE n,n0;
         Obj t;    Obj t;
   
         struct canvas *can;    struct canvas *can;
         MAT m;    MAT m;
         pointer **mb;    pointer **mb;
         double **tabe, *px, *px1, *px2;    double **tabe, *px, *px1, *px2;
         Q one;    Q one;
         int width, height, ix, iy;    int width, height, ix, iy;
         int id;    int id;
   
         STOQ(-2,m2); STOQ(2,p2);    STOQ(-2,m2); STOQ(2,p2);
         STOQ(1,one);    STOQ(1,one);
         MKNODE(n,p2,0); MKNODE(defrange,m2,n);    MKNODE(n,p2,0); MKNODE(defrange,m2,n);
         poly = 0; vl = 0; geom = 0; ri = 0;    poly = 0; vl = 0; geom = 0; ri = 0;
         v[0] = v[1] = 0;    v[0] = v[1] = 0;
         for ( ; arg; arg = NEXT(arg) ){    for ( ; arg; arg = NEXT(arg) ){
                 switch ( OID(BDY(arg)) ) {      switch ( OID(BDY(arg)) ) {
                 case O_P:      case O_P:
                         poly = (P)BDY(arg);        poly = (P)BDY(arg);
                         get_vars_recursive((Obj)poly,&vl);        get_vars_recursive((Obj)poly,&vl);
                         for(vl0=vl,i=0;vl0;vl0=NEXT(vl0)){        for(vl0=vl,i=0;vl0;vl0=NEXT(vl0)){
                                 if(vl0->v->attr==(pointer)V_IND){          if(vl0->v->attr==(pointer)V_IND){
                                         if(i>=2){            if(i>=2){
                                                 error("ifplot : invalid argument");              error("ifplot : invalid argument");
                                         } else {            } else {
                                                 v[i++]=vl0->v;              v[i++]=vl0->v;
                                         }            }
                                 }          }
                         }        }
                         break;        break;
                 case O_LIST:      case O_LIST:
                         list = (LIST)BDY(arg);        list = (LIST)BDY(arg);
                         if ( OID(BDY(BDY(list))) == O_P )        if ( OID(BDY(BDY(list))) == O_P )
                                 if ( ri > 1 )          if ( ri > 1 )
                                         error("ifplot : invalid argument");            error("ifplot : invalid argument");
                                 else          else
                                         range[ri++] = list;            range[ri++] = list;
                         else        else
                                 geom = list;          geom = list;
                         break;        break;
                 default:      default:
                         error("ifplot : invalid argument"); break;        error("ifplot : invalid argument"); break;
                 }      }
         }    }
         if ( !poly ) error("ifplot : invalid argument");    if ( !poly ) error("ifplot : invalid argument");
         switch ( ri ) {    switch ( ri ) {
                 case 0:      case 0:
                         if ( !v[1] ) error("ifplot : please specify all variables");        if ( !v[1] ) error("ifplot : please specify all variables");
                         MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n);        MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n);
                         MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n);        MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n);
                         break;        break;
                 case 1:      case 1:
                         if ( !v[1] ) error("ifplot : please specify all variables");        if ( !v[1] ) error("ifplot : please specify all variables");
                         av[0] = VR((P)BDY(BDY(range[0])));        av[0] = VR((P)BDY(BDY(range[0])));
                         if ( v[0] == av[0] ) {        if ( v[0] == av[0] ) {
                                 xrange = range[0];          xrange = range[0];
                                 MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n);          MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n);
                         } else if ( v[1] == av[0] ) {        } else if ( v[1] == av[0] ) {
                                 MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n);          MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n);
                                 yrange = range[0];          yrange = range[0];
                         } else        } else
                                 error("ifplot : invalid argument");          error("ifplot : invalid argument");
                         break;        break;
                 case 2:      case 2:
                         av[0] = VR((P)BDY(BDY(range[0])));        av[0] = VR((P)BDY(BDY(range[0])));
                         av[1] = VR((P)BDY(BDY(range[1])));        av[1] = VR((P)BDY(BDY(range[1])));
                         if ( ((v[0] == av[0]) && (!v[1] || v[1] == av[1])) ||        if ( ((v[0] == av[0]) && (!v[1] || v[1] == av[1])) ||
                                  ((v[0] == av[1]) && (!v[1] || v[1] == av[0])) ) {           ((v[0] == av[1]) && (!v[1] || v[1] == av[0])) ) {
                                         xrange = range[0]; yrange = range[1];            xrange = range[0]; yrange = range[1];
                         } else error("ifplot : invalid argument");        } else error("ifplot : invalid argument");
                         break;        break;
                 default:      default:
                         error("ifplot : cannot happen"); break;        error("ifplot : cannot happen"); break;
         }    }
         can = canvas[id = search_canvas()];    can = canvas[id = search_canvas()];
         if ( !geom ) {    if ( !geom ) {
                 width = 300;      width = 300;
                 height = 300;      height = 300;
                 can->width = 300;      can->width = 300;
                 can->height = 300;      can->height = 300;
         } else {    } else {
                 can->width = QTOS((Q)BDY(BDY(geom)));      can->width = QTOS((Q)BDY(BDY(geom)));
                 can->height = QTOS((Q)BDY(NEXT(BDY(geom))));      can->height = QTOS((Q)BDY(NEXT(BDY(geom))));
                 width = can->width;      width = can->width;
                 height = can->height;      height = can->height;
         }    }
         if ( xrange ) {    if ( xrange ) {
                 n = BDY(xrange); can->vx = VR((P)BDY(n)); n = NEXT(n);      n = BDY(xrange); can->vx = VR((P)BDY(n)); n = NEXT(n);
                 can->qxmin = (Q)BDY(n); n = NEXT(n); can->qxmax = (Q)BDY(n);      can->qxmin = (Q)BDY(n); n = NEXT(n); can->qxmax = (Q)BDY(n);
                 can->xmin = ToReal(can->qxmin); can->xmax = ToReal(can->qxmax);      can->xmin = ToReal(can->qxmin); can->xmax = ToReal(can->qxmax);
         }    }
         if ( yrange ) {    if ( yrange ) {
                 n = BDY(yrange); can->vy = VR((P)BDY(n)); n = NEXT(n);      n = BDY(yrange); can->vy = VR((P)BDY(n)); n = NEXT(n);
                 can->qymin = (Q)BDY(n); n = NEXT(n); can->qymax = (Q)BDY(n);      can->qymin = (Q)BDY(n); n = NEXT(n); can->qymax = (Q)BDY(n);
                 can->ymin = ToReal(can->qymin); can->ymax = ToReal(can->qymax);      can->ymin = ToReal(can->qymin); can->ymax = ToReal(can->qymax);
         }    }
         can->wname = "ifcheck";    can->wname = "ifcheck";
         can->formula = poly;    can->formula = poly;
         tabe = (double **)ALLOCA((width+1)*sizeof(double *));    tabe = (double **)ALLOCA((width+1)*sizeof(double *));
         for ( i = 0; i <= width; i++ )    for ( i = 0; i <= width; i++ )
                 tabe[i] = (double *)ALLOCA((height+1)*sizeof(double));      tabe[i] = (double *)ALLOCA((height+1)*sizeof(double));
         for(i=0;i<=width;i++)for(j=0;j<=height;j++)tabe[i][j]=0;    for(i=0;i<=width;i++)for(j=0;j<=height;j++)tabe[i][j]=0;
         ccalc(tabe,can,0);    ccalc(tabe,can,0);
         MKMAT(m,width,height);    MKMAT(m,width,height);
         mb = BDY(m);    mb = BDY(m);
         for( ix=0; ix<width; ix++ ){    for( ix=0; ix<width; ix++ ){
                 for( iy=0; iy<height; iy++){      for( iy=0; iy<height; iy++){
                         if ( tabe[ix][iy] >= 0 ){        if ( tabe[ix][iy] >= 0 ){
                                 if ( (tabe[ix+1][iy] <= 0) ||          if ( (tabe[ix+1][iy] <= 0) ||
                                         (tabe[ix][iy+1] <= 0 ) ||            (tabe[ix][iy+1] <= 0 ) ||
                                         (tabe[ix+1][iy+1] <= 0 ) ) mb[ix][iy] = (Obj)one;            (tabe[ix+1][iy+1] <= 0 ) ) mb[ix][iy] = (Obj)one;
                         } else {        } else {
                                 if ( (tabe[ix+1][iy] >= 0 ) ||          if ( (tabe[ix+1][iy] >= 0 ) ||
                                         ( tabe[ix][iy+1] >= 0 ) ||            ( tabe[ix][iy+1] >= 0 ) ||
                                         ( tabe[ix+1][iy+1] >= 0 )) mb[ix][iy] = (Obj)one;            ( tabe[ix+1][iy+1] >= 0 )) mb[ix][iy] = (Obj)one;
                         }        }
                 }      }
         }    }
         *rp = (Obj)m;    *rp = (Obj)m;
 }  }
   
 void ccalc(double **tab,struct canvas *can,int nox)  void ccalc(double **tab,struct canvas *can,int nox)
 {  {
          double x,y,xmin,ymin,xstep,ystep;     double x,y,xmin,ymin,xstep,ystep;
          int ix,iy;     int ix,iy;
          Real r,rx,ry;     Real r,rx,ry;
          Obj fr,g;     Obj fr,g;
          int w,h;     int w,h;
          V vx,vy;     V vx,vy;
          Obj t,s;     Obj t,s;
   
          MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);     MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
          vx = can->vx;     vx = can->vx;
          vy = can->vy;     vy = can->vy;
          w = can->width; h = can->height;     w = can->width; h = can->height;
          xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;     xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;
          ymin = can->ymin; ystep = (can->ymax-can->ymin)/h;     ymin = can->ymin; ystep = (can->ymax-can->ymin)/h;
          MKReal(1.0,rx); MKReal(1.0,ry);     MKReal(1.0,rx); MKReal(1.0,ry);
          for( ix = 0, x = xmin; ix < w+1 ; ix++, x += xstep ) {     for( ix = 0, x = xmin; ix < w+1 ; ix++, x += xstep ) {
             BDY(rx) = x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);        BDY(rx) = x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);
             devalr(CO,t,&g);        devalr(CO,t,&g);
             for( iy = 0, y = ymin; iy < h+1 ; iy++, y += ystep ) {        for( iy = 0, y = ymin; iy < h+1 ; iy++, y += ystep ) {
                BDY(ry) = y;           BDY(ry) = y;
                substr(CO,0,g,vy,y?(Obj)ry:0,&t);           substr(CO,0,g,vy,y?(Obj)ry:0,&t);
                devalr(CO,t,&s);           devalr(CO,t,&s);
                tab[ix][iy] = ToReal(s);           tab[ix][iy] = ToReal(s);
             }        }
          }     }
 }  }
 /* end plot time check */  /* end plot time check */
   
 static void  static void
 Pitvversion(Q *rp)  Pitvversion(NODE arg, Q *rp)
 {  {
         STOQ(ASIR_VERSION, *rp);    STOQ(INT_ASIR_VERSION, *rp);
 }  }
   
 extern int      bigfloat;  extern int  bigfloat;
   
 static void  static void
 Pitv(NODE arg, Obj *rp)  Pitv(NODE arg, Obj *rp)
 {  {
         Num     a, i, s;    Num  a, i, s;
         Itv     c;    Itv  c;
         double  inf, sup;    double  inf, sup;
   
 #if 1  #if 1
         if ( bigfloat )    if ( bigfloat )
                 Pitvbf(arg, rp);      Pitvbf(arg, rp);
         else    else
                 Pitvd(arg,rp);      Pitvd(arg,rp);
 #else  #else
         asir_assert(ARG0(arg),O_N,"itv");    asir_assert(ARG0(arg),O_N,"itv");
         if ( argc(arg) > 1 ) {    if ( argc(arg) > 1 ) {
                 asir_assert(ARG1(arg),O_N,"itv");      asir_assert(ARG1(arg),O_N,"itv");
                 istoitv((Num)ARG0(arg),(Num)ARG1(arg),&c);      istoitv((Num)ARG0(arg),(Num)ARG1(arg),&c);
         } else {    } else {
                 a = (Num)ARG0(arg);      a = (Num)ARG0(arg);
                 if ( ! a ) {      if ( ! a ) {
                         *rp = 0;        *rp = 0;
                         return;        return;
                 }      }
                 else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat) {      else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat) {
                         *rp = (Obj)a;        *rp = (Obj)a;
                         return;        return;
                 }      }
                 else if ( NID(a) == N_IntervalDouble ) {      else if ( NID(a) == N_IntervalDouble ) {
                         inf = INF((IntervalDouble)a);        inf = INF((IntervalDouble)a);
                         sup = SUP((IntervalDouble)a);        sup = SUP((IntervalDouble)a);
                         double2bf(inf, (BF *)&i);        double2bf(inf, (BF *)&i);
                         double2bf(sup, (BF *)&s);        double2bf(sup, (BF *)&s);
                         istoitv(i,s,&c);        istoitv(i,s,&c);
                 }      }
                 else istoitv(a,a,&c);      else istoitv(a,a,&c);
         }    }
         if ( NID( c ) == N_IntervalBigFloat )    if ( NID( c ) == N_IntervalBigFloat )
                 addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp);      addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp);
         else *rp = (Obj)c;    else *rp = (Obj)c;
 #endif  #endif
 }  }
   
 static void  static void
 Pitvbf(NODE arg, Obj *rp)  Pitvbf(NODE arg, Obj *rp)
 {  {
         Num     a, i, s;    Num  a, i, s;
         Itv     c;    IntervalBigFloat  c;
         BF      ii,ss;    Num  ii,ss;
         double  inf, sup;    Real di, ds;
     double  inf, sup;
     int current_roundmode;
   
         asir_assert(ARG0(arg),O_N,"intvalbf");    asir_assert(ARG0(arg),O_N,"intvalbf");
         a = (Num)ARG0(arg);    a = (Num)ARG0(arg);
         if ( argc(arg) > 1 ) {    if ( argc(arg) > 1 ) {
                 asir_assert(ARG1(arg),O_N,"intvalbf");      asir_assert(ARG1(arg),O_N,"intvalbf");
                 i = (Num)ARG0(arg);  
                 s = (Num)ARG1(arg);      i = (Num)ARG0(arg);
                 ToBf(i, &ii);      s = (Num)ARG1(arg);
                 ToBf(s, &ss);      current_roundmode = mpfr_roundmode;
                 istoitv((Num)ii,(Num)ss,&c);      mpfr_roundmode = MPFR_RNDD;
         } else {      ii = tobf(i, DEFAULTPREC);
                 if ( ! a ) {      mpfr_roundmode = MPFR_RNDU;
                         *rp = 0;      ss = tobf(s, DEFAULTPREC);
                         return;      istoitv(ii,ss,(Itv *)&c);
                 }  //    MKIntervalBigFloat((BF)ii,(BF)ss,c);
                 else if ( NID(a) == N_IP ) {  //    ToBf(s, &ss);
                         itvtois((Itv)a, &i, &s);      mpfr_roundmode = current_roundmode;
                         ToBf(i, &ii);    } else {
                         ToBf(s, &ss);      if ( ! a ) {
                         istoitv((Num)ii,(Num)ss,&c);        *rp = 0;
                 }        return;
                 else if ( NID(a) == N_IntervalBigFloat) {      }
                         *rp = (Obj)a;      else if ( NID(a) == N_IP ) {
                         return;        itvtois((Itv)a, &i, &s);
                 }        current_roundmode = mpfr_roundmode;
                 else if ( NID(a) == N_IntervalDouble ) {        mpfr_roundmode = MPFR_RNDD;
                         inf = INF((IntervalDouble)a);        ii = tobf(i, DEFAULTPREC);
                         sup = SUP((IntervalDouble)a);        mpfr_roundmode = MPFR_RNDU;
                         double2bf(inf, (BF *)&i);        ss = tobf(s, DEFAULTPREC);
                         double2bf(sup, (BF *)&s);        istoitv(ii,ss,(Itv *)&c);
                         istoitv(i,s,&c);  //      MKIntervalBigFloat((BF)ii,(BF)ss,c);
                 }        mpfr_roundmode = current_roundmode;
                 else {      }
                         ToBf(a, (BF *)&i);      else if ( NID(a) == N_IntervalBigFloat) {
                         istoitv(i,i,&c);        *rp = (Obj)a;
                 }        return;
         }      }
         if ( c && OID( c ) == O_N && NID( c ) == N_IntervalBigFloat )      else if ( NID(a) == N_IntervalDouble ) {
                 addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp);        inf = INF((IntervalDouble)a);
         else *rp = (Obj)c;        sup = SUP((IntervalDouble)a);
         current_roundmode = mpfr_roundmode;
         //double2bf(inf, (BF *)&i);
         //double2bf(sup, (BF *)&s);
         mpfr_roundmode = MPFR_RNDD;
         MKReal(inf,di);
         ii = tobf((Num)di, DEFAULTPREC);
         mpfr_roundmode = MPFR_RNDU;
         MKReal(sup,ds);
         ss = tobf((Num)ds, DEFAULTPREC);
         istoitv(ii,ss,(Itv *)&c);
   //      MKIntervalBigFloat((BF)ii,(BF)ss,c);
         mpfr_roundmode = current_roundmode;
       }
       else {
         current_roundmode = mpfr_roundmode;
         mpfr_roundmode = MPFR_RNDD;
         ii = tobf(a, DEFAULTPREC);
         mpfr_roundmode = MPFR_RNDU;
         ss = tobf(a, DEFAULTPREC);
         //ToBf(a, (BF *)&i);
         istoitv(ii,ss,(Itv *)&c);
   //      MKIntervalBigFloat((BF)ii,(BF)ss,c);
         mpfr_roundmode = current_roundmode;
       }
     }
   //  if ( c && OID( c ) == O_N && NID( c ) == N_IntervalBigFloat )
   //    addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp);
   //  else *rp = (Obj)c;
     *rp = (Obj)c;
 }  }
   
 static void  static void
 Pitvd(NODE arg, Obj *rp)  Pitvd(NODE arg, Obj *rp)
 {  {
         double  inf, sup;    double  inf, sup;
         Num     a, a0, a1, t;    Num  a, a0, a1, t;
         Itv     ia;    Itv  ia;
         IntervalDouble  d;    IntervalDouble  d;
   
         asir_assert(ARG0(arg),O_N,"intvald");    asir_assert(ARG0(arg),O_N,"intvald");
         a0 = (Num)ARG0(arg);    a0 = (Num)ARG0(arg);
         if ( argc(arg) > 1 ) {    if ( argc(arg) > 1 ) {
                 asir_assert(ARG1(arg),O_N,"intvald");      asir_assert(ARG1(arg),O_N,"intvald");
                 a1 = (Num)ARG1(arg);      a1 = (Num)ARG1(arg);
         } else {    } else {
                 if ( a0 && OID(a0)==O_N && NID(a0)==N_IntervalDouble ) {      if ( a0 && OID(a0)==O_N && NID(a0)==N_IntervalDouble ) {
                         inf = INF((IntervalDouble)a0);        inf = INF((IntervalDouble)a0);
                         sup = SUP((IntervalDouble)a0);        sup = SUP((IntervalDouble)a0);
                         MKIntervalDouble(inf,sup,d);        MKIntervalDouble(inf,sup,d);
                         *rp = (Obj)d;        *rp = (Obj)d;
                         return;        return;
                 }      }
                 a1 = (Num)ARG0(arg);      a1 = (Num)ARG0(arg);
         }    }
         if ( compnum(0,a0,a1) > 0 ) {    if ( compnum(0,a0,a1) > 0 ) {
                 t = a0; a0 = a1; a1 = t;      t = a0; a0 = a1; a1 = t;
         }    }
         inf = ToRealDown(a0);    inf = toRealDown(a0);
         sup = ToRealUp(a1);    sup = toRealUp(a1);
         MKIntervalDouble(inf,sup,d);    MKIntervalDouble(inf,sup,d);
         *rp = (Obj)d;    *rp = (Obj)d;
 }  }
   
 static void  static void
 Pinf(NODE arg, Obj *rp)  Pinf(NODE arg, Obj *rp)
 {  {
         Num     a, i, s;    Num  a, i, s;
         Real    r;    Real  r;
         double  d;    double  d;
   
         a = (Num)ARG0(arg);    a = (Num)ARG0(arg);
         if ( ! a ) {    if ( ! a ) {
                 *rp = 0;      *rp = 0;
         } else if  ( OID(a) == O_N ) {    } else if  ( OID(a) == O_N ) {
                 switch ( NID(a) ) {      switch ( NID(a) ) {
                         case N_IntervalDouble:        case N_IntervalDouble:
                                 d = INF((IntervalDouble)a);          d = INF((IntervalDouble)a);
                                 MKReal(d, r);          MKReal(d, r);
                                 *rp = (Obj)r;          *rp = (Obj)r;
                                 break;          break;
                         case N_IP:        case N_IP:
                         case N_IntervalBigFloat:        case N_IntervalBigFloat:
                         case N_IntervalQuad:        case N_IntervalQuad:
                                 itvtois((Itv)ARG0(arg),&i,&s);          itvtois((Itv)ARG0(arg),&i,&s);
                                 *rp = (Obj)i;          *rp = (Obj)i;
                                 break;          break;
                         default:        default:
                                 *rp = (Obj)a;          *rp = (Obj)a;
                                 break;          break;
                 }      }
         } else {    } else {
                 *rp = (Obj)a;      *rp = (Obj)a;
         }    }
 }  }
   
 static void  static void
 Psup(NODE arg, Obj *rp)  Psup(NODE arg, Obj *rp)
 {  {
         Num     a, i, s;    Num  a, i, s;
         Real    r;    Real  r;
         double  d;    double  d;
   
         a = (Num)ARG0(arg);    a = (Num)ARG0(arg);
         if ( ! a ) {    if ( ! a ) {
                 *rp = 0;      *rp = 0;
         } else if  ( OID(a) == O_N ) {    } else if  ( OID(a) == O_N ) {
                 switch ( NID(a) ) {      switch ( NID(a) ) {
                         case N_IntervalDouble:        case N_IntervalDouble:
                                 d = SUP((IntervalDouble)a);          d = SUP((IntervalDouble)a);
                                 MKReal(d, r);          MKReal(d, r);
                                 *rp = (Obj)r;          *rp = (Obj)r;
                                 break;          break;
                         case N_IP:        case N_IP:
                         case N_IntervalBigFloat:        case N_IntervalBigFloat:
                         case N_IntervalQuad:        case N_IntervalQuad:
                                 itvtois((Itv)ARG0(arg),&i,&s);          itvtois((Itv)ARG0(arg),&i,&s);
                                 *rp = (Obj)s;          *rp = (Obj)s;
                                 break;          break;
                         default:        default:
                                 *rp = (Obj)a;          *rp = (Obj)a;
                                 break;          break;
                 }      }
         } else {    } else {
                         *rp = (Obj)a;        *rp = (Obj)a;
         }    }
 }  }
   
 static void  static void
 Pmid(NODE arg, Obj *rp)  Pmid(NODE arg, Obj *rp)
 {  {
         Num     a, s;    Num  a, s;
         Real    r;    Real  r;
         double  d;    double  d;
   
         a = (Num)ARG0(arg);    a = (Num)ARG0(arg);
         if ( ! a ) {    if ( ! a ) {
                 *rp = 0;      *rp = 0;
         } else switch (OID(a)) {    } else switch (OID(a)) {
                 case O_N:      case O_N:
                         if ( NID(a) == N_IntervalDouble ) {        if ( NID(a) == N_IntervalDouble ) {
                                 d = ( INF((IntervalDouble)a)+SUP((IntervalDouble)a) ) / 2.0;          d = ( INF((IntervalDouble)a)+SUP((IntervalDouble)a) ) / 2.0;
                                 MKReal(d, r);          MKReal(d, r);
                                 *rp = (Obj)r;          *rp = (Obj)r;
                         } else if ( NID(a) == N_IntervalQuad ) {        } else if ( NID(a) == N_IntervalQuad ) {
                                 error("mid: not supported operation");          error("mid: not supported operation");
                                 *rp = 0;          *rp = 0;
                         } else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat ) {        } else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat ) {
                                 miditvp((Itv)ARG0(arg),&s);          miditvp((Itv)ARG0(arg),&s);
                                 *rp = (Obj)s;          *rp = (Obj)s;
                         } else {        } else {
                                 *rp = (Obj)a;          *rp = (Obj)a;
                         }        }
                         break;        break;
 #if 0  #if 0
                 case O_P:      case O_P:
                 case O_R:      case O_R:
                 case O_LIST:      case O_LIST:
                 case O_VECT:      case O_VECT:
                 case O_MAT:      case O_MAT:
 #endif  #endif
                 default:      default:
                         *rp = (Obj)a;        *rp = (Obj)a;
                         break;        break;
         }    }
 }  }
   
 static void  static void
 Pcup(NODE arg, Obj *rp)  Pcup(NODE arg, Obj *rp)
 {  {
         Itv     s;    Itv  s;
         Num     a, b;    Num  a, b;
   
         asir_assert(ARG0(arg),O_N,"cup");    asir_assert(ARG0(arg),O_N,"cup");
         asir_assert(ARG1(arg),O_N,"cup");    asir_assert(ARG1(arg),O_N,"cup");
         a = (Num)ARG0(arg);    a = (Num)ARG0(arg);
         b = (Num)ARG1(arg);    b = (Num)ARG1(arg);
         if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {    if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {
                 cupitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp);      cupitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp);
         } else {    } else {
                 cupitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);      cupitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);
                 *rp = (Obj)s;      *rp = (Obj)s;
         }    }
 }  }
   
 static void  static void
 Pcap(NODE arg, Obj *rp)  Pcap(NODE arg, Obj *rp)
 {  {
         Itv     s;    Itv  s;
         Num     a, b;    Num  a, b;
   
         asir_assert(ARG0(arg),O_N,"cap");    asir_assert(ARG0(arg),O_N,"cap");
         asir_assert(ARG1(arg),O_N,"cap");    asir_assert(ARG1(arg),O_N,"cap");
         a = (Num)ARG0(arg);    a = (Num)ARG0(arg);
         b = (Num)ARG1(arg);    b = (Num)ARG1(arg);
         if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {    if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {
                 capitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp);      capitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp);
         } else {    } else {
                 capitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);      capitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);
                 *rp = (Obj)s;      *rp = (Obj)s;
         }    }
 }  }
   
 static void  static void
Line 509  Pwidth(arg,rp)
Line 610  Pwidth(arg,rp)
 NODE arg;  NODE arg;
 Obj *rp;  Obj *rp;
 {  {
         Num     s;    Num  s;
         Num     a;    Num  a;
   
         asir_assert(ARG0(arg),O_N,"width");    asir_assert(ARG0(arg),O_N,"width");
         a = (Num)ARG0(arg);    a = (Num)ARG0(arg);
         if ( ! a ) {    if ( ! a ) {
                 *rp = 0;      *rp = 0;
         } else if ( NID(a) == N_IntervalDouble ) {    } else if ( NID(a) == N_IntervalDouble ) {
                 widthitvd((IntervalDouble)a, (Num *)rp);      widthitvd((IntervalDouble)a, (Num *)rp);
         } else {    } else {
                 widthitvp((Itv)ARG0(arg),&s);      widthitvp((Itv)ARG0(arg),&s);
                 *rp = (Obj)s;      *rp = (Obj)s;
         }    }
 }  }
   
 static void  static void
Line 529  Pabsitv(arg,rp)
Line 630  Pabsitv(arg,rp)
 NODE arg;  NODE arg;
 Obj *rp;  Obj *rp;
 {  {
         Num     s;    Num  s;
         Num     a, b;    Num  a, b;
   
         asir_assert(ARG0(arg),O_N,"absitv");    asir_assert(ARG0(arg),O_N,"absitv");
         a = (Num)ARG0(arg);    a = (Num)ARG0(arg);
         if ( ! a ) {    if ( ! a ) {
                 *rp = 0;      *rp = 0;
         } else if ( NID(a) == N_IntervalDouble ) {    } else if ( NID(a) == N_IntervalDouble ) {
                 absitvd((IntervalDouble)a, (Num *)rp);      absitvd((IntervalDouble)a, (Num *)rp);
         } else {    } else {
                 absitvp((Itv)ARG0(arg),&s);      absitvp((Itv)ARG0(arg),&s);
                 *rp = (Obj)s;      *rp = (Obj)s;
         }    }
 }  }
   
 static void  static void
Line 549  Pdistance(arg,rp)
Line 650  Pdistance(arg,rp)
 NODE arg;  NODE arg;
 Obj *rp;  Obj *rp;
 {  {
         Num     s;    Num  s;
         Num     a, b;    Num  a, b;
   
         asir_assert(ARG0(arg),O_N,"distance");    asir_assert(ARG0(arg),O_N,"distance");
         asir_assert(ARG1(arg),O_N,"distance");    asir_assert(ARG1(arg),O_N,"distance");
         a = (Num)ARG0(arg);    a = (Num)ARG0(arg);
         b = (Num)ARG1(arg);    b = (Num)ARG1(arg);
         if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {    if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {
                 distanceitvd((IntervalDouble)a, (IntervalDouble)b, (Num *)rp);      distanceitvd((IntervalDouble)a, (IntervalDouble)b, (Num *)rp);
         } else {    } else {
                 distanceitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);      distanceitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);
                 *rp = (Obj)s;      *rp = (Obj)s;
         }    }
 }  }
   
 static void  static void
Line 569  Pinitv(arg,rp)
Line 670  Pinitv(arg,rp)
 NODE arg;  NODE arg;
 Obj *rp;  Obj *rp;
 {  {
         int     s;    int  s;
         Q       q;    Q  q;
   
         asir_assert(ARG0(arg),O_N,"intval");    asir_assert(ARG0(arg),O_N,"intval");
         asir_assert(ARG1(arg),O_N,"intval");    asir_assert(ARG1(arg),O_N,"intval");
         if ( ! ARG1(arg) ) {    if ( ! ARG1(arg) ) {
                 if ( ! ARG0(arg) ) s = 1;      if ( ! ARG0(arg) ) s = 1;
                 else s = 0;      else s = 0;
         }    }
         else if ( NID(ARG1(arg)) == N_IntervalDouble ) {    else if ( NID(ARG1(arg)) == N_IntervalDouble ) {
                 s = initvd((Num)ARG0(arg),(IntervalDouble)ARG1(arg));      s = initvd((Num)ARG0(arg),(IntervalDouble)ARG1(arg));
   
         } else if ( NID(ARG1(arg)) == N_IP || NID(ARG1(arg)) == N_IntervalBigFloat ) {    } else if ( NID(ARG1(arg)) == N_IP || NID(ARG1(arg)) == N_IntervalBigFloat ) {
                 if ( ! ARG0(arg) ) s = initvp((Num)ARG0(arg),(Itv)ARG1(arg));      if ( ! ARG0(arg) ) s = initvp((Num)ARG0(arg),(Itv)ARG1(arg));
                 else if ( NID(ARG0(arg)) == N_IP ) {      else if ( NID(ARG0(arg)) == N_IP ) {
                         s = itvinitvp((Itv)ARG0(arg),(Itv)ARG1(arg));        s = itvinitvp((Itv)ARG0(arg),(Itv)ARG1(arg));
                 } else {      } else {
                         s = initvp((Num)ARG0(arg),(Itv)ARG1(arg));        s = initvp((Num)ARG0(arg),(Itv)ARG1(arg));
                 }      }
         } else {    } else {
                 s = ! compnum(0,(Num)ARG0(arg),(Num)ARG1(arg));      s = ! compnum(0,(Num)ARG0(arg),(Num)ARG1(arg));
         }    }
         STOQ(s,q);    STOQ(s,q);
         *rp = (Obj)q;    *rp = (Obj)q;
 }  }
   
 static void  static void
Line 600  Pdisjitv(arg,rp)
Line 701  Pdisjitv(arg,rp)
 NODE arg;  NODE arg;
 Obj *rp;  Obj *rp;
 {  {
         Itv     s;    Itv  s;
   
         asir_assert(ARG0(arg),O_N,"disjitv");    asir_assert(ARG0(arg),O_N,"disjitv");
         asir_assert(ARG1(arg),O_N,"disjitv");    asir_assert(ARG1(arg),O_N,"disjitv");
         error("disjitv: not implemented yet");    error("disjitv: not implemented yet");
         if ( ! s ) *rp = 0;    if ( ! s ) *rp = 0;
         else *rp = (Obj)ONE;    else *rp = (Obj)ONE;
 }  }
   
   static void
   PzeroRewriteMode(NODE arg, Obj *rp)
   {
     Q  a, r;
   
     STOQ(zerorewrite,r);
     *rp = (Obj)r;
   
     if (arg) {
       a = (Q)ARG0(arg);
       if(!a) {
         zerorewrite = 0;
       } else if ( (NUM(a)&&INT(a)) ){
         zerorewrite = 1;
       }
     }
   }
   
   static void
   PzeroRewriteCountClear(NODE arg, Obj *rp)
   {
     Q  a, r;
   
     STOQ(zerorewriteCount,r);
     *rp = (Obj)r;
   
     if (arg) {
       a = (Q)ARG0(arg);
       if(a &&(NUM(a)&&INT(a))){
         zerorewriteCount = 0;
       }
     }
   }
   
   static void
   PzeroRewriteCount(NODE arg, Obj *rp)
   {
     Q  r;
   
     STOQ(zerorewriteCount,r);
     *rp = (Obj)r;
   }
   
   
 #endif  #endif
 extern int      printmode;  extern int  printmode;
   
 static void     pprintmode( void )  static void  pprintmode( void )
 {  {
         switch (printmode) {    switch (printmode) {
 #if defined(INTERVAL)  #if defined(INTERVAL)
                 case MID_PRINTF_E:      case MID_PRINTF_E:
                         fprintf(stderr,"Interval printing mode is a mitpoint type.\n");        fprintf(stderr,"Interval printing mode is a mitpoint type.\n");
 #endif  #endif
                 case PRINTF_E:      case PRINTF_E:
                         fprintf(stderr,"Printf's double printing mode is \"%%.16e\".\n");        fprintf(stderr,"Printf's double printing mode is \"%%.16e\".\n");
                         break;        break;
 #if defined(INTERVAL)  #if defined(INTERVAL)
                 case MID_PRINTF_G:      case MID_PRINTF_G:
                         fprintf(stderr,"Interval printing mode is a mitpoint type.\n");        fprintf(stderr,"Interval printing mode is a mitpoint type.\n");
 #endif  #endif
                 default:      default:
                 case PRINTF_G:      case PRINTF_G:
                         fprintf(stderr,"Printf's double printing mode is \"%%g\".\n");        fprintf(stderr,"Printf's double printing mode is \"%%g\".\n");
                         break;        break;
         }    }
 }  }
   
 static void  static void
 Pprintmode(NODE arg, Obj *rp)  Pprintmode(NODE arg, Obj *rp)
 {  {
         int     l;    int  l;
         Q       a, r;    Q  a, r;
   
         a = (Q)ARG0(arg);    a = (Q)ARG0(arg);
         if(!a||(NUM(a)&&INT(a))){    if(!a||(NUM(a)&&INT(a))){
                 l=QTOS(a);      l=QTOS(a);
                 if ( l < 0 ) l = 0;      if ( l < 0 ) l = 0;
 #if defined(INTERVAL)  #if defined(INTERVAL)
                 else if ( l > MID_PRINTF_E ) l = 0;      else if ( l > MID_PRINTF_E ) l = 0;
 #else  #else
                 else if ( l > PRINTF_E ) l = 0;      else if ( l > PRINTF_E ) l = 0;
 #endif  #endif
                 STOQ(printmode,r);      STOQ(printmode,r);
                 *rp = (Obj)r;      *rp = (Obj)r;
                 printmode = l;      printmode = l;
                 pprintmode();      pprintmode();
     } else {
       *rp = 0;
     }
   }
   
   #if defined(INTERVAL)
   void
   Ppi_itvd(NODE arg, Obj *rp)
   {
           double inf, sup;
           IntervalDouble c;
       FPMINUSINF
           sscanf("3.1415926535897932384626433832795028841971693993751", "%lf", &inf);
       FPPLUSINF
       sscanf("3.1415926535897932384626433832795028841971693993752", "%lf", &sup);
       FPNEAREST
       MKIntervalDouble(inf,sup,c);
           *rp = (Obj)c;
   }
   void
   Pe_itvd(NODE arg, Obj *rp)
   {
           double inf, sup;
           IntervalDouble c;
       FPMINUSINF
           sscanf( "2.7182818284590452353602874713526624977572470936999", "%lf", &inf);
       FPPLUSINF
       sscanf( "2.7182818284590452353602874713526624977572470937000", "%lf", &sup);
       FPNEAREST
       MKIntervalDouble(inf,sup,c);
           *rp = (Obj)c;
   }
   void
   Pln2_itvd(NODE arg, Obj *rp)
   {
           double inf, sup;
           IntervalDouble c;
       FPMINUSINF
           sscanf( "0.69314718055994530941723212145817656807550013436025", "%lf", &inf);
       FPPLUSINF
       sscanf( "0.69314718055994530941723212145817656807550013436026", "%lf", &sup);
       FPNEAREST
       MKIntervalDouble(inf,sup,c);
           *rp = (Obj)c;
   }
   
   void mpfi_func(NODE arg, int (*mpfi_f)(), int prec, Obj *rp)
   {
     Num a, ii, ss;
     Itv c;
     BF inf, sup;
     int arg1prec;
     mpfi_t mpitv, rv;
   
   
   /*
     if ( argc(arg) == 2 ) {
       prec = QTOS((Q)ARG1(arg))*3.32193;
       if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
       else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
     } else {
       prec = 0;
       prec = mpfr_get_default_prec();
     }
   */
     if ( prec > 0 ) arg1prec = prec;
     else arg1prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : mpfr_get_default_prec();
     a = ARG0(arg);
     itvtois((Itv)a, &ii, &ss);
   
     inf = (BF)tobf(ii, arg1prec);
     sup = (BF)tobf(ss, arg1prec);
   
     mpfi_init2(rv,arg1prec);
     mpfi_init2(mpitv,arg1prec);
     mpfr_set(&(mpitv->left),  BDY(inf), MPFR_RNDD);
     mpfr_set(&(mpitv->right), BDY(sup), MPFR_RNDU);
   
     (*mpfi_f)(rv, mpitv);
     //mpfi_sin(rv, mpitv);
   
     MPFRTOBF(&(rv->left), inf);
     MPFRTOBF(&(rv->right), sup);
   
     if ( !cmpbf((Num)inf,0) ) inf = 0;
     if ( !cmpbf((Num)sup,0) ) sup = 0;
   
     if ( inf || sup ) {
           istoitv((Num)inf, (Num)sup, &c);
     } else {
       c = 0;
     }
     *rp = (Obj)c;
     //mpfi_clear(rv);
     mpfi_clear(mpitv);
   }
   
   void mpfi_func_d(NODE arg, int (*mpfi_f)(), Obj *rp)
   {
     Obj bfv;
     Num ii, ss;
     IntervalDouble c;
     double inf, sup;
   
     mpfi_func(arg, mpfi_f, 53, &bfv);
     itvtois((Itv)bfv, &ii, &ss);
     inf = toRealDown(ii);
     sup = toRealUp(ss);
     MKIntervalDouble(inf,sup,c);
     *rp = (Obj)c;
   }
   
   
   void
   Psinitvd(NODE arg, Obj *rp)
   {
     Obj bfv;
     Num ii,ss;
     IntervalDouble c;
     double  ai,as,mas, bi,bs;
     double  inf,sup;
   
   #if 1
     mpfi_func(arg, mpfi_sin, 53, &bfv);
     itvtois((Itv)bfv, &ii, &ss);
     inf = toRealDown(ii);
     sup = toRealUp(ss);
     MKIntervalDouble(inf,sup,c);
     *rp = (Obj)c;
   #else
       a = ARG0(arg);
       Num2double(a,&ai,&as);
       FPMINUSINF
           inf = sin(ai);
       FPPLUSINF
           sup = sin(as);
       FPNEAREST
       MKIntervalDouble(inf,sup,c);
           *rp = (Obj)c;
   #endif
   }
   
   void
   Psinitv(NODE arg, Obj *rp)
   {
     //Num a;
     //Itv c;
     //BF inf, sup;
     //int prec;
     //BF r,re,im;
     //mpfi_t mpitv, rv;
   
   #if 1
     mpfi_func(arg, mpfi_sin, 0, rp);
   #else
     prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : mpfr_get_default_prec();
     a = ARG0(arg);
     itvtois((Itv)a, (Num *)&inf, (Num *)&sup);
   
     mpfi_init2(rv,prec);
     mpfi_init2(mpitv,prec);
     mpfr_set(&(mpitv->left), inf->body, MPFR_RNDD);
     mpfr_set(&(mpitv->right), BDY(sup), MPFR_RNDU);
   
     //(*mpfi_f)(rv, mpitv);
     mpfi_sin(rv, mpitv);
   
     MPFRTOBF(&(rv->left), inf);
     MPFRTOBF(&(rv->right), sup);
     istoitv((Num)inf, (Num)sup, &c);
     *rp = (Obj)c;
     mpfi_clear(rv);
     mpfi_clear(mpitv);
   #endif
   }
   
   
   
   
   //void evalitvr(VL, Obj, int, int, Obj *);
   
   static void
   Pevalitv(NODE arg, Obj *rp)
   {
     int prec;
   
     asir_assert(ARG0(arg),O_R,"evalitv");
     if ( argc(arg) == 2 ) {
       long int mpfr_prec_max = MPFR_PREC_MAX;
       prec = QTOS((Q)ARG1(arg))*3.32193;
       if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
       //else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
       else if ( prec > mpfr_prec_max ) prec = mpfr_prec_max;
     } else
       prec = 0;
     evalitvr(CO,(Obj)ARG0(arg),prec,EvalIntervalBigFloat,rp);
   }
   
   static void
   Pevalitvd(NODE arg, Obj *rp)
   {
     asir_assert(ARG0(arg),O_R,"evalitvd");
     evalitvr(CO,(Obj)ARG0(arg),53,EvalIntervalDouble,rp);
   }
   
   // in parse/puref.c
   void instoobj(PFINS ins,Obj *rp);
   
   // in this
   void evalitvr(VL, Obj, int, int, Obj *);
   void evalitvp(VL,   P, int, int, P *);
   void evalitvv(VL,   V, int, int, Obj *);
   void evalitvins(PFINS, int, int, Obj *);
   
   
   
   void evalitvr(VL vl,Obj a,int prec, int type, Obj *c)
   {
     Obj nm,dn;
   
     if ( !a )
       *c = 0;
     else {
       switch ( OID(a) ) {
         case O_N:
           toInterval((Num)a, prec, type, (Num *)c);
                   break;
         case O_P:
           evalitvp(vl,(P)a,prec,type,(P *)c);
                   break;
         case O_R:
           evalitvp(vl,NM((R)a),prec,type,(P *)&nm);
                   evalitvp(vl,DN((R)a),prec,type,(P *)&dn);
           divr(vl,nm,dn,c);
           break;
         default:
           error("evalr : not implemented"); break;
       }
     }
   }
   
   void evalitvp(VL vl,P p,int prec, int type, P *pr)
   {
     P t;
     DCP dc,dcr0,dcr;
     Obj u;
   
     if ( !p || NUM(p) ) {
           toInterval((Num)p, prec, type, (Num *)pr);
       //*pr = p;
     } else {
       for ( dcr0 = 0, dc = DC((P)p); dc; dc = NEXT(dc) ) {
         evalitvp(vl,COEF(dc),prec,type, &t);
         if ( t ) {
           NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
         }
       }
       if ( !dcr0 ) {
         *pr = 0; return;
       } else {
         NEXT(dcr) = 0; MKP(VR(p),dcr0,t);
       }
       if ( NUM(t) ) {
         //*pr = t;
           toInterval((Num)t, prec, type, (Num *)pr);
                   return;
           } else if ( (VR(t) != VR(p)) || ((vid)VR(p)->attr != V_PF) ) {
         *pr = t; return;
       } else {
           evalitvv(vl,VR(p),prec,type,&u);
                   substr(vl,1,(Obj)t,VR(p),u, (Obj *)&t);
                   if ( t && NUM(t) ) {
                   toInterval((Num)t, prec, type, (Num *)pr);
                   } else {
                           *pr = t;
                   }
       }
     }
   }
   
   void evalitvv(VL vl,V v,int prec, int type, Obj *rp)
   {
     PFINS ins,tins;
     PFAD ad,tad;
     PF pf;
     P t;
     int i;
   
     if ( (vid)v->attr != V_PF ) {
       MKV(v,t); *rp = (Obj)t;
     } else {
       ins = (PFINS)v->priv; ad = ins->ad; pf = ins->pf;
       tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
       tins->pf = pf;
       for ( i = 0, tad = tins->ad; i < pf->argc; i++ ) {
         tad[i].d = ad[i].d;
                   evalitvr(vl,ad[i].arg,prec,type,&tad[i].arg);
       }
       evalitvins(tins,prec,type,rp);
     }
   }
   
   void evalitvins(PFINS ins,int prec, int type, Obj *rp)
   {
     PF pf;
     PFINS tins;
     PFAD ad,tad;
     int i;
     Q q;
     V v;
     P x;
     NODE n0,n;
   
   
     pf = ins->pf; ad = ins->ad;
     tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD));
     tins->pf = pf; tad = tins->ad;
     for ( i = 0; i < pf->argc; i++ ) {
       //tad[i].d = ad[i].d; evalr(CO,ad[i].arg,prec,&tad[i].arg);
       tad[i].d = ad[i].d; evalitvr(CO,ad[i].arg,prec,type,&tad[i].arg);
      }
     for ( i = 0; i < pf->argc; i++ )
       if ( tad[i].d || (tad[i].arg && !NUM(tad[i].arg)) ) break;
     if ( (i != pf->argc) || !pf->intervalfunc[type] ) { ///////////////////////////
       instoobj(tins,rp);
     } else {
           int IsCPLX = 0;
       for ( n0 = 0, i = 0; i < pf->argc; i++ ) {
         NEXTNODE(n0,n); BDY(n) = (pointer)tad[i].arg;
                   if (tad[i].arg && NID(tad[i].arg) == N_C) IsCPLX = 1;
       }
       if ( prec ) {
         NEXTNODE(n0,n); STOQ(prec,q); BDY(n) = (pointer)q;
       }
       if ( n0 ) NEXT(n) = 0;
   
   
           if ( IsCPLX ) {
           instoobj(tins,rp);
         } else {          } else {
                 *rp = 0;          (*pf->intervalfunc[type])(n0,rp);
         }          }
     }
 }  }
   
   
   
   void Pitvbf_pi(NODE arg, Obj *rp)
   {
     BF inf, sup;
     IntervalBigFloat c;
     mpfi_t rv;
     int prec;
   
     prec = (arg) ? QTOS((Q)ARG0(arg)) : 0; //mpfr_get_default_prec();
     prec ? mpfi_init2(rv,prec) : mpfi_init(rv);
   
     mpfi_const_pi(rv);
   
     MPFRTOBF(&(rv->left), inf);
     MPFRTOBF(&(rv->right), sup);
   
     MKIntervalBigFloat(inf,sup,c);
   
     *rp = (Obj)c;
   }
   
   void Pitvd_pi(NODE arg, Obj *rp)
   {
     BF bfinf, bfsup;
     IntervalDouble c;
     mpfi_t rv;
     double inf, sup;
   
     mpfi_init2(rv,53);
   
     mpfi_const_pi(rv);
   
     MPFRTOBF(&(rv->left),  bfinf);
     MPFRTOBF(&(rv->right), bfsup);
   
     inf = toRealDown((Num)bfinf);
     sup = toRealUp((Num)bfsup);
   
     MKIntervalDouble(inf,sup,c);
   
     *rp = (Obj)c;
   }
   
   void Pitvbf_e(NODE arg,Obj *rp)
   {
     BF inf, sup;
     IntervalBigFloat c;
     mpfi_t rv;
     mpfi_t one;
     int prec;
   
     prec = (arg) ? QTOS((Q)ARG0(arg)) : 0; //mpfr_get_default_prec();
     prec ? mpfi_init2(rv,prec) : mpfi_init(rv);
   
     mpfi_init(one);
     mpfi_set_ui(one,1);
     mpfi_exp(rv,one);
   
     MPFRTOBF(&(rv->left), inf);
     MPFRTOBF(&(rv->right), sup);
   
     MKIntervalBigFloat(inf,sup,c);
     //istoitv((Num)inf, (Num)sup, &c);
     *rp = (Obj)c;
     mpfi_clear(one);
   }
   
   void Pitvd_e(NODE arg, Obj *rp)
   {
     BF bfinf, bfsup;
     IntervalDouble c;
     mpfi_t rv;
     mpfi_t one;
     double inf, sup;
   
     mpfi_init2(rv,53);
   
     mpfi_init2(one, 53);
     mpfi_set_ui(one,1);
     mpfi_exp(rv,one);
   
   
     MPFRTOBF(&(rv->left),  bfinf);
     MPFRTOBF(&(rv->right), bfsup);
   
     inf = toRealDown((Num)bfinf);
     sup = toRealUp((Num)bfsup);
   
     MKIntervalDouble(inf,sup,c);
   
     *rp = (Obj)c;
   }
   
   void (*pi_itv_ft[])() =         {Pitvd_pi,              0,      Pitvbf_pi};
   void (*e_itv_ft[])() =          {Pitvd_e,               0,      Pitvbf_e};
   //void (*sin_itv_ft[])() =      {Psinitvd,      0,      Psinitv};
   void (*sin_itv_ft[])() =        {Pitvd_sin,             0,      Pitvbf_sin};
   void (*cos_itv_ft[])() =        {Pitvbf_cos,    0,      Pitvbf_cos};
   void (*tan_itv_ft[])() =        {Pitvd_tan,             0,      Pitvbf_tan};
   void (*asin_itv_ft[])() =       {Pitvd_asin,    0,      Pitvbf_asin};
   void (*acos_itv_ft[])() =       {Pitvd_acos,    0,      Pitvbf_acos};
   void (*atan_itv_ft[])() =       {Pitvd_atan,    0,      Pitvbf_atan};
   void (*sinh_itv_ft[])() =       {Pitvd_sinh,    0,      Pitvbf_sinh};
   void (*cosh_itv_ft[])() =       {Pitvd_cosh,    0,      Pitvbf_cosh};
   void (*tanh_itv_ft[])() =       {Pitvd_tanh,    0,      Pitvbf_tanh};
   void (*asinh_itv_ft[])() =      {Pitvd_asinh,   0,      Pitvbf_asinh};
   void (*acosh_itv_ft[])() =      {Pitvd_acosh,   0,      Pitvbf_acosh};
   void (*atanh_itv_ft[])() =      {Pitvd_atanh,   0,      Pitvbf_atanh};
   void (*exp_itv_ft[])() =        {Pitvd_exp,             0,      Pitvbf_exp};
   void (*log_itv_ft[])() =        {Pitvd_log,             0,      Pitvbf_log};
   void (*abs_itv_ft[])() =        {0};
   void (*pow_itv_ft[])() =        {Pitvbf_pow,    0,      Pitvbf_pow};
   //void (*pow_itv_ft[])() =      {0,     0,      0};
   
   
   void Pitvbf_sin(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_sin, 0, rp);
   }
   
   void Pitvbf_cos(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_cos, 0, rp);
   }
   
   void Pitvbf_tan(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_tan, 0, rp);
   }
   
   void Pitvbf_asin(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_asin, 0, rp);
   }
   
   void Pitvbf_acos(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_acos, 0, rp);
   }
   
   void Pitvbf_atan(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_atan, 0, rp);
   }
   
   void Pitvbf_sinh(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_sinh, 0, rp);
   }
   
   void Pitvbf_cosh(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_cosh, 0, rp);
   }
   
   void Pitvbf_tanh(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_tanh, 0, rp);
   }
   
   void Pitvbf_asinh(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_asinh, 0, rp);
   }
   
   void Pitvbf_acosh(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_acosh, 0, rp);
   }
   
   void Pitvbf_atanh(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_atanh, 0, rp);
   }
   
   void Pitvbf_exp(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_exp, 0, rp);
   }
   
   void Pitvbf_log(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_log, 0, rp);
   }
   
   void Pitvbf_abs(NODE arg,Obj *rp)
   {
     mpfi_func(arg, mpfi_abs, 0, rp);
   }
   
   void Pitvd_sin(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_sin, rp);
   }
   
   void Pitvd_cos(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_cos, rp);
   }
   
   void Pitvd_tan(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_tan, rp);
   }
   
   void Pitvd_asin(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_asin, rp);
   }
   
   void Pitvd_acos(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_acos, rp);
   }
   
   void Pitvd_atan(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_atan, rp);
   }
   
   void Pitvd_sinh(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_sinh, rp);
   }
   
   void Pitvd_cosh(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_cosh, rp);
   }
   
   void Pitvd_tanh(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_tanh, rp);
   }
   
   void Pitvd_asinh(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_asinh, rp);
   }
   
   void Pitvd_acosh(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_acosh, rp);
   }
   
   void Pitvd_atanh(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_atanh, rp);
   }
   
   void Pitvd_exp(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_exp, rp);
   }
   
   void Pitvd_log(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_log, rp);
   }
   
   void Pitvd_abs(NODE arg,Obj *rp)
   {
     mpfi_func_d(arg, mpfi_abs, rp);
   }
   
   /*
   void mp_factorial(NODE arg,Num *rp)
   {
     struct oNODE arg0;
     Num a,a1;
   
     a = (Num)ARG0(arg);
     if ( !a ) *rp = (Num)ONE;
     else if ( INT(a) ) Pfac(arg,rp);
     else {
       addnum(0,a,(Num)ONE,&a1);
       arg0.body = (pointer)a1;
       arg0.next = arg->next;
       Pmpfr_gamma(&arg0,rp);
     }
   }
   */
   
   void Pitvbf_pow(NODE arg,Num *rp)
   {
     Num a,e;
     BF r,re,im;
     C c;
     mpc_t mpc,a1,e1;
     int prec;
   
     prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : mpfr_get_default_prec();
     a = ARG0(arg);
     e = ARG1(arg);
     if ( !e ) {
       *rp = (Num)ONE;
     } else if ( !a ) {
       *rp = 0;
     } else if ( NID(a) == N_C || NID(e) == N_C ) {
       error("itv_pow() : not support args");
       *rp = 0;
     } else {
       Num ii, ss;
       Itv c;
       BF inf, sup;
       mpfi_t a_val, e_val, rv;
   
       mpfi_init2(rv,prec);
   
       itvtois((Itv)a, &ii, &ss);
       inf = (BF)tobf(ii, prec);
       sup = (BF)tobf(ss, prec);
       mpfi_init2(a_val,prec);
       mpfr_set(&(a_val->left),  BDY(inf), MPFR_RNDD);
       mpfr_set(&(a_val->right), BDY(sup), MPFR_RNDU);
   
       itvtois((Itv)e, &ii, &ss);
       inf = (BF)tobf(ii, prec);
       sup = (BF)tobf(ss, prec);
       mpfi_init2(e_val,prec);
       mpfr_set(&(e_val->left),  BDY(inf), MPFR_RNDD);
       mpfr_set(&(e_val->right), BDY(sup), MPFR_RNDU);
   
   
       mpfi_log(rv, a_val);
       mpfi_mul(a_val, rv, e_val);
       mpfi_exp(rv, a_val);
   
       MPFRTOBF(&(rv->left), inf);
       MPFRTOBF(&(rv->right), sup);
   
       if ( !cmpbf((Num)inf,0) ) inf = 0;
       if ( !cmpbf((Num)sup,0) ) sup = 0;
   
       if ( inf || sup ) {
         istoitv((Num)inf, (Num)sup, &c);
       } else {
         c = 0;
       }
       *rp = (Num)c;
     //mpfi_clear(rv);
       mpfi_clear(a_val);
       mpfi_clear(e_val);
     }
   }
   
   #endif
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.13
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>