[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.1 and 1.5

version 1.1, 2000/12/22 10:03:28 version 1.5, 2015/08/08 14:19:41
Line 1 
Line 1 
 /*  /*
  * $OpenXM: $   * $OpenXM: OpenXM_contrib2/asir2000/engine/d-itv.c,v 1.4 2009/03/27 14:42:29 ohara 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 PARI  #if defined(PARI)
 #include "genpari.h"  #include "genpari.h"
 #endif  #endif
   
Line 28  void printbinint(int d)
Line 28  void printbinint(int d)
                 }                  }
         }          }
         fprintf(stderr,"\n");          fprintf(stderr,"\n");
   #if defined(__MINGW32__) || defined(__MINGW64__)
           fflush(stderr);
   #endif
 }  }
 #endif  #endif
   
Line 48  double NatToRealUp(N a, int *expo)
Line 51  double NatToRealUp(N a, int *expo)
   
 #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__)
           fflush(stderr);
   #endif
         printbinint(p[j]);          printbinint(p[j]);
 #endif  #endif
                 kk = (1<< (BSH - tail)) - 1;                  kk = (1<< (BSH - tail)) - 1;
Line 231  double ToRealDown(Num a)
Line 237  double ToRealDown(Num a)
                 case N_IP:                  case N_IP:
                         inf = ToRealDown(INF((Itv)a));                          inf = ToRealDown(INF((Itv)a));
                         break;                          break;
                 case N_ID:                  case N_IntervalDouble:
                         inf = INF((ItvD)a); break;                          inf = INF((IntervalDouble)a); break;
                 case N_A:                  case N_A:
                 default:                  default:
                         inf = 0.0;                          inf = 0.0;
Line 256  double ToRealUp(Num a)
Line 262  double ToRealUp(Num a)
                         sup = PARI2doubleUp((BF)a); break;                          sup = PARI2doubleUp((BF)a); break;
                 case N_IP:                  case N_IP:
                         sup = ToRealUp(SUP((Itv)a)); break;                          sup = ToRealUp(SUP((Itv)a)); break;
                 case N_ID:                  case N_IntervalDouble:
                         sup = SUP((ItvD)a); break;                          sup = SUP((IntervalDouble)a); break;
                 case N_A:                  case N_A:
                 default:                  default:
                         sup = 0.0;                          sup = 0.0;
Line 287  void Num2double(Num a, double *inf, double *sup)
Line 293  void Num2double(Num a, double *inf, double *sup)
                         *inf = ToRealDown(INF((Itv)a));                          *inf = ToRealDown(INF((Itv)a));
                         *sup = ToRealUp(SUP((Itv)a));                          *sup = ToRealUp(SUP((Itv)a));
                         break;                          break;
                 case N_ID:                  case N_IntervalDouble:
                         *inf = INF((ItvD)a);                          *inf = INF((IntervalDouble)a);
                         *sup = SUP((ItvD)a);                          *sup = SUP((IntervalDouble)a);
                         break;                          break;
                 case N_A:                  case N_A:
                 default:                  default:
Line 301  void Num2double(Num a, double *inf, double *sup)
Line 307  void Num2double(Num a, double *inf, double *sup)
 }  }
   
   
 void additvd(Num a, Num b, ItvD *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 = (ItvD)b;                  *c = (IntervalDouble)b;
         } else if ( !b ) {          } else if ( !b ) {
                 *c = (ItvD)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);
Line 323  void additvd(Num a, Num b, ItvD *c)
Line 329  void additvd(Num a, Num b, ItvD *c)
                 inf = ai + bi;                  inf = ai + bi;
                 sup = mas - bs;         /*  as + bs = ( - ( (-as) - bs ) ) */                  sup = mas - bs;         /*  as + bs = ( - ( (-as) - bs ) ) */
                 FPNEAREST                  FPNEAREST
                 MKItvD(inf,(-sup),*c);                  MKIntervalDouble(inf,(-sup),*c);
         }          }
 }  }
   
 void subitvd(Num a, Num b, ItvD *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 = (ItvD)b;                  *c = (IntervalDouble)b;
         } else if ( !b ) {          } else if ( !b ) {
                 *c = (ItvD)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);
