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