Annotation of OpenXM_contrib2/asir2018/plot/calc.c, Revision 1.2
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: *
1.2 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2018/plot/calc.c,v 1.1 2018/09/19 05:45:08 noro Exp $
1.1 noro 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:
1.2 ! noro 113: subq(can->qxmax,can->qxmin,&dx); STOZ(can->width,w); divq(dx,(Q)w,&xstep);
! 114: subq(can->qymax,can->qymin,&dy); STOZ(can->height,h); divq(dy,(Q)h,&ystep);
1.1 noro 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);
1.2 ! noro 123: STOZ(iy,q1); substp(CO,g,can->vy,(P)q1,(P *)&t); ptozp((P)t,1,&c,&g1);
1.1 noro 124: for(ix=0;ix<can->width;ix++){
1.2 ! noro 125: STOZ(ix,q1);substp(CO,g1,can->vx,(P)q1,(P *)&t);
1.1 noro 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;
1.2 ! noro 153: subq(can->qxmax,can->qxmin,&dx); STOZ(can->width,w); divq(dx,(Q)w,&xstep);
! 154: subq(can->qymax,can->qymin,&dy); STOZ(can->height,h); divq(dy,(Q)h,&ystep);
1.1 noro 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;
1.2 ! noro 162: subq(can->qxmax,can->qxmin,&dx); STOZ(can->width,w); divq(dx,(Q)w,&xstep);
! 163: subq(can->qymax,can->qymin,&dy); STOZ(can->height,h); divq(dy,(Q)h,&ystep);
1.1 noro 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++){
1.2 ! noro 171: STOZ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
1.1 noro 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 {
1.2 ! noro 179: STOZ(iy,q1);substp(CO,g1,can->vy,(P)q1,&t);
1.1 noro 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++){
1.2 ! noro 189: STOZ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
1.1 noro 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 {
1.2 ! noro 198: STOZ(ix,q1);substp(CO,g1,can->vx,(P)q1,&t);
1.1 noro 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));
1.2 ! noro 219: for(d=ZTOS(DEG(dc)),dc=NEXT(dc);dc;d=ZTOS(DEG(dc)),dc=NEXT(dc)){
! 220: t=t*pwrreal0(r,(d-ZTOS(DEG(dc))))+BDY((Real)COEF(dc));
1.1 noro 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:
1.2 ! noro 236: subq(can->qxmax,can->qxmin,&dx); STOZ(can->width,w); divq(dx,(Q)w,&xstep);
! 237: subq(can->qymax,can->qymin,&dy); STOZ(can->height,h); divq(dy,(Q)h,&ystep);
1.1 noro 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);
1.2 ! noro 247: STOZ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
1.1 noro 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);
1.2 ! noro 261: STOZ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
1.1 noro 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++){
1.2 ! noro 320: STOZ(i,q);usubstqp(f,(Q)q,&t);
1.1 noro 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--){
1.2 ! noro 328: STOZ(j,q); usubstqp(f,(Q)q,&t);
1.1 noro 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:
1.2 ! noro 397: STOZ(can->prec,prec); arg = mknode(1,prec); Psetprec(arg,&t);
1.1 noro 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);
1.2 ! noro 403: STOZ(can->width,w);
1.1 noro 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));
1.2 ! noro 425: STOZ(can->height-1,h1);
1.1 noro 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>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to calc.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / plot