Return to itvnum.c CVS log

Up to [local] / OpenXM_contrib2 / asir2018 / builtin

fix interval arithmetic

/*
 * $OpenXM: OpenXM_contrib2/asir2018/builtin/itvnum.c,v 1.4 2019/10/17 03:03:12 kondoh Exp $
 */

#include "ca.h"
#include "parse.h"
#include "version.h"
#if !defined(ANDROID)
#include "../plot/ifplot.h"
#endif

// in engine/bf.c
Num tobf(Num,int);

#if defined(INTERVAL)
static void Pitv(NODE, Obj *);
static void Pitvd(NODE, Obj *);
static void Pitvbf(NODE, Obj *);
static void Pinf(NODE, Obj *);
static void Psup(NODE, Obj *);
static void Pmid(NODE, Obj *);
static void Pabsitv(NODE, Obj *);
static void Pdisjitv(NODE, Obj *);
static void Pinitv(NODE, Obj *);
static void Pcup(NODE, Obj *);
static void Pcap(NODE, Obj *);
static void Pwidth(NODE, Obj *);
static void Pdistance(NODE, Obj *);
static void Pitvversion(NODE, Q *);
static void PzeroRewriteMode(NODE, Obj *);
static void PzeroRewriteCountClear(NODE, Obj *);
static void PzeroRewriteCount(NODE, Obj *);
//void miditvp(Itv,Num *);
//void absitvp(Itv,Num *);
//int initvd(Num,IntervalDouble);
//int initvp(Num,Itv);
//int itvinitvp(Itv,Itv);
#endif
static void Pprintmode(NODE, Obj *);

/* plot time check func */
static void ccalc(double **, struct canvas *, int);
static void Pifcheck(NODE, Obj *);

#if  defined(__osf__) && 0
int  end;
#endif

struct ftab interval_tab[] = {
  {"printmode",Pprintmode,1},
#if defined(INTERVAL)
  {"itvd",Pitvd,-2},
  {"intvald",Pitvd,-2},
  {"itv",Pitv,-2},
  {"intval",Pitv,-2},
  {"itvbf",Pitvbf,-2},
  {"intvalbf",Pitvbf,-2},
  {"inf",Pinf,1},
  {"sup",Psup,1},
  {"absintval",Pabsitv,1},
  {"disintval",Pdisjitv,2},
  {"inintval",Pinitv,2},
  {"cup",Pcup,2},
  {"cap",Pcap,2},
  {"mid",Pmid,1},
  {"width",Pwidth,1},
  {"diam",Pwidth,1},
  {"distance",Pdistance,2},
  {"iversion",Pitvversion,-1},
  {"intvalversion",Pitvversion,-1},
  {"zerorewritemode",PzeroRewriteMode,-1},
  {"zeroRewriteMode",PzeroRewriteMode,-1},
  {"zeroRewriteCountClear",PzeroRewriteCountClear,-1},
  {"zeroRewriteCount",PzeroRewriteCount,-1},
/* plot time check */
  {"ifcheck",Pifcheck,-7},
#endif
  {0,0,0},
};

extern int mpfr_roundmode;

#if defined(INTERVAL)

