![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to f-itv.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / engine|
Return to f-itv.c CVS log | Up to [local] / OpenXM_contrib2 / asir2018 / engine |
1.1 noro 1: /*
1.5 ! kondoh 2: * $OpenXM: OpenXM_contrib2/asir2018/engine/f-itv.c,v 1.4 2019/10/17 03:03:12 kondoh Exp $
1.1 noro 3: */
4: #if defined(INTERVAL)
5: #include "ca.h"
6: #include "base.h"
1.3 kondoh 7: #if 0
1.1 noro 8: #if defined(PARI)
9: #include "genpari.h"
10: #include "itv-pari.h"
11: long get_pariprec();
12: #endif
13:
14: void ToBf(Num a, BF *rp)
15: {
16: GEN pa, pb, pc;
17: BF bn, bd;
18: BF c;
19: long ltop, lbot;
20:
21: if ( ! a ) {
22: *rp = 0;
23: }
24: else switch ( NID(a) ) {
25: case N_Q:
26: ltop = avma;
27: ritopa_i(NM((Q)a), SGN((Q)a), &pa);
28: pb = cgetr(get_pariprec());
29: mpaff(pa, pb);
30: if ( INT((Q)a) ) {
31: lbot = avma;
32: pb = gerepile(ltop, lbot, pb);
33: patori(pb, &c);
34: cgiv(pb);
35: *rp = c;
36: } else {
37: patori(pb, &bn);
38: ritopa_i(DN((Q)a), 1, &pa);
39: pb = cgetr(get_pariprec());
40: mpaff(pa, pb);
41: lbot = avma;
42: pb = gerepile(ltop, lbot, pb);
43: patori(pb, &bd);
44: cgiv(pb);
45: divbf((Num)bn,(Num)bd, (Num *)&c);
46: *rp = c;
47: }
48: break;
49: case N_R:
50: pb = dbltor(BDY((Real)a));
51: patori(pb, &c);
52: cgiv(pb);
53: *rp = c;
54: break;
55: case N_B:
56: *rp = (BF)a;
57: break;
58: default:
59: error("ToBf : not implemented");
60: break;
61: }
62: }
63:
64: void double2bf(double d, BF *rp)
65: {
66: GEN p;
67:
68: p = dbltor(d);
69: patori(p, rp);
70: cgiv(p);
71: }
72:
73: double bf2double(BF a)
74: {
75: GEN p;
76:
77: ritopa(a, &p);
78: return rtodbl(p);
79: }
80:
81: void getulp(BF a, BF *au)
82: {
83: GEN b, c;
84: long lb, i;
85:
86: ritopa(a, &b);
87: lb = lg(b);
88: c = cgetr(lb);
89: c[1] = b[1];
90: for (i=2;i<lb-1;i++) c[i] = 0;
91: c[i] = 1;
92: setsigne(c, 1);
93: patori(c,au);
94: cgiv(c);
95: cgiv(b);
96: }
97:
98: void addulp(IntervalBigFloat a, IntervalBigFloat *c)
99: {
100: Num ai, as, aiu, asu, inf, sup;
101:
102: itvtois((Itv)a,&ai,&as);
103: if ( ai ) {
104: getulp((BF)ai, (BF *)&aiu);
105: subbf(ai,aiu,&inf);
106: } else inf = ai;
107: if ( as ) {
108: getulp((BF)as, (BF *)&asu);
109: addbf(as,asu,&sup);
110: } else sup = as;
111: istoitv(inf,sup, (Itv *)c);
112: }
1.3 kondoh 113: #endif
114:
115: // in builtin/bfaux.c
116: extern int mpfr_roundmode;
117:
118: // in engine/bf.c
119: Num tobf(Num,int);
120:
121: #define BFPREC(a) (((BF)(a))->body->_mpfr_prec)
1.1 noro 122:
1.5 ! kondoh 123: double mpfr2dblDown(mpfr_t a)
! 124: {
! 125: return mpfr_get_d(a,MPFR_RNDD);
! 126: }
! 127:
! 128: double mpfr2dblUp(mpfr_t a)
! 129: {
! 130: return mpfr_get_d(a,MPFR_RNDU);
! 131: }
! 132:
! 133:
! 134: void toInterval(Num a, int prec, int type, Num *rp)
! 135: {
! 136: if ( ! a ) {
! 137: *rp = 0;
! 138: } else if (type == EvalIntervalDouble) {
! 139: if (NID(a)==N_C) {
! 140: double inf,sup;
! 141: C z;
! 142: IntervalDouble re, im;
! 143:
! 144: if ( ! ((C)a)->r ) {
! 145: re = 0;
! 146: } else {
! 147: inf = toRealDown(((C)a)->r);
! 148: sup = toRealUp(((C)a)->r);
! 149: MKIntervalDouble(inf,sup,re);
! 150: }
! 151: if ( ! ((C)a)->i ) {
! 152: im = 0;
! 153: } else {
! 154: inf = toRealDown(((C)a)->i);
! 155: sup = toRealUp(((C)a)->i);
! 156: MKIntervalDouble(inf,sup,im);
! 157: }
! 158: if ( !re && !im )
! 159: z = 0;
! 160: else
! 161: reimtocplx((Num)re,(Num)im,(Num *)&z);
! 162: *rp = (Num)z;
! 163: } else {
! 164: double inf,sup;
! 165: IntervalDouble c;
! 166:
! 167: inf = toRealDown(a);
! 168: sup = toRealUp(a);
! 169:
! 170: MKIntervalDouble(inf,sup,c);
! 171: *rp = (Num) c;
! 172: }
! 173: } else if (type == EvalIntervalBigFloat) {
! 174: if (NID(a)==N_C) {
! 175: Num ai,as;
! 176: Num inf,sup;
! 177: C z;
! 178: IntervalBigFloat re, im;
! 179: int current_roundmode;
! 180:
! 181: current_roundmode = mpfr_roundmode;
! 182:
! 183: if ( ! ((C)a)->r )
! 184: re = 0;
! 185: else {
! 186: itvtois((Itv)((C)a)->r,&ai,&as);
! 187: mpfr_roundmode = MPFR_RNDD;
! 188: inf = tobf(ai, prec);
! 189: mpfr_roundmode = MPFR_RNDU;
! 190: sup = tobf(as, prec);
! 191: istoitv(inf,sup,(Itv *)&re);
! 192: }
! 193:
! 194: if ( ! ((C)a)->i )
! 195: im = 0;
! 196: else {
! 197: itvtois((Itv)((C)a)->i,&ai,&as);
! 198: mpfr_roundmode = MPFR_RNDD;
! 199: inf = tobf(ai, prec);
! 200: mpfr_roundmode = MPFR_RNDU;
! 201: sup = tobf(as, prec);
! 202: istoitv(inf,sup,(Itv *)&im);
! 203: }
! 204:
! 205: mpfr_roundmode = current_roundmode;
! 206: reimtocplx((Num)re,(Num)im,(Num *)&z);
! 207: *rp = (Num)z;
! 208: } else {
! 209: Num ai,as;
! 210: Num inf,sup;
! 211: IntervalBigFloat c;
! 212: int current_roundmode;
! 213:
! 214: itvtois((Itv)a,&ai,&as);
! 215:
! 216: current_roundmode = mpfr_roundmode;
! 217: mpfr_roundmode = MPFR_RNDD;
! 218: inf = tobf(ai, prec);
! 219: mpfr_roundmode = MPFR_RNDU;
! 220: sup = tobf(as, prec);
! 221: istoitv(inf,sup,(Itv *)&c);
! 222: mpfr_roundmode = current_roundmode;
! 223: *rp = (Num) c;
! 224: }
! 225: } else {
! 226: error("toInterval: not supported types.");
! 227: *rp = 0;
! 228: }
! 229: }
! 230:
1.3 kondoh 231:
232: void additvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *rp)
1.1 noro 233: {
1.3 kondoh 234: Num ai,as,bi,bs;
235: IntervalBigFloat c;
236: //,mas,mbs;
237: //,tmp;
1.1 noro 238: Num inf,sup;
1.3 kondoh 239: //GEN pa, pb, z;
240: //long ltop, lbot;
241: int current_roundmode;
1.1 noro 242:
243: if ( !a )
1.3 kondoh 244: *rp = b;
1.1 noro 245: else if ( !b )
1.3 kondoh 246: *rp = a;
1.1 noro 247: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
1.3 kondoh 248: addnum(0,(Num)a,(Num)b,(Num *)rp);
1.1 noro 249: else {
1.3 kondoh 250: itvtois((Itv)a,&ai,&as);
251: itvtois((Itv)b,&bi,&bs);
252:
253: current_roundmode = mpfr_roundmode;
254:
255: mpfr_roundmode = MPFR_RNDD;
256: ai = tobf(ai, DEFAULTPREC);
257: bi = tobf(bi, DEFAULTPREC);
258: //addnum(0,ai,bi,&inf);
259: addbf(ai,bi,&inf);
260:
261: mpfr_roundmode = MPFR_RNDU;
262: as = tobf(as, DEFAULTPREC);
263: bs = tobf(bs, DEFAULTPREC);
264: //addnum(0,as,bs,&sup);
265: addbf(as,bs,&sup);
266:
267: istoitv(inf,sup,(Itv *)&c);
268: mpfr_roundmode = current_roundmode;
269: //MKIntervalBigFloat((BF)inf, (BF)sup, c);
270: *rp = c;
1.1 noro 271: #if 0
272: printexpr(CO, ai);
273: printexpr(CO, as);
274: printexpr(CO, bi);
275: printexpr(CO, bs);
276: #endif
277: return;
278: }
279: }
280:
1.3 kondoh 281: void subitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *rp)
1.1 noro 282: {
1.3 kondoh 283: Num ai,as,bi,bs;
284: IntervalBigFloat c;
285: //,mas, mbs;
286: Num inf,sup;
287: //,tmp;
288: //GEN pa, pb, z;
289: //long ltop, lbot;
290: int current_roundmode;
1.1 noro 291:
292: if ( !a )
1.3 kondoh 293: chsgnnum((Num)b,(Num *)rp);
1.1 noro 294: else if ( !b )
1.3 kondoh 295: *rp = a;
296: else if ( (NID(a) < N_IntervalBigFloat) && (NID(b) < N_IntervalBigFloat ) )
297: subnum(0,(Num)a,(Num)b,(Num *)rp);
1.1 noro 298: else {
1.3 kondoh 299: itvtois((Itv)a,&ai,&as);
300: itvtois((Itv)b,&bi,&bs);
301:
302: current_roundmode = mpfr_roundmode;
303:
304: mpfr_roundmode = MPFR_RNDD;
305: ai = tobf(ai, DEFAULTPREC);
306: bi = tobf(bi, DEFAULTPREC);
307:
308: mpfr_roundmode = MPFR_RNDU;
309: as = tobf(as, DEFAULTPREC);
310: bs = tobf(bs, DEFAULTPREC);
311: //subnum(0,as,bi,&sup);
312: subbf(as,bi,&sup);
313:
314: mpfr_roundmode = MPFR_RNDD;
315: //subnum(0,ai,bs,&inf);
316: subbf(ai,bs,&inf);
317:
318: istoitv(inf,sup,(Itv *)&c);
319: mpfr_roundmode = current_roundmode;
320: //MKIntervalBigFloat((BF)inf, (BF)sup, c);
321: *rp = c;
322:
323: //addulp((IntervalBigFloat)tmp, c);
1.1 noro 324: }
325: }
326:
327: void mulitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
328: {
329: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t,tmp;
330: Num inf, sup;
331: int ba,bb;
1.3 kondoh 332: int current_roundmode;
1.1 noro 333:
334: if ( !a || !b )
335: *c = 0;
336: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
337: mulnum(0,(Num)a,(Num)b,(Num *)c);
338: else {
1.3 kondoh 339: itvtois((Itv)a,&ai,&as);
340: itvtois((Itv)b,&bi,&bs);
341:
342: current_roundmode = mpfr_roundmode;
343:
344: mpfr_roundmode = MPFR_RNDU;
345: as = tobf(as, DEFAULTPREC);
346: bs = tobf(bs, DEFAULTPREC);
347:
348: mpfr_roundmode = MPFR_RNDD;
349: ai = tobf(ai, DEFAULTPREC);
350: bi = tobf(bi, DEFAULTPREC);
351:
352: // itvtois((Itv)a,&inf,&sup);
353: // ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
354: // itvtois((Itv)b,&inf,&sup);
355: // ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
1.1 noro 356:
357: /* MUST check if ai, as, bi, bs are bigfloat. */
358: chsgnnum(as,&a2);
359: chsgnnum(bs,&b2);
360:
361:
362: if ( ba = (compnum(0,0,a2) > 0) ) {
363: a1 = ai;
364: } else {
365: a1 = a2;
366: a2 = ai;
367: }
368: if ( bb = (compnum(0,0,b2) > 0) ) {
369: b1 = bi;
370: } else {
371: b1 = b2;
372: b2 = bi;
373: }
374: mulnum(0,a2,b2,&t);
375: subnum(0,0,t,&c2);
376: if ( compnum(0,0,b1) > 0 ) {
377: mulnum(0,a2,b1,&t);
378: subnum(0,0,t,&c1);
379: if ( compnum(0,0,a1) > 0 ) {
380: mulnum(0,a1,b2,&t);
381: subnum(0,0,t,&p);
382: if ( compnum(0,c1,p) > 0 ) c1 = p;
383: mulnum(0,a1,b1,&t);
384: subnum(0,0,t,&p);
385: if ( compnum(0,c2,p) > 0 ) c2 = p;
386: }
387: } else {
388: if ( compnum(0,0,a1) > 0 ) {
389: mulnum(0,a1,b2,&t);
390: subnum(0,0,t,&c1);
391: } else {
392: mulnum(0,a1,b1,&c1);
393: }
394: }
395: if ( ba == bb ) {
396: subnum(0,0,c2,&t);
1.3 kondoh 397: istoitv(c1,t,(Itv *)c);
1.1 noro 398: } else {
399: subnum(0,0,c1,&t);
1.3 kondoh 400: istoitv(c2,t,(Itv *)c);
1.1 noro 401: }
1.3 kondoh 402: //addulp((IntervalBigFloat)tmp, c);
1.1 noro 403: }
1.3 kondoh 404: mpfr_roundmode = current_roundmode;
1.1 noro 405: }
406:
407: int initvf(Num a, Itv b)
408: {
409: Num inf, sup;
410:
411: itvtois(b, &inf, &sup);
412: if ( compnum(0,inf,a) > 0 ) return 0;
413: else if ( compnum(0,a,sup) > 0 ) return 0;
414: else return 1;
415: }
416:
417: int itvinitvf(Itv a, Itv b)
418: {
419: Num ai,as,bi,bs;
420:
421: itvtois(a, &ai, &as);
422: itvtois(b, &bi, &bs);
423: if ( compnum(0,ai,bi) >= 0 && compnum(0,bs,as) >= 0 ) return 1;
424: else return 0;
425: }
426:
427: void divitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
428: {
429: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t,tmp;
430: Num inf, sup;
431: int ba,bb;
1.3 kondoh 432: int current_roundmode;
1.1 noro 433:
434: if ( !b )
435: error("divitv : division by 0");
436: else if ( !a )
437: *c = 0;
438: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
439: divnum(0,(Num)a,(Num)b,(Num *)c);
440: else {
1.3 kondoh 441: itvtois((Itv)a,&ai,&as);
442: itvtois((Itv)b,&bi,&bs);
443:
444: current_roundmode = mpfr_roundmode;
445:
446: mpfr_roundmode = MPFR_RNDU;
447: as = tobf(as, DEFAULTPREC);
448: bs = tobf(bs, DEFAULTPREC);
449:
450: mpfr_roundmode = MPFR_RNDD;
451: ai = tobf(ai, DEFAULTPREC);
452: bi = tobf(bi, DEFAULTPREC);
453:
454: // itvtois((Itv)a,&inf,&sup);
455: // ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
456: // itvtois((Itv)b,&inf,&sup);
457: // ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
1.1 noro 458: /* MUST check if ai, as, bi, bs are bigfloat. */
459: chsgnnum(as,&a2);
460: chsgnnum(bs,&b2);
461: if ( ba = (compnum(0,0,a2) > 0) ) {
462: a1 = ai;
463: } else {
464: a1 = a2;
465: a2 = ai;
466: }
467: if ( bb = (compnum(0,0,b2) > 0) ) {
468: b1 = bi;
469: } else {
470: b1 = b2;
471: b2 = bi;
472: }
473: if ( compnum(0,0,b1) >= 0 ) {
474: error("divitv : division by interval including 0.");
475: } else {
476: divnum(0,a2,b1,&c2);
477: if ( compnum(0,0,a1) > 0 ) {
478: divnum(0,a1,b1,&c1);
479: } else {
480: divnum(0,a1,b2,&t);
481: subnum(0,0,t,&c1);
482: }
483: }
484: if ( ba == bb ) {
485: subnum(0,0,c2,&t);
1.3 kondoh 486: istoitv(c1,t,(Itv *)c);
1.1 noro 487: } else {
488: subnum(0,0,c1,&t);
1.3 kondoh 489: istoitv(c2,t,(Itv *)c);
1.1 noro 490: }
1.3 kondoh 491: //addulp((IntervalBigFloat)tmp, c);
1.1 noro 492: }
1.3 kondoh 493: mpfr_roundmode = current_roundmode;
1.1 noro 494: }
495:
496: void pwritvf(Itv a, Num e, Itv *c)
497: {
1.4 kondoh 498: long ei;
1.1 noro 499: Itv t;
500:
501: if ( !e )
502: *c = (Itv)ONE;
503: else if ( !a )
504: *c = 0;
505: else if ( NID(a) <= N_B )
506: pwrnum(0,(Num)a,e,(Num *)c);
507: else if ( !INT(e) ) {
508: #if defined(PARI) && 0
509: gpui_ri((Obj)a,(Obj)c,(Obj *)c);
510: #else
1.3 kondoh 511: error("pwritvf() : can't calculate a fractional power");
1.1 noro 512: #endif
513: } else {
1.4 kondoh 514: //ei = QTOS((Q)e);
515: ei = mpz_get_si(BDY((Q)e));
1.3 kondoh 516: if (ei<0) ei = -ei;
1.1 noro 517: pwritv0f(a,ei,&t);
1.4 kondoh 518: // if ( SGN((Q)e) < 0 )
519: if ( sgnq((Q)e) < 0 )
1.3 kondoh 520: divitvf((IntervalBigFloat)ONE,(IntervalBigFloat)t,(IntervalBigFloat *)c); /* bug ?? */
1.1 noro 521: else
522: *c = t;
523: }
524: }
525:
1.4 kondoh 526: void pwritv0f(Itv a, long e, Itv *c)
1.1 noro 527: {
528: Num inf, sup;
529: Num ai,Xmin,Xmax;
530: IntervalBigFloat tmp;
531: Q ne;
1.3 kondoh 532: int current_roundmode;
1.1 noro 533:
534: if ( e == 1 )
535: *c = a;
536: else {
1.4 kondoh 537: STOZ(e,ne);
1.1 noro 538: if ( !(e%2) ) {
539: if ( initvp(0,a) ) {
540: Xmin = 0;
541: chsgnnum(INF(a),&ai);
542: if ( compnum(0,ai,SUP(a)) > 0 ) {
543: Xmax = ai;
544: } else {
545: Xmax = SUP(a);
546: }
547: } else {
548: if ( compnum(0,INF(a),0) > 0 ) {
549: Xmin = INF(a);
550: Xmax = SUP(a);
551: } else {
552: Xmin = SUP(a);
553: Xmax = INF(a);
554: }
555: }
556: } else {
557: Xmin = INF(a);
558: Xmax = SUP(a);
559: }
1.3 kondoh 560:
561: current_roundmode = mpfr_roundmode;
1.1 noro 562: if ( ! Xmin ) inf = 0;
1.3 kondoh 563: else {
564: mpfr_roundmode = MPFR_RNDD;
565: pwrbf(Xmin,(Num)ne,&inf);
566: }
1.1 noro 567: if ( ! Xmax ) sup = 0;
1.3 kondoh 568: else {
569: mpfr_roundmode = MPFR_RNDU;
570: pwrbf(Xmax,(Num)ne,&sup);
571: }
1.1 noro 572: istoitv(inf,sup,(Itv *)&tmp);
1.3 kondoh 573: mpfr_roundmode = current_roundmode;
574: *c = (Itv)tmp;
575: // addulp(tmp, (IntervalBigFloat *)c);
1.1 noro 576: }
577: }
578:
1.3 kondoh 579: void chsgnitvf(IntervalBigFloat a, IntervalBigFloat *rp)
1.1 noro 580: {
581: Num inf,sup;
1.3 kondoh 582: IntervalBigFloat c;
583: int current_roundmode;
1.1 noro 584:
585: if ( !a )
1.3 kondoh 586: *rp = 0;
587: else if ( NID(a) < N_IntervalBigFloat )
588: chsgnnum((Num)a,(Num *)rp);
1.1 noro 589: else {
1.3 kondoh 590: current_roundmode = mpfr_roundmode;
591:
592: mpfr_roundmode = MPFR_RNDU;
593: chsgnnum((Num)INF(a),&sup);
594: mpfr_roundmode = MPFR_RNDD;
595: chsgnnum((Num)SUP(a),&inf);
596: //MKIntervalBigFloat((BF)inf,(BF)sup,c);
597: istoitv(inf,sup,(Itv *)&c);
598: *rp = c;
599: mpfr_roundmode = current_roundmode;
1.1 noro 600: }
601: }
602:
603: int cmpitvf(Itv a, Itv b)
604: {
605: Num ai,as,bi,bs;
606: int s,t;
607:
608: if ( !a ) {
609: if ( !b || (NID(b)<=N_B) ) {
610: return compnum(0,(Num)a,(Num)b);
611: } else {
612: itvtois(b,&bi,&bs);
613: if ( compnum(0,(Num)a,bs) > 0 ) return 1;
614: else if ( compnum(0,bi,(Num)a) > 0 ) return -1;
615: else return 0;
616: }
617: } else if ( !b ) {
618: if ( !a || (NID(a)<=N_B) ) {
619: return compnum(0,(Num)a,(Num)b);
620: } else {
621: itvtois(a,&ai,&as);
622: if ( compnum(0,ai,(Num)b) > 0 ) return 1;
623: else if ( compnum(0,(Num)b,as) > 0 ) return -1;
624: else return 0;
625: }
626: } else {
627: itvtois(a,&ai,&as);
628: itvtois(b,&bi,&bs);
629: s = compnum(0,ai,bs) ;
630: t = compnum(0,bi,as) ;
631: if ( s > 0 ) return 1;
632: else if ( t > 0 ) return -1;
633: else return 0;
634: }
635: }
636:
637: void miditvf(Itv a, Num *b)
638: {
639: Num ai,as;
640: Num t;
641:
642: if ( ! a ) *b = 0;
643: else if ( (NID(a) <= N_B) )
644: *b = (Num)a;
645: else {
1.4 kondoh 646: //STOQ(2,TWO);
1.1 noro 647: itvtois(a,&ai,&as);
648: addbf(ai,as,&t);
649: divbf(t,(Num)TWO,b);
650: }
651: }
652:
653: void cupitvf(Itv a, Itv b, Itv *c)
654: {
655: Num ai,as,bi,bs;
656: Num inf,sup;
657:
658: itvtois(a,&ai,&as);
659: itvtois(b,&bi,&bs);
660: if ( compnum(0,ai,bi) > 0 ) inf = bi;
661: else inf = ai;
662: if ( compnum(0,as,bs) > 0 ) sup = as;
663: else sup = bs;
664: istoitv(inf,sup,c);
665: }
666:
667: void capitvf(Itv a, Itv b, Itv *c)
668: {
669: Num ai,as,bi,bs;
670: Num inf,sup;
671:
672: itvtois(a,&ai,&as);
673: itvtois(b,&bi,&bs);
674: if ( compnum(0,ai,bi) > 0 ) inf = ai;
675: else inf = bi;
676: if ( compnum(0,as,bs) > 0 ) sup = bs;
677: else sup = as;
678: if ( compnum(0,inf,sup) > 0 ) *c = 0;
679: else istoitv(inf,sup,c);
680: }
681:
682:
683: void widthitvf(Itv a, Num *b)
684: {
685: Num ai,as;
686:
687: if ( ! a ) *b = 0;
688: else if ( (NID(a) <= N_B) )
689: *b = (Num)a;
690: else {
691: itvtois(a,&ai,&as);
692: subnum(0,as,ai,b);
693: }
694: }
695:
696: void absitvf(Itv a, Num *b)
697: {
698: Num ai,as,bi,bs;
699:
700: if ( ! a ) *b = 0;
701: else if ( (NID(a) <= N_B) ) {
702: if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
703: else *b = (Num)a;
704: } else {
705: itvtois(a,&ai,&as);
706: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
707: else bi = ai;
708: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
709: else bs = as;
710: if ( compnum(0,bi,bs) > 0 ) *b = bi;
711: else *b = bs;
712: }
713: }
714:
715: void distanceitvf(Itv a, Itv b, Num *c)
716: {
717: Num ai,as,bi,bs;
718: Num di,ds;
719: Itv d;
720:
721: itvtois(a,&ai,&as);
722: itvtois(b,&bi,&bs);
723: subnum(0,ai,bi,&di);
724: subnum(0,as,bs,&ds);
725: istoitv(di,ds,&d);
726: absitvf(d,c);
727: }
728:
729: #endif