Annotation of OpenXM_contrib2/asir2000/engine/f-itv.c, Revision 1.10
1.1 saito 1: /*
1.10 ! kondoh 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/f-itv.c,v 1.9 2019/06/04 07:11:22 kondoh Exp $
1.1 saito 3: */
4: #if defined(INTERVAL)
5: #include "ca.h"
6: #include "base.h"
1.9 kondoh 7: #if 0
1.3 ohara 8: #if defined(PARI)
1.1 saito 9: #include "genpari.h"
10: #include "itv-pari.h"
1.7 ohara 11: long get_pariprec();
1.1 saito 12: #endif
13:
14: void ToBf(Num a, BF *rp)
15: {
1.8 noro 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: }
1.1 saito 62: }
63:
64: void double2bf(double d, BF *rp)
65: {
1.8 noro 66: GEN p;
1.1 saito 67:
1.8 noro 68: p = dbltor(d);
69: patori(p, rp);
70: cgiv(p);
1.1 saito 71: }
72:
1.8 noro 73: double bf2double(BF a)
1.1 saito 74: {
1.8 noro 75: GEN p;
1.1 saito 76:
1.8 noro 77: ritopa(a, &p);
78: return rtodbl(p);
1.1 saito 79: }
80:
81: void getulp(BF a, BF *au)
82: {
1.8 noro 83: GEN b, c;
84: long lb, i;
1.1 saito 85:
1.8 noro 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);
1.1 saito 96: }
97:
1.2 kondoh 98: void addulp(IntervalBigFloat a, IntervalBigFloat *c)
1.1 saito 99: {
1.8 noro 100: Num ai, as, aiu, asu, inf, sup;
1.1 saito 101:
1.8 noro 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);
1.1 saito 112: }
1.9 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 saito 122:
1.10 ! 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.9 kondoh 231:
232: void additvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *rp)
1.1 saito 233: {
1.9 kondoh 234: Num ai,as,bi,bs;
235: IntervalBigFloat c;
236: //,mas,mbs;
237: //,tmp;
1.8 noro 238: Num inf,sup;
1.9 kondoh 239: //GEN pa, pb, z;
240: //long ltop, lbot;
241: int current_roundmode;
1.8 noro 242:
243: if ( !a )
1.9 kondoh 244: *rp = b;
1.8 noro 245: else if ( !b )
1.9 kondoh 246: *rp = a;
1.8 noro 247: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
1.9 kondoh 248: addnum(0,(Num)a,(Num)b,(Num *)rp);
1.8 noro 249: else {
1.9 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 saito 271: #if 0
272: printexpr(CO, ai);
273: printexpr(CO, as);
274: printexpr(CO, bi);
275: printexpr(CO, bs);
276: #endif
1.8 noro 277: return;
278: }
1.1 saito 279: }
280:
1.9 kondoh 281: void subitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *rp)
1.1 saito 282: {
1.9 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.8 noro 291:
292: if ( !a )
1.9 kondoh 293: chsgnnum((Num)b,(Num *)rp);
1.8 noro 294: else if ( !b )
1.9 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.8 noro 298: else {
1.9 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.8 noro 324: }
1.1 saito 325: }
326:
1.2 kondoh 327: void mulitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 328: {
1.8 noro 329: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t,tmp;
330: Num inf, sup;
331: int ba,bb;
1.9 kondoh 332: int current_roundmode;
1.8 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.9 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.8 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.9 kondoh 397: istoitv(c1,t,(Itv *)c);
1.8 noro 398: } else {
399: subnum(0,0,c1,&t);
1.9 kondoh 400: istoitv(c2,t,(Itv *)c);
1.8 noro 401: }
1.9 kondoh 402: //addulp((IntervalBigFloat)tmp, c);
1.8 noro 403: }
1.9 kondoh 404: mpfr_roundmode = current_roundmode;
1.1 saito 405: }
406:
407: int initvf(Num a, Itv b)
408: {
1.8 noro 409: Num inf, sup;
1.1 saito 410:
1.8 noro 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;
1.1 saito 415: }
416:
417: int itvinitvf(Itv a, Itv b)
418: {
1.8 noro 419: Num ai,as,bi,bs;
1.1 saito 420:
1.8 noro 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;
1.1 saito 425: }
426:
1.2 kondoh 427: void divitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 428: {
1.8 noro 429: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t,tmp;
430: Num inf, sup;
431: int ba,bb;
1.9 kondoh 432: int current_roundmode;
1.8 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.9 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 saito 458: /* MUST check if ai, as, bi, bs are bigfloat. */
1.8 noro 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.9 kondoh 486: istoitv(c1,t,(Itv *)c);
1.8 noro 487: } else {
488: subnum(0,0,c1,&t);
1.9 kondoh 489: istoitv(c2,t,(Itv *)c);
1.8 noro 490: }
1.9 kondoh 491: //addulp((IntervalBigFloat)tmp, c);
1.8 noro 492: }
1.9 kondoh 493: mpfr_roundmode = current_roundmode;
1.1 saito 494: }
495:
496: void pwritvf(Itv a, Num e, Itv *c)
497: {
1.8 noro 498: int ei;
499: Itv t;
1.1 saito 500:
1.8 noro 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) ) {
1.3 ohara 508: #if defined(PARI) && 0
1.8 noro 509: gpui_ri((Obj)a,(Obj)c,(Obj *)c);
1.1 saito 510: #else
1.9 kondoh 511: error("pwritvf() : can't calculate a fractional power");
1.1 saito 512: #endif
1.8 noro 513: } else {
514: ei = QTOS((Q)e);
1.9 kondoh 515: if (ei<0) ei = -ei;
1.8 noro 516: pwritv0f(a,ei,&t);
517: if ( SGN((Q)e) < 0 )
1.9 kondoh 518: divitvf((IntervalBigFloat)ONE,(IntervalBigFloat)t,(IntervalBigFloat *)c); /* bug ?? */
1.8 noro 519: else
520: *c = t;
521: }
1.1 saito 522: }
523:
524: void pwritv0f(Itv a, int e, Itv *c)
525: {
1.8 noro 526: Num inf, sup;
527: Num ai,Xmin,Xmax;
528: IntervalBigFloat tmp;
529: Q ne;
1.9 kondoh 530: int current_roundmode;
1.8 noro 531:
532: if ( e == 1 )
533: *c = a;
534: else {
535: STOQ(e,ne);
536: if ( !(e%2) ) {
537: if ( initvp(0,a) ) {
538: Xmin = 0;
539: chsgnnum(INF(a),&ai);
540: if ( compnum(0,ai,SUP(a)) > 0 ) {
541: Xmax = ai;
542: } else {
543: Xmax = SUP(a);
544: }
545: } else {
546: if ( compnum(0,INF(a),0) > 0 ) {
547: Xmin = INF(a);
548: Xmax = SUP(a);
549: } else {
550: Xmin = SUP(a);
551: Xmax = INF(a);
552: }
553: }
554: } else {
555: Xmin = INF(a);
556: Xmax = SUP(a);
557: }
1.9 kondoh 558:
559: current_roundmode = mpfr_roundmode;
1.8 noro 560: if ( ! Xmin ) inf = 0;
1.9 kondoh 561: else {
562: mpfr_roundmode = MPFR_RNDD;
563: pwrbf(Xmin,(Num)ne,&inf);
564: }
1.8 noro 565: if ( ! Xmax ) sup = 0;
1.9 kondoh 566: else {
567: mpfr_roundmode = MPFR_RNDU;
568: pwrbf(Xmax,(Num)ne,&sup);
569: }
1.8 noro 570: istoitv(inf,sup,(Itv *)&tmp);
1.9 kondoh 571: mpfr_roundmode = current_roundmode;
572: *c = (Itv)tmp;
573: // addulp(tmp, (IntervalBigFloat *)c);
1.8 noro 574: }
1.1 saito 575: }
576:
1.9 kondoh 577: void chsgnitvf(IntervalBigFloat a, IntervalBigFloat *rp)
1.1 saito 578: {
1.8 noro 579: Num inf,sup;
1.9 kondoh 580: IntervalBigFloat c;
581: int current_roundmode;
1.1 saito 582:
1.8 noro 583: if ( !a )
1.9 kondoh 584: *rp = 0;
585: else if ( NID(a) < N_IntervalBigFloat )
586: chsgnnum((Num)a,(Num *)rp);
1.8 noro 587: else {
1.9 kondoh 588: current_roundmode = mpfr_roundmode;
589:
590: mpfr_roundmode = MPFR_RNDU;
591: chsgnnum((Num)INF(a),&sup);
592: mpfr_roundmode = MPFR_RNDD;
593: chsgnnum((Num)SUP(a),&inf);
594: //MKIntervalBigFloat((BF)inf,(BF)sup,c);
595: istoitv(inf,sup,(Itv *)&c);
596: *rp = c;
597: mpfr_roundmode = current_roundmode;
1.8 noro 598: }
1.1 saito 599: }
600:
601: int cmpitvf(Itv a, Itv b)
602: {
1.8 noro 603: Num ai,as,bi,bs;
604: int s,t;
1.1 saito 605:
1.8 noro 606: if ( !a ) {
607: if ( !b || (NID(b)<=N_B) ) {
608: return compnum(0,(Num)a,(Num)b);
609: } else {
610: itvtois(b,&bi,&bs);
611: if ( compnum(0,(Num)a,bs) > 0 ) return 1;
612: else if ( compnum(0,bi,(Num)a) > 0 ) return -1;
613: else return 0;
614: }
615: } else if ( !b ) {
616: if ( !a || (NID(a)<=N_B) ) {
617: return compnum(0,(Num)a,(Num)b);
618: } else {
619: itvtois(a,&ai,&as);
620: if ( compnum(0,ai,(Num)b) > 0 ) return 1;
621: else if ( compnum(0,(Num)b,as) > 0 ) return -1;
622: else return 0;
623: }
624: } else {
625: itvtois(a,&ai,&as);
626: itvtois(b,&bi,&bs);
627: s = compnum(0,ai,bs) ;
628: t = compnum(0,bi,as) ;
629: if ( s > 0 ) return 1;
630: else if ( t > 0 ) return -1;
631: else return 0;
632: }
1.1 saito 633: }
634:
635: void miditvf(Itv a, Num *b)
636: {
1.8 noro 637: Num ai,as;
638: Num t;
1.1 saito 639:
1.8 noro 640: if ( ! a ) *b = 0;
641: else if ( (NID(a) <= N_B) )
642: *b = (Num)a;
643: else {
644: STOQ(2,TWO);
645: itvtois(a,&ai,&as);
646: addbf(ai,as,&t);
647: divbf(t,(Num)TWO,b);
648: }
1.1 saito 649: }
650:
651: void cupitvf(Itv a, Itv b, Itv *c)
652: {
1.8 noro 653: Num ai,as,bi,bs;
654: Num inf,sup;
1.1 saito 655:
1.8 noro 656: itvtois(a,&ai,&as);
657: itvtois(b,&bi,&bs);
658: if ( compnum(0,ai,bi) > 0 ) inf = bi;
659: else inf = ai;
660: if ( compnum(0,as,bs) > 0 ) sup = as;
661: else sup = bs;
662: istoitv(inf,sup,c);
1.1 saito 663: }
664:
665: void capitvf(Itv a, Itv b, Itv *c)
666: {
1.8 noro 667: Num ai,as,bi,bs;
668: Num inf,sup;
1.1 saito 669:
1.8 noro 670: itvtois(a,&ai,&as);
671: itvtois(b,&bi,&bs);
672: if ( compnum(0,ai,bi) > 0 ) inf = ai;
673: else inf = bi;
674: if ( compnum(0,as,bs) > 0 ) sup = bs;
675: else sup = as;
676: if ( compnum(0,inf,sup) > 0 ) *c = 0;
677: else istoitv(inf,sup,c);
1.1 saito 678: }
679:
680:
681: void widthitvf(Itv a, Num *b)
682: {
1.8 noro 683: Num ai,as;
1.1 saito 684:
1.8 noro 685: if ( ! a ) *b = 0;
686: else if ( (NID(a) <= N_B) )
687: *b = (Num)a;
688: else {
689: itvtois(a,&ai,&as);
690: subnum(0,as,ai,b);
691: }
1.1 saito 692: }
693:
694: void absitvf(Itv a, Num *b)
695: {
1.8 noro 696: Num ai,as,bi,bs;
1.1 saito 697:
1.8 noro 698: if ( ! a ) *b = 0;
699: else if ( (NID(a) <= N_B) ) {
700: if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
701: else *b = (Num)a;
702: } else {
703: itvtois(a,&ai,&as);
704: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
705: else bi = ai;
706: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
707: else bs = as;
708: if ( compnum(0,bi,bs) > 0 ) *b = bi;
709: else *b = bs;
710: }
1.1 saito 711: }
712:
713: void distanceitvf(Itv a, Itv b, Num *c)
714: {
1.8 noro 715: Num ai,as,bi,bs;
716: Num di,ds;
717: Itv d;
718:
719: itvtois(a,&ai,&as);
720: itvtois(b,&bi,&bs);
721: subnum(0,ai,bi,&di);
722: subnum(0,as,bs,&ds);
723: istoitv(di,ds,&d);
724: absitvf(d,c);
1.1 saito 725: }
726:
727: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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 / asir2000 / engine