/* plot time check */
static void
Pifcheck(NODE arg, Obj *rp)
{
  Z m2,p2,s_id;
  NODE defrange;
  LIST xrange,yrange,range[2],list,geom;
  VL vl,vl0;
  V v[2],av[2];
  int ri,i,j,sign;
  P poly;
  P var;
  NODE n,n0;
  Obj t;

  struct canvas *can;
  MAT m;
  pointer **mb;
  double **tabe, *px, *px1, *px2;
  Z one;
  int width, height, ix, iy;
  int id;

  STOZ(-2,m2); STOZ(2,p2);
  STOZ(1,one);
  MKNODE(n,p2,0); MKNODE(defrange,m2,n);
  poly = 0; vl = 0; geom = 0; ri = 0;
  v[0] = v[1] = 0;
  for ( ; arg; arg = NEXT(arg) ){
    switch ( OID(BDY(arg)) ) {
    case O_P:
      poly = (P)BDY(arg);
      get_vars_recursive((Obj)poly,&vl);
      for(vl0=vl,i=0;vl0;vl0=NEXT(vl0)){
        if(vl0->v->attr==(pointer)V_IND){
          if(i>=2){
            error("ifplot : invalid argument");
          } else {
            v[i++]=vl0->v;
          }
        }
      }
      break;
    case O_LIST:
      list = (LIST)BDY(arg);
      if ( OID(BDY(BDY(list))) == O_P )
        if ( ri > 1 )
          error("ifplot : invalid argument");
        else
          range[ri++] = list;
      else
        geom = list;
      break;
    default:
      error("ifplot : invalid argument"); break;
    }
  }
  if ( !poly ) error("ifplot : invalid argument");
  switch ( ri ) {
    case 0:
      if ( !v[1] ) error("ifplot : please specify all variables");
      MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n);
      MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n);
      break;
    case 1:
      if ( !v[1] ) error("ifplot : please specify all variables");
      av[0] = VR((P)BDY(BDY(range[0])));
      if ( v[0] == av[0] ) {
        xrange = range[0];
        MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n);
      } else if ( v[1] == av[0] ) {
        MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n);
        yrange = range[0];
      } else
        error("ifplot : invalid argument");
      break;
    case 2:
      av[0] = VR((P)BDY(BDY(range[0])));
      av[1] = VR((P)BDY(BDY(range[1])));
      if ( ((v[0] == av[0]) && (!v[1] || v[1] == av[1])) ||
         ((v[0] == av[1]) && (!v[1] || v[1] == av[0])) ) {
          xrange = range[0]; yrange = range[1];
      } else error("ifplot : invalid argument");
      break;
    default:
      error("ifplot : cannot happen"); break;
  }
  can = canvas[id = search_canvas()];
  if ( !geom ) {
    width = 300;
    height = 300;
    can->width = 300;
    can->height = 300;
  } else {
    can->width = ZTOS((Z)BDY(BDY(geom)));
    can->height = ZTOS((Z)BDY(NEXT(BDY(geom))));
    width = can->width;
    height = can->height;
  }
  if ( xrange ) {
    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->xmin = ToReal(can->qxmin); can->xmax = ToReal(can->qxmax); 
  }
  if ( yrange ) {
    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->ymin = ToReal(can->qymin); can->ymax = ToReal(can->qymax); 
  }
  can->wname = "ifcheck";
  can->formula = poly;
  tabe = (double **)ALLOCA((width+1)*sizeof(double *));
  for ( i = 0; i <= width; i++ )
    tabe[i] = (double *)ALLOCA((height+1)*sizeof(double));
  for(i=0;i<=width;i++)for(j=0;j<=height;j++)tabe[i][j]=0;
  ccalc(tabe,can,0);
  MKMAT(m,width,height);
  mb = BDY(m);
  for( ix=0; ix<width; ix++ ){
    for( iy=0; iy<height; iy++){
      if ( tabe[ix][iy] >= 0 ){
        if ( (tabe[ix+1][iy] <= 0) ||
          (tabe[ix][iy+1] <= 0 ) ||
          (tabe[ix+1][iy+1] <= 0 ) ) mb[ix][iy] = (Obj)one;
      } else {
        if ( (tabe[ix+1][iy] >= 0 ) ||
          ( tabe[ix][iy+1] >= 0 ) ||
          ( tabe[ix+1][iy+1] >= 0 )) mb[ix][iy] = (Obj)one;
      }
    }
  }
  *rp = (Obj)m;
}

void ccalc(double **tab,struct canvas *can,int nox)
{
   double x,y,xmin,ymin,xstep,ystep;
   int ix,iy;
   Real r,rx,ry;
   Obj fr,g;
   int w,h;
   V vx,vy;
   Obj t,s;

   MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
   vx = can->vx;
   vy = can->vy;
   w = can->width; h = can->height; 
   xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;
   ymin = can->ymin; ystep = (can->ymax-can->ymin)/h;
   MKReal(1.0,rx); MKReal(1.0,ry);
   for( ix = 0, x = xmin; ix < w+1 ; ix++, x += xstep ) {
      BDY(rx) = x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);
      devalr(CO,t,&g);
      for( iy = 0, y = ymin; iy < h+1 ; iy++, y += ystep ) {
         BDY(ry) = y;
         substr(CO,0,g,vy,y?(Obj)ry:0,&t);
         devalr(CO,t,&s);
         tab[ix][iy] = ToReal(s);
      }
   }
}
/* end plot time check */

