![[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: }