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