static void
Pitvversion(NODE arg, Q *rp)
{
	Z r;
  STOZ(INT_ASIR_VERSION, r);
	*rp = (Q)r;
}

extern int  bigfloat;

static void
Pitv(NODE arg, Obj *rp)
{
  Num  a, i, s;
  Itv  c;
  double  inf, sup;

#if 1
  if ( bigfloat )
    Pitvbf(arg, rp);
  else
    Pitvd(arg,rp);
#else
  asir_assert(ARG0(arg),O_N,"itv");
  if ( argc(arg) > 1 ) {
    asir_assert(ARG1(arg),O_N,"itv");
    istoitv((Num)ARG0(arg),(Num)ARG1(arg),&c);
  } else {
    a = (Num)ARG0(arg);
    if ( ! a ) {
      *rp = 0;
      return;
    }
    else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat) {
      *rp = (Obj)a;
      return;
    }
    else if ( NID(a) == N_IntervalDouble ) {
      inf = INF((IntervalDouble)a);
      sup = SUP((IntervalDouble)a);
      double2bf(inf, (BF *)&i);
      double2bf(sup, (BF *)&s);
      istoitv(i,s,&c);
    }
    else istoitv(a,a,&c);
  }
  if ( NID( c ) == N_IntervalBigFloat )
    addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp);
  else *rp = (Obj)c;
#endif
}