Line 347  void subitvd(Num a, Num b, ItvD *c)
Line 353  void subitvd(Num a, Num b, ItvD *c)
                 inf = ai - bs;                  inf = ai - bs;
                 sup = ( bi - as );      /* as - bi = ( - ( bi - as ) ) */                  sup = ( bi - as );      /* as - bi = ( - ( bi - as ) ) */
                 FPNEAREST                  FPNEAREST
                 MKItvD(inf,(-sup),*c);                  MKIntervalDouble(inf,(-sup),*c);
         }          }
 }  }
   
 void    mulitvd(Num a, Num b, ItvD *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;
Line 405  void mulitvd(Num a, Num b, ItvD *c)
Line 411  void mulitvd(Num a, Num b, ItvD *c)
                         sup = - c1;                          sup = - c1;
                 }                  }
                 FPNEAREST                  FPNEAREST
                 MKItvD(inf,sup,*c);                  MKIntervalDouble(inf,sup,*c);
         }          }
 }  }
   
 int     initvd(Num a, ItvD b)  int     initvd(Num a, IntervalDouble b)
 {  {
         Real inf, sup;          Real inf, sup;
   
         if ( NID(b) == N_ID ) {          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;
Line 422  int     initvd(Num a, ItvD b)
Line 428  int     initvd(Num a, ItvD b)
         else return 1;          else return 1;
 }  }
   
 void    divitvd(Num a, Num b, ItvD *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;
Line 474  void divitvd(Num a, Num b, ItvD *c)
Line 480  void divitvd(Num a, Num b, ItvD *c)
                         sup = - c1;                          sup = - c1;
                 }                  }
                 FPNEAREST                  FPNEAREST
                 MKItvD(inf,sup,*c);                  MKIntervalDouble(inf,sup,*c);
         }          }
 }  }
   
 void    pwritvd(Num a, Num e, ItvD *c)  void    pwritvd(Num a, Num e, IntervalDouble *c)
 {  {
         int     ei;          int     ei;
         ItvD    t;          IntervalDouble  t;
   
         if ( !e ) {          if ( !e ) {
                 *c = (ItvD)ONE;                  *c = (IntervalDouble)ONE;
         }  else if ( !a ) {          }  else if ( !a ) {
                 *c = 0;                  *c = 0;
 #if     0  #if     0
Line 492  void pwritvd(Num a, Num e, ItvD *c)
Line 498  void pwritvd(Num a, Num e, ItvD *c)
                 pwrnum(0,a,e,c);                  pwrnum(0,a,e,c);
 #endif  #endif
         } else if ( !INT(e) ) {          } else if ( !INT(e) ) {
 #if 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("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((ItvD)a,ei,&t);                  pwritv0d((IntervalDouble)a,ei,&t);
                 if ( SGN((Q)e) < 0 )                  if ( SGN((Q)e) < 0 )
                         divnum(0,(Num)ONE,(Num)t,(Num *)c);                          divnum(0,(Num)ONE,(Num)t,(Num *)c);
                 else                  else
Line 512  void pwritvd(Num a, Num e, ItvD *c)
Line 513  void pwritvd(Num a, Num e, ItvD *c)
         }          }
 }  }
   
 void    pwritv0d(ItvD a, int e, ItvD *c)  void    pwritv0d(IntervalDouble a, int e, IntervalDouble *c)
 {  {
         double inf, sup;          double inf, sup;
         double t, Xmin, Xmax;          double t, Xmin, Xmax;
Line 545  void pwritv0d(ItvD a, int e, ItvD *c)
Line 546  void pwritv0d(ItvD a, int e, ItvD *c)
                 FPMINUSINF                  FPMINUSINF
                 inf = (Xmin!=0.0)?pwrreal0(Xmin,e):0.0;                  inf = (Xmin!=0.0)?pwrreal0(Xmin,e):0.0;
                 FPNEAREST                  FPNEAREST
                 MKItvD(inf,sup,*c);                  MKIntervalDouble(inf,sup,*c);
         }          }
 }  }
   
 void    chsgnitvd(ItvD a, ItvD *c)  void    chsgnitvd(IntervalDouble a, IntervalDouble *c)
 {  {
         double inf,sup;          double inf,sup;
   
Line 562  void chsgnitvd(ItvD a, ItvD *c)
Line 563  void chsgnitvd(ItvD a, ItvD *c)
         else {          else {
                 inf = - SUP(a);                  inf = - SUP(a);
                 sup = - INF(a);                  sup = - INF(a);
                 MKItvD(inf,sup,*c);                  MKIntervalDouble(inf,sup,*c);
         }          }
 }  }
   
 int     cmpitvd(ItvD a, ItvD b)  int     cmpitvd(IntervalDouble a, IntervalDouble b)
   /*    a > b:  1       */
   /*    a = b:  0       */
   /*    a < b: -1       */
 {  {
         double  ai,as,bi,bs;          double  ai,as,bi,bs;
   
Line 601  int cmpitvd(ItvD a, ItvD b)
Line 605  int cmpitvd(ItvD a, ItvD b)
         }          }
 }  }
   
 void    miditvd(ItvD a, Num *b)  void    miditvd(IntervalDouble a, Num *b)
 {  {
         double  t;          double  t;
         Real    rp;          Real    rp;
Line 616  void miditvd(ItvD a, Num *b)
Line 620  void miditvd(ItvD a, Num *b)
         }          }
 }  }
   
 void    cupitvd(ItvD a, ItvD b, ItvD *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;
Line 627  void cupitvd(ItvD a, ItvD b, ItvD *c)
Line 631  void cupitvd(ItvD a, ItvD b, ItvD *c)
         bs = SUP(b);          bs = SUP(b);
         inf = MIN(ai,bi);          inf = MIN(ai,bi);
         sup = MAX(as,bs);          sup = MAX(as,bs);
         MKItvD(inf,sup,*c);          MKIntervalDouble(inf,sup,*c);
 }  }
   
 void    capitvd(ItvD a, ItvD b, ItvD *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;
Line 642  void capitvd(ItvD a, ItvD b, ItvD *c)
Line 646  void capitvd(ItvD a, ItvD b, ItvD *c)
         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 MKItvD(inf,sup,*c);          else MKIntervalDouble(inf,sup,*c);
 }  }
   
   
 void    widthitvd(ItvD a, Num *b)  void    widthitvd(IntervalDouble a, Num *b)
 {  {
         double  t;          double  t;
         Real    rp;          Real    rp;
Line 659  void widthitvd(ItvD a, Num *b)
Line 663  void widthitvd(ItvD a, Num *b)
         }          }
 }  }
   
 void    absitvd(ItvD a, Num *b)  void    absitvd(IntervalDouble a, Num *b)
 {  {
         double  ai,as,bi,bs;          double  ai,as,bi,bs;
         double  t;          double  t;
Line 677  void absitvd(ItvD a, Num *b)
Line 681  void absitvd(ItvD a, Num *b)
         }          }
 }  }
   
 void    distanceitvd(ItvD a, ItvD 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;
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.5
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>