[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.2 and 1.8

version 1.2, 2000/08/21 08:31:50 version 1.8, 2011/08/10 04:51:58
Line 23 
Line 23 
  * shall be made on your publication or presentation in any form of the   * shall be made on your publication or presentation in any form of the
  * results obtained by use of the SOFTWARE.   * results obtained by use of the SOFTWARE.
  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by   * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
  * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification   * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
  * for such modification or the source code of the modified part of the   * for such modification or the source code of the modified part of the
  * SOFTWARE.   * SOFTWARE.
  *   *
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.1.1.1 1999/12/03 07:39:12 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/plot/calc.c,v 1.7 2003/02/14 22:29:19 ohara Exp $
 */  */
 #include "ca.h"  #include "ca.h"
   #include "parse.h"
 #include "ifplot.h"  #include "ifplot.h"
 #include <math.h>  #include <math.h>
 #if PARI  #if defined(PARI)
 #include "genpari.h"  #include "genpari.h"
 #endif  #endif
   
 double usubstrp(P,double);  #ifndef MAXSHORT
   #define MAXSHORT ((short)0x7fff)
   #endif
   
 void calc(tab,can)  void calc(double **tab,struct canvas *can,int nox)
 double **tab;  
 struct canvas *can;  
 {  {
         double x,y,xmin,ymin,xstep,ystep;          double x,y,xmin,ymin,xstep,ystep;
         int ix,iy;          int ix,iy;
Line 68  struct canvas *can;
Line 69  struct canvas *can;
         V vx,vy;          V vx,vy;
         Obj t,s;          Obj t,s;
   
         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;          vx = can->vx;
         vy = can->vy;          vy = can->vy;
Line 77  struct canvas *can;
Line 78  struct canvas *can;
         ymin = can->ymin; ystep = (can->ymax-can->ymin)/h;          ymin = can->ymin; ystep = (can->ymax-can->ymin)/h;
         MKReal(1.0,rx); MKReal(1.0,ry); /* dummy real */          MKReal(1.0,rx); MKReal(1.0,ry); /* dummy real */
         for( ix = 0, x = xmin; ix < w ; ix++, x += xstep ) {          for( ix = 0, x = xmin; ix < w ; ix++, x += xstep ) {
 #if 0                  BDY(rx) = x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);
                 MKReal(x,r); substp(CO,fr,vx,(P)r,&g);                  devalr(CO,t,&g);
                 marker(can,DIR_X,ix);                  if ( !nox ) marker(can,DIR_X,ix);
                 for( iy = 0, y = ymin; iy < h ; iy++, y += ystep )  
                         tab[ix][iy] = usubstrp(g,y);  
 #endif  
                 BDY(rx) = x; substr(CO,0,fr,vx,x?(P)rx:0,&t); devalr(CO,t,&g);  
                 marker(can,DIR_X,ix);  
                 for( iy = 0, y = ymin; iy < h ; iy++, y += ystep ) {                  for( iy = 0, y = ymin; iy < h ; iy++, y += ystep ) {
                         BDY(ry) = y;                          BDY(ry) = y;
                         substr(CO,0,g,vy,y?(P)ry:0,&t); devalr(CO,t,&s);                          substr(CO,0,g,vy,y?(Obj)ry:0,&t);
                           devalr(CO,t,&s);
                         tab[ix][iy] = ToReal(s);                          tab[ix][iy] = ToReal(s);
                 }                  }
         }          }
 }  }
   
 double usubstrp(p,r)  double usubstrp(P p,double r)
 P p;  
 double r;  
 {  {
         double t;          double t;
         DCP dc;          DCP dc;
Line 118  double r;
Line 113  double r;
         return t;          return t;
 }  }
   
 void qcalc(tab,can)  void qcalc(char **tab,struct canvas *can)
 char **tab;  
 struct canvas *can;  
 {  {
         Q dx,dy,w,h,xstep,ystep,c,q1,q2;          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;
         DCP dc;  
         int ix,iy;          int ix,iy;
         int *a,*pa;          int *a,*pa;
         char *px;  
         VECT ss;          VECT ss;
   
   
         subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);          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);          subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
         MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);          MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
