Annotation of OpenXM_contrib2/asir2000/builtin/itvnum.c, Revision 1.6
1.1 saito 1: /*
1.6 ! saito 2: * $OpenXM: OpenXM_contrib2/asir2000/builtin/itvnum.c,v 1.5 2005/07/14 22:46:03 kondoh Exp $
1.1 saito 3: */
4:
5: #include "ca.h"
6: #include "parse.h"
7: #include "version.h"
1.6 ! saito 8: #include "../plot/ifplot.h"
1.1 saito 9:
10: #if defined(INTERVAL)
11:
12: static void Pitv(NODE, Obj *);
13: static void Pitvd(NODE, Obj *);
14: static void Pitvbf(NODE, Obj *);
15: static void Pinf(NODE, Obj *);
16: static void Psup(NODE, Obj *);
17: static void Pmid(NODE, Obj *);
18: static void Pabsitv(NODE, Obj *);
19: static void Pdisjitv(NODE, Obj *);
20: static void Pinitv(NODE, Obj *);
21: static void Pcup(NODE, Obj *);
22: static void Pcap(NODE, Obj *);
23: static void Pwidth(NODE, Obj *);
24: static void Pdistance(NODE, Obj *);
1.3 saito 25: static void Pitvversion(Q *);
1.1 saito 26: #endif
27: static void Pprintmode(NODE, Obj *);
28:
1.6 ! saito 29: /* plot time check func */
! 30: static void ccalc(double **, struct canvas *, int);
! 31: static void Pifcheck(NODE, Obj *);
! 32:
1.1 saito 33: #if defined(__osf__) && 0
34: int end;
35: #endif
36:
37: struct ftab interval_tab[] = {
38: {"printmode",Pprintmode,1},
39: #if defined(INTERVAL)
40: {"itvd",Pitvd,-2},
41: {"intvald",Pitvd,-2},
42: {"itv",Pitv,-2},
43: {"intval",Pitv,-2},
44: {"itvbf",Pitvbf,-2},
45: {"intvalbf",Pitvbf,-2},
46: {"inf",Pinf,1},
47: {"sup",Psup,1},
48: {"absintval",Pabsitv,1},
49: {"disintval",Pdisjitv,2},
50: {"inintval",Pinitv,2},
51: {"cup",Pcup,2},
52: {"cap",Pcap,2},
53: {"mid",Pmid,1},
54: {"width",Pwidth,1},
55: {"diam",Pwidth,1},
56: {"distance",Pdistance,2},
57: {"iversion",Pitvversion,0},
1.6 ! saito 58: /* plot time check */
! 59: {"ifcheck",Pifcheck,-7},
1.1 saito 60: #endif
61: {0,0,0},
62: };
63:
64: #if defined(INTERVAL)
1.6 ! saito 65:
! 66: /* plot time check */
! 67: static void
! 68: Pifcheck(NODE arg, Obj *rp)
! 69: {
! 70: Q m2,p2,s_id;
! 71: NODE defrange;
! 72: LIST xrange,yrange,range[2],list,geom;
! 73: VL vl,vl0;
! 74: V v[2],av[2];
! 75: int ri,i,j,sign;
! 76: P poly;
! 77: P var;
! 78: NODE n,n0;
! 79: Obj t;
! 80:
! 81: struct canvas *can;
! 82: MAT m;
! 83: pointer **mb;
! 84: double **tabe, *px, *px1, *px2;
! 85: Q one;
! 86: int width, height, ix, iy;
! 87: int id;
! 88:
! 89: STOQ(-2,m2); STOQ(2,p2);
! 90: STOQ(1,one);
! 91: MKNODE(n,p2,0); MKNODE(defrange,m2,n);
! 92: poly = 0; vl = 0; geom = 0; ri = 0;
! 93: v[0] = v[1] = 0;
! 94: for ( ; arg; arg = NEXT(arg) ){
! 95: switch ( OID(BDY(arg)) ) {
! 96: case O_P:
! 97: poly = (P)BDY(arg);
! 98: get_vars_recursive((Obj)poly,&vl);
! 99: for ( vl0 = vl, i = 0; vl0; vl0 = NEXT(vl0) )
! 100: if ( vl0->v->attr == (pointer)V_IND )
! 101: if ( i >= 2 )
! 102: error("ifplot : invalid argument");
! 103: else
! 104: v[i++] = vl0->v;
! 105: break;
! 106: case O_LIST:
! 107: list = (LIST)BDY(arg);
! 108: if ( OID(BDY(BDY(list))) == O_P )
! 109: if ( ri > 1 )
! 110: error("ifplot : invalid argument");
! 111: else
! 112: range[ri++] = list;
! 113: else
! 114: geom = list;
! 115: break;
! 116: default:
! 117: error("ifplot : invalid argument"); break;
! 118: }
! 119: }
! 120: if ( !poly ) error("ifplot : invalid argument");
! 121: switch ( ri ) {
! 122: case 0:
! 123: if ( !v[1] ) error("ifplot : please specify all variables");
! 124: MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n);
! 125: MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n);
! 126: break;
! 127: case 1:
! 128: if ( !v[1] ) error("ifplot : please specify all variables");
! 129: av[0] = VR((P)BDY(BDY(range[0])));
! 130: if ( v[0] == av[0] ) {
! 131: xrange = range[0];
! 132: MKV(v[1],var); MKNODE(n,var,defrange); MKLIST(yrange,n);
! 133: } else if ( v[1] == av[0] ) {
! 134: MKV(v[0],var); MKNODE(n,var,defrange); MKLIST(xrange,n);
! 135: yrange = range[0];
! 136: } else
! 137: error("ifplot : invalid argument");
! 138: break;
! 139: case 2:
! 140: av[0] = VR((P)BDY(BDY(range[0])));
! 141: av[1] = VR((P)BDY(BDY(range[1])));
! 142: if ( ((v[0] == av[0]) && (!v[1] || v[1] == av[1])) ||
! 143: ((v[0] == av[1]) && (!v[1] || v[1] == av[0])) ) {
! 144: xrange = range[0]; yrange = range[1];
! 145: } else error("ifplot : invalid argument");
! 146: break;
! 147: default:
! 148: error("ifplot : cannot happen"); break;
! 149: }
! 150: can = canvas[id = search_canvas()];
! 151: if ( !geom ) {
! 152: width = 300;
! 153: height = 300;
! 154: can->width = 300;
! 155: can->height = 300;
! 156: } else {
! 157: can->width = QTOS((Q)BDY(BDY(geom)));
! 158: can->height = QTOS((Q)BDY(NEXT(BDY(geom))));
! 159: width = can->width;
! 160: height = can->height;
! 161: }
! 162: if ( xrange ) {
! 163: n = BDY(xrange); can->vx = VR((P)BDY(n)); n = NEXT(n);
! 164: can->qxmin = (Q)BDY(n); n = NEXT(n); can->qxmax = (Q)BDY(n);
! 165: can->xmin = ToReal(can->qxmin); can->xmax = ToReal(can->qxmax);
! 166: }
! 167: if ( yrange ) {
! 168: n = BDY(yrange); can->vy = VR((P)BDY(n)); n = NEXT(n);
! 169: can->qymin = (Q)BDY(n); n = NEXT(n); can->qymax = (Q)BDY(n);
! 170: can->ymin = ToReal(can->qymin); can->ymax = ToReal(can->qymax);
! 171: }
! 172: can->wname = "ifcheck";
! 173: can->formula = poly;
! 174: tabe = (double **)ALLOCA((width+1)*sizeof(double *));
! 175: for ( i = 0; i <= width; i++ )
! 176: tabe[i] = (double *)ALLOCA((height+1)*sizeof(double));
! 177: for(i=0;i<=width;i++)for(j=0;j<=height;j++)tabe[i][j]=0;
! 178: ccalc(tabe,can,0);
! 179: MKMAT(m,width,height);
! 180: mb = BDY(m);
! 181: for( ix=0; ix<width; ix++ ){
! 182: for( iy=0; iy<height; iy++){
! 183: if ( tabe[ix][iy] >= 0 ){
! 184: if ( (tabe[ix+1][iy] <= 0) ||
! 185: (tabe[ix][iy+1] <= 0 ) ||
! 186: (tabe[ix+1][iy+1] <= 0 ) ) mb[ix][iy] = (Obj)one;
! 187: } else {
! 188: if ( (tabe[ix+1][iy] >= 0 ) ||
! 189: ( tabe[ix][iy+1] >= 0 ) ||
! 190: ( tabe[ix+1][iy+1] >= 0 )) mb[ix][iy] = (Obj)one;
! 191: }
! 192: }
! 193: }
! 194: *rp = (Obj)m;
! 195: }
! 196:
! 197: void ccalc(double **tab,struct canvas *can,int nox)
! 198: {
! 199: double x,y,xmin,ymin,xstep,ystep;
! 200: int ix,iy;
! 201: Real r,rx,ry;
! 202: Obj fr,g;
! 203: int w,h;
! 204: V vx,vy;
! 205: Obj t,s;
! 206:
! 207: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
! 208: vx = can->vx;
! 209: vy = can->vy;
! 210: w = can->width; h = can->height;
! 211: xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;
! 212: ymin = can->ymin; ystep = (can->ymax-can->ymin)/h;
! 213: MKReal(1.0,rx); MKReal(1.0,ry);
! 214: for( ix = 0, x = xmin; ix < w+1 ; ix++, x += xstep ) {
! 215: BDY(rx) = x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);
! 216: devalr(CO,t,&g);
! 217: for( iy = 0, y = ymin; iy < h+1 ; iy++, y += ystep ) {
! 218: BDY(ry) = y;
! 219: substr(CO,0,g,vy,y?(Obj)ry:0,&t);
! 220: devalr(CO,t,&s);
! 221: tab[ix][iy] = ToReal(s);
! 222: }
! 223: }
! 224: }
! 225: /* end plot time check */
! 226:
1.1 saito 227: static void
1.4 saito 228: Pitvversion(Q *rp)
1.1 saito 229: {
1.3 saito 230: STOQ(ASIR_VERSION, *rp);
1.1 saito 231: }
232:
233: extern int bigfloat;
234:
235: static void
236: Pitv(NODE arg, Obj *rp)
237: {
238: Num a, i, s;
239: Itv c;
240: double inf, sup;
241:
242: #if 1
243: if ( bigfloat )
244: Pitvbf(arg, rp);
245: else
246: Pitvd(arg,rp);
247: #else
248: asir_assert(ARG0(arg),O_N,"itv");
249: if ( argc(arg) > 1 ) {
250: asir_assert(ARG1(arg),O_N,"itv");
251: istoitv((Num)ARG0(arg),(Num)ARG1(arg),&c);
252: } else {
253: a = (Num)ARG0(arg);
254: if ( ! a ) {
255: *rp = 0;
256: return;
257: }
1.2 kondoh 258: else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat) {
1.1 saito 259: *rp = (Obj)a;
260: return;
261: }
1.2 kondoh 262: else if ( NID(a) == N_IntervalDouble ) {
263: inf = INF((IntervalDouble)a);
264: sup = SUP((IntervalDouble)a);
1.1 saito 265: double2bf(inf, (BF *)&i);
266: double2bf(sup, (BF *)&s);
267: istoitv(i,s,&c);
268: }
269: else istoitv(a,a,&c);
270: }
1.2 kondoh 271: if ( NID( c ) == N_IntervalBigFloat )
272: addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp);
1.1 saito 273: else *rp = (Obj)c;
274: #endif
275: }
276:
277: static void
278: Pitvbf(NODE arg, Obj *rp)
279: {
280: Num a, i, s;
281: Itv c;
282: BF ii,ss;
283: double inf, sup;
284:
285: asir_assert(ARG0(arg),O_N,"intvalbf");
286: a = (Num)ARG0(arg);
287: if ( argc(arg) > 1 ) {
288: asir_assert(ARG1(arg),O_N,"intvalbf");
289: i = (Num)ARG0(arg);
290: s = (Num)ARG1(arg);
291: ToBf(i, &ii);
292: ToBf(s, &ss);
293: istoitv((Num)ii,(Num)ss,&c);
294: } else {
295: if ( ! a ) {
296: *rp = 0;
297: return;
298: }
299: else if ( NID(a) == N_IP ) {
300: itvtois((Itv)a, &i, &s);
301: ToBf(i, &ii);
302: ToBf(s, &ss);
303: istoitv((Num)ii,(Num)ss,&c);
304: }
1.2 kondoh 305: else if ( NID(a) == N_IntervalBigFloat) {
1.1 saito 306: *rp = (Obj)a;
307: return;
308: }
1.2 kondoh 309: else if ( NID(a) == N_IntervalDouble ) {
310: inf = INF((IntervalDouble)a);
311: sup = SUP((IntervalDouble)a);
1.1 saito 312: double2bf(inf, (BF *)&i);
313: double2bf(sup, (BF *)&s);
314: istoitv(i,s,&c);
315: }
316: else {
317: ToBf(a, (BF *)&i);
318: istoitv(i,i,&c);
319: }
320: }
1.2 kondoh 321: if ( c && OID( c ) == O_N && NID( c ) == N_IntervalBigFloat )
322: addulp((IntervalBigFloat)c, (IntervalBigFloat *)rp);
1.1 saito 323: else *rp = (Obj)c;
324: }
325:
326: static void
327: Pitvd(NODE arg, Obj *rp)
328: {
329: double inf, sup;
330: Num a, a0, a1, t;
331: Itv ia;
1.2 kondoh 332: IntervalDouble d;
1.1 saito 333:
334: asir_assert(ARG0(arg),O_N,"intvald");
335: a0 = (Num)ARG0(arg);
336: if ( argc(arg) > 1 ) {
337: asir_assert(ARG1(arg),O_N,"intvald");
338: a1 = (Num)ARG1(arg);
339: } else {
1.2 kondoh 340: if ( a0 && OID(a0)==O_N && NID(a0)==N_IntervalDouble ) {
341: inf = INF((IntervalDouble)a0);
342: sup = SUP((IntervalDouble)a0);
343: MKIntervalDouble(inf,sup,d);
1.1 saito 344: *rp = (Obj)d;
345: return;
346: }
347: a1 = (Num)ARG0(arg);
348: }
349: if ( compnum(0,a0,a1) > 0 ) {
350: t = a0; a0 = a1; a1 = t;
351: }
352: inf = ToRealDown(a0);
353: sup = ToRealUp(a1);
1.2 kondoh 354: MKIntervalDouble(inf,sup,d);
1.1 saito 355: *rp = (Obj)d;
356: }
357:
358: static void
359: Pinf(NODE arg, Obj *rp)
360: {
361: Num a, i, s;
362: Real r;
363: double d;
364:
365: a = (Num)ARG0(arg);
366: if ( ! a ) {
367: *rp = 0;
368: } else if ( OID(a) == O_N ) {
369: switch ( NID(a) ) {
1.2 kondoh 370: case N_IntervalDouble:
371: d = INF((IntervalDouble)a);
1.1 saito 372: MKReal(d, r);
373: *rp = (Obj)r;
374: break;
375: case N_IP:
1.2 kondoh 376: case N_IntervalBigFloat:
377: case N_IntervalQuad:
1.1 saito 378: itvtois((Itv)ARG0(arg),&i,&s);
379: *rp = (Obj)i;
380: break;
1.5 kondoh 381: default:
1.1 saito 382: *rp = (Obj)a;
383: break;
384: }
385: } else {
386: *rp = (Obj)a;
387: }
388: }
389:
390: static void
391: Psup(NODE arg, Obj *rp)
392: {
393: Num a, i, s;
394: Real r;
395: double d;
396:
397: a = (Num)ARG0(arg);
398: if ( ! a ) {
399: *rp = 0;
400: } else if ( OID(a) == O_N ) {
401: switch ( NID(a) ) {
1.2 kondoh 402: case N_IntervalDouble:
403: d = SUP((IntervalDouble)a);
1.1 saito 404: MKReal(d, r);
405: *rp = (Obj)r;
406: break;
407: case N_IP:
1.2 kondoh 408: case N_IntervalBigFloat:
409: case N_IntervalQuad:
1.1 saito 410: itvtois((Itv)ARG0(arg),&i,&s);
411: *rp = (Obj)s;
412: break;
1.5 kondoh 413: default:
1.1 saito 414: *rp = (Obj)a;
415: break;
416: }
417: } else {
418: *rp = (Obj)a;
419: }
420: }
421:
422: static void
423: Pmid(NODE arg, Obj *rp)
424: {
425: Num a, s;
426: Real r;
427: double d;
428:
429: a = (Num)ARG0(arg);
430: if ( ! a ) {
431: *rp = 0;
432: } else switch (OID(a)) {
433: case O_N:
1.2 kondoh 434: if ( NID(a) == N_IntervalDouble ) {
435: d = ( INF((IntervalDouble)a)+SUP((IntervalDouble)a) ) / 2.0;
1.1 saito 436: MKReal(d, r);
437: *rp = (Obj)r;
1.2 kondoh 438: } else if ( NID(a) == N_IntervalQuad ) {
1.1 saito 439: error("mid: not supported operation");
440: *rp = 0;
1.2 kondoh 441: } else if ( NID(a) == N_IP || NID(a) == N_IntervalBigFloat ) {
1.1 saito 442: miditvp((Itv)ARG0(arg),&s);
443: *rp = (Obj)s;
444: } else {
445: *rp = (Obj)a;
446: }
447: break;
448: #if 0
449: case O_P:
450: case O_R:
451: case O_LIST:
452: case O_VECT:
453: case O_MAT:
454: #endif
1.5 kondoh 455: default:
1.1 saito 456: *rp = (Obj)a;
457: break;
458: }
459: }
460:
461: static void
462: Pcup(NODE arg, Obj *rp)
463: {
464: Itv s;
465: Num a, b;
466:
467: asir_assert(ARG0(arg),O_N,"cup");
468: asir_assert(ARG1(arg),O_N,"cup");
469: a = (Num)ARG0(arg);
470: b = (Num)ARG1(arg);
1.2 kondoh 471: if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {
472: cupitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp);
1.1 saito 473: } else {
474: cupitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);
475: *rp = (Obj)s;
476: }
477: }
478:
479: static void
480: Pcap(NODE arg, Obj *rp)
481: {
482: Itv s;
483: Num a, b;
484:
485: asir_assert(ARG0(arg),O_N,"cap");
486: asir_assert(ARG1(arg),O_N,"cap");
487: a = (Num)ARG0(arg);
488: b = (Num)ARG1(arg);
1.2 kondoh 489: if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {
490: capitvd((IntervalDouble)a, (IntervalDouble)b, (IntervalDouble *)rp);
1.1 saito 491: } else {
492: capitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);
493: *rp = (Obj)s;
494: }
495: }
496:
497: static void
498: Pwidth(arg,rp)
499: NODE arg;
500: Obj *rp;
501: {
502: Num s;
503: Num a;
504:
505: asir_assert(ARG0(arg),O_N,"width");
506: a = (Num)ARG0(arg);
507: if ( ! a ) {
508: *rp = 0;
1.2 kondoh 509: } else if ( NID(a) == N_IntervalDouble ) {
510: widthitvd((IntervalDouble)a, (Num *)rp);
1.1 saito 511: } else {
512: widthitvp((Itv)ARG0(arg),&s);
513: *rp = (Obj)s;
514: }
515: }
516:
517: static void
518: Pabsitv(arg,rp)
519: NODE arg;
520: Obj *rp;
521: {
522: Num s;
523: Num a, b;
524:
525: asir_assert(ARG0(arg),O_N,"absitv");
526: a = (Num)ARG0(arg);
527: if ( ! a ) {
528: *rp = 0;
1.2 kondoh 529: } else if ( NID(a) == N_IntervalDouble ) {
530: absitvd((IntervalDouble)a, (Num *)rp);
1.1 saito 531: } else {
532: absitvp((Itv)ARG0(arg),&s);
533: *rp = (Obj)s;
534: }
535: }
536:
537: static void
538: Pdistance(arg,rp)
539: NODE arg;
540: Obj *rp;
541: {
542: Num s;
543: Num a, b;
544:
545: asir_assert(ARG0(arg),O_N,"distance");
546: asir_assert(ARG1(arg),O_N,"distance");
547: a = (Num)ARG0(arg);
548: b = (Num)ARG1(arg);
1.2 kondoh 549: if ( a && NID(a) == N_IntervalDouble && b && NID(b) == N_IntervalDouble ) {
550: distanceitvd((IntervalDouble)a, (IntervalDouble)b, (Num *)rp);
1.1 saito 551: } else {
552: distanceitvp((Itv)ARG0(arg),(Itv)ARG1(arg),&s);
553: *rp = (Obj)s;
554: }
555: }
556:
557: static void
558: Pinitv(arg,rp)
559: NODE arg;
560: Obj *rp;
561: {
562: int s;
563: Q q;
564:
565: asir_assert(ARG0(arg),O_N,"intval");
566: asir_assert(ARG1(arg),O_N,"intval");
567: if ( ! ARG1(arg) ) {
568: if ( ! ARG0(arg) ) s = 1;
569: else s = 0;
570: }
1.2 kondoh 571: else if ( NID(ARG1(arg)) == N_IntervalDouble ) {
572: s = initvd((Num)ARG0(arg),(IntervalDouble)ARG1(arg));
1.1 saito 573:
1.2 kondoh 574: } else if ( NID(ARG1(arg)) == N_IP || NID(ARG1(arg)) == N_IntervalBigFloat ) {
1.1 saito 575: if ( ! ARG0(arg) ) s = initvp((Num)ARG0(arg),(Itv)ARG1(arg));
576: else if ( NID(ARG0(arg)) == N_IP ) {
577: s = itvinitvp((Itv)ARG0(arg),(Itv)ARG1(arg));
578: } else {
579: s = initvp((Num)ARG0(arg),(Itv)ARG1(arg));
580: }
581: } else {
582: s = ! compnum(0,(Num)ARG0(arg),(Num)ARG1(arg));
583: }
584: STOQ(s,q);
585: *rp = (Obj)q;
586: }
587:
588: static void
589: Pdisjitv(arg,rp)
590: NODE arg;
591: Obj *rp;
592: {
593: Itv s;
594:
595: asir_assert(ARG0(arg),O_N,"disjitv");
596: asir_assert(ARG1(arg),O_N,"disjitv");
597: error("disjitv: not implemented yet");
598: if ( ! s ) *rp = 0;
599: else *rp = (Obj)ONE;
600: }
601:
602: #endif
603: extern int printmode;
604:
605: static void pprintmode( void )
606: {
607: switch (printmode) {
608: #if defined(INTERVAL)
609: case MID_PRINTF_E:
610: fprintf(stderr,"Interval printing mode is a mitpoint type.\n");
611: #endif
612: case PRINTF_E:
613: fprintf(stderr,"Printf's double printing mode is \"%%.16e\".\n");
614: break;
615: #if defined(INTERVAL)
616: case MID_PRINTF_G:
617: fprintf(stderr,"Interval printing mode is a mitpoint type.\n");
618: #endif
619: default:
620: case PRINTF_G:
621: fprintf(stderr,"Printf's double printing mode is \"%%g\".\n");
622: break;
623: }
624: }
625:
626: static void
627: Pprintmode(NODE arg, Obj *rp)
628: {
629: int l;
630: Q a, r;
631:
632: a = (Q)ARG0(arg);
633: if ( !a || NUM(a) && INT(a) ) {
634: l = QTOS(a);
635: if ( l < 0 ) l = 0;
636: #if defined(INTERVAL)
637: else if ( l > MID_PRINTF_E ) l = 0;
638: #else
639: else if ( l > PRINTF_E ) l = 0;
640: #endif
641: STOQ(printmode,r);
642: *rp = (Obj)r;
643: printmode = l;
644: pprintmode();
645: } else {
646: *rp = 0;
647: }
648: }
649:
650:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>