Annotation of OpenXM_contrib2/asir2000/plot/calc.c, Revision 1.7
1.2 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
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 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.7 ! ohara 48: * $OpenXM: OpenXM_contrib2/asir2000/plot/calc.c,v 1.6 2002/08/02 08:59:47 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
1.5 noro 51: #include "parse.h"
1.1 noro 52: #include "ifplot.h"
53: #include <math.h>
1.7 ! ohara 54: #if defined(PARI)
1.1 noro 55: #include "genpari.h"
56: #endif
57:
1.6 noro 58: #ifndef MAXSHORT
59: #define MAXSHORT ((short)0x7fff)
60: #endif
61:
1.5 noro 62: void calc(double **tab,struct canvas *can,int nox)
1.1 noro 63: {
64: double x,y,xmin,ymin,xstep,ystep;
65: int ix,iy;
66: Real r,rx,ry;
67: Obj fr,g;
68: int w,h;
69: V vx,vy;
70: Obj t,s;
71:
1.4 noro 72: if ( !nox ) initmarker(can,"Evaluating...");
1.1 noro 73: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
74: vx = can->vx;
75: vy = can->vy;
76: w = can->width; h = can->height;
77: xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;
78: ymin = can->ymin; ystep = (can->ymax-can->ymin)/h;
79: MKReal(1.0,rx); MKReal(1.0,ry); /* dummy real */
80: for( ix = 0, x = xmin; ix < w ; ix++, x += xstep ) {
81: #if 0
82: MKReal(x,r); substp(CO,fr,vx,(P)r,&g);
1.4 noro 83: if ( !nox ) marker(can,DIR_X,ix);
1.1 noro 84: for( iy = 0, y = ymin; iy < h ; iy++, y += ystep )
85: tab[ix][iy] = usubstrp(g,y);
86: #endif
1.5 noro 87: BDY(rx) = x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);
88: devalr(CO,t,&g);
1.4 noro 89: if ( !nox ) marker(can,DIR_X,ix);
1.1 noro 90: for( iy = 0, y = ymin; iy < h ; iy++, y += ystep ) {
91: BDY(ry) = y;
1.5 noro 92: substr(CO,0,g,vy,y?(Obj)ry:0,&t);
93: devalr(CO,t,&s);
1.1 noro 94: tab[ix][iy] = ToReal(s);
95: }
96: }
97: }
98:
1.5 noro 99: double usubstrp(P p,double r)
1.1 noro 100: {
101: double t;
102: DCP dc;
103: int d;
104: double pwrreal0();
105:
106: if ( !p )
107: t = 0.0;
108: else if ( NUM(p) )
109: t = BDY((Real)p);
110: else {
111: dc = DC(p); t = BDY((Real)COEF(dc));
112: for ( d = QTOS(DEG(dc)), dc = NEXT(dc); dc;
113: d = QTOS(DEG(dc)), dc = NEXT(dc) ) {
114: t = t*pwrreal0(r,(d-QTOS(DEG(dc))))+BDY((Real)COEF(dc));
115: }
116: if ( d )
117: t *= pwrreal0(r,d);
118: }
119: return t;
120: }
121:
1.5 noro 122: void qcalc(char **tab,struct canvas *can)
1.1 noro 123: {
1.5 noro 124: Q dx,dy,w,h,xstep,ystep,c,q1;
1.1 noro 125: P g,g1,f1,f2,x,y,t,s;
126: int ix,iy;
127: int *a,*pa;
128: VECT ss;
129:
130:
131: subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
132: subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
133: MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
134: addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
135: MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
136: addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
137: ptozp(f2,1,&c,&g);
138: a = (int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
139: initmarker(can,"Horizontal scan...");
140: for( ix = 0; ix < can->width; ix++ ) {
141: marker(can,DIR_X,ix);
142: STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
143: if ( !g1 )
144: for ( iy = 0; iy < can->height; iy++ )
145: tab[ix][iy] = 1;
146: else if ( !NUM(g1) ) {
1.5 noro 147: sturmseq(CO,g1,&ss); seproot(ss,0,can->height,a);
1.1 noro 148: for ( iy = 0, pa = a; iy < can->height; iy++, pa++ )
149: if ( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
150: tab[ix][iy] = 1;
151: }
152: }
153: initmarker(can,"Vertical scan...");
154: for( iy = 0; iy < can->height; iy++ ) {
155: marker(can,DIR_Y,iy);
156: STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
157: if ( !g1 )
158: for ( ix = 0; ix < can->width; ix++ )
159: tab[ix][iy] = 1;
160: else if ( !NUM(g1) ) {
1.5 noro 161: sturmseq(CO,g1,&ss); seproot(ss,0,can->width,a);
1.1 noro 162: for ( ix = 0, pa = a; ix < can->width; ix++, pa++ )
163: if ( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
164: tab[ix][iy] = 1;
165: }
166: }
167: }
168:
1.5 noro 169: void sturmseq(VL vl,P p,VECT *rp)
1.1 noro 170: {
171: P g1,g2,q,r,s;
172: P *t;
173: V v;
174: VECT ret;
175: int i,j;
176: Q a,b,c,d,h,l,m,x;
177:
178: v = VR(p); t = (P *)ALLOCA((deg(v,p)+1)*sizeof(P));
179: g1 = t[0] = p; diffp(vl,p,v,(P *)&a); ptozp((P)a,1,&c,&g2); t[1] = g2;
180: for ( i = 1, h = ONE, x = ONE; ; ) {
181: if ( NUM(g2) )
182: break;
183: subq(DEG(DC(g1)),DEG(DC(g2)),&d);
184: l = (Q)LC(g2);
185: if ( SGN(l) < 0 ) {
186: chsgnq(l,&a); l = a;
187: }
188: addq(d,ONE,&a); pwrq(l,a,&b); mulp(vl,(P)b,g1,(P *)&a);
189: divsrp(vl,(P)a,g2,&q,&r);
190: if ( !r )
191: break;
192: chsgnp(r,&s); r = s; i++;
193: if ( NUM(r) ) {
194: t[i] = r; break;
195: }
196: pwrq(h,d,&m); g1 = g2;
197: mulq(m,x,&a); divsp(vl,r,(P)a,&g2); t[i] = g2;
198: x = (Q)LC(g1);
199: if ( SGN(x) < 0 ) {
200: chsgnq(x,&a); x = a;
201: }
202: pwrq(x,d,&a); mulq(a,h,&b); divq(b,m,&h);
203: }
204: MKVECT(ret,i+1);
205: for ( j = 0; j <= i; j++ )
206: ret->body[j] = (pointer)t[j];
207: *rp = ret;
208: }
209:
1.5 noro 210: void seproot(VECT s,int min,int max,int *ar)
1.1 noro 211: {
212: P f;
213: P *ss;
214: Q q,t;
215: int i,j,k;
216:
217: ss = (P *)s->body; f = ss[0];
218: for ( i = min; i <= max; i++ ) {
219: STOQ(i,q); usubstqp(f,q,&t);
220: if ( !t )
221: ar[i] = -1;
222: else {
223: ar[i] = numch(s,q,t); break;
224: }
225: }
226: if ( i > max )
227: return;
228: for ( j = max; j >= min; j-- ) {
229: STOQ(j,q); usubstqp(f,q,&t);
230: if ( !t )
231: ar[j] = -1;
232: else {
233: if ( i != j )
234: ar[j] = numch(s,q,t);
235: break;
236: }
237: }
238: if ( j <= i+1 )
239: return;
240: if ( ar[i] == ar[j] ) {
241: for ( k = i+1; k < j; k++ )
242: ar[k] = ar[i];
243: return;
244: }
245: k = (i+j)/2;
246: seproot(s,i,k,ar);
247: seproot(s,k,j,ar);
248: }
249:
1.5 noro 250: int numch(VECT s,Q n,Q a0)
1.1 noro 251: {
252: int len,i,c;
253: Q a;
254: P *ss;
255:
256: len = s->len; ss = (P *)s->body;
257: for ( i = 1, c = 0; i < len; i++ ) {
258: usubstqp(ss[i],n,&a);
259: if ( a ) {
260: if ( (SGN(a)>0 && SGN(a0)<0) || (SGN(a)<0 && SGN(a0)>0) )
261: c++;
262: a0 = a;
263: }
264: }
265: return c;
266: }
267:
1.5 noro 268: void usubstqp(P p,Q r,Q *v)
1.1 noro 269: {
270: Q d,d1,a,b,t;
271: DCP dc;
272:
273: if ( !p )
274: *v = 0;
275: else if ( NUM(p) )
276: *v = (Q)p;
277: else {
278: dc = DC(p); t = (Q)COEF(dc);
279: for ( d = DEG(dc), dc = NEXT(dc); dc;
280: d = DEG(dc), dc = NEXT(dc) ) {
281: subq(d,DEG(dc),&d1); pwrq(r,d1,&a);
282: mulq(t,a,&b); addq(b,(Q)COEF(dc),&t);
283: }
284: if ( d ) {
285: pwrq(r,d,&a); mulq(t,a,&b); t = b;
286: }
287: *v = t;
288: }
289: }
290:
1.5 noro 291: void plotcalc(struct canvas *can)
1.1 noro 292: {
293: double x,xmin,xstep,ymax,ymin,dy;
294: int ix;
295: Real r;
296: Obj fr;
297: double usubstrp();
298: int w,h;
299: double *tab;
300: POINT *pa;
301: Real rx;
302: Obj t,s;
303:
304: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
305: w = can->width; h = can->height;
306: xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;
307: tab = (double *)ALLOCA(w*sizeof(double));
308: MKReal(1,rx); /* dummy real number */
309: for( ix = 0, x = xmin; ix < w ; ix++, x += xstep ) {
310: /* full substitution */
311: BDY(rx) = x;
1.5 noro 312: substr(CO,0,fr,can->vx,x?(Obj)rx:0,&s);
313: devalr(CO,(Obj)s,&t);
1.1 noro 314: if ( t && (OID(t)!=O_N || NID((Num)t)!=N_R) )
315: error("plotcalc : invalid evaluation");
316: tab[ix] = ToReal((Num)t);
317: #if 0
318: tab[ix] = usubstrp(fr,x);
319: #endif
320: }
321: if ( can->ymax == can->ymin ) {
322: for ( ymax = ymin = tab[0], ix = 1; ix < w; ix++ ) {
323: if ( tab[ix] > ymax )
324: ymax = tab[ix];
325: if ( tab[ix] < ymin )
326: ymin = tab[ix];
327: }
328: can->ymax = ymax; can->ymin = ymin;
329: } else {
330: ymax = can->ymax; ymin = can->ymin;
331: }
332: dy = ymax-ymin;
333: can->pa = (struct pa *)MALLOC(sizeof(struct pa));
334: can->pa[0].length = w;
335: can->pa[0].pos = pa = (POINT *)MALLOC(w*sizeof(POINT));
336: for ( ix = 0; ix < w; ix++ ) {
337: double t;
338:
339: XC(pa[ix]) = ix;
340: t = (h - 1)*(ymax - tab[ix])/dy;
341: if ( t > MAXSHORT )
342: YC(pa[ix]) = MAXSHORT;
343: else if ( t < -MAXSHORT )
344: YC(pa[ix]) = -MAXSHORT;
345: else
1.5 noro 346: YC(pa[ix]) = (long)t;
1.6 noro 347: }
348: }
349:
350: void polarplotcalc(struct canvas *can)
351: {
352: double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep;
353: int i,nstep;
354: double usubstrp();
355: int w,h;
356: double tr,p;
357: double *tabx,*taby;
358: POINT *pa;
359: Real r;
360: Obj fr,t,s;
361:
362: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
363: w = can->width; h = can->height; nstep = can->nzstep;
364: pmin = can->zmin; pstep = (can->zmax-can->zmin)/nstep;
365: tabx = (double *)ALLOCA(nstep*sizeof(double));
366: taby = (double *)ALLOCA(nstep*sizeof(double));
367: MKReal(1,r); /* dummy real number */
368: for( i = 0, p = pmin; i < nstep ; i++, p += pstep ) {
369: /* full substitution */
370: BDY(r) = p;
371: substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);
372: devalr(CO,(Obj)s,&t);
373: if ( t && (OID(t)!=O_N || NID((Num)t)!=N_R) )
374: error("polarplotcalc : invalid evaluation");
375: tr = ToReal((Num)t);
376: tabx[i] = tr*cos(p);
377: taby[i] = tr*sin(p);
378: }
379: xmax = xmin = tabx[0];
380: ymax = ymin = taby[0];
381: for ( i = 1; i < nstep; i++ ) {
382: if ( tabx[i] > xmax ) xmax = tabx[i];
383: if ( tabx[i] < xmin ) xmin = tabx[i];
384: if ( taby[i] > ymax ) ymax = taby[i];
385: if ( taby[i] < ymin ) ymin = taby[i];
386: }
387: can->xmax = xmax; can->xmin = xmin;
388: can->ymax = ymax; can->ymin = ymin;
389: dx = xmax-xmin;
390: dy = ymax-ymin;
391: can->pa = (struct pa *)MALLOC(sizeof(struct pa));
392: can->pa[0].length = nstep;
393: can->pa[0].pos = pa = (POINT *)MALLOC(w*sizeof(POINT));
394: for ( i = 0; i < nstep; i++ ) {
395: XC(pa[i]) = (w-1)*(tabx[i]-xmin)/dx;
396: YC(pa[i]) = (h-1)*(ymax-taby[i])/dy;
1.1 noro 397: }
398: }
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 / asir2000 / plot