static void
Pitvbf(NODE arg, Obj *rp)
{
  Num  a, i, s;
  IntervalBigFloat  c;
  Num  ii,ss;
  Real di, ds;
  double  inf, sup;
  int current_roundmode;

  asir_assert(ARG0(arg),O_N,"intvalbf");
  a = (Num)ARG0(arg);
  if ( argc(arg) > 1 ) {
    asir_assert(ARG1(arg),O_N,"intvalbf");

    i = (Num)ARG0(arg);
    s = (Num)ARG1(arg);
    current_roundmode = mpfr_roundmode;
    mpfr_roundmode = MPFR_RNDD;
    ii = tobf(i, DEFAULTPREC);
    mpfr_roundmode = MPFR_RNDU;
    ss = tobf(s, DEFAULTPREC);
    istoitv(ii,ss,(Itv *)&c);
//    MKIntervalBigFloat((BF)ii,(BF)ss,c);
//    ToBf(s, &ss);
    mpfr_roundmode = current_roundmode;
  } else {
    if ( ! a ) {
      *rp = 0;
      return;
    }
    else if ( NID(a) == N_IP ) {
      itvtois((Itv)a, &i, &s);
      current_roundmode = mpfr_roundmode;
      mpfr_roundmode = MPFR_RNDD;
      ii = tobf(i, DEFAULTPREC);
      mpfr_roundmode = MPFR_RNDU;
      ss = tobf(s, DEFAULTPREC);
      istoitv(ii,ss,(Itv *)&c);
//      MKIntervalBigFloat((BF)ii,(BF)ss,c);
      mpfr_roundmode = current_roundmode;
    }
    else if ( NID(a) == N_IntervalBigFloat) {
      *rp = (Obj)a;
      return;
    }
    else if ( NID(a) == N_IntervalDouble ) {
      inf = INF((IntervalDouble)a);
      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
Pitvd(NODE arg, Obj *rp)
{
  double  inf, sup;
  Num  a, a0, a1, t;
  Itv  ia;
  IntervalDouble  d;

  asir_assert(ARG0(arg),O_N,"intvald");
  a0 = (Num)ARG0(arg);
  if ( argc(arg) > 1 ) {
    asir_assert(ARG1(arg),O_N,"intvald");
    a1 = (Num)ARG1(arg);
  } else {
    if ( a0 && OID(a0)==O_N && NID(a0)==N_IntervalDouble ) {
      inf = INF((IntervalDouble)a0);
      sup = SUP((IntervalDouble)a0);
      MKIntervalDouble(inf,sup,d);
      *rp = (Obj)d;
      return;
    }
    a1 = (Num)ARG0(arg);
  }
  if ( compnum(0,a0,a1) > 0 ) {
    t = a0; a0 = a1; a1 = t;
  }
  inf = ToRealDown(a0);
  sup = ToRealUp(a1);
  MKIntervalDouble(inf,sup,d);
  *rp = (Obj)d;
}

static void
Pinf(NODE arg, Obj *rp)
{
  Num  a, i, s;
  Real  r;
  double  d;

  a = (Num)ARG0(arg);
  if ( ! a ) {
    *rp = 0;
  } else if  ( OID(a) == O_N ) {
    switch ( NID(a) ) {
      case N_IntervalDouble:
        d = INF((IntervalDouble)a);
        MKReal(d, r);
        *rp = (Obj)r;
        break;
      case N_IP:
      case N_IntervalBigFloat:
      case N_IntervalQuad:
        itvtois((Itv)ARG0(arg),&i,&s);
        *rp = (Obj)i;
        break;
      default:
        *rp = (Obj)a;
        break;
    }
  } else {
    *rp = (Obj)a;
  }
}

static void
Psup(NODE arg, Obj *rp)
{
  Num  a, i, s;
  Real  r;
  double  d;

  a = (Num)ARG0(arg);
  if ( ! a ) {
    *rp = 0;
  } else if  ( OID(a) == O_N ) {
    switch ( NID(a) ) {
      case N_IntervalDouble:
        d = SUP((IntervalDouble)a);
        MKReal(d, r);
        *rp = (Obj)r;
        break;
      case N_IP:
      case N_IntervalBigFloat:
      case N_IntervalQuad:
        itvtois((Itv)ARG0(arg),&i,&s);
        *rp = (Obj)s;
        break;
      default:
        *rp = (Obj)a;
        break;
    }
  } else {
      *rp = (Obj)a;
  }
}

static void
Pmid(NODE arg, Obj *rp)
{
  Num  a, s;
  Real  r;
  double  d;

  a = (Num)ARG0(arg);
  if ( ! a ) {
    *rp = 0;
  } else switch (OID(a)) {
    case O_N:
      if ( NID(a) == N_IntervalDouble ) {
        d = ( INF((IntervalDouble)a)+SUP((IntervalDouble)a) ) / 2.0;
        MKReal(d, r);
        *rp = (Obj)r;
      } else if ( NID(a) == N_IntervalQuad ) {
        error("mid: not supported operation");
        *rp = 0;
      } else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat ) {
        miditvp((Itv)ARG0(arg),&s);
        *rp = (Obj)s;
      } else {
        *rp = (Obj)a;
      }
      break;
#if 0
    case O_P:
    case O_R:
    case O_LIST:
    case O_VECT:
    case O_MAT:
#endif
    default:
      *rp = (Obj)a;
      break;
  }
}

static void
Pcup(NODE arg, Obj *rp)
{
  Itv  s;
  Num  a, b;

  asir_assert(ARG0(arg),O_N,"cup");
  asir_assert(ARG1(arg),O_N,"cup");
  a = (Num)ARG0(arg);
  b = (Num)ARG1(arg);
  if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {
    cupitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp);
  } else {
    cupitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);
    *rp = (Obj)s;
  }
}

static void
Pcap(NODE arg, Obj *rp)
{
  Itv  s;
  Num  a, b;

  asir_assert(ARG0(arg),O_N,"cap");
  asir_assert(ARG1(arg),O_N,"cap");
  a = (Num)ARG0(arg);
  b = (Num)ARG1(arg);
  if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {
    capitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp);
  } else {
    capitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);
    *rp = (Obj)s;
  }
}

static void
Pwidth(arg,rp)
NODE arg;
Obj *rp;
{
  Num  s;
  Num  a;

  asir_assert(ARG0(arg),O_N,"width");
  a = (Num)ARG0(arg);
  if ( ! a ) {
    *rp = 0;
  } else if ( NID(a) == N_IntervalDouble ) {
    widthitvd((IntervalDouble)a, (Num *)rp);
  } else {
    widthitvp((Itv)ARG0(arg),&s);
    *rp = (Obj)s;
  }
}