Line 147  struct canvas *can;
Line 137  struct canvas *can;
                         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) ) {
                         strum(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;
Line 161  struct canvas *can;
Line 151  struct canvas *can;
                         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) ) {
                         strum(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;
Line 169  struct canvas *can;
Line 159  struct canvas *can;
         }          }
 }  }
   
 strum(vl,p,rp)  void sturmseq(VL vl,P p,VECT *rp)
 VL vl;  
 P p;  
 VECT *rp;  
 {  {
         P g1,g2,q,r,s;          P g1,g2,q,r,s;
         P *t;          P *t;
Line 213  VECT *rp;
Line 200  VECT *rp;
         *rp = ret;          *rp = ret;
 }  }
   
 seproot(s,min,max,ar)  void seproot(VECT s,int min,int max,int *ar)
 VECT s;  
 int min,max;  
 int *ar;  
 {  {
         P f;          P f;
         P *ss;          P *ss;
Line 256  int *ar;
Line 240  int *ar;
         seproot(s,k,j,ar);          seproot(s,k,j,ar);
 }  }
   
 numch(s,n,a0)  int numch(VECT s,Q n,Q a0)
 VECT s;  
 Q n,a0;  
 {  {
         int len,i,c;          int len,i,c;
         Q a;          Q a;
Line 276  Q n,a0;
Line 258  Q n,a0;
         return c;          return c;
 }  }
   
 usubstqp(p,r,v)  void usubstqp(P p,Q r,Q *v)
 P p;  
 Q r;  
 Q *v;  
 {  {
         Q d,d1,a,b,t;          Q d,d1,a,b,t;
         DCP dc;          DCP dc;
Line 302  Q *v;
Line 281  Q *v;
         }          }
 }  }
   
 void plotcalc(can)  void plotcalc(struct canvas *can)
 struct canvas *can;  
 {  {
         double x,xmin,xstep,ymax,ymin,dy;          double x,xmin,xstep,ymax,ymin,dy;
         int ix;          int ix;
Line 324  struct canvas *can;
Line 302  struct canvas *can;
         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?(P)rx:0,&s); devalr(CO,(Obj)s,&t);                  substr(CO,0,fr,can->vx,x?(Obj)rx:0,&s);
                   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 0  
                 tab[ix] = usubstrp(fr,x);  
 #endif  
         }          }
         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++ ) {
Line 348  struct canvas *can;
Line 324  struct canvas *can;
         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++ ) {
 #ifndef MAXSHORT  
 #define MAXSHORT ((short)0x7fff)  
 #endif  
                 double t;                  double t;
   
                 XC(pa[ix]) = ix;                  XC(pa[ix]) = ix;
Line 360  struct canvas *can;
Line 333  struct canvas *can;
                 else if ( t < -MAXSHORT )                  else if ( t < -MAXSHORT )
                         YC(pa[ix]) = -MAXSHORT;                          YC(pa[ix]) = -MAXSHORT;
                 else                  else
                         YC(pa[ix]) = t;                          YC(pa[ix]) = (long)t;
         }          }
 }  }
   
   void polarplotcalc(struct canvas *can)
   {
           double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep;
           int i,nstep;
           double usubstrp();
           int w,h;
           double tr,p;
           double *tabx,*taby;
           POINT *pa;
           Real r;
           Obj fr,t,s;
   
           MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
           w = can->width; h = can->height; nstep = can->nzstep;
           pmin = can->zmin; pstep = (can->zmax-can->zmin)/nstep;
           tabx = (double *)ALLOCA(nstep*sizeof(double));
           taby = (double *)ALLOCA(nstep*sizeof(double));
           MKReal(1,r); /* dummy real number */
           for( i = 0, p = pmin; i < nstep ; i++, p += pstep ) {
                   /* full substitution */
                   BDY(r) = p;
                   substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);
                   devalr(CO,(Obj)s,&t);
                   if ( t && (OID(t)!=O_N || NID((Num)t)!=N_R) )
                           error("polarplotcalc : invalid evaluation");
                   tr = ToReal((Num)t);
                   tabx[i] = tr*cos(p);
                   taby[i] = tr*sin(p);
           }
           xmax = xmin = tabx[0];
           ymax = ymin = taby[0];
           for ( i = 1; i < nstep; i++ ) {
                   if ( tabx[i] > xmax ) xmax = tabx[i];
                   if ( tabx[i] < xmin ) xmin = tabx[i];
                   if ( taby[i] > ymax ) ymax = taby[i];
                   if ( taby[i] < ymin ) ymin = taby[i];
           }
           can->xmax = xmax; can->xmin = xmin;
           can->ymax = ymax; can->ymin = ymin;
           dx = xmax-xmin;
           dy = ymax-ymin;
           can->pa = (struct pa *)MALLOC(sizeof(struct pa));
           can->pa[0].length = nstep;
           can->pa[0].pos = pa = (POINT *)MALLOC(w*sizeof(POINT));
           for ( i = 0; i < nstep; i++ ) {
                   XC(pa[i]) = (w-1)*(tabx[i]-xmin)/dx;
                   YC(pa[i]) = (h-1)*(ymax-taby[i])/dy;
           }
   }
   
   #if defined(INTERVAL)
   #define NORMAL 1
   #define TRANSFER 1
   #define RECURSION 1
   #define RECTRANS 1
   
   void ineqncalc(double **tab,struct canvas *can,int nox)
   {
           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);
           vx = can->vx; vy = can->vy;
           w = can->width; h = can->height;
           xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;
           ymax = can->ymax; ystep = (can->ymax-can->ymin)/h;
           MKReal(1.0,rx); MKReal(1.0,ry); /* dummy real */
   
           for( iy = 0, y = ymax; iy <= h ; iy++, y -= ystep ) {
                   BDY(ry) = y; substr(CO,0,fr,vy,y?(Obj)ry:0,&t);
                   devalr(CO,t,&g);
                   if ( !nox ) marker(can,DIR_Y,iy);
                   for( ix = 0, x = xmin; ix <= w ; ix++, x += xstep ) {
                           BDY(rx) = x;
                           substr(CO,0,g,vx,x?(Obj)rx:0,&t);
                           devalr(CO,t,&s);
                           tab[iy][ix] = ToReal(s);
                   }
           }
   }
   
   void reccalc(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 rxwid, rywid, rxlw, rylw, rzero;
           Itv itx, ity;
           Obj fr1;
           P py, px, tmp;
   
           double **tbl , xs, ys;
           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);
           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;
   
           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);
                   substr(CO,0,(Obj)fr,vy,(Obj)ity,&t);
                   for( ix = 0, x = can->xmin; ix < w ; ix++, x += xstep ) {
                           MKReal(x,rx);
                           MKReal(x+xstep,rx1);
                           istoitv((Num)rx,(Num)rx1,&itx);
                           substr(CO,0,(Obj)t,vx,(Obj)itx,&g);
                           if (compnum(0,0,(Num)g)) mask[iy][ix] = -1;
                           else mask[iy][ix] = 0;
                   }
           }
   }
   #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
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.8
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>