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