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