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

Diff for /OpenXM_contrib2/asir2000/engine/d-itv.c between version 1.5 and 1.9

version 1.5, 2015/08/08 14:19:41 version 1.9, 2019/06/04 07:11:22
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2000/engine/d-itv.c,v 1.4 2009/03/27 14:42:29 ohara Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/d-itv.c,v 1.8 2018/03/29 01:32:51 noro Exp $
 */  */
 #if defined(INTERVAL)  #if defined(INTERVAL)
 #include <float.h>  #include <float.h>
 #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"
 #endif  #endif
   #endif
   
 #if defined(ITVDEBUG)  #if defined(ITVDEBUG)
 void printbinint(int d)  void printbinint(int d)
 {  {
         int     i, j, mask;    int  i, j, mask;
         union {    union {
                 int     x;      int  x;
                 char    c[sizeof(int)];      char  c[sizeof(int)];
         } a;    } a;
   
         a.x = d;    a.x = d;
         for(i=sizeof(int)-1;i>=0;i--) {    for(i=sizeof(int)-1;i>=0;i--) {
                 mask = 0x80;      mask = 0x80;
                 for(j=0;j<8;j++) {      for(j=0;j<8;j++) {
                         if (a.c[i] & mask) fprintf(stderr,"1");        if (a.c[i] & mask) fprintf(stderr,"1");
                         else fprintf(stderr,"0");        else fprintf(stderr,"0");
                         mask >>= 1;        mask >>= 1;
                 }      }
         }    }
         fprintf(stderr,"\n");    fprintf(stderr,"\n");
 #if defined(__MINGW32__) || defined(__MINGW64__)  
         fflush(stderr);  
 #endif  
 }  }
 #endif  #endif
   
 double NatToRealUp(N a, int *expo)  double NatToRealUp(N a, int *expo)
 {  {
         int *p;    int *p;
         int l,all,i,j,s,tb,top,tail;    int l,all,i,j,s,tb,top,tail;
         unsigned int t,m[2];    unsigned int t,m[2];
         N       b,c;    N  b,c;
         int     kk, carry, rem;    int  kk, carry, rem;
   
         p = BD(a); l = PL(a) - 1;    p = BD(a); l = PL(a) - 1;
         for ( top = 0, t = p[l]; t; t >>= 1, top++ );    for ( top = 0, t = p[l]; t; t >>= 1, top++ );
         all = top + BSH*l; tail = (53-top)%BSH; j = l-(53-top)/BSH-1;    all = top + BSH*l; tail = (53-top)%BSH; j = l-(53-top)/BSH-1;
   
   
         if ( j >= 0 ) {    if ( j >= 0 ) {
   
 #if defined(ITVDEBUG)  #if defined(ITVDEBUG)
         fprintf(stderr," %d : tail = %d\n", j, tail);    fprintf(stderr," %d : tail = %d\n", j, tail);
 #if defined(__MINGW32__) || defined(__MINGW64__)    printbinint(p[j]);
         fflush(stderr);  
 #endif  #endif
         printbinint(p[j]);      kk = (1<< (BSH - tail)) - 1;
 #endif      rem = p[j] & kk;
                 kk = (1<< (BSH - tail)) - 1;      if ( !rem )
                 rem = p[j] & kk;        for (i=0;i<j;i++)
                 if ( !rem )          if ( p[j] ) {
                         for (i=0;i<j;i++)            carry = 1;
                                 if ( p[j] ) {            break;
                                         carry = 1;          }
                                         break;      else carry = 1;
                                 }      if ( carry ) {
                 else carry = 1;        b = NALLOC(j+1+1);
                 if ( carry ) {        PL(b) = j+1;
                         b = NALLOC(j+1+1);        p = BD(b);
                         PL(b) = j+1;        for(i=0;i<j;i++) p[i] = 0;
                         p = BD(b);        p[j]=(1<< (BSH - tail));
                         for(i=0;i<j;i++) p[i] = 0;  
                         p[j]=(1<< (BSH - tail));  
   
                         addn(a,b,&c);        addn(a,b,&c);
   
                         p = BD(c); l = PL(c) - 1;        p = BD(c); l = PL(c) - 1;
                         for ( top = 0, t = p[l]; t; t >>= 1, top++ );        for ( top = 0, t = p[l]; t; t >>= 1, top++ );
                         all = top + BSH*l; tail = (53-top)%BSH; i = l-(53-top)/BSH-1;        all = top + BSH*l; tail = (53-top)%BSH; i = l-(53-top)/BSH-1;
                 } else i = j;      } else i = j;
         } else i = j;    } else i = j;
   
   
         m[1] = i < 0 ? 0 : p[i]>>(BSH-tail);    m[1] = i < 0 ? 0 : p[i]>>(BSH-tail);
         for ( j = 1, i++, tb = tail; i <= l; i++ ) {    for ( j = 1, i++, tb = tail; i <= l; i++ ) {
                 s = 32-tb; t = i < 0 ? 0 : p[i];      s = 32-tb; t = i < 0 ? 0 : p[i];
                 if ( BSH > s ) {      if ( BSH > s ) {
                         m[j] |= ((t&((1<<s)-1))<<tb);        m[j] |= ((t&((1<<s)-1))<<tb);
                         if ( !j )        if ( !j )
                                 break;          break;
                         else {        else {
                                 j--; m[j] = t>>s; tb = BSH-s;          j--; m[j] = t>>s; tb = BSH-s;
                         }        }
                 } else {      } else {
                         m[j] |= (t<<tb); tb += BSH;        m[j] |= (t<<tb); tb += BSH;
                 }      }
         }    }
         s = (all-1)+1023;    s = (all-1)+1023;
         m[0] = (m[0]&((1<<20)-1))|(MIN(2046,s)<<20); *expo = MAX(s-2046,0);    m[0] = (m[0]&((1<<20)-1))|(MIN(2046,s)<<20); *expo = MAX(s-2046,0);
 #ifdef vax  #ifdef vax
         t = m[0]; m[0] = m[1]; m[1] = t; itod(m);    t = m[0]; m[0] = m[1]; m[1] = t; itod(m);
 #endif  #endif
 #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__)  #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID)
         t = m[0]; m[0] = m[1]; m[1] = t;    t = m[0]; m[0] = m[1]; m[1] = t;
 #endif  #endif
         return *((double *)m);    return *((double *)m);
 }  }
   
 static double   Q2doubleDown(Q a)  static double  Q2doubleDown(Q a)
 {  {
         double nm,dn,t,snm,rd;    double nm,dn,t,snm,rd;
         int enm,edn,e;    int enm,edn,e;
         unsigned int *p,s;    unsigned int *p,s;
   
         nm = (SGN(a)>0)? NatToReal(NM(a),&enm): NatToRealUp(NM(a),&enm);    nm = (SGN(a)>0)? NatToReal(NM(a),&enm): NatToRealUp(NM(a),&enm);
   
         if ( INT(a) )    if ( INT(a) )
                 if ( enm >= 1 )      if ( enm >= 1 )
                         error("Q2doubleDown : Overflow");        error("Q2doubleDown : Overflow");
                 else      else
                         return SGN(a)>0 ? nm : -nm;        return SGN(a)>0 ? nm : -nm;
         else {    else {
                 dn = (SGN(a)>0)? NatToReal(DN(a),&edn): NatToRealUp(DN(a),&edn);      dn = (SGN(a)>0)? NatToReal(DN(a),&edn): NatToRealUp(DN(a),&edn);
   
                 if ( ((e = enm - edn) >= 1024) || (e <= -1023) )      if ( ((e = enm - edn) >= 1024) || (e <= -1023) )
                         error("Q2doubleDown : Overflow"); /* XXX */        error("Q2doubleDown : Overflow"); /* XXX */
                 else {      else {
                         t = 0.0; p = (unsigned int *)&t;        /* Machine */        t = 0.0; p = (unsigned int *)&t;   /* Machine */
                         *p = ((e+1023)<<20);        *p = ((e+1023)<<20);
 #ifdef vax  #ifdef vax
                         s = p[0]; p[0] = p[1]; p[1] = s; itod(p);        s = p[0]; p[0] = p[1]; p[1] = s; itod(p);
 #endif  #endif
 #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__)  #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) | defined(ANDROID)
                         s = p[0]; p[0] = p[1]; p[1] = s;        s = p[0]; p[0] = p[1]; p[1] = s;
 #endif  #endif
                         FPMINUSINF        FPMINUSINF
                         snm = (SGN(a)>0) ? nm : -nm;        snm = (SGN(a)>0) ? nm : -nm;
                         rd  = snm / dn * t;        rd  = snm / dn * t;
                         FPNEAREST        FPNEAREST
                         return rd;        return rd;
                 }      }
         }    }
 }  }
   
   
 static double   Q2doubleUp(Q a)  static double  Q2doubleUp(Q a)
 {  {
         double nm,dn,t,snm,rd;    double nm,dn,t,snm,rd;
         int enm,edn,e;    int enm,edn,e;
         unsigned int *p,s;    unsigned int *p,s;
   
         nm = (SGN(a)>0)? NatToRealUp(NM(a),&enm): NatToReal(NM(a),&enm);    nm = (SGN(a)>0)? NatToRealUp(NM(a),&enm): NatToReal(NM(a),&enm);
   
         if ( INT(a) )    if ( INT(a) )
                 if ( enm >= 1 )      if ( enm >= 1 )
                         error("Q2doubleUp : Overflow");        error("Q2doubleUp : Overflow");
                 else      else
                         return SGN(a)>0 ? nm : -nm;        return SGN(a)>0 ? nm : -nm;
         else {    else {
                 dn = (SGN(a)>0)? NatToRealUp(DN(a),&edn): NatToReal(DN(a),&edn);      dn = (SGN(a)>0)? NatToRealUp(DN(a),&edn): NatToReal(DN(a),&edn);
   
                 if ( ((e = enm - edn) >= 1024) || (e <= -1023) )      if ( ((e = enm - edn) >= 1024) || (e <= -1023) )
                         error("Q2doubleUp : Overflow"); /* XXX */        error("Q2doubleUp : Overflow"); /* XXX */
                 else {      else {
                         t = 0.0; p = (unsigned int *)&t;        /* Machine */        t = 0.0; p = (unsigned int *)&t;   /* Machine */
                         *p = ((e+1023)<<20);        *p = ((e+1023)<<20);
 #ifdef vax  #ifdef vax
                         s = p[0]; p[0] = p[1]; p[1] = s; itod(p);        s = p[0]; p[0] = p[1]; p[1] = s; itod(p);
 #endif  #endif
 #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__)  #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID)
                         s = p[0]; p[0] = p[1]; p[1] = s;        s = p[0]; p[0] = p[1]; p[1] = s;
 #endif  #endif
 #if 0  #if 0
                         FPPLUSINF        FPPLUSINF
                         snm = (SGN(a)>0) ? nm : -nm;        snm = (SGN(a)>0) ? nm : -nm;
                         rd  = snm / dn * t;        rd  = snm / dn * t;
 #endif  #endif
   
                         FPMINUSINF        FPMINUSINF
                         snm = (SGN(a)>0) ? -nm : nm;        snm = (SGN(a)>0) ? -nm : nm;
                         rd  = snm / dn * t;        rd  = snm / dn * t;
                         FPNEAREST        FPNEAREST
                         return (-rd);        return (-rd);
                 }      }
         }    }
 }  }
   
 static double   PARI2doubleDown(BF a)  #if 0
   static double  PARI2doubleDown(BF a)
 {  {
         GEN     p;          //GEN     p;
         double  d;          Num     p;
     double  d;
   
         ritopa(a, &p);          ritopa(a, &p);
         d = rtodbl(p);          d = rtodbl(p);
         cgiv(p);    cgiv(p);
         return d;    return d;
 }  }
   
 static double   PARI2doubleUp(BF a)  static double  PARI2doubleUp(BF a)
 {  {
         return PARI2doubleDown(a);    return PARI2doubleDown(a);
 }  }
   #endif
   
 double  subulpd(double d)  double  subulpd(double d)
 {  {
         double inf;    double inf;
   
         FPMINUSINF    FPMINUSINF
         inf = d - DBL_MIN;    inf = d - DBL_MIN;
         FPNEAREST    FPNEAREST
         return inf;    return inf;
 }  }
   
 double  addulpd(double d)  double  addulpd(double d)
 {  {
         double md, sup;    double md, sup;
   
         md = -d;    md = -d;
         FPMINUSINF    FPMINUSINF
         sup = md - DBL_MIN;    sup = md - DBL_MIN;
         FPNEAREST    FPNEAREST
         return (-sup);    return (-sup);
 }  }
   
 double  ToRealDown(Num a)  double  ToRealDown(Num a)
 {  {
         double inf;    double inf;
   
         if ( ! a ) return 0.0;    if ( ! a ) return 0.0;
         switch ( NID(a) ) {    switch ( NID(a) ) {
                 case N_Q:      case N_Q:
                         inf = Q2doubleDown((Q)a); break;        inf = Q2doubleDown((Q)a); break;
                 case N_R:      case N_R:
                         inf = subulpd(BDY((Real)a)); break;        inf = subulpd(BDY((Real)a)); break;
                 case N_B:      case N_B:
                         inf = PARI2doubleDown((BF)a); break;        //inf = PARI2doubleDown((BF)a); break;
                 case N_IP:          inf = 0;
                         inf = ToRealDown(INF((Itv)a));        error("ToRealDown: not supported operands.");
                         break;          break;
                 case N_IntervalDouble:      case N_IP:
                         inf = INF((IntervalDouble)a); break;        inf = ToRealDown(INF((Itv)a));
                 case N_A:        break;
                 default:      case N_IntervalDouble:
                         inf = 0.0;        inf = INF((IntervalDouble)a); break;
                         error("ToRealDown: not supported operands.");      case N_A:
                         break;      default:
         }        inf = 0.0;
         return inf;        error("ToRealDown: not supported operands.");
         break;
     }
     return inf;
 }  }
   
 double  ToRealUp(Num a)  double  ToRealUp(Num a)
 {  {
         double sup;    double sup;
   
         if ( ! a ) return 0.0;    if ( ! a ) return 0.0;
         switch ( NID(a) ) {    switch ( NID(a) ) {
                 case N_Q:      case N_Q:
                         sup = Q2doubleUp((Q)a); break;        sup = Q2doubleUp((Q)a); break;
                 case N_R:      case N_R:
                         sup = addulpd(BDY((Real)a)); break;        sup = addulpd(BDY((Real)a)); break;
                 case N_B:      case N_B:
                         sup = PARI2doubleUp((BF)a); break;        //sup = PARI2doubleUp((BF)a); break;
                 case N_IP:        sup = 0;
                         sup = ToRealUp(SUP((Itv)a)); break;        error("ToRealUp: not supported operands.");
                 case N_IntervalDouble:        break;
                         sup = SUP((IntervalDouble)a); break;      case N_IP:
                 case N_A:        sup = ToRealUp(SUP((Itv)a)); break;
                 default:      case N_IntervalDouble:
                         sup = 0.0;        sup = SUP((IntervalDouble)a); break;
                         error("ToRealUp: not supported operands.");      case N_A:
                         break;      default:
         }        sup = 0.0;
         return sup;        error("ToRealUp: not supported operands.");
         break;
     }
     return sup;
 }  }
   
   
 void    Num2double(Num a, double *inf, double *sup)  void  Num2double(Num a, double *inf, double *sup)
 {  {
         switch ( NID(a) ) {    switch ( NID(a) ) {
                 case N_Q:      case N_Q:
                         *inf = Q2doubleDown((Q)a);        *inf = Q2doubleDown((Q)a);
                         *sup = Q2doubleUp((Q)a);        *sup = Q2doubleUp((Q)a);
                         break;        break;
                 case N_R:      case N_R:
                         *inf = BDY((Real)a);        *inf = BDY((Real)a);
                         *sup = BDY((Real)a);        *sup = BDY((Real)a);
                         break;        break;
                 case N_B:      case N_B:
                         *inf = PARI2doubleDown((BF)a);        //*inf = PARI2doubleDown((BF)a);
                         *sup = PARI2doubleUp((BF)a);        //*sup = PARI2doubleUp((BF)a);
                         break;        *inf = mpfr_get_d(BDY((BF)a), MPFR_RNDD);
                 case N_IP:        *sup = mpfr_get_d(BDY((BF)a), MPFR_RNDU);
                         *inf = ToRealDown(INF((Itv)a));        break;
                         *sup = ToRealUp(SUP((Itv)a));      case N_IP:
                         break;        *inf = ToRealDown(INF((Itv)a));
                 case N_IntervalDouble:        *sup = ToRealUp(SUP((Itv)a));
                         *inf = INF((IntervalDouble)a);        break;
                         *sup = SUP((IntervalDouble)a);      case N_IntervalDouble:
                         break;        *inf = INF((IntervalDouble)a);
                 case N_A:        *sup = SUP((IntervalDouble)a);
                 default:        break;
                         *inf = 0.0;      case N_A:
                         *sup = 0.0;      default:
                         error("Num2double: not supported operands.");        *inf = 0.0;
                         break;        *sup = 0.0;
         }        error("Num2double: not supported operands.");
         break;
     }
 }  }
   
   
 void additvd(Num a, Num b, IntervalDouble *c)  void additvd(Num a, Num b, IntervalDouble *c)
 {  {
         double  ai,as,mas, bi,bs;    double  ai,as,mas, bi,bs;
         double  inf,sup;    double  inf,sup;
   
         if ( !a ) {    if ( !a ) {
                 *c = (IntervalDouble)b;      *c = (IntervalDouble)b;
         } else if ( !b ) {    } else if ( !b ) {
                 *c = (IntervalDouble)a;      *c = (IntervalDouble)a;
 #if     0  #if  0
         } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {    } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
                 addnum(0,a,b,c);      addnum(0,a,b,c);
 #endif  #endif
         } else {    } else {
   
                 Num2double((Num)a,&ai,&as);      Num2double((Num)a,&ai,&as);
                 Num2double((Num)b,&bi,&bs);      Num2double((Num)b,&bi,&bs);
                 mas = -as;      mas = -as;
                 FPMINUSINF      FPMINUSINF
                 inf = ai + bi;      inf = ai + bi;
                 sup = mas - bs;         /*  as + bs = ( - ( (-as) - bs ) ) */      sup = mas - bs;    /*  as + bs = ( - ( (-as) - bs ) ) */
                 FPNEAREST      FPNEAREST
                 MKIntervalDouble(inf,(-sup),*c);      MKIntervalDouble(inf,(-sup),*c);
         }    }
 }  }
   
 void subitvd(Num a, Num b, IntervalDouble *c)  void subitvd(Num a, Num b, IntervalDouble *c)
 {  {
         double  ai,as,mas, bi,bs;    double  ai,as,mas, bi,bs;
         double  inf,sup;    double  inf,sup;
   
         if ( !a ) {    if ( !a ) {
                 *c = (IntervalDouble)b;      *c = (IntervalDouble)b;
         } else if ( !b ) {    } else if ( !b ) {
                 *c = (IntervalDouble)a;      *c = (IntervalDouble)a;
 #if     0  #if  0
         } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {    } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
                 addnum(0,a,b,c);      addnum(0,a,b,c);
 #endif  #endif
         } else {    } else {
                 Num2double(a,&ai,&as);      Num2double(a,&ai,&as);
                 Num2double(b,&bi,&bs);      Num2double(b,&bi,&bs);
                 FPMINUSINF      FPMINUSINF
                 inf = ai - bs;      inf = ai - bs;
                 sup = ( bi - as );      /* as - bi = ( - ( bi - as ) ) */      sup = ( bi - as );  /* as - bi = ( - ( bi - as ) ) */
                 FPNEAREST      FPNEAREST
                 MKIntervalDouble(inf,(-sup),*c);      MKIntervalDouble(inf,(-sup),*c);
         }    }
 }  }
   
 void    mulitvd(Num a, Num b, IntervalDouble *c)  void  mulitvd(Num a, Num b, IntervalDouble *c)
 {  {
         double  ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p;    double  ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p;
         double  inf, sup;    double  inf, sup;
         int     ba,bb;    int  ba,bb;
   
         if ( !a || !b ) {    if ( !a || !b ) {
                 *c = 0;      *c = 0;
 #if     0  #if  0
         } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {    } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
                 mulnum(0,a,b,c);      mulnum(0,a,b,c);
 #endif  #endif
         } else {    } else {
                 Num2double(a,&ai,&as);      Num2double(a,&ai,&as);
                 Num2double(b,&bi,&bs);      Num2double(b,&bi,&bs);
   
                 FPMINUSINF      FPMINUSINF
   
                 a2 = -as;      a2 = -as;
                 b2 = -bs;      b2 = -bs;
   
                 if ( ba = ( a2 < 0.0 ) ) {      if ( ba = ( a2 < 0.0 ) ) {
                         a1 = ai;        a1 = ai;
                 } else {      } else {
                         a1 = a2;        a1 = a2;
                         a2 = ai;        a2 = ai;
                 }      }
                 if ( bb = ( b2 < 0.0 ) ) {      if ( bb = ( b2 < 0.0 ) ) {
                         b1 = bi;        b1 = bi;
                 } else {      } else {
                         b1 = b2;        b1 = b2;
                         b2 = bi;        b2 = bi;
                 }      }
   
                 c2 = - a2 * b2;      c2 = - a2 * b2;
                 if ( b1 < 0.0 ) {      if ( b1 < 0.0 ) {
                         c1 = - a2 * b1;        c1 = - a2 * b1;
                         if ( a1 < 0.0 ) {        if ( a1 < 0.0 ) {
                                 p = - a1 * b2;          p = - a1 * b2;
                                 if ( p < c1 ) c1 = p;          if ( p < c1 ) c1 = p;
                                 p = - a1 * b1;          p = - a1 * b1;
                                 if ( p < c2 ) c2 = p;          if ( p < c2 ) c2 = p;
                         }        }
                 } else {      } else {
                         c1 = (a1<0.0)?( - a1 * b2 ):( a1 * b1 );        c1 = (a1<0.0)?( - a1 * b2 ):( a1 * b1 );
                 }      }
                 if ( ba == bb ) {      if ( ba == bb ) {
                         inf = c1;        inf = c1;
                         sup = - c2;        sup = - c2;
                 } else {      } else {
                         inf = c2;        inf = c2;
                         sup = - c1;        sup = - c1;
                 }      }
                 FPNEAREST      FPNEAREST
                 MKIntervalDouble(inf,sup,*c);      MKIntervalDouble(inf,sup,*c);
         }    }
 }  }
   
 int     initvd(Num a, IntervalDouble b)  int     initvd(Num a, IntervalDouble b)
 {  {
         Real inf, sup;    Real inf, sup;
   
         if ( NID(b) == N_IntervalDouble ) {    if ( NID(b) == N_IntervalDouble ) {
                 MKReal(INF(b), inf);      MKReal(INF(b), inf);
                 MKReal(SUP(b), sup);      MKReal(SUP(b), sup);
         } else return 0;    } else return 0;
         if ( compnum(0,(Num)inf,a) > 0 ) return 0;    if ( compnum(0,(Num)inf,a) > 0 ) return 0;
         else if ( compnum(0,a,(Num)sup) > 0 ) return 0;    else if ( compnum(0,a,(Num)sup) > 0 ) return 0;
         else return 1;    else return 1;
 }  }
   
 void    divitvd(Num a, Num b, IntervalDouble *c)  void  divitvd(Num a, Num b, IntervalDouble *c)
 {  {
         double  ai,as,bi,bs,a1,a2,b1,b2,c1,c2;    double  ai,as,bi,bs,a1,a2,b1,b2,c1,c2;
         double  inf, sup;    double  inf, sup;
         int     ba,bb;    int  ba,bb;
   
         if ( !b ) {    if ( !b ) {
                 *c = 0;      *c = 0;
                 error("divitvd : division by 0");      error("divitvd : division by 0");
         } else if ( !a ) {    } else if ( !a ) {
                 *c = 0;      *c = 0;
 #if     0  #if  0
         } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {    } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
                 divnum(0,a,b,c);      divnum(0,a,b,c);
 #endif  #endif
         } else {    } else {
                 Num2double(a,&ai,&as);      Num2double(a,&ai,&as);
                 Num2double(b,&bi,&bs);      Num2double(b,&bi,&bs);
   
                 FPMINUSINF      FPMINUSINF
   
                 a2 = -as;      a2 = -as;
                 b2 = -bs;      b2 = -bs;
   
                 if ( ba = ( a2 < 0.0 ) ) {      if ( ba = ( a2 < 0.0 ) ) {
                         a1 = ai;        a1 = ai;
                 } else {      } else {
                         a1 = a2;        a1 = a2;
                         a2 = ai;        a2 = ai;
                 }      }
                 if ( bb = ( b2 < 0.0 ) ) {      if ( bb = ( b2 < 0.0 ) ) {
                         b1 = bi;        b1 = bi;
                 } else {      } else {
                         b1 = b2;        b1 = b2;
                         b2 = bi;        b2 = bi;
                 }      }
   
                 if ( b1 <= 0.0 ) {      if ( b1 <= 0.0 ) {
                         *c = 0;        *c = 0;
                         error("divitvd : division by 0 interval");        error("divitvd : division by 0 interval");
                 } else {      } else {
                         c2 = a2 / b1;        c2 = a2 / b1;
                         c1 = (a1<0.0)?( a1 / b1 ):( - a1 / b2 );        c1 = (a1<0.0)?( a1 / b1 ):( - a1 / b2 );
                 }      }
                 if ( ba == bb ) {      if ( ba == bb ) {
                         inf = c1;        inf = c1;
                         sup = - c2;        sup = - c2;
                 } else {      } else {
                         inf = c2;        inf = c2;
                         sup = - c1;        sup = - c1;
                 }      }
                 FPNEAREST      FPNEAREST
                 MKIntervalDouble(inf,sup,*c);      MKIntervalDouble(inf,sup,*c);
         }    }
 }  }
   
 void    pwritvd(Num a, Num e, IntervalDouble *c)  void  pwritvd(Num a, Num e, IntervalDouble *c)
 {  {
         int     ei;    int  ei;
         IntervalDouble  t;    IntervalDouble  t;
   
         if ( !e ) {    if ( !e ) {
                 *c = (IntervalDouble)ONE;      *c = (IntervalDouble)ONE;
         }  else if ( !a ) {    }  else if ( !a ) {
                 *c = 0;      *c = 0;
 #if     0  #if  0
         } else if ( NID(a) <= N_IP ) {    } else if ( NID(a) <= N_IP ) {
                 pwrnum(0,a,e,c);      pwrnum(0,a,e,c);
 #endif  #endif
         } else if ( !INT(e) ) {    } else if ( !INT(e) ) {
 #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("pwritvd : can't calculate a fractional power");      error("pwritvd : can't calculate a fractional power");
 #endif  #endif
         } else {    } else {
                 ei = QTOS((Q)e);      ei = QTOS((Q)e);
                 pwritv0d((IntervalDouble)a,ei,&t);      if (ei<0) ei = -ei;
                 if ( SGN((Q)e) < 0 )      pwritv0d((IntervalDouble)a,ei,&t);
                         divnum(0,(Num)ONE,(Num)t,(Num *)c);      if ( SGN((Q)e) < 0 )
                 else        divnum(0,(Num)ONE,(Num)t,(Num *)c);
                         *c = t;      else
         }        *c = t;
     }
 }  }
   
 void    pwritv0d(IntervalDouble a, int e, IntervalDouble *c)  void  pwritv0d(IntervalDouble a, int e, IntervalDouble *c)
 {  {
         double inf, sup;    double inf, sup;
         double t, Xmin, Xmax;    double t, Xmin, Xmax;
         int i;    int i;
   
         if ( e == 1 )    if ( e == 1 )
                 *c = a;      *c = a;
         else {    else {
                 if ( !(e%2) ) {      if ( !(e%2) ) {
                         if ( initvd(0,a) ) {        if ( initvd(0,a) ) {
                                 Xmin = 0.0;          Xmin = 0.0;
                                 t = - INF(a);          t = - INF(a);
                                 Xmax = SUP(a);          Xmax = SUP(a);
                                 if ( t > Xmax ) Xmax = t;          if ( t > Xmax ) Xmax = t;
                         } else {        } else {
                                 if ( INF(a) > 0.0 ) {          if ( 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);
                 }      }
                 FPPLUSINF      FPPLUSINF
                 sup = (Xmax!=0.0)?pwrreal0(Xmax,e):0.0;      sup = (Xmax!=0.0)?pwrreal0(Xmax,e):0.0;
                 FPMINUSINF      FPMINUSINF
                 inf = (Xmin!=0.0)?pwrreal0(Xmin,e):0.0;      inf = (Xmin!=0.0)?pwrreal0(Xmin,e):0.0;
                 FPNEAREST      FPNEAREST
                 MKIntervalDouble(inf,sup,*c);      MKIntervalDouble(inf,sup,*c);
         }    }
 }  }
   
 void    chsgnitvd(IntervalDouble a, IntervalDouble *c)  void  chsgnitvd(IntervalDouble a, IntervalDouble *c)
 {  {
         double inf,sup;    double inf,sup;
   
         if ( !a )    if ( !a )
                 *c = 0;      *c = 0;
 #if     0  #if  0
         else if ( NID(a) <= N_IP )    else if ( NID(a) <= N_IP )
                 chsgnnum(a,c);      chsgnnum(a,c);
 #endif  #endif
         else {    else {
                 inf = - SUP(a);      inf = - SUP(a);
                 sup = - INF(a);      sup = - INF(a);
                 MKIntervalDouble(inf,sup,*c);      MKIntervalDouble(inf,sup,*c);
         }    }
 }  }
   
 int     cmpitvd(IntervalDouble a, IntervalDouble b)  int  cmpitvd(IntervalDouble a, IntervalDouble b)
 /*    a > b:  1       */  /*    a > b:  1       */
 /*    a = b:  0       */  /*    a = b:  0       */
 /*    a < b: -1       */  /*    a < b: -1       */
 {  {
         double  ai,as,bi,bs;    double  ai,as,bi,bs;
   
         if ( !a ) {    if ( !a ) {
                 if ( !b || (NID(b)<=N_IP) ) {      if ( !b || (NID(b)<=N_IP) ) {
                         return compnum(0,(Num)a,(Num)b);        return compnum(0,(Num)a,(Num)b);
                 } else {      } else {
                         bi = INF(b);        bi = INF(b);
                         bs = SUP(b);        bs = SUP(b);
                         if ( bi > 0.0 ) return -1;        if ( bi > 0.0 ) return -1;
                         else if ( bs < 0.0 ) return 1;        else if ( bs < 0.0 ) return 1;
                         else return 0;        else return 0;
                 }      }
         } else if ( !b ) {    } else if ( !b ) {
                 if ( !a || (NID(a)<=N_IP) ) {      if ( !a || (NID(a)<=N_IP) ) {
                         return compnum(0,(Num)a,(Num)b);        return compnum(0,(Num)a,(Num)b);
                 } else {      } else {
                         ai = INF(a);        ai = INF(a);
                         as = SUP(a);        as = SUP(a);
                         if ( ai > 0.0 ) return 1;        if ( ai > 0.0 ) return 1;
                         else if ( as < 0.0 ) return -1;        else if ( as < 0.0 ) return -1;
                         else return 0;        else return 0;
                 }      }
         } else {    } else {
                 bi = INF(b);      bi = INF(b);
                 bs = SUP(b);      bs = SUP(b);
                 ai = INF(a);      ai = INF(a);
                 as = SUP(a);      as = SUP(a);
                 if ( ai > bs ) return 1;      if ( ai > bs ) return 1;
                 else if ( bi > as ) return -1;      else if ( bi > as ) return -1;
                 else return 0;      else return 0;
         }    }
 }  }
   
 void    miditvd(IntervalDouble a, Num *b)  void  miditvd(IntervalDouble a, Num *b)
 {  {
         double  t;    double  t;
         Real    rp;    Real  rp;
   
         if ( ! a ) *b = 0;    if ( ! a ) *b = 0;
         else if ( (NID(a) <= N_IP) )    else if ( (NID(a) <= N_IP) )
                 *b = (Num)a;      *b = (Num)a;
         else {    else {
                 t = (INF(a) + SUP(a))/2.0;      t = (INF(a) + SUP(a))/2.0;
                 MKReal(t,rp);      MKReal(t,rp);
                 *b = (Num)rp;      *b = (Num)rp;
         }    }
 }  }
   
 void    cupitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)  void  cupitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
 {  {
         double  ai,as,bi,bs;    double  ai,as,bi,bs;
         double  inf,sup;    double  inf,sup;
   
         ai = INF(a);    ai = INF(a);
         as = SUP(a);    as = SUP(a);
         bi = INF(b);    bi = INF(b);
         bs = SUP(b);    bs = SUP(b);
         inf = MIN(ai,bi);    inf = MIN(ai,bi);
         sup = MAX(as,bs);    sup = MAX(as,bs);
         MKIntervalDouble(inf,sup,*c);    MKIntervalDouble(inf,sup,*c);
 }  }
   
 void    capitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)  void  capitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
 {  {
         double  ai,as,bi,bs;    double  ai,as,bi,bs;
         double  inf,sup;    double  inf,sup;
   
         ai = INF(a);    ai = INF(a);
         as = SUP(a);    as = SUP(a);
         bi = INF(b);    bi = INF(b);
         bs = SUP(b);    bs = SUP(b);
         inf = MAX(ai,bi);    inf = MAX(ai,bi);
         sup = MIN(as,bs);    sup = MIN(as,bs);
         if ( inf > sup ) *c = 0;    if ( inf > sup ) *c = 0;
         else MKIntervalDouble(inf,sup,*c);    else MKIntervalDouble(inf,sup,*c);
 }  }
   
   
 void    widthitvd(IntervalDouble a, Num *b)  void  widthitvd(IntervalDouble a, Num *b)
 {  {
         double  t;    double  t;
         Real    rp;    Real  rp;
   
         if ( ! a ) *b = 0;    if ( ! a ) *b = 0;
         else {    else {
                 t = SUP(a) - INF(a);      t = SUP(a) - INF(a);
                 MKReal(t,rp);      MKReal(t,rp);
                 *b = (Num)rp;      *b = (Num)rp;
         }    }
 }  }
   
 void    absitvd(IntervalDouble a, Num *b)  void  absitvd(IntervalDouble a, Num *b)
 {  {
         double  ai,as,bi,bs;    double  ai,as,bi,bs;
         double  t;    double  t;
         Real    rp;    Real  rp;
   
         if ( ! a ) *b = 0;    if ( ! a ) *b = 0;
         else {    else {
                 ai = INF(a);      ai = INF(a);
                 as = SUP(a);      as = SUP(a);
                 if (ai < 0) bi = -ai;   else    bi = ai;      if (ai < 0) bi = -ai;  else  bi = ai;
                 if (as < 0) bs = -as;   else    bs = as;      if (as < 0) bs = -as;  else  bs = as;
                 if ( bi > bs )  t = bi; else    t = bs;      if ( bi > bs )  t = bi; else  t = bs;
                 MKReal(t,rp);      MKReal(t,rp);
                 *b = (Num)rp;      *b = (Num)rp;
         }    }
 }  }
   
 void    distanceitvd(IntervalDouble a, IntervalDouble b, Num *c)  void  distanceitvd(IntervalDouble a, IntervalDouble b, Num *c)
 {  {
         double  ai,as,bi,bs;    double  ai,as,bi,bs;
         double  di,ds;    double  di,ds;
         double  t;    double  t;
         Real    rp;    Real  rp;
   
         ai = INF(a);    ai = INF(a);
         as = SUP(a);    as = SUP(a);
         bi = INF(b);    bi = INF(b);
         bs = SUP(b);    bs = SUP(b);
         di = bi - ai;    di = bi - ai;
         ds = bs - as;    ds = bs - as;
   
         if (di < 0) di = -di;    if (di < 0) di = -di;
         if (ds < 0) ds = -ds;    if (ds < 0) ds = -ds;
         if (di > ds)    t = di; else    t = ds;    if (di > ds)  t = di; else  t = ds;
         MKReal(t,rp);    MKReal(t,rp);
         *c = (Num)rp;    *c = (Num)rp;
 }  }
   
 #endif  #endif
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.9
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>