static void
Pabsitv(arg,rp)
NODE arg;
Obj *rp;
{
  Num  s;
  Num  a, b;

  asir_assert(ARG0(arg),O_N,"absitv");
  a = (Num)ARG0(arg);
  if ( ! a ) {
    *rp = 0;
  } else if ( NID(a) == N_IntervalDouble ) {
    absitvd((IntervalDouble)a, (Num *)rp);
  } else {
    absitvp((Itv)ARG0(arg),&s);
    *rp = (Obj)s;
  }
}

static void
Pdistance(arg,rp)
NODE arg;
Obj *rp;
{
  Num  s;
  Num  a, b;

  asir_assert(ARG0(arg),O_N,"distance");
  asir_assert(ARG1(arg),O_N,"distance");
  a = (Num)ARG0(arg);
  b = (Num)ARG1(arg);
  if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {
    distanceitvd((IntervalDouble)a, (IntervalDouble)b, (Num *)rp);
  } else {
    distanceitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);
    *rp = (Obj)s;
  }
}

static void
Pinitv(NODE arg, Obj *rp)
{
  int  s;
  Z  q;

  asir_assert(ARG0(arg),O_N,"intval");
  asir_assert(ARG1(arg),O_N,"intval");
  if ( ! ARG1(arg) ) {
    if ( ! ARG0(arg) ) s = 1;
    else s = 0;
  }
  else if ( NID(ARG1(arg)) == N_IntervalDouble ) {
    s = initvd((Num)ARG0(arg),(IntervalDouble)ARG1(arg));

  } else if ( NID(ARG1(arg)) == N_IP || NID(ARG1(arg)) == N_IntervalBigFloat ) {
    if ( ! ARG0(arg) ) s = initvp((Num)ARG0(arg),(Itv)ARG1(arg));
    else if ( NID(ARG0(arg)) == N_IP ) {
      s = itvinitvp((Itv)ARG0(arg),(Itv)ARG1(arg));
    } else {
      s = initvp((Num)ARG0(arg),(Itv)ARG1(arg));
    }
  } else {
    s = ! compnum(0,(Num)ARG0(arg),(Num)ARG1(arg));
  }
  STOZ(s,q);
  *rp = (Obj)q;
}

static void
Pdisjitv(arg,rp)
NODE arg;
Obj *rp;
{
  Itv  s;

  asir_assert(ARG0(arg),O_N,"disjitv");
  asir_assert(ARG1(arg),O_N,"disjitv");
  error("disjitv: not implemented yet");
  if ( ! s ) *rp = 0;
  else *rp = (Obj)ONE;
}

static void
PzeroRewriteMode(NODE arg, Obj *rp)
{
  Q  a;
  Z r;

  STOZ(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;
  Z  r;

  STOZ(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)
{
  Z  r;

  STOZ(zerorewriteCount,r);
  *rp = (Obj)r;
}
  

#endif
extern int  printmode;

static void  pprintmode( void )
{
  switch (printmode) {
#if defined(INTERVAL)
    case MID_PRINTF_E:
      fprintf(stderr,"Interval printing mode is a mitpoint type.\n");
#endif
    case PRINTF_E:
      fprintf(stderr,"Printf's double printing mode is \"%%.16e\".\n");
      break;
#if defined(INTERVAL)
    case MID_PRINTF_G:
      fprintf(stderr,"Interval printing mode is a mitpoint type.\n");
#endif
    default:
    case PRINTF_G:
      fprintf(stderr,"Printf's double printing mode is \"%%g\".\n");
      break;
  }
}

static void
Pprintmode(NODE arg, Obj *rp)
{
  int  l;
  Z  a, r;

  a = (Z)ARG0(arg);
  if(!a||(NUM(a)&&INT(a))){
    l=ZTOS(a);
    if ( l < 0 ) l = 0;
#if defined(INTERVAL)
    else if ( l > MID_PRINTF_E ) l = 0;
#else
    else if ( l > PRINTF_E ) l = 0;
#endif
    STOZ(printmode,r);
    *rp = (Obj)r;
    printmode = l;
    pprintmode();
  } else {
    *rp = 0;
  }
}