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

Diff for /OpenXM_contrib2/asir2000/plot/calc.c between version 1.8 and 1.9

version 1.8, 2011/08/10 04:51:58 version 1.9, 2013/12/19 06:04:09
Line 45 
Line 45 
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *   *
  * $OpenXM: OpenXM_contrib2/asir2000/plot/calc.c,v 1.7 2003/02/14 22:29:19 ohara Exp $   * $OpenXM: OpenXM_contrib2/asir2000/plot/calc.c,v 1.8 2011/08/10 04:51:58 saito Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "parse.h"  #include "parse.h"
Line 59 
Line 59 
 #define MAXSHORT ((short)0x7fff)  #define MAXSHORT ((short)0x7fff)
 #endif  #endif
   
 void calc(double **tab,struct canvas *can,int nox)  void calc(double **tab,struct canvas *can,int nox){
 {          //memory_plot,MODE_IFPLOTD,MODE_INEQND,MODE_INEQNANDD,MODE_INEQNORD
         double x,y,xmin,ymin,xstep,ystep;          //MODE_INEQNXORD,conplotmainD
           double x,y,xstep,ystep;
         int ix,iy;          int ix,iy;
         Real r,rx,ry;          Real r,rx,ry;
         Obj fr,g;          Obj fr,g,t,s;
         int w,h;  
         V vx,vy;  
         Obj t,s;  
   
         if ( !nox ) initmarker(can,"Evaluating...");          if(!nox)initmarker(can,"Evaluating...");
         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;          xstep=(can->xmax-can->xmin)/can->width;
         vy = can->vy;          ystep=(can->ymax-can->ymin)/can->height;
         w = can->width; h = can->height;          MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
         xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;          BDY(rx)=can->xmin;
         ymin = can->ymin; ystep = (can->ymax-can->ymin)/h;          substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&t); devalr(CO,t,&g);
         MKReal(1.0,rx); MKReal(1.0,ry); /* dummy real */          BDY(ry)=can->ymin;
         for( ix = 0, x = xmin; ix < w ; ix++, x += xstep ) {          substr(CO,0,g,can->vy,can->ymin?(Obj)ry:0,&t); devalr(CO,t,&s);
                 BDY(rx) = x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);          can->vmax=can->vmin=ToReal(s);
           for(ix=0,x=can->xmin; ix<can->width; ix++,x+=xstep){
                   BDY(rx)=x; substr(CO,0,fr,can->vx,x?(Obj)rx:0,&t);
                 devalr(CO,t,&g);                  devalr(CO,t,&g);
                 if ( !nox ) marker(can,DIR_X,ix);                  if(!nox)marker(can,DIR_X,ix);
                 for( iy = 0, y = ymin; iy < h ; iy++, y += ystep ) {                  for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
                         BDY(ry) = y;                          BDY(ry)=y;
                         substr(CO,0,g,vy,y?(Obj)ry:0,&t);                          substr(CO,0,g,can->vy,y?(Obj)ry:0,&t);
                         devalr(CO,t,&s);                          devalr(CO,t,&s);
                         tab[ix][iy] = ToReal(s);                          tab[ix][iy]=ToReal(s);
                           if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
                           if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
                 }                  }
         }          }
 }  }
   
 double usubstrp(P p,double r)  void calcq(double **tab,struct canvas *can,int nox){
 {          //MODE_IFPLOTQ,MODE_INEQNQ,MODE_INEQNANDQ,MODE_INEQNORQ,MODE_INEQNXORQ
         double t;          //plotoverD
           Q dx,dy,xstep,ystep,q1,w,h,c;
           P g,g1,f1,f2,x,y;
           int ix,iy;
           Obj fr,gm,t,s;
           Real r,rx,ry;
   
           MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
           MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
           BDY(rx)=can->xmin;
           substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&t); devalr(CO,t,&gm);
           BDY(ry)=can->ymin;
           substr(CO,0,gm,can->vy,can->ymin?(Obj)ry:0,&t); devalr(CO,t,&s);
           can->vmax=can->vmin=ToReal(s);
   
           subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
           subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
           MKV(can->vx,x); mulp(CO,(P)xstep,x,(P *)&t);
           addp(CO,(P)can->qxmin,(P)t,(P *)&s); substp(CO,can->formula,can->vx,(P)s,&f1);
           MKV(can->vy,y); mulp(CO,(P)ystep,y,(P *)&t);
           addp(CO,(P)can->qymin,(P)t,(P *)&s); substp(CO,f1,can->vy,(P)s,&f2);
           ptozp(f2,1,&c,&g);
           if(!nox) initmarker(can,"Evaluating...");
           for(iy=0;iy<can->height;iy++){
                   marker(can,DIR_Y,iy);
                   STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,(P *)&t); ptozp((P)t,1,&c,&g1);
                   for(ix=0;ix<can->width;ix++){
                           STOQ(ix,q1);substp(CO,g1,can->vx,(P)q1,(P *)&t);
                           devalr(CO,t,&s);
                           tab[ix][iy]=ToReal(s);
                           if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
                           if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
                   }
           }
   }
   
   void calcb(double **tab,struct canvas *can,int nox){
           //MODE_IFPLOTB,MODE_INEQNB,MODE_INEQNANDB,MODE_INEQNORB,MODE_INEQNXORB
           Q dx,dy,xstep,ystep,q1,w,h,c;
           P g,g1,f1,f2,x,y,t,s;
           int ix,iy,*a,*pa;
           VECT ss;
           Obj fr,gm,tm,sm;
           Real r,rx,ry;
   
           MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
           MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
           BDY(rx)=can->xmin;
           substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&tm); devalr(CO,tm,&gm);
           BDY(ry)=can->ymin;
           substr(CO,0,gm,can->vy,can->ymin?(Obj)ry:0,&tm); devalr(CO,tm,&sm);
           can->vmax=can->vmin=ToReal(sm);
   
           for(iy=0;iy<can->height;iy++)for(ix=0;ix<can->width;ix++)tab[ix][iy]=1.0;
           subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
           subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
           MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
           addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
           MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
           addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
           ptozp(f2,1,&c,&g);
           a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
           for(iy=0;iy<can->height;iy++)for(ix=0;ix<can->width;ix++)tab[ix][iy]=1.0;
           subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
           subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
           MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
           addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
           MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
           addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
           ptozp(f2,1,&c,&g);
           a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
           for(ix=0;ix<can->width;ix++){
                   STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
                   if(!g1)for(iy=0;iy<can->height;iy++)tab[ix][iy]=0.0;
                   else if(!NUM(g1)){
                           sturmseq(CO,g1,&ss);
                           seproot(ss,0,can->width,a);
                           for(iy=0,pa=a;iy<can->height;iy++,pa++){
                                   if(*pa<0||(*(pa+1)>=0&&(*pa>*(pa+1))))tab[ix][iy]=0.0;
                                   else {
                                           STOQ(iy,q1);substp(CO,g1,can->vy,(P)q1,&t);
                                           devalr(CO,(Obj)t,(Obj *)&s);
                                           tab[ix][iy]=ToReal(s);
                                           if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
                                           if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
                                   }
                           }
                   }
           }
           for(iy=0;iy<can->height;iy++){
                   STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
                   if(!g1) for(ix=0;ix<can->width;ix++)tab[ix][iy]=0.0;
                   else if(!NUM(g1)){
                           sturmseq(CO,g1,&ss);
                           seproot(ss,0,can->height,a);
                           for(ix=0,pa=a;ix<can->width;ix++,pa++){
                                   if(tab[ix][iy]!=0.0){
                                           if(*pa<0||(*(pa+1)>=0&&(*pa>*(pa+1))))tab[ix][iy]=0.0;
                                           else {
                                                   STOQ(ix,q1);substp(CO,g1,can->vx,(P)q1,&t);
                                                   devalr(CO,(Obj)t,(Obj *)&s);
                                                   tab[ix][iy]=ToReal(s);
                                                   if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
                                                   if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
                                           }
                                   }
                           }
                   }
           }
   }
   
   double usubstrp(P p,double r){
         DCP dc;          DCP dc;
         int d;          int d;
         double pwrreal0();          double t,pwrreal0();
   
         if ( !p )          if(!p) t=0.0;
                 t = 0.0;          else if(NUM(p))t=BDY((Real)p);
         else if ( NUM(p) )  
                 t = BDY((Real)p);  
         else {          else {
                 dc = DC(p); t = BDY((Real)COEF(dc));                  dc=DC(p); t=BDY((Real)COEF(dc));
                 for ( d = QTOS(DEG(dc)), dc = NEXT(dc); dc;                  for(d=QTOS(DEG(dc)),dc=NEXT(dc);dc;d=QTOS(DEG(dc)),dc=NEXT(dc)){
                         d = QTOS(DEG(dc)), dc = NEXT(dc) ) {                          t=t*pwrreal0(r,(d-QTOS(DEG(dc))))+BDY((Real)COEF(dc));
                         t = t*pwrreal0(r,(d-QTOS(DEG(dc))))+BDY((Real)COEF(dc));  
                 }                  }
                 if ( d )                  if(d)t*=pwrreal0(r,d);
                         t *= pwrreal0(r,d);  
         }          }
         return t;          return t;
 }  }
   
 void qcalc(char **tab,struct canvas *can)  void qcalc(char **tab,struct canvas *can){
 {          //qifplotmain(Old type)
         Q dx,dy,w,h,xstep,ystep,c,q1;          Q dx,dy,w,h,xstep,ystep,c,q1;
         P g,g1,f1,f2,x,y,t,s;          P g,g1,f1,f2,x,y,t,s;
         int ix,iy;          int ix,iy;
Line 128  void qcalc(char **tab,struct canvas *can)
Line 237  void qcalc(char **tab,struct canvas *can)
         MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);          MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
         addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);          addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
         ptozp(f2,1,&c,&g);          ptozp(f2,1,&c,&g);
         a = (int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));          a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
         initmarker(can,"Horizontal scan...");          initmarker(can,"Horizontal scan...");
         for( ix = 0; ix < can->width; ix++ ) {          for( ix=0; ix < can->width; ix++ ){
                 marker(can,DIR_X,ix);                  marker(can,DIR_X,ix);
                 STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);                  STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
                 if ( !g1 )                  if( !g1 )
                         for ( iy = 0; iy < can->height; iy++ )                          for(iy=0; iy < can->height; iy++ )
                                 tab[ix][iy] = 1;                                  tab[ix][iy]=1;
                 else if ( !NUM(g1) ) {                  else if( !NUM(g1) ){
                         sturmseq(CO,g1,&ss); seproot(ss,0,can->height,a);                          sturmseq(CO,g1,&ss); seproot(ss,0,can->height,a);
                         for ( iy = 0, pa = a; iy < can->height; iy++, pa++ )                          for(iy=0, pa=a; iy < can->height; iy++, pa++ )
                                 if ( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )                                  if( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
                                         tab[ix][iy] = 1;                                          tab[ix][iy]=1;
                 }                  }
         }          }
         initmarker(can,"Vertical scan...");          initmarker(can,"Vertical scan...");
         for( iy = 0; iy < can->height; iy++ ) {          for( iy=0; iy < can->height; iy++ ){
                 marker(can,DIR_Y,iy);                  marker(can,DIR_Y,iy);
                 STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);                  STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
                 if ( !g1 )                  if( !g1 )
                         for ( ix = 0; ix < can->width; ix++ )                          for(ix=0; ix < can->width; ix++ )
                                 tab[ix][iy] = 1;                                  tab[ix][iy]=1;
                 else if ( !NUM(g1) ) {                  else if( !NUM(g1) ){
                         sturmseq(CO,g1,&ss); seproot(ss,0,can->width,a);                          sturmseq(CO,g1,&ss); seproot(ss,0,can->width,a);
                         for ( ix = 0, pa = a; ix < can->width; ix++, pa++ )                          for(ix=0, pa=a; ix < can->width; ix++, pa++ )
                                 if ( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )                                  if( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
                                         tab[ix][iy] = 1;                                          tab[ix][iy]=1;
                 }                  }
         }          }
 }  }
   
 void sturmseq(VL vl,P p,VECT *rp)  void sturmseq(VL vl,P p,VECT *rp){
 {          P g1,g2,q,r,s,*t;
         P g1,g2,q,r,s;  
         P *t;  
         V v;          V v;
         VECT ret;          VECT ret;
         int i,j;          int i,j;
         Q a,b,c,d,h,l,m,x;          Q a,b,c,d,h,l,m,x;
   
         v = VR(p); t = (P *)ALLOCA((deg(v,p)+1)*sizeof(P));          v=VR(p);t=(P *)ALLOCA((deg(v,p)+1)*sizeof(P));
         g1 = t[0] = p; diffp(vl,p,v,(P *)&a); ptozp((P)a,1,&c,&g2); t[1] = g2;          g1=t[0]=p;diffp(vl,p,v,(P *)&a);ptozp((P)a,1,&c,&g2);t[1]=g2;
         for ( i = 1, h = ONE, x = ONE; ; ) {          for(i=1,h=ONE,x=ONE;;){
                 if ( NUM(g2) )                  if(NUM(g2)) break;
                         break;  
                 subq(DEG(DC(g1)),DEG(DC(g2)),&d);                  subq(DEG(DC(g1)),DEG(DC(g2)),&d);
                 l = (Q)LC(g2);                  l=(Q)LC(g2);
                 if ( SGN(l) < 0 ) {                  if(SGN(l)<0){
                         chsgnq(l,&a); l = a;                          chsgnq(l,&a);l=a;
                 }                  }
                 addq(d,ONE,&a); pwrq(l,a,&b); mulp(vl,(P)b,g1,(P *)&a);                  addq(d,ONE,&a);pwrq(l,a,&b);mulp(vl,(P)b,g1,(P *)&a);
                 divsrp(vl,(P)a,g2,&q,&r);                  divsrp(vl,(P)a,g2,&q,&r);
                 if ( !r )                  if(!r) break;
                         break;                  chsgnp(r,&s);r=s;i++;
                 chsgnp(r,&s); r = s; i++;                  if(NUM(r)){
                 if ( NUM(r) ) {                          t[i]=r;break;
                         t[i] = r; break;  
                 }                  }
                 pwrq(h,d,&m); g1 = g2;                  pwrq(h,d,&m);g1=g2;
                 mulq(m,x,&a); divsp(vl,r,(P)a,&g2); t[i] = g2;                  mulq(m,x,&a);divsp(vl,r,(P)a,&g2);t[i]=g2;
                 x = (Q)LC(g1);                  x=(Q)LC(g1);
                 if ( SGN(x) < 0 ) {                  if(SGN(x)<0){
                         chsgnq(x,&a); x = a;                          chsgnq(x,&a);x=a;
                 }                  }
                 pwrq(x,d,&a); mulq(a,h,&b); divq(b,m,&h);                  pwrq(x,d,&a);mulq(a,h,&b);divq(b,m,&h);
         }          }
         MKVECT(ret,i+1);          MKVECT(ret,i+1);
         for ( j = 0; j <= i; j++ )          for(j=0;j<=i;j++)
                 ret->body[j] = (pointer)t[j];                  ret->body[j]=(pointer)t[j];
         *rp = ret;          *rp=ret;
 }  }
   
 void seproot(VECT s,int min,int max,int *ar)  void seproot(VECT s,int min,int max,int *ar){
 {          P f,*ss;
         P f;  
         P *ss;  
         Q q,t;          Q q,t;
         int i,j,k;          int i,j,k;
   
         ss = (P *)s->body; f = ss[0];          ss=(P *)s->body;f=ss[0];
         for ( i = min; i <= max; i++ ) {          for(i=min;i<=max;i++){
                 STOQ(i,q); usubstqp(f,q,&t);                  STOQ(i,q);usubstqp(f,q,&t);
                 if ( !t )                  if(!t)ar[i]=-1;
                         ar[i] = -1;  
                 else {                  else {
                         ar[i] = numch(s,q,t); break;                          ar[i]=numch(s,q,t);break;
                 }                  }
         }          }
         if ( i > max )          if(i>max) return;
                 return;          for(j=max;j>= min;j--){
         for ( j = max; j >= min; j-- ) {  
                 STOQ(j,q); usubstqp(f,q,&t);                  STOQ(j,q); usubstqp(f,q,&t);
                 if ( !t )                  if(!t)ar[j]=-1;
                         ar[j] = -1;  
                 else {                  else {
                         if ( i != j )                          if(i!=j)ar[j]=numch(s,q,t);
                                 ar[j] = numch(s,q,t);  
                         break;                          break;
                 }                  }
         }          }
         if ( j <= i+1 )          if(j<=i+1) return;
           if(ar[i]==ar[j]){
                   for(k=i+1;k<j;k++)ar[k]=ar[i];
                 return;                  return;
         if ( ar[i] == ar[j] ) {  
                 for ( k = i+1; k < j; k++ )  
                         ar[k] = ar[i];  
                 return;  
         }          }
         k = (i+j)/2;          k=(i+j)/2;
         seproot(s,i,k,ar);          seproot(s,i,k,ar);
         seproot(s,k,j,ar);          seproot(s,k,j,ar);
 }  }
   
 int numch(VECT s,Q n,Q a0)  int numch(VECT s,Q n,Q a0){
 {  
         int len,i,c;          int len,i,c;
         Q a;          Q a;
         P *ss;          P *ss;
   
         len = s->len; ss = (P *)s->body;          len=s->len;ss=(P *)s->body;
         for ( i = 1, c = 0; i < len; i++ ) {          for(i=1,c=0;i<len;i++){
                 usubstqp(ss[i],n,&a);                  usubstqp(ss[i],n,&a);
                 if ( a ) {                  if(a){
                         if ( (SGN(a)>0 && SGN(a0)<0) || (SGN(a)<0 && SGN(a0)>0) )                          if((SGN(a)>0 && SGN(a0)<0)||(SGN(a)<0&&SGN(a0)>0))c++;
                                 c++;                          a0=a;
                         a0 = a;  
                 }                  }
         }          }
         return c;          return c;
 }  }
   
 void usubstqp(P p,Q r,Q *v)  void usubstqp(P p,Q r,Q *v){
 {  
         Q d,d1,a,b,t;          Q d,d1,a,b,t;
         DCP dc;          DCP dc;
   
         if ( !p )          if(!p)
                 *v = 0;                  *v=0;
         else if ( NUM(p) )          else if(NUM(p))*v=(Q)p;
                 *v = (Q)p;  
         else {          else {
                 dc = DC(p); t = (Q)COEF(dc);                  dc=DC(p);t=(Q)COEF(dc);
                 for ( d = DEG(dc), dc = NEXT(dc); dc;                  for(d=DEG(dc),dc=NEXT(dc);dc;d=DEG(dc),dc=NEXT(dc)){
                         d = DEG(dc), dc = NEXT(dc) ) {                          subq(d,DEG(dc),&d1);pwrq(r,d1,&a);
                         subq(d,DEG(dc),&d1); pwrq(r,d1,&a);                          mulq(t,a,&b);addq(b,(Q)COEF(dc),&t);
                         mulq(t,a,&b); addq(b,(Q)COEF(dc),&t);  
                 }                  }
                 if ( d ) {                  if(d){
                         pwrq(r,d,&a); mulq(t,a,&b); t = b;                          pwrq(r,d,&a);mulq(t,a,&b);t=b;
                 }                  }
                 *v = t;                  *v=t;
         }          }
 }  }
   
 void plotcalc(struct canvas *can)  void plotcalc(struct canvas *can){
 {          //plot,memory_plot,plotover,plot_resize
         double x,xmin,xstep,ymax,ymin,dy;          double x,xmin,xstep,ymax,ymin,dy,*tab,usubstrp();
         int ix;          int ix,w,h;
         Real r;          Real r,rx;
         Obj fr;          Obj fr,t,s;
         double usubstrp();  
         int w,h;  
         double *tab;  
         POINT *pa;          POINT *pa;
         Real rx;  
         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);
         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;
         tab = (double *)ALLOCA(w*sizeof(double));          tab=(double *)ALLOCA(w*sizeof(double));
         MKReal(1,rx); /* dummy real number */          MKReal(1,rx); // dummy real number
         for( ix = 0, x = xmin; ix < w ; ix++, x += xstep ) {          for(ix=0,x=xmin;ix<w;ix++,x+=xstep){
                 /* full substitution */                  // full substitution
                 BDY(rx) = x;                  BDY(rx)=x;
                 substr(CO,0,fr,can->vx,x?(Obj)rx:0,&s);                  substr(CO,0,fr,can->vx,x?(Obj)rx:0,&s);
                 devalr(CO,(Obj)s,&t);                  devalr(CO,(Obj)s,&t);
                 if ( t && (OID(t)!=O_N || NID((Num)t)!=N_R) )                  if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
                         error("plotcalc : invalid evaluation");                          error("plotcalc : invalid evaluation");
                 tab[ix] = ToReal((Num)t);                  tab[ix]=ToReal((Num)t);
         }          }
         if ( can->ymax == can->ymin ) {          if(can->ymax==can->ymin){
                 for ( ymax = ymin = tab[0], ix = 1; ix < w; ix++ ) {                  for(ymax=ymin=tab[0],ix=1;ix<w;ix++){
                         if ( tab[ix] > ymax )                          if(tab[ix]>ymax)ymax=tab[ix];
                                 ymax = tab[ix];                          if(tab[ix]<ymin)ymin=tab[ix];
                         if ( tab[ix] < ymin )  
                                 ymin = tab[ix];  
                 }                  }
                 can->ymax = ymax; can->ymin = ymin;                  can->ymax=ymax;can->ymin=ymin;
         } else {          } else {
                 ymax = can->ymax; ymin = can->ymin;                  ymax=can->ymax;ymin=can->ymin;
         }          }
         dy = ymax-ymin;          dy=ymax-ymin;
         can->pa = (struct pa *)MALLOC(sizeof(struct pa));          can->pa=(struct pa *)MALLOC(sizeof(struct pa));
         can->pa[0].length = w;          can->pa[0].length=w;
         can->pa[0].pos = pa = (POINT *)MALLOC(w*sizeof(POINT));          can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
         for ( ix = 0; ix < w; ix++ ) {          for(ix=0;ix<w;ix++){
                 double t;                  double t;
                   XC(pa[ix])=ix;
                 XC(pa[ix]) = ix;                  t=(h-1)*(ymax-tab[ix])/dy;
                 t = (h - 1)*(ymax - tab[ix])/dy;                  if(t>MAXSHORT)YC(pa[ix])=MAXSHORT;
                 if ( t > MAXSHORT )                  else if(t<-MAXSHORT)YC(pa[ix])=-MAXSHORT;
                         YC(pa[ix]) = MAXSHORT;                  else YC(pa[ix])=(long)t;
                 else if ( t < -MAXSHORT )  
                         YC(pa[ix]) = -MAXSHORT;  
                 else  
                         YC(pa[ix]) = (long)t;  
         }          }
 }  }
   
 void polarplotcalc(struct canvas *can)  void polarplotcalc(struct canvas *can){
 {          //polarplotNG
         double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep;          double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep,tr,p, *tabx,*taby;
         int i,nstep;  
         double usubstrp();          double usubstrp();
         int w,h;          int i,nstep,w,h;
         double tr,p;  
         double *tabx,*taby;  
         POINT *pa;          POINT *pa;
         Real r;          Real r;
         Obj fr,t,s;          Obj fr,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);
         w = can->width; h = can->height; nstep = can->nzstep;          w=can->width; h=can->height; nstep=can->nzstep;
         pmin = can->zmin; pstep = (can->zmax-can->zmin)/nstep;          pmin=can->zmin; pstep=(can->zmax-can->zmin)/nstep;
         tabx = (double *)ALLOCA(nstep*sizeof(double));          tabx=(double *)ALLOCA(nstep*sizeof(double));
         taby = (double *)ALLOCA(nstep*sizeof(double));          taby=(double *)ALLOCA(nstep*sizeof(double));
         MKReal(1,r); /* dummy real number */          MKReal(1,r); // dummy real number
         for( i = 0, p = pmin; i < nstep ; i++, p += pstep ) {          for(i=0,p=pmin;i<nstep;i++,p+= pstep){
                 /* full substitution */                  // full substitution
                 BDY(r) = p;                  BDY(r)=p;
                 substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);                  substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);
                 devalr(CO,(Obj)s,&t);                  devalr(CO,(Obj)s,&t);
                 if ( t && (OID(t)!=O_N || NID((Num)t)!=N_R) )                  if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
                         error("polarplotcalc : invalid evaluation");                          error("polarplotcalc : invalid evaluation");
                 tr = ToReal((Num)t);                  tr=ToReal((Num)t);
                 tabx[i] = tr*cos(p);                  tabx[i]=tr*cos(p);
                 taby[i] = tr*sin(p);                  taby[i]=tr*sin(p);
         }          }
         xmax = xmin = tabx[0];          xmax=xmin=tabx[0];
         ymax = ymin = taby[0];          ymax=ymin=taby[0];
         for ( i = 1; i < nstep; i++ ) {          for(i=1;i<nstep;i++){
                 if ( tabx[i] > xmax ) xmax = tabx[i];                  if(tabx[i]>xmax)xmax=tabx[i];
                 if ( tabx[i] < xmin ) xmin = tabx[i];                  if(tabx[i]<xmin)xmin=tabx[i];
                 if ( taby[i] > ymax ) ymax = taby[i];                  if(taby[i]>ymax)ymax=taby[i];
                 if ( taby[i] < ymin ) ymin = taby[i];                  if(taby[i]<ymin)ymin=taby[i];
         }          }
         can->xmax = xmax; can->xmin = xmin;          can->xmax=xmax;can->xmin=xmin;
         can->ymax = ymax; can->ymin = ymin;          can->ymax=ymax;can->ymin=ymin;
         dx = xmax-xmin;          dx=xmax-xmin;
         dy = ymax-ymin;          dy=ymax-ymin;
         can->pa = (struct pa *)MALLOC(sizeof(struct pa));          can->pa=(struct pa *)MALLOC(sizeof(struct pa));
         can->pa[0].length = nstep;          can->pa[0].length=nstep;
         can->pa[0].pos = pa = (POINT *)MALLOC(w*sizeof(POINT));          can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
         for ( i = 0; i < nstep; i++ ) {          for(i=0;i<nstep;i++){
                 XC(pa[i]) = (w-1)*(tabx[i]-xmin)/dx;                  XC(pa[i])=(w-1)*(tabx[i]-xmin)/dx;
                 YC(pa[i]) = (h-1)*(ymax-taby[i])/dy;                  YC(pa[i])=(h-1)*(ymax-taby[i])/dy;
         }          }
 }  }
   
 #if defined(INTERVAL)  /*
 #define NORMAL 1  void ineqncalc(double **tab,struct canvas *can,int nox){
 #define TRANSFER 1          double x,y,xmin,ymin,xstep,ystep;
 #define RECURSION 1          int ix,iy,w,h;
 #define RECTRANS 1          Real r,rx,ry;
           Obj fr,g,t,s;
           V vx,vy;
   
 void ineqncalc(double **tab,struct canvas *can,int nox)          if(!nox) initmarker(can,"Evaluating...");
 {  
         double x, y, xmin, ymax, xstep, ystep;  
         int ix, iy;  
         Real r, rx, ry;  
         Obj fr, g, t, s;  
         int w, h;  
         V vx, vy;  
   
         if ( !nox ) initmarker(can,"Evaluating...");  
         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; vy = can->vy;          vx=can->vx;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;
         ymax = can->ymax; ystep = (can->ymax-can->ymin)/h;          ymin=can->ymin;ystep=(can->ymin-can->ymin)/h;
         MKReal(1.0,rx); MKReal(1.0,ry); /* dummy real */          MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
   
         for( iy = 0, y = ymax; iy <= h ; iy++, y -= ystep ) {          for(ix=0,x=xmin;ix<=w;ix++,x+=xstep){
                 BDY(ry) = y; substr(CO,0,fr,vy,y?(Obj)ry:0,&t);                  BDY(rx)=x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);
                 devalr(CO,t,&g);                  devalr(CO,t,&g);
                 if ( !nox ) marker(can,DIR_Y,iy);                  if(!nox) marker(can,DIR_X,ix);
                 for( ix = 0, x = xmin; ix <= w ; ix++, x += xstep ) {                  for(iy=0,y=ymin;iy<=h;iy++,y+=ystep){
                         BDY(rx) = x;                          BDY(ry)=y;
                         substr(CO,0,g,vx,x?(Obj)rx:0,&t);                          substr(CO,0,g,vy,y?(Obj)ry:0,&t);
                         devalr(CO,t,&s);                          devalr(CO,t,&s);
                         tab[iy][ix] = ToReal(s);                          tab[ix][iy]=ToReal(s);
                 }                  }
         }          }
 }  }
   */
   
 void reccalc(P fr, V vx, V vy, int ixlw, int ixhg, int iylw, int iyhg,  #if defined(INTERVAL)
         double *xval, double *yval, int w, int **mask, int itvsize)  void itvcalc(double **mask, struct canvas *can, int nox){
 {          //MODE_ITVIFPLOT
         int i, j;          double x,y,xstep,ystep,dx,dy,wx,wy;
         int ixmd, iymd;          int idv,ix,iy,idx,idy;
         double inf, sup, dt;          Itv ity,itx,ddx,ddy;
         Real rxwid, rywid, rxlw, rylw, rzero;          Real r,rx,ry,rx1,ry1,rdx,rdy,rdx1,rdy1;
         Itv itx, ity;          V vx,vy;
         Obj fr1;          Obj fr,g,t,s;
         P py, px, tmp;  
   
         double **tbl , xs, ys;          idv=can->division;
         int ywidth, xwidth;  
         Real rx, ry;  
         Obj g, t, s;  
   
         if ( (ixlw > ixhg) || (iylw > iyhg) ) return;  
         NEWItvP(itx); NEWItvP(ity);  
         MKReal(xval[ixhg + 1] - xval[ixlw], rxwid);  
         MKReal(yval[iyhg + 1] - yval[iylw], rywid);  
         MKReal(0.0,rzero);  
         istoitv((Num)rzero, (Num)rxwid, &itx);  
         istoitv((Num)rzero, (Num)rywid, &ity);  
   
         MKReal(xval[ixlw],rxlw);MKReal(yval[iylw],rylw);  
         MKV(vx,tmp); addp(CO,(P)tmp,(P)rxlw,&px);  
         MKV(vy,tmp); addp(CO,(P)tmp,(P)rylw,&py);  
         substp(CO,(P)fr,(V)vx,(P)px,(P*)&fr1);  
         substp(CO,(P)fr1,(V)vy,(P)py,(P*)&fr1);  
         substr(CO,0,(Obj)fr1,vx,(Obj)itx,&fr1);  
         substr(CO,0,(Obj)fr1,vy,(Obj)ity,&fr1);  
         Num2double((Num)fr1,&inf,&sup);  
   
         if(inf > sup){dt=inf;inf=sup;sup=dt;}  
         if ( inf <= 0.0 && sup >= 0.0 ) {  
                 if ( (ixhg - ixlw <= itvsize) && (iyhg - iylw <= itvsize) ){  
                         ywidth = iyhg - iylw + 2; xwidth = ixhg - ixlw + 2;  
                         MKReal(1.0, rx); MKReal(1.0, ry);  
                         tbl = (double **)ALLOCA((ywidth)*sizeof(double *));  
                         for (i=0;i<ywidth;i++)  
                                 tbl[i] = (double *)ALLOCA((xwidth)*sizeof(double));  
                         for (i = 0; i < ywidth; i++){  
                                 BDY(ry) = yval[iylw+i];  
                                 substr(CO,0,(Obj)fr,vy,(Obj)ry,&t);  
                                 devalr(CO, t, &g);  
                                 for (j = 0; j < xwidth; j++){  
                                         BDY(rx) = xval[ixlw+j];  
                                         substr(CO,0,g,vx,(Obj)rx,&t);  
                                         devalr(CO, t, &s);  
                                         tbl[i][j] = ToReal(s);  
                                 }  
                         }  
                         for ( i = 0; i < ywidth - 1; i++){  
                                 for ( j = 0; j < xwidth - 1; j++){  
                                         if ( tbl[i][j] >= 0.0 ){  
                                                 if ( (tbl[i+1][j] <= 0.0) ||  
                                                         (tbl[i+1][j+1] <= 0.0) ||  
                                                         (tbl[i][j+1]   <= 0.0) ) mask[w-(i+iylw)][j+ixlw] = 0;  
                                                 else mask[w-(i+iylw)][j+ixlw] = 1;  
                                         } else {  
                                                 if ( (tbl[i+1][j] >= 0.0) ||  
                                                         (tbl[i+1][j+1] >= 0.0) ||  
                                                         (tbl[i][j+1]   >= 0.0) ) mask[w-(i+iylw)][j+ixlw] = 0;  
                                                 else mask[w-(i+iylw)][j+ixlw] = 1;  
                                         }  
                                 }  
                         }  
                 } else {  
                         ixmd = (ixhg + ixlw)/2; iymd = (iyhg + iylw)/2;  
                         reccalc((P)fr,vx,vy,ixlw,ixmd,iylw,iymd,xval,yval,w,mask,itvsize);  
                         reccalc((P)fr,vx,vy,ixmd+1,ixhg,iylw,iymd,xval,yval,w,mask,itvsize);  
                         reccalc((P)fr,vx,vy,ixlw,ixmd,iymd+1,iyhg,xval,yval,w,mask,itvsize);  
                         reccalc((P)fr,vx,vy,ixmd+1,ixhg,iymd+1,iyhg,xval,yval,w,mask,itvsize);  
                 }  
         } else {  
                 for (i=w-iyhg; i> w-iylw; i--)  
                         for (j=ixlw; j <= ixhg; j++) mask[i][j] = -1;  
         }  
 }  
   
 void itvcalc(int **mask, struct canvas *can, int nox, int itvsize)  
 {  
         double xstep, ystep, xv, yv, *xval, *yval;  
         int imx, imy, i, w, h;  
         Real r;  
         Obj fr;  
         V vx, vy;  
   
         if ( !nox ) initmarker(can,"Evaluating...");  
         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; vy = can->vy;          vx=can->vx; vy=can->vy;
         w = can->width; h = can->height;          xstep=(can->xmax-can->xmin)/can->width;
         xstep = (can->xmax - can->xmin)/w;          ystep=(can->ymax-can->ymin)/can->height;
         ystep = (can->ymax - can->ymin)/h;          if(idv!=0){
                   wx=xstep/can->division;
         xval = (double *)ALLOCA((w+1)*sizeof(double));                  wy=ystep/can->division;
         yval = (double *)ALLOCA((h+1)*sizeof(double));          }
         for (i=0, xv=can->xmin; i<= w; i++, xv += xstep) xval[i] = xv;          MKReal(can->ymin,ry1);
         for (i=0, yv=can->ymin; i<= h; i++, yv += ystep) yval[i] = yv;          for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
         imx = w/2; imy = h/2;                  ry=ry1;
                   MKReal(y+ystep,ry1);
         reccalc((P)fr,vx,vy,0,imx,0,imy,xval,yval,w-1,mask,itvsize);  
         reccalc((P)fr,vx,vy,imx+1,w-1,0,imy,xval,yval,w-1,mask,itvsize);  
         reccalc((P)fr,vx,vy,0,imx,imy+1,h-1,xval,yval,w-1,mask,itvsize);  
         reccalc((P)fr,vx,vy,imx+1,w-1,imy+1,h-1,xval,yval,w-1,mask,itvsize);  
 }  
   
 #if NORMAL  
 void itvcalc1(int **mask, struct canvas *can, int nox)  
 {  
         double x, y, xstep, ystep;  
         int ix, iy, w, h;  
         Itv ity, itx;  
         Real r, rx, ry, rx1, ry1;  
         Obj fr, g, t;  
         V vx, vy;  
   
         if ( !nox ) initmarker(can,"Evaluating...");  
         MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);  
         vx = can->vx; vy = can->vy;  
         w = can->width; h = can->height;  
         xstep = (can->xmax - can->xmin)/w;  
         ystep = (can->ymax - can->ymin)/h;  
   
         for( iy = 0, y = can->ymax; iy < h ; iy++, y -= ystep ) {  
                 MKReal(y, ry);  
                 MKReal(y-ystep,ry1);  
                 istoitv((Num)(ry1),(Num)ry,&ity);                  istoitv((Num)(ry1),(Num)ry,&ity);
                 substr(CO,0,(Obj)fr,vy,(Obj)ity,&t);                  substr(CO,0,(Obj)fr,vy,(Obj)ity,&t);
                 for( ix = 0, x = can->xmin; ix < w ; ix++, x += xstep ) {                  MKReal(can->xmin,rx1);
                         MKReal(x,rx);                  for(ix=0,x=can->xmin; ix<can->width; ix++,x+=xstep){
                           rx=rx1;
                         MKReal(x+xstep,rx1);                          MKReal(x+xstep,rx1);
                         istoitv((Num)rx,(Num)rx1,&itx);                          istoitv((Num)(rx1),(Num)rx,&itx);
                         substr(CO,0,(Obj)t,vx,(Obj)itx,&g);                          substr(CO,0,(Obj)fr,vx,(Obj)itx,&t);
                         if (compnum(0,0,(Num)g)) mask[iy][ix] = -1;                          MKReal(can->ymin,ry1);
                         else mask[iy][ix] = 0;                          for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
                                   ry=ry1;
                                   MKReal(y+ystep,ry1);
                                   istoitv((Num)ry,(Num)ry1,&ity);
                                   substr(CO,0,(Obj)t,vy,(Obj)ity,&g);
                                   if(compnum(0,0,(Num)g))mask[ix][iy]=-1;
                                   else {
                                           mask[ix][iy]=0;
   /*
                                           if(idv==0) mask[ix][iy]=0;
                                           else {
                                                   MKReal(y,rdy1);
                                                   for(idy=0,dy=y;idy<idv;dy+=wy,idy++){
                                                           rdy=rdy1;
                                                           MKReal(dy+wy,rdy1);
                                                           istoitv((Num)rdy,(Num)rdy1,&ddy);
                                                           substr(CO,0,(Obj)fr,vy,(Obj)ddy,&t);
                                                           MKReal(x,rdx1);
                                                           for(idx=0,dx=x;idx<idx;dx+=wx,idx++){
                                                                   rdx=rdx1;
                                                                   MKReal(dx+wx,rdx1);
                                                                   istoitv((Num)rdx,(Num)rdx1,&ddx);
                                                                   substr(CO,0,(Obj)t,vx,(Obj)ddx,&g);
                                                                   if(!compnum(0,0,(Num)g)){
                                                                           mask[ix][iy]=0;
                                                                           break;
                                                                   }
                                                           }
                                                           if(mask[ix][iy]==0)break;
                                                   }
                                           }
   */
                                   }
                           }
                 }                  }
         }          }
 }  }
 #endif  
   
 #if TRANSFER  
 void itvcalc2(int **mask, struct canvas *can, int nox)  
 {  
         double x, y, xstep, ystep;  
         int ix, iy, w, h;  
         Itv ity, itx;  
         Real r, rx, ry, rzero;  
         Obj fr, g, t, s;  
         V vx, vy;  
         P mp, fr2;  
         Real qx,qy;  
   
         if ( !nox ) initmarker(can,"Evaluating...");  
         MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);  
         vx = can->vx; vy = can->vy;  
         w = can->width; h = can->height;  
         xstep = (can->xmax - can->xmin)/w;  
         ystep = (can->ymax - can->ymin)/h;  
   
         MKReal(0.0, rzero);  
         MKReal(ystep,ry);  
         istoitv((Num)rzero,(Num)ry,&ity);  
         MKReal(xstep,rx);  
         istoitv((Num)rzero,(Num)rx,&itx);  
         for( iy = 0, y = can->ymin; iy < h ; iy++, y += ystep ) {  
                 MKReal(y,qy);  
                 MKV(vy,mp);subp(CO,mp,(P)qy,&mp);  
                 substp(CO,(P)fr,(V)vy,(P)mp,&fr2);  
                 substr(CO,0,(Obj)fr2,vy,(Obj)ity,&t);  
                 for( ix = 0, x = can->xmin; ix < w ; ix++, x += xstep ) {  
                         MKReal(x,qx);  
                         MKV(vx,mp);addp(CO,mp,(P)qx,&mp);  
                         substp(CO,(P)t,(V)vx,(P)mp,(P*)&s);  
                         substr(CO,0,(Obj)s,vx,(Obj)itx,&g);  
                         if (compnum(0,0,(Num)g)) mask[iy][ix] = -1;  
                         else mask[iy][ix] = 0;  
                 }  
         }  
 }  
 #endif  
   
 #if RECURSION  
 void reccalc3(P fr, V vx, V vy, int ixlw, int ixhg, int iylw, int iyhg,  
         double *xval, double *yval,int w, int **mask, int itvsize)  
 {  
         int i, j, ixmd, iymd;  
         double inf, sup, dt;  
         Real rxlw, rxhg, rylw, ryhg;  
         Itv itx, ity;  
         Obj fr1;  
   
         NEWItvP(itx); NEWItvP(ity);  
   
         if ( (ixlw > ixhg) || (iylw > iyhg) ) return;  
   
         MKReal(xval[ixlw], rxlw); MKReal(xval[ixhg + 1], rxhg);  
         MKReal(yval[iylw], rylw); MKReal(yval[iyhg + 1], ryhg);  
         istoitv((Num)rxlw, (Num)rxhg, &itx);  
         istoitv((Num)rylw, (Num)ryhg, &ity);  
         substr(CO, 0, (Obj)fr, vy, (Obj)ity, &fr1);  
         substr(CO, 0, (Obj)fr1, vx, (Obj)itx, &fr1);  
         Num2double((Num)fr1, &inf, &sup);  
   
         if( itvsize <= 0 ) itvsize = 1;  
         if(inf > sup){dt = inf; inf = sup; sup = dt;}  
         if ( inf <= 0.0 && sup >= 0.0 ) {  
                 if ( (ixhg - ixlw < itvsize) || (iyhg - iylw < itvsize) ){  
                         for(i = iylw; i <=iyhg; i++)  
                                 for(j = ixlw; j <=ixhg; j++) mask[w-i][j] = 0;  
                 } else {  
                         ixmd = (ixhg + ixlw)/2; iymd = (iyhg + iylw)/2;  
                         reccalc3((P)fr,vx,vy,ixlw,ixmd,iylw,iymd,xval,yval,w,mask,itvsize);  
                         reccalc3((P)fr,vx,vy,ixmd+1,ixhg,iylw,iymd,xval,yval,w,mask,itvsize);  
                         reccalc3((P)fr,vx,vy,ixlw,ixmd,iymd+1,iyhg,xval,yval,w,mask,itvsize);  
                         reccalc3((P)fr,vx,vy,ixmd+1,ixhg,iymd+1,iyhg,xval,yval,w,mask,itvsize);  
                 }  
         } else {  
                 for (i = iylw; i<= iyhg; i++)  
                         for (j = ixlw; j <= ixhg; j++) mask[w-i][j] = -1;  
         }  
 }  
   
 void itvcalc3(int **mask, struct canvas *can, int nox, int itvsize)  
 {  
         double xstep, ystep, xv, yv, *xval, *yval;  
         int i, imx, imy, w, h;  
         Real r;  
         Obj fr;  
         V vx, vy;  
   
         if ( !nox ) initmarker(can,"Evaluating...");  
         MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);  
         vx = can->vx; vy = can->vy;  
         w = can->width; h = can->height;  
         xstep = (can->xmax - can->xmin)/w;  
         ystep = (can->ymax - can->ymin)/h;  
   
         xval = (double *)ALLOCA((w+1)*sizeof(double));  
         yval = (double *)ALLOCA((h+1)*sizeof(double));  
         for (i = 0, xv = can->xmin; i <= w; i++, xv += xstep) xval[i] = xv;  
         for (i = 0, yv = can->ymin; i <= h; i++, yv += ystep) yval[i] = yv;  
         imx = w/2; imy = h/2;  
   
         reccalc3((P)fr,vx,vy,0,imx,0,imy,xval,yval,w-1,mask,itvsize);  
         reccalc3((P)fr,vx,vy,imx+1,w-1,0,imy,xval,yval,w-1,mask,itvsize);  
         reccalc3((P)fr,vx,vy,0,imx,imy+1,h-1,xval,yval,w-1,mask,itvsize);  
         reccalc3((P)fr,vx,vy,imx+1,w-1,imy+1,h-1,xval,yval,w-1,mask,itvsize);  
 }  
 #endif  
   
 #if RECTRANS  
 void reccalc4(P fr, V vx, V vy, int ixlw, int ixhg, int iylw, int iyhg,  
         double *xval, double *yval, int w, int **mask, int itvsize)  
 {  
         int i, j;  
         int ixmd, iymd;  
         double inf, sup, dt;  
         Real rxlw, rylw, rzero, rxwid, rywid;  
         Itv itx, ity;  
         Obj fr1;  
         P py, px, tmp;  
   
         if ( (ixlw > ixhg) || (iylw > iyhg)) return;  
         NEWItvP(itx); NEWItvP(ity);  
   
         MKReal(0.0,rzero);  
         MKReal(xval[ixhg + 1] - xval[ixlw], rxwid);  
         MKReal(yval[iyhg + 1] - yval[iylw], rywid);  
         istoitv((Num)rzero, (Num)rxwid, &itx);  
         istoitv((Num)rzero, (Num)rywid, &ity);  
   
         MKReal(xval[ixlw], rxlw); MKReal(yval[iylw], rylw);  
         MKV(vx,tmp);addp(CO, (P)tmp, (P)rxlw, &px);  
         MKV(vy,tmp);addp(CO, (P)tmp, (P)rylw, &py);  
         substp(CO, (P)fr, (V)vx, (P)px, (P*)&fr1);  
         substp(CO, (P)fr1, (V)vy, (P)py, (P*)&fr1);  
         substr(CO, 0, (Obj)fr1, vx, (Obj)itx, &fr1);  
         substr(CO, 0, (Obj)fr1, vy, (Obj)ity, &fr1);  
         Num2double((Num)fr1, &inf, &sup);  
   
         if(inf > sup){dt=inf;inf=sup;sup=dt;}  
         if ( inf <= 0.0 && sup >= 0.0 ) {  
                 if ( (ixhg - ixlw <= itvsize) && (iyhg - iylw <= itvsize) ){  
                         for(i = iylw; i <= iyhg; i++)  
                                 for(j = ixlw; j <= ixhg; j++) mask[w-i][j]=0;  
                 } else {  
                         ixmd = (ixhg + ixlw)/2; iymd = (iyhg + iylw)/2;  
   
                         reccalc4((P)fr,vx,vy,ixlw,ixmd,iylw,iymd,xval,yval,w,mask,itvsize);  
                         reccalc4((P)fr,vx,vy,ixmd+1,ixhg,iylw,iymd,xval,yval,w,mask,itvsize);  
                         reccalc4((P)fr,vx,vy,ixlw,ixmd,iymd+1,iyhg,xval,yval,w,mask,itvsize);  
                         reccalc4((P)fr,vx,vy,ixmd+1,ixhg,iymd+1,iyhg,xval,yval,w,mask,itvsize);  
                 }  
         } else {  
                 for (i = iylw; i<= iyhg; i++)  
                         for (j = ixlw; j <= ixhg; j++) mask[w-i][j] = -1;  
         }  
 }  
   
 void itvcalc4(int **mask, struct canvas *can, int nox, int itvsize)  
 {  
         double xstep, ystep, xv, yv, *xval, *yval;  
         int i, imx, imy, w, h;  
         Real r;  
         Obj fr;  
         V vx, vy;  
   
         if ( !nox ) initmarker(can,"Evaluating...");  
         MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);  
         vx = can->vx; vy = can->vy;  
         w = can->width; h = can->height;  
         xstep = (can->xmax - can->xmin)/w;  
         ystep = (can->ymax - can->ymin)/h;  
   
         xval = (double *)ALLOCA((w+1)*sizeof(double));  
         yval = (double *)ALLOCA((h+1)*sizeof(double));  
         for (i=0, xv=can->xmin; i <= w; i++, xv += xstep) xval[i] = xv;  
         for (i=0, yv=can->ymin; i <= h; i++, yv += ystep) yval[i] = yv;  
         imx = w/2; imy = h/2;  
   
         reccalc4((P)fr,vx,vy,0,imx,0,imy,xval,yval,w-1,mask,itvsize);  
         reccalc4((P)fr,vx,vy,imx+1,w-1,0,imy,xval,yval,w-1,mask,itvsize);  
         reccalc4((P)fr,vx,vy,0,imx,imy+1,h-1,xval,yval,w-1,mask,itvsize);  
         reccalc4((P)fr,vx,vy,imx+1,w-1,imy+1,h-1,xval,yval,w-1,mask,itvsize);  
 }  
 #endif  
 #endif  #endif
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.9
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>