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

Annotation of OpenXM_contrib2/asir2018/plot/calc.c, Revision 1.1

1.1     ! noro        1: /*
        !             2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
        !             3:  * All rights reserved.
        !             4:  *
        !             5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
        !             6:  * non-exclusive and royalty-free license to use, copy, modify and
        !             7:  * redistribute, solely for non-commercial and non-profit purposes, the
        !             8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
        !             9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
        !            10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
        !            11:  * third party developer retains all rights, including but not limited to
        !            12:  * copyrights, in and to the SOFTWARE.
        !            13:  *
        !            14:  * (1) FLL does not grant you a license in any way for commercial
        !            15:  * purposes. You may use the SOFTWARE only for non-commercial and
        !            16:  * non-profit purposes only, such as academic, research and internal
        !            17:  * business use.
        !            18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
        !            19:  * international copyright treaties. If you make copies of the SOFTWARE,
        !            20:  * with or without modification, as permitted hereunder, you shall affix
        !            21:  * to all such copies of the SOFTWARE the above copyright notice.
        !            22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
        !            23:  * shall be made on your publication or presentation in any form of the
        !            24:  * results obtained by use of the SOFTWARE.
        !            25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
        !            26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
        !            27:  * for such modification or the source code of the modified part of the
        !            28:  * SOFTWARE.
        !            29:  *
        !            30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
        !            31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
        !            32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
        !            33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
        !            34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
        !            35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
        !            36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
        !            37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
        !            38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
        !            39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
        !            40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
        !            41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
        !            42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
        !            43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
        !            44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
        !            45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
        !            46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
        !            47:  *
        !            48:  * $OpenXM$
        !            49: */
        !            50: #include "ca.h"
        !            51: #include "parse.h"
        !            52: #include "ifplot.h"
        !            53: #include <math.h>
        !            54: #if defined(PARI)
        !            55: #include "genpari.h"
        !            56: #endif
        !            57:
        !            58: #ifndef MAXSHORT
        !            59: #define MAXSHORT ((short)0x7fff)
        !            60: #endif
        !            61:
        !            62: void calc(double **tab,struct canvas *can,int nox){
        !            63:   //memory_plot,IFPLOTD,INEQND,INEQNANDD,INEQNORD
        !            64:   //INEQNXORD,conplotmainD
        !            65:   double x,y,xstep,ystep;
        !            66:   int ix,iy;
        !            67:   Real r,rx,ry;
        !            68:   Obj fr,g,t,s;
        !            69:
        !            70:   if(!nox)initmarker(can,"Evaluating...");
        !            71:   todouble((Obj)can->formula,(Obj *)&fr);
        !            72:   xstep=(can->xmax-can->xmin)/can->width;
        !            73:   ystep=(can->ymax-can->ymin)/can->height;
        !            74:   MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
        !            75:   BDY(rx)=can->xmin;
        !            76:   substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&t); devalr(CO,t,&g);
        !            77:   BDY(ry)=can->ymin;
        !            78:   substr(CO,0,g,can->vy,can->ymin?(Obj)ry:0,&t); devalr(CO,t,&s);
        !            79:   can->vmax=can->vmin=ToReal(s);
        !            80:   for(ix=0,x=can->xmin; ix<can->width; ix++,x+=xstep){
        !            81:     BDY(rx)=x; substr(CO,0,fr,can->vx,x?(Obj)rx:0,&t);
        !            82:     devalr(CO,t,&g);
        !            83:     if(!nox)marker(can,DIR_X,ix);
        !            84:     for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
        !            85:       BDY(ry)=y;
        !            86:       substr(CO,0,g,can->vy,y?(Obj)ry:0,&t);
        !            87:       devalr(CO,t,&s);
        !            88:       tab[ix][iy]=ToReal(s);
        !            89:       if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
        !            90:       if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
        !            91:     }
        !            92:   }
        !            93: }
        !            94:
        !            95: void calcq(double **tab,struct canvas *can,int nox){
        !            96:   //IFPLOTQ,INEQNQ,INEQNANDQ,INEQNORQ,INEQNXORQ
        !            97:   //plotoverD
        !            98:   Z w,h,q1;
        !            99:   Q dx,dy,xstep,ystep,c;
        !           100:   P g,g1,f1,f2,x,y;
        !           101:   int ix,iy;
        !           102:   Obj fr,gm,t,s;
        !           103:   Real r,rx,ry;
        !           104:
        !           105:   todouble((Obj)can->formula,&fr);
        !           106:   MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
        !           107:   BDY(rx)=can->xmin;
        !           108:   substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&t); devalr(CO,t,&gm);
        !           109:   BDY(ry)=can->ymin;
        !           110:   substr(CO,0,gm,can->vy,can->ymin?(Obj)ry:0,&t); devalr(CO,t,&s);
        !           111:   can->vmax=can->vmin=ToReal(s);
        !           112:
        !           113:   subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,(Q)w,&xstep);
        !           114:   subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,(Q)h,&ystep);
        !           115:   MKV(can->vx,x); mulp(CO,(P)xstep,x,(P *)&t);
        !           116:   addp(CO,(P)can->qxmin,(P)t,(P *)&s); substp(CO,can->formula,can->vx,(P)s,&f1);
        !           117:   MKV(can->vy,y); mulp(CO,(P)ystep,y,(P *)&t);
        !           118:   addp(CO,(P)can->qymin,(P)t,(P *)&s); substp(CO,f1,can->vy,(P)s,&f2);
        !           119:   ptozp(f2,1,&c,&g);
        !           120:   if(!nox) initmarker(can,"Evaluating...");
        !           121:   for(iy=0;iy<can->height;iy++){
        !           122:     marker(can,DIR_Y,iy);
        !           123:     STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,(P *)&t); ptozp((P)t,1,&c,&g1);
        !           124:     for(ix=0;ix<can->width;ix++){
        !           125:       STOQ(ix,q1);substp(CO,g1,can->vx,(P)q1,(P *)&t);
        !           126:       devalr(CO,t,&s);
        !           127:       tab[ix][iy]=ToReal(s);
        !           128:       if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
        !           129:       if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
        !           130:     }
        !           131:   }
        !           132: }
        !           133:
        !           134: void calcb(double **tab,struct canvas *can,int nox){
        !           135:   //IFPLOTB,INEQNB,INEQNANDB,INEQNORB,INEQNXORB
        !           136:   Z w,h,q1;
        !           137:   Q dx,dy,xstep,ystep,c;
        !           138:   P g,g1,f1,f2,x,y,t,s;
        !           139:   int ix,iy,*a,*pa;
        !           140:   VECT ss;
        !           141:   Obj fr,gm,tm,sm;
        !           142:   Real r,rx,ry;
        !           143:
        !           144:   todouble((Obj)can->formula,&fr);
        !           145:   MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
        !           146:   BDY(rx)=can->xmin;
        !           147:   substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&tm); devalr(CO,tm,&gm);
        !           148:   BDY(ry)=can->ymin;
        !           149:   substr(CO,0,gm,can->vy,can->ymin?(Obj)ry:0,&tm); devalr(CO,tm,&sm);
        !           150:   can->vmax=can->vmin=ToReal(sm);
        !           151:
        !           152:   for(iy=0;iy<can->height;iy++)for(ix=0;ix<can->width;ix++)tab[ix][iy]=1.0;
        !           153:   subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,(Q)w,&xstep);
        !           154:   subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,(Q)h,&ystep);
        !           155:   MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
        !           156:   addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
        !           157:   MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
        !           158:   addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
        !           159:   ptozp(f2,1,&c,&g);
        !           160:   a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
        !           161:   for(iy=0;iy<can->height;iy++)for(ix=0;ix<can->width;ix++)tab[ix][iy]=1.0;
        !           162:   subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,(Q)w,&xstep);
        !           163:   subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,(Q)h,&ystep);
        !           164:   MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
        !           165:   addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
        !           166:   MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
        !           167:   addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
        !           168:   ptozp(f2,1,&c,&g);
        !           169:   a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
        !           170:   for(ix=0;ix<can->width;ix++){
        !           171:     STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
        !           172:     if(!g1)for(iy=0;iy<can->height;iy++)tab[ix][iy]=0.0;
        !           173:     else if(!NUM(g1)){
        !           174:       sturmseq(CO,g1,&ss);
        !           175:       seproot(ss,0,can->width,a);
        !           176:       for(iy=0,pa=a;iy<can->height;iy++,pa++){
        !           177:         if(*pa<0||(*(pa+1)>=0&&(*pa>*(pa+1))))tab[ix][iy]=0.0;
        !           178:         else {
        !           179:           STOQ(iy,q1);substp(CO,g1,can->vy,(P)q1,&t);
        !           180:           devalr(CO,(Obj)t,(Obj *)&s);
        !           181:           tab[ix][iy]=ToReal(s);
        !           182:           if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
        !           183:           if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
        !           184:         }
        !           185:       }
        !           186:     }
        !           187:   }
        !           188:   for(iy=0;iy<can->height;iy++){
        !           189:     STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
        !           190:     if(!g1) for(ix=0;ix<can->width;ix++)tab[ix][iy]=0.0;
        !           191:     else if(!NUM(g1)){
        !           192:       sturmseq(CO,g1,&ss);
        !           193:       seproot(ss,0,can->height,a);
        !           194:       for(ix=0,pa=a;ix<can->width;ix++,pa++){
        !           195:         if(tab[ix][iy]!=0.0){
        !           196:           if(*pa<0||(*(pa+1)>=0&&(*pa>*(pa+1))))tab[ix][iy]=0.0;
        !           197:           else {
        !           198:             STOQ(ix,q1);substp(CO,g1,can->vx,(P)q1,&t);
        !           199:             devalr(CO,(Obj)t,(Obj *)&s);
        !           200:             tab[ix][iy]=ToReal(s);
        !           201:             if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
        !           202:             if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
        !           203:           }
        !           204:         }
        !           205:       }
        !           206:     }
        !           207:   }
        !           208: }
        !           209:
        !           210: double usubstrp(P p,double r){
        !           211:   DCP dc;
        !           212:   int d;
        !           213:   double t,pwrreal0();
        !           214:
        !           215:   if(!p) t=0.0;
        !           216:   else if(NUM(p))t=BDY((Real)p);
        !           217:   else {
        !           218:     dc=DC(p); t=BDY((Real)COEF(dc));
        !           219:     for(d=QTOS(DEG(dc)),dc=NEXT(dc);dc;d=QTOS(DEG(dc)),dc=NEXT(dc)){
        !           220:       t=t*pwrreal0(r,(d-QTOS(DEG(dc))))+BDY((Real)COEF(dc));
        !           221:     }
        !           222:     if(d)t*=pwrreal0(r,d);
        !           223:   }
        !           224:   return t;
        !           225: }
        !           226:
        !           227: void qcalc(char **tab,struct canvas *can){
        !           228:   //qifplotmain(Old type)
        !           229:   Z w,h,q1;
        !           230:   Q dx,dy,xstep,ystep,c;
        !           231:   P g,g1,f1,f2,x,y,t,s;
        !           232:   int ix,iy;
        !           233:   int *a,*pa;
        !           234:   VECT ss;
        !           235:
        !           236:   subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,(Q)w,&xstep);
        !           237:   subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,(Q)h,&ystep);
        !           238:   MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
        !           239:   addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
        !           240:   MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
        !           241:   addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
        !           242:   ptozp(f2,1,&c,&g);
        !           243:   a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
        !           244:   initmarker(can,"Horizontal scan...");
        !           245:   for( ix=0; ix < can->width; ix++ ){
        !           246:     marker(can,DIR_X,ix);
        !           247:     STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
        !           248:     if( !g1 )
        !           249:       for(iy=0; iy < can->height; iy++ )
        !           250:         tab[ix][iy]=1;
        !           251:     else if( !NUM(g1) ){
        !           252:       sturmseq(CO,g1,&ss); seproot(ss,0,can->height,a);
        !           253:       for(iy=0, pa=a; iy < can->height; iy++, pa++ )
        !           254:         if( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
        !           255:           tab[ix][iy]=1;
        !           256:     }
        !           257:   }
        !           258:   initmarker(can,"Vertical scan...");
        !           259:   for( iy=0; iy < can->height; iy++ ){
        !           260:     marker(can,DIR_Y,iy);
        !           261:     STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
        !           262:     if( !g1 )
        !           263:       for(ix=0; ix < can->width; ix++ )
        !           264:         tab[ix][iy]=1;
        !           265:     else if( !NUM(g1) ){
        !           266:       sturmseq(CO,g1,&ss); seproot(ss,0,can->width,a);
        !           267:       for(ix=0, pa=a; ix < can->width; ix++, pa++ )
        !           268:         if( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
        !           269:           tab[ix][iy]=1;
        !           270:     }
        !           271:   }
        !           272: }
        !           273:
        !           274: void sturmseq(VL vl,P p,VECT *rp){
        !           275:   P g1,g2,q,r,s,*t;
        !           276:   V v;
        !           277:   VECT ret;
        !           278:   int i,j;
        !           279:   Q b,c,l,m,h,x;
        !           280:   Z d,a;
        !           281:
        !           282:   v=VR(p);t=(P *)ALLOCA((deg(v,p)+1)*sizeof(P));
        !           283:   g1=t[0]=p;diffp(vl,p,v,(P *)&a);ptozp((P)a,1,&c,&g2);t[1]=g2;
        !           284:   for(i=1,h=(Q)ONE,x=(Q)ONE;;){
        !           285:     if(NUM(g2)) break;
        !           286:     subz(DEG(DC(g1)),DEG(DC(g2)),&d);
        !           287:     l=(Q)LC(g2);
        !           288:     if(sgnq(l)<0){
        !           289:       chsgnq(l,&b);l=b;
        !           290:     }
        !           291:     addz(d,ONE,&a);pwrq(l,(Q)a,&b);mulp(vl,(P)b,g1,&s);
        !           292:     divsrp(vl,s,g2,&q,&r);
        !           293:     if(!r) break;
        !           294:     chsgnp(r,&s);r=s;i++;
        !           295:     if(NUM(r)){
        !           296:       t[i]=r;break;
        !           297:     }
        !           298:     pwrq(h,(Q)d,&m);g1=g2;
        !           299:     mulq(m,x,(Q *)&a);divsp(vl,r,(P)a,&g2);t[i]=g2;
        !           300:     x=(Q)LC(g1);
        !           301:     if(sgnq(x)<0){
        !           302:       chsgnq(x,&b);x=b;
        !           303:     }
        !           304:     pwrq(x,(Q)d,(Q *)&a);mulq((Q)a,h,&b);divq(b,m,&h);
        !           305:   }
        !           306:   MKVECT(ret,i+1);
        !           307:   for(j=0;j<=i;j++)
        !           308:     ret->body[j]=(pointer)t[j];
        !           309:   *rp=ret;
        !           310: }
        !           311:
        !           312: void seproot(VECT s,int min,int max,int *ar){
        !           313:   P f,*ss;
        !           314:   Q t;
        !           315:   Z q;
        !           316:   int i,j,k;
        !           317:
        !           318:   ss=(P *)s->body;f=ss[0];
        !           319:   for(i=min;i<=max;i++){
        !           320:     STOQ(i,q);usubstqp(f,(Q)q,&t);
        !           321:     if(!t)ar[i]=-1;
        !           322:     else {
        !           323:       ar[i]=numch(s,(Q)q,t);break;
        !           324:     }
        !           325:   }
        !           326:   if(i>max) return;
        !           327:   for(j=max;j>= min;j--){
        !           328:     STOQ(j,q); usubstqp(f,(Q)q,&t);
        !           329:     if(!t)ar[j]=-1;
        !           330:     else {
        !           331:       if(i!=j)ar[j]=numch(s,(Q)q,t);
        !           332:       break;
        !           333:     }
        !           334:   }
        !           335:   if(j<=i+1) return;
        !           336:   if(ar[i]==ar[j]){
        !           337:     for(k=i+1;k<j;k++)ar[k]=ar[i];
        !           338:     return;
        !           339:   }
        !           340:   k=(i+j)/2;
        !           341:   seproot(s,i,k,ar);
        !           342:   seproot(s,k,j,ar);
        !           343: }
        !           344:
        !           345: int numch(VECT s,Q n,Q a0){
        !           346:   int len,i,c;
        !           347:   Q a;
        !           348:   P *ss;
        !           349:
        !           350:   len=s->len;ss=(P *)s->body;
        !           351:   for(i=1,c=0;i<len;i++){
        !           352:     usubstqp(ss[i],n,&a);
        !           353:     if(a){
        !           354:       if((sgnq(a)>0 && sgnq(a0)<0)||(sgnq(a)<0&&sgnq(a0)>0))c++;
        !           355:       a0=a;
        !           356:     }
        !           357:   }
        !           358:   return c;
        !           359: }
        !           360:
        !           361: void usubstqp(P p,Q r,Q *v){
        !           362:   Q a,b,t;
        !           363:   Z d,d1;
        !           364:   DCP dc;
        !           365:
        !           366:   if(!p)
        !           367:     *v=0;
        !           368:   else if(NUM(p))*v=(Q)p;
        !           369:   else {
        !           370:     dc=DC(p);t=(Q)COEF(dc);
        !           371:     for(d=DEG(dc),dc=NEXT(dc);dc;d=DEG(dc),dc=NEXT(dc)){
        !           372:       subz(d,DEG(dc),&d1);pwrq(r,(Q)d1,&a);
        !           373:       mulq(t,a,&b);addq(b,(Q)COEF(dc),&t);
        !           374:     }
        !           375:     if(d){
        !           376:       pwrq(r,(Q)d,&a);mulq(t,a,&b);t=b;
        !           377:     }
        !           378:     *v=t;
        !           379:   }
        !           380: }
        !           381:
        !           382: Num tobf(Num,int);
        !           383: void Psetprec(NODE arg,Obj *rp);
        !           384:
        !           385: void plotcalcbf(struct canvas *can){
        !           386:   Obj fr,s,t;
        !           387:   Num xmin,xmax,ymin,ymax,xstep;
        !           388:   Num u,v,ha,dx,dy,x;
        !           389:   Num *tab;
        !           390:   Real r;
        !           391:   int ix;
        !           392:   POINT *pa;
        !           393:   double rr;
        !           394:   Z prec,w,h1;
        !           395:   NODE arg;
        !           396:
        !           397:   STOQ(can->prec,prec); arg = mknode(1,prec); Psetprec(arg,&t);
        !           398:   evalr(CO,(Obj)can->formula,can->prec,&fr);
        !           399:   MKReal(can->xmin,r); xmin = tobf((Num)r,can->prec);
        !           400:   MKReal(can->xmax,r); xmax = tobf((Num)r,can->prec);
        !           401:   MKReal(can->ymin,r); ymin = tobf((Num)r,can->prec);
        !           402:   MKReal(can->ymax,r); ymax = tobf((Num)r,can->prec);
        !           403:   STOQ(can->width,w);
        !           404:   subbf(xmax,xmin,&dx); divbf(dx,(Num)w,&xstep);
        !           405:   tab=(Num *)MALLOC(can->width*sizeof(Num));
        !           406:   for(ix=0,x=xmin;ix<can->width;ix++){
        !           407:     substr(CO,0,fr,can->vx,(Obj)x,(Obj *)&s);
        !           408:     evalr(CO,(Obj)s,can->prec,&t);
        !           409:     if(t&&(OID(t)!=O_N))
        !           410:       error("plotcalcbf : invalid evaluation");
        !           411:     tab[ix]=(Num)t;
        !           412:     addbf(x,xstep,&u); x = u;
        !           413:   }
        !           414:   if(!cmpbf(ymax,ymin)){
        !           415:     for(ymax=ymin=tab[0],ix=1;ix<can->width;ix++){
        !           416:       if(cmpbf(tab[ix],ymax)>0)ymax=tab[ix];
        !           417:       if(cmpbf(tab[ix],ymin)<0)ymin=tab[ix];
        !           418:     }
        !           419:     can->ymax=ToReal(ymax);can->ymin=ToReal(ymin);
        !           420:   }
        !           421:   subbf(ymax,ymin,&dy);
        !           422:   can->pa=(struct pa *)MALLOC(sizeof(struct pa));
        !           423:   can->pa[0].length=can->width;
        !           424:   can->pa[0].pos=pa=(POINT *)MALLOC(can->width*sizeof(POINT));
        !           425:   STOQ(can->height-1,h1);
        !           426:   for(ix=0;ix<can->width;ix++){
        !           427:     XC(pa[ix])=ix;
        !           428:     subbf(ymax,tab[ix],&u); divbf(u,dy,&v); mulbf(v,(Num)h1,&u);
        !           429:     rr = ToReal(u);
        !           430:     if(rr>MAXSHORT)YC(pa[ix])=MAXSHORT;
        !           431:     else if(rr<-MAXSHORT)YC(pa[ix])=-MAXSHORT;
        !           432:     else YC(pa[ix])=(long)rr;
        !           433:   }
        !           434: }
        !           435:
        !           436: void plotcalc(struct canvas *can){
        !           437:   //plot,memory_plot,plotover,plot_resize
        !           438:   double x,xmin,xstep,ymax,ymin,dy,*tab,usubstrp();
        !           439:   int ix,w,h;
        !           440:   Real r,rx;
        !           441:   Obj fr,t,s;
        !           442:   POINT *pa;
        !           443:
        !           444:   if ( can->prec ) {
        !           445:     plotcalcbf(can);
        !           446:     return;
        !           447:   }
        !           448:   todouble((Obj)can->formula,&fr);
        !           449:   w=can->width;h=can->height;
        !           450:   xmin=can->xmin;xstep=(can->xmax-can->xmin)/w;
        !           451:   tab=(double *)ALLOCA(w*sizeof(double));
        !           452:   MKReal(1,rx); // dummy real number
        !           453:   for(ix=0,x=xmin;ix<w;ix++,x+=xstep){
        !           454:     // full substitution
        !           455:     BDY(rx)=x;
        !           456:     substr(CO,0,fr,can->vx,x?(Obj)rx:0,&s);
        !           457:     devalr(CO,(Obj)s,&t);
        !           458:     if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
        !           459:       error("plotcalc : invalid evaluation");
        !           460:     tab[ix]=ToReal((Num)t);
        !           461:   }
        !           462:   if(can->ymax==can->ymin){
        !           463:     for(ymax=ymin=tab[0],ix=1;ix<w;ix++){
        !           464:       if(tab[ix]>ymax)ymax=tab[ix];
        !           465:       if(tab[ix]<ymin)ymin=tab[ix];
        !           466:     }
        !           467:     can->ymax=ymax;can->ymin=ymin;
        !           468:   } else {
        !           469:     ymax=can->ymax;ymin=can->ymin;
        !           470:   }
        !           471:   dy=ymax-ymin;
        !           472:   can->pa=(struct pa *)MALLOC(sizeof(struct pa));
        !           473:   can->pa[0].length=w;
        !           474:   can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
        !           475:   for(ix=0;ix<w;ix++){
        !           476:     double t;
        !           477:     XC(pa[ix])=ix;
        !           478:     t=(h-1)*(ymax-tab[ix])/dy;
        !           479:     if(t>MAXSHORT)YC(pa[ix])=MAXSHORT;
        !           480:     else if(t<-MAXSHORT)YC(pa[ix])=-MAXSHORT;
        !           481:     else YC(pa[ix])=(long)t;
        !           482:   }
        !           483: }
        !           484:
        !           485: void polarcalc(struct canvas *can){
        !           486:   double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep,tr,p,*tabx,*taby;
        !           487:   double usubstrp();
        !           488:   int i,nstep,w,h;
        !           489:   POINT *pa;
        !           490:   Real r;
        !           491:   Obj fr,t,s;
        !           492:
        !           493:   todouble((Obj)can->formula,&fr);
        !           494:   w=can->width; h=can->height; nstep=can->nzstep;
        !           495:   pmin=can->zmin; pstep=(can->zmax-can->zmin)/nstep;
        !           496:   tabx=(double *)ALLOCA(nstep*sizeof(double));
        !           497:   taby=(double *)ALLOCA(nstep*sizeof(double));
        !           498:   MKReal(1,r); // dummy real number
        !           499:
        !           500:   for(i=0,p=pmin;i<nstep;i++,p+= pstep){
        !           501:     // full substitution
        !           502:     BDY(r)=p;
        !           503:     substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);
        !           504:     devalr(CO,(Obj)s,&t);
        !           505:     if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
        !           506:       error("polarcalc : invalid evaluation");
        !           507:     tr=ToReal((Num)t);
        !           508:     tabx[i]=tr*cos(p);
        !           509:     taby[i]=tr*sin(p);
        !           510:   }
        !           511:   xmax=xmin=tabx[0];
        !           512:   ymax=ymin=taby[0];
        !           513:   for(i=1;i<nstep;i++){
        !           514:     if(tabx[i]>xmax)xmax=tabx[i];
        !           515:     if(tabx[i]<xmin)xmin=tabx[i];
        !           516:     if(taby[i]>ymax)ymax=taby[i];
        !           517:     if(taby[i]<ymin)ymin=taby[i];
        !           518:   }
        !           519:   can->xmax=xmax;can->xmin=xmin;
        !           520:   can->ymax=ymax;can->ymin=ymin;
        !           521:   dx=xmax-xmin;
        !           522:   dy=ymax-ymin;
        !           523:   can->pa=(struct pa *)MALLOC(sizeof(struct pa));
        !           524:   can->pa[0].length=nstep;
        !           525:   can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
        !           526:   for(i=0;i<nstep;i++){
        !           527:     XC(pa[i])=(w-1)*(tabx[i]-xmin)/dx;
        !           528:     YC(pa[i])=(h-1)*(ymax-taby[i])/dy;
        !           529:   }
        !           530: }
        !           531:
        !           532: void polarcalcNG(struct canvas *can){
        !           533:   //polarplotNG
        !           534:   double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep,tr,p, *tabx,*taby;
        !           535:   double usubstrp();
        !           536:   int i,ix,iy,nstep,w,h;
        !           537:   POINT *pa;
        !           538:   Real r;
        !           539:   Obj fr,t,s;
        !           540:
        !           541:   todouble((Obj)can->formula,&fr);
        !           542:   w=can->width; h=can->height; nstep=can->nzstep;
        !           543:   pmin=can->zmin; pstep=(can->zmax-can->zmin)/nstep;
        !           544:   tabx=(double *)ALLOCA(nstep*sizeof(double));
        !           545:   taby=(double *)ALLOCA(nstep*sizeof(double));
        !           546:   MKReal(1,r); // dummy real number
        !           547:
        !           548:   for(i=0,p=pmin;i<nstep;i++,p+= pstep){
        !           549:     // full substitution
        !           550:     BDY(r)=p;
        !           551:     substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);
        !           552:     devalr(CO,(Obj)s,&t);
        !           553:     if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
        !           554:       error("polarcalc : invalid evaluation");
        !           555:     tr=ToReal((Num)t);
        !           556:     tabx[i]=tr*cos(p);
        !           557:     taby[i]=tr*sin(p);
        !           558:     if(i==0){
        !           559:       xmax=xmin=tabx[0];
        !           560:       ymax=ymin=taby[0];
        !           561:     } else {
        !           562:       if(tabx[i]>xmax)xmax=tabx[i];
        !           563:       if(tabx[i]<xmin)xmin=tabx[i];
        !           564:       if(taby[i]>ymax)ymax=taby[i];
        !           565:       if(taby[i]<ymin)ymin=taby[i];
        !           566:     }
        !           567:   }
        !           568:   can->xmax=xmax;can->xmin=xmin;
        !           569:   can->ymax=ymax;can->ymin=ymin;
        !           570:   dx=xmax-xmin;
        !           571:   dy=ymax-ymin;
        !           572:   can->pa=(struct pa *)MALLOC(sizeof(struct pa));
        !           573:   can->pa[0].length=nstep;
        !           574:   can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
        !           575:   for(i=0;i<nstep;i++){
        !           576:     XC(pa[i])=(w-1)*(tabx[i]-xmin)/dx;
        !           577:     YC(pa[i])=(h-1)*(ymax-taby[i])/dy;
        !           578:   }
        !           579: }
        !           580:
        !           581: /*
        !           582: void ineqncalc(double **tab,struct canvas *can,int nox){
        !           583:   double x,y,xmin,ymin,xstep,ystep;
        !           584:   int ix,iy,w,h;
        !           585:   Real r,rx,ry;
        !           586:   Obj fr,g,t,s;
        !           587:   V vx,vy;
        !           588:
        !           589:   if(!nox) initmarker(can,"Evaluating...");
        !           590:   todouble((Obj)can->formula,&fr);
        !           591:   vx=can->vx;vy=can->vy;
        !           592:   w=can->width;h=can->height;
        !           593:   xmin=can->xmin;xstep=(can->xmax-can->xmin)/w;
        !           594:   ymin=can->ymin;ystep=(can->ymin-can->ymin)/h;
        !           595:   MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
        !           596:
        !           597:   for(ix=0,x=xmin;ix<=w;ix++,x+=xstep){
        !           598:     BDY(rx)=x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);
        !           599:     devalr(CO,t,&g);
        !           600:     if(!nox) marker(can,DIR_X,ix);
        !           601:     for(iy=0,y=ymin;iy<=h;iy++,y+=ystep){
        !           602:       BDY(ry)=y;
        !           603:       substr(CO,0,g,vy,y?(Obj)ry:0,&t);
        !           604:       devalr(CO,t,&s);
        !           605:       tab[ix][iy]=ToReal(s);
        !           606:     }
        !           607:   }
        !           608: }
        !           609: */
        !           610:
        !           611: #if defined(INTERVAL)
        !           612: void itvcalc(double **mask, struct canvas *can, int nox){
        !           613:   //ITVIFPLOT
        !           614:   double x,y,xstep,ystep,dx,dy,wx,wy;
        !           615:   int idv,ix,iy,idx,idy;
        !           616:   Itv ity,itx,ddx,ddy;
        !           617:   Real r,rx,ry,rx1,ry1,rdx,rdy,rdx1,rdy1;
        !           618:   V vx,vy;
        !           619:   Obj fr,g,t,s;
        !           620:
        !           621:   idv=can->division;
        !           622:   todouble((Obj)can->formula,&fr);
        !           623:   vx=can->vx; vy=can->vy;
        !           624:   xstep=(can->xmax-can->xmin)/can->width;
        !           625:   ystep=(can->ymax-can->ymin)/can->height;
        !           626:   if(idv!=0){
        !           627:     wx=xstep/can->division;
        !           628:     wy=ystep/can->division;
        !           629:   }
        !           630:   MKReal(can->ymin,ry1);
        !           631:   for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
        !           632:     ry=ry1;
        !           633:     MKReal(y+ystep,ry1);
        !           634:     istoitv((Num)(ry1),(Num)ry,&ity);
        !           635:     substr(CO,0,(Obj)fr,vy,(Obj)ity,&t);
        !           636:     MKReal(can->xmin,rx1);
        !           637:     for(ix=0,x=can->xmin; ix<can->width; ix++,x+=xstep){
        !           638:       rx=rx1;
        !           639:       MKReal(x+xstep,rx1);
        !           640:       istoitv((Num)(rx1),(Num)rx,&itx);
        !           641:       substr(CO,0,(Obj)fr,vx,(Obj)itx,&t);
        !           642:       MKReal(can->ymin,ry1);
        !           643:       for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
        !           644:         ry=ry1;
        !           645:         MKReal(y+ystep,ry1);
        !           646:         istoitv((Num)ry,(Num)ry1,&ity);
        !           647:         substr(CO,0,(Obj)t,vy,(Obj)ity,&g);
        !           648:         if(compnum(0,0,(Num)g))mask[ix][iy]=-1;
        !           649:         else {
        !           650:           mask[ix][iy]=0;
        !           651: /*
        !           652:           if(idv==0) mask[ix][iy]=0;
        !           653:           else {
        !           654:             MKReal(y,rdy1);
        !           655:             for(idy=0,dy=y;idy<idv;dy+=wy,idy++){
        !           656:               rdy=rdy1;
        !           657:               MKReal(dy+wy,rdy1);
        !           658:               istoitv((Num)rdy,(Num)rdy1,&ddy);
        !           659:               substr(CO,0,(Obj)fr,vy,(Obj)ddy,&t);
        !           660:               MKReal(x,rdx1);
        !           661:               for(idx=0,dx=x;idx<idx;dx+=wx,idx++){
        !           662:                 rdx=rdx1;
        !           663:                 MKReal(dx+wx,rdx1);
        !           664:                 istoitv((Num)rdx,(Num)rdx1,&ddx);
        !           665:                 substr(CO,0,(Obj)t,vx,(Obj)ddx,&g);
        !           666:                 if(!compnum(0,0,(Num)g)){
        !           667:                   mask[ix][iy]=0;
        !           668:                   break;
        !           669:                 }
        !           670:               }
        !           671:               if(mask[ix][iy]==0)break;
        !           672:             }
        !           673:           }
        !           674: */
        !           675:         }
        !           676:       }
        !           677:     }
        !           678:   }
        !           679: }
        !           680: #endif

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>