Annotation of OpenXM_contrib2/asir2000/engine/d-itv.c, Revision 1.10
1.1 saito 1: /*
1.10 ! kondoh 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/d-itv.c,v 1.9 2019/06/04 07:11:22 kondoh Exp $
1.1 saito 3: */
4: #if defined(INTERVAL)
5: #include <float.h>
6: #include "ca.h"
7: #include "base.h"
1.9 kondoh 8: #if 0
1.3 ohara 9: #if defined(PARI)
1.1 saito 10: #include "genpari.h"
11: #endif
1.9 kondoh 12: #endif
1.1 saito 13:
14: #if defined(ITVDEBUG)
15: void printbinint(int d)
16: {
1.8 noro 17: int i, j, mask;
18: union {
19: int x;
20: char c[sizeof(int)];
21: } a;
22:
23: a.x = d;
24: for(i=sizeof(int)-1;i>=0;i--) {
25: mask = 0x80;
26: for(j=0;j<8;j++) {
27: if (a.c[i] & mask) fprintf(stderr,"1");
28: else fprintf(stderr,"0");
29: mask >>= 1;
30: }
31: }
32: fprintf(stderr,"\n");
1.1 saito 33: }
34: #endif
35:
36: double NatToRealUp(N a, int *expo)
37: {
1.8 noro 38: int *p;
39: int l,all,i,j,s,tb,top,tail;
40: unsigned int t,m[2];
41: N b,c;
42: int kk, carry, rem;
43:
44: p = BD(a); l = PL(a) - 1;
45: for ( top = 0, t = p[l]; t; t >>= 1, top++ );
46: all = top + BSH*l; tail = (53-top)%BSH; j = l-(53-top)/BSH-1;
1.1 saito 47:
48:
1.8 noro 49: if ( j >= 0 ) {
1.1 saito 50:
51: #if defined(ITVDEBUG)
1.8 noro 52: fprintf(stderr," %d : tail = %d\n", j, tail);
53: printbinint(p[j]);
1.1 saito 54: #endif
1.8 noro 55: kk = (1<< (BSH - tail)) - 1;
56: rem = p[j] & kk;
57: if ( !rem )
58: for (i=0;i<j;i++)
59: if ( p[j] ) {
60: carry = 1;
61: break;
62: }
63: else carry = 1;
64: if ( carry ) {
65: b = NALLOC(j+1+1);
66: PL(b) = j+1;
67: p = BD(b);
68: for(i=0;i<j;i++) p[i] = 0;
69: p[j]=(1<< (BSH - tail));
70:
71: addn(a,b,&c);
72:
73: p = BD(c); l = PL(c) - 1;
74: for ( top = 0, t = p[l]; t; t >>= 1, top++ );
75: all = top + BSH*l; tail = (53-top)%BSH; i = l-(53-top)/BSH-1;
76: } else i = j;
77: } else i = j;
78:
79:
80: m[1] = i < 0 ? 0 : p[i]>>(BSH-tail);
81: for ( j = 1, i++, tb = tail; i <= l; i++ ) {
82: s = 32-tb; t = i < 0 ? 0 : p[i];
83: if ( BSH > s ) {
84: m[j] |= ((t&((1<<s)-1))<<tb);
85: if ( !j )
86: break;
87: else {
88: j--; m[j] = t>>s; tb = BSH-s;
89: }
90: } else {
91: m[j] |= (t<<tb); tb += BSH;
92: }
93: }
94: s = (all-1)+1023;
95: m[0] = (m[0]&((1<<20)-1))|(MIN(2046,s)<<20); *expo = MAX(s-2046,0);
1.1 saito 96: #ifdef vax
1.8 noro 97: t = m[0]; m[0] = m[1]; m[1] = t; itod(m);
1.1 saito 98: #endif
1.7 ohara 99: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID)
1.8 noro 100: t = m[0]; m[0] = m[1]; m[1] = t;
1.1 saito 101: #endif
1.8 noro 102: return *((double *)m);
1.1 saito 103: }
104:
1.8 noro 105: static double Q2doubleDown(Q a)
1.1 saito 106: {
1.8 noro 107: double nm,dn,t,snm,rd;
108: int enm,edn,e;
109: unsigned int *p,s;
1.1 saito 110:
1.8 noro 111: nm = (SGN(a)>0)? NatToReal(NM(a),&enm): NatToRealUp(NM(a),&enm);
1.1 saito 112:
1.8 noro 113: if ( INT(a) )
114: if ( enm >= 1 )
115: error("Q2doubleDown : Overflow");
116: else
117: return SGN(a)>0 ? nm : -nm;
118: else {
119: dn = (SGN(a)>0)? NatToReal(DN(a),&edn): NatToRealUp(DN(a),&edn);
1.1 saito 120:
1.8 noro 121: if ( ((e = enm - edn) >= 1024) || (e <= -1023) )
122: error("Q2doubleDown : Overflow"); /* XXX */
123: else {
124: t = 0.0; p = (unsigned int *)&t; /* Machine */
125: *p = ((e+1023)<<20);
1.1 saito 126: #ifdef vax
1.8 noro 127: s = p[0]; p[0] = p[1]; p[1] = s; itod(p);
1.1 saito 128: #endif
1.7 ohara 129: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) | defined(ANDROID)
1.8 noro 130: s = p[0]; p[0] = p[1]; p[1] = s;
1.1 saito 131: #endif
1.8 noro 132: FPMINUSINF
133: snm = (SGN(a)>0) ? nm : -nm;
134: rd = snm / dn * t;
135: FPNEAREST
136: return rd;
137: }
138: }
1.1 saito 139: }
140:
141:
1.8 noro 142: static double Q2doubleUp(Q a)
143: {
144: double nm,dn,t,snm,rd;
145: int enm,edn,e;
146: unsigned int *p,s;
147:
148: nm = (SGN(a)>0)? NatToRealUp(NM(a),&enm): NatToReal(NM(a),&enm);
149:
150: if ( INT(a) )
151: if ( enm >= 1 )
152: error("Q2doubleUp : Overflow");
153: else
154: return SGN(a)>0 ? nm : -nm;
155: else {
156: dn = (SGN(a)>0)? NatToRealUp(DN(a),&edn): NatToReal(DN(a),&edn);
157:
158: if ( ((e = enm - edn) >= 1024) || (e <= -1023) )
159: error("Q2doubleUp : Overflow"); /* XXX */
160: else {
161: t = 0.0; p = (unsigned int *)&t; /* Machine */
162: *p = ((e+1023)<<20);
1.1 saito 163: #ifdef vax
1.8 noro 164: s = p[0]; p[0] = p[1]; p[1] = s; itod(p);
1.1 saito 165: #endif
1.7 ohara 166: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID)
1.8 noro 167: s = p[0]; p[0] = p[1]; p[1] = s;
1.1 saito 168: #endif
169: #if 0
1.8 noro 170: FPPLUSINF
171: snm = (SGN(a)>0) ? nm : -nm;
172: rd = snm / dn * t;
1.1 saito 173: #endif
174:
1.8 noro 175: FPMINUSINF
176: snm = (SGN(a)>0) ? -nm : nm;
177: rd = snm / dn * t;
178: FPNEAREST
179: return (-rd);
180: }
181: }
1.1 saito 182: }
183:
1.9 kondoh 184: #if 0
1.8 noro 185: static double PARI2doubleDown(BF a)
1.1 saito 186: {
1.9 kondoh 187: //GEN p;
188: Num p;
1.8 noro 189: double d;
1.1 saito 190:
191: ritopa(a, &p);
192: d = rtodbl(p);
1.8 noro 193: cgiv(p);
194: return d;
1.1 saito 195: }
196:
1.8 noro 197: static double PARI2doubleUp(BF a)
1.1 saito 198: {
1.8 noro 199: return PARI2doubleDown(a);
1.1 saito 200: }
1.9 kondoh 201: #endif
1.1 saito 202:
1.8 noro 203: double subulpd(double d)
204: {
205: double inf;
206:
207: FPMINUSINF
208: inf = d - DBL_MIN;
209: FPNEAREST
210: return inf;
211: }
212:
213: double addulpd(double d)
214: {
215: double md, sup;
216:
217: md = -d;
218: FPMINUSINF
219: sup = md - DBL_MIN;
220: FPNEAREST
221: return (-sup);
222: }
223:
1.10 ! kondoh 224: double toRealDown(Num a)
1.8 noro 225: {
226: double inf;
227:
228: if ( ! a ) return 0.0;
229: switch ( NID(a) ) {
230: case N_Q:
231: inf = Q2doubleDown((Q)a); break;
232: case N_R:
233: inf = subulpd(BDY((Real)a)); break;
234: case N_B:
1.9 kondoh 235: //inf = PARI2doubleDown((BF)a); break;
1.10 ! kondoh 236: inf = mpfr_get_d(BDY((BF)a),MPFR_RNDD);;
! 237: //error("toRealDown: not supported operands.");
! 238: break;
1.8 noro 239: case N_IP:
1.10 ! kondoh 240: inf = toRealDown(INF((Itv)a));
1.8 noro 241: break;
242: case N_IntervalDouble:
243: inf = INF((IntervalDouble)a); break;
1.10 ! kondoh 244: case N_IntervalBigFloat:
! 245: inf = mpfr_get_d(BDY((BF)INF((IntervalBigFloat)a)),MPFR_RNDD);;
! 246: break;
1.8 noro 247: case N_A:
248: default:
249: inf = 0.0;
1.10 ! kondoh 250: error("toRealDown: not supported operands.");
1.8 noro 251: break;
252: }
253: return inf;
254: }
255:
1.10 ! kondoh 256: double toRealUp(Num a)
1.8 noro 257: {
258: double sup;
259:
260: if ( ! a ) return 0.0;
261: switch ( NID(a) ) {
262: case N_Q:
263: sup = Q2doubleUp((Q)a); break;
264: case N_R:
265: sup = addulpd(BDY((Real)a)); break;
266: case N_B:
1.9 kondoh 267: //sup = PARI2doubleUp((BF)a); break;
1.10 ! kondoh 268: sup = mpfr_get_d(BDY((BF)a),MPFR_RNDU);;
! 269: //error("toRealUp: not supported operands.");
1.9 kondoh 270: break;
1.8 noro 271: case N_IP:
1.10 ! kondoh 272: sup = toRealUp(SUP((Itv)a)); break;
1.8 noro 273: case N_IntervalDouble:
274: sup = SUP((IntervalDouble)a); break;
1.10 ! kondoh 275: case N_IntervalBigFloat:
! 276: sup = mpfr_get_d(BDY((BF)SUP((IntervalBigFloat)a)),MPFR_RNDU);;
! 277: break;
1.8 noro 278: case N_A:
279: default:
280: sup = 0.0;
1.10 ! kondoh 281: error("toRealUp: not supported operands.");
1.8 noro 282: break;
283: }
284: return sup;
285: }
286:
287:
288: void Num2double(Num a, double *inf, double *sup)
289: {
1.10 ! kondoh 290: *inf = 0.0;
! 291: *sup = 0.0;
! 292: if ( a && NUM(a) ) {
! 293: switch ( NID(a) ) {
! 294: case N_Q:
! 295: *inf = Q2doubleDown((Q)a);
! 296: *sup = Q2doubleUp((Q)a);
! 297: break;
! 298: case N_R:
! 299: *inf = BDY((Real)a);
! 300: *sup = BDY((Real)a);
! 301: break;
! 302: case N_B:
! 303: //*inf = PARI2doubleDown((BF)a);
! 304: //*sup = PARI2doubleUp((BF)a);
! 305: *inf = mpfr_get_d(BDY((BF)a), MPFR_RNDD);
! 306: *sup = mpfr_get_d(BDY((BF)a), MPFR_RNDU);
! 307: break;
! 308: case N_IP:
! 309: *inf = toRealDown(INF((Itv)a));
! 310: *sup = toRealUp(SUP((Itv)a));
! 311: break;
! 312: case N_IntervalDouble:
! 313: *inf = INF((IntervalDouble)a);
! 314: *sup = SUP((IntervalDouble)a);
! 315: break;
! 316: case N_A:
! 317: default:
! 318: *inf = 0.0;
! 319: *sup = 0.0;
! 320: error("Num2double: not supported operands.");
! 321: break;
! 322: }
1.8 noro 323: }
1.1 saito 324: }
325:
326:
1.2 kondoh 327: void additvd(Num a, Num b, IntervalDouble *c)
1.1 saito 328: {
1.8 noro 329: double ai,as,mas, bi,bs;
330: double inf,sup;
1.1 saito 331:
1.8 noro 332: if ( !a ) {
333: *c = (IntervalDouble)b;
334: } else if ( !b ) {
335: *c = (IntervalDouble)a;
336: #if 0
337: } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
338: addnum(0,a,b,c);
339: #endif
340: } else {
341:
342: Num2double((Num)a,&ai,&as);
343: Num2double((Num)b,&bi,&bs);
344: mas = -as;
345: FPMINUSINF
346: inf = ai + bi;
347: sup = mas - bs; /* as + bs = ( - ( (-as) - bs ) ) */
348: FPNEAREST
349: MKIntervalDouble(inf,(-sup),*c);
350: }
1.1 saito 351: }
352:
1.2 kondoh 353: void subitvd(Num a, Num b, IntervalDouble *c)
1.1 saito 354: {
1.8 noro 355: double ai,as,mas, bi,bs;
356: double inf,sup;
1.1 saito 357:
1.8 noro 358: if ( !a ) {
359: *c = (IntervalDouble)b;
360: } else if ( !b ) {
361: *c = (IntervalDouble)a;
362: #if 0
363: } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
364: addnum(0,a,b,c);
365: #endif
366: } else {
367: Num2double(a,&ai,&as);
368: Num2double(b,&bi,&bs);
369: FPMINUSINF
370: inf = ai - bs;
371: sup = ( bi - as ); /* as - bi = ( - ( bi - as ) ) */
372: FPNEAREST
373: MKIntervalDouble(inf,(-sup),*c);
374: }
375: }
376:
377: void mulitvd(Num a, Num b, IntervalDouble *c)
378: {
379: double ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p;
380: double inf, sup;
381: int ba,bb;
382:
383: if ( !a || !b ) {
384: *c = 0;
385: #if 0
386: } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
387: mulnum(0,a,b,c);
388: #endif
389: } else {
390: Num2double(a,&ai,&as);
391: Num2double(b,&bi,&bs);
392:
393: FPMINUSINF
394:
395: a2 = -as;
396: b2 = -bs;
397:
398: if ( ba = ( a2 < 0.0 ) ) {
399: a1 = ai;
400: } else {
401: a1 = a2;
402: a2 = ai;
403: }
404: if ( bb = ( b2 < 0.0 ) ) {
405: b1 = bi;
406: } else {
407: b1 = b2;
408: b2 = bi;
409: }
410:
411: c2 = - a2 * b2;
412: if ( b1 < 0.0 ) {
413: c1 = - a2 * b1;
414: if ( a1 < 0.0 ) {
415: p = - a1 * b2;
416: if ( p < c1 ) c1 = p;
417: p = - a1 * b1;
418: if ( p < c2 ) c2 = p;
419: }
420: } else {
421: c1 = (a1<0.0)?( - a1 * b2 ):( a1 * b1 );
422: }
423: if ( ba == bb ) {
424: inf = c1;
425: sup = - c2;
426: } else {
427: inf = c2;
428: sup = - c1;
429: }
430: FPNEAREST
431: MKIntervalDouble(inf,sup,*c);
432: }
1.1 saito 433: }
434:
1.2 kondoh 435: int initvd(Num a, IntervalDouble b)
1.1 saito 436: {
1.8 noro 437: Real inf, sup;
1.1 saito 438:
1.8 noro 439: if ( NID(b) == N_IntervalDouble ) {
440: MKReal(INF(b), inf);
441: MKReal(SUP(b), sup);
442: } else return 0;
443: if ( compnum(0,(Num)inf,a) > 0 ) return 0;
444: else if ( compnum(0,a,(Num)sup) > 0 ) return 0;
445: else return 1;
446: }
447:
448: void divitvd(Num a, Num b, IntervalDouble *c)
449: {
450: double ai,as,bi,bs,a1,a2,b1,b2,c1,c2;
451: double inf, sup;
452: int ba,bb;
453:
454: if ( !b ) {
455: *c = 0;
456: error("divitvd : division by 0");
457: } else if ( !a ) {
458: *c = 0;
459: #if 0
460: } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
461: divnum(0,a,b,c);
462: #endif
463: } else {
464: Num2double(a,&ai,&as);
465: Num2double(b,&bi,&bs);
466:
467: FPMINUSINF
468:
469: a2 = -as;
470: b2 = -bs;
471:
472: if ( ba = ( a2 < 0.0 ) ) {
473: a1 = ai;
474: } else {
475: a1 = a2;
476: a2 = ai;
477: }
478: if ( bb = ( b2 < 0.0 ) ) {
479: b1 = bi;
480: } else {
481: b1 = b2;
482: b2 = bi;
483: }
484:
485: if ( b1 <= 0.0 ) {
486: *c = 0;
487: error("divitvd : division by 0 interval");
488: } else {
489: c2 = a2 / b1;
490: c1 = (a1<0.0)?( a1 / b1 ):( - a1 / b2 );
491: }
492: if ( ba == bb ) {
493: inf = c1;
494: sup = - c2;
495: } else {
496: inf = c2;
497: sup = - c1;
498: }
499: FPNEAREST
500: MKIntervalDouble(inf,sup,*c);
501: }
502: }
503:
504: void pwritvd(Num a, Num e, IntervalDouble *c)
505: {
506: int ei;
507: IntervalDouble t;
508:
509: if ( !e ) {
510: *c = (IntervalDouble)ONE;
511: } else if ( !a ) {
512: *c = 0;
513: #if 0
514: } else if ( NID(a) <= N_IP ) {
515: pwrnum(0,a,e,c);
1.1 saito 516: #endif
1.8 noro 517: } else if ( !INT(e) ) {
1.3 ohara 518: #if defined(PARI) && 0
1.8 noro 519: gpui_ri((Obj)a,(Obj)c,(Obj *)c);
1.1 saito 520: #else
1.8 noro 521: error("pwritvd : can't calculate a fractional power");
1.1 saito 522: #endif
1.8 noro 523: } else {
524: ei = QTOS((Q)e);
1.9 kondoh 525: if (ei<0) ei = -ei;
1.8 noro 526: pwritv0d((IntervalDouble)a,ei,&t);
527: if ( SGN((Q)e) < 0 )
528: divnum(0,(Num)ONE,(Num)t,(Num *)c);
529: else
530: *c = t;
531: }
532: }
533:
534: void pwritv0d(IntervalDouble a, int e, IntervalDouble *c)
535: {
536: double inf, sup;
537: double t, Xmin, Xmax;
538: int i;
539:
540: if ( e == 1 )
541: *c = a;
542: else {
543: if ( !(e%2) ) {
544: if ( initvd(0,a) ) {
545: Xmin = 0.0;
546: t = - INF(a);
547: Xmax = SUP(a);
548: if ( t > Xmax ) Xmax = t;
549: } else {
550: if ( INF(a) > 0.0 ) {
551: Xmin = INF(a);
552: Xmax = SUP(a);
553: } else {
554: Xmin = - SUP(a);
555: Xmax = - INF(a);
556: }
557: }
558: } else {
559: Xmin = INF(a);
560: Xmax = SUP(a);
561: }
562: FPPLUSINF
563: sup = (Xmax!=0.0)?pwrreal0(Xmax,e):0.0;
564: FPMINUSINF
565: inf = (Xmin!=0.0)?pwrreal0(Xmin,e):0.0;
566: FPNEAREST
567: MKIntervalDouble(inf,sup,*c);
568: }
569: }
570:
571: void chsgnitvd(IntervalDouble a, IntervalDouble *c)
572: {
573: double inf,sup;
574:
575: if ( !a )
576: *c = 0;
577: #if 0
578: else if ( NID(a) <= N_IP )
579: chsgnnum(a,c);
580: #endif
581: else {
582: inf = - SUP(a);
583: sup = - INF(a);
584: MKIntervalDouble(inf,sup,*c);
585: }
1.1 saito 586: }
587:
1.8 noro 588: int cmpitvd(IntervalDouble a, IntervalDouble b)
1.2 kondoh 589: /* a > b: 1 */
590: /* a = b: 0 */
591: /* a < b: -1 */
1.1 saito 592: {
1.8 noro 593: double ai,as,bi,bs;
1.1 saito 594:
1.8 noro 595: if ( !a ) {
596: if ( !b || (NID(b)<=N_IP) ) {
597: return compnum(0,(Num)a,(Num)b);
598: } else {
599: bi = INF(b);
600: bs = SUP(b);
601: if ( bi > 0.0 ) return -1;
602: else if ( bs < 0.0 ) return 1;
603: else return 0;
604: }
605: } else if ( !b ) {
606: if ( !a || (NID(a)<=N_IP) ) {
607: return compnum(0,(Num)a,(Num)b);
608: } else {
609: ai = INF(a);
610: as = SUP(a);
611: if ( ai > 0.0 ) return 1;
612: else if ( as < 0.0 ) return -1;
613: else return 0;
614: }
615: } else {
616: bi = INF(b);
617: bs = SUP(b);
618: ai = INF(a);
619: as = SUP(a);
620: if ( ai > bs ) return 1;
621: else if ( bi > as ) return -1;
622: else return 0;
623: }
624: }
625:
626: void miditvd(IntervalDouble a, Num *b)
627: {
628: double t;
629: Real rp;
630:
631: if ( ! a ) *b = 0;
632: else if ( (NID(a) <= N_IP) )
633: *b = (Num)a;
634: else {
635: t = (INF(a) + SUP(a))/2.0;
636: MKReal(t,rp);
637: *b = (Num)rp;
638: }
639: }
640:
641: void cupitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
642: {
643: double ai,as,bi,bs;
644: double inf,sup;
645:
646: ai = INF(a);
647: as = SUP(a);
648: bi = INF(b);
649: bs = SUP(b);
650: inf = MIN(ai,bi);
651: sup = MAX(as,bs);
652: MKIntervalDouble(inf,sup,*c);
653: }
654:
655: void capitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
656: {
657: double ai,as,bi,bs;
658: double inf,sup;
659:
660: ai = INF(a);
661: as = SUP(a);
662: bi = INF(b);
663: bs = SUP(b);
664: inf = MAX(ai,bi);
665: sup = MIN(as,bs);
666: if ( inf > sup ) *c = 0;
667: else MKIntervalDouble(inf,sup,*c);
668: }
669:
670:
671: void widthitvd(IntervalDouble a, Num *b)
672: {
673: double t;
674: Real rp;
675:
676: if ( ! a ) *b = 0;
677: else {
678: t = SUP(a) - INF(a);
679: MKReal(t,rp);
680: *b = (Num)rp;
681: }
682: }
683:
684: void absitvd(IntervalDouble a, Num *b)
685: {
686: double ai,as,bi,bs;
687: double t;
688: Real rp;
689:
690: if ( ! a ) *b = 0;
691: else {
692: ai = INF(a);
693: as = SUP(a);
694: if (ai < 0) bi = -ai; else bi = ai;
695: if (as < 0) bs = -as; else bs = as;
696: if ( bi > bs ) t = bi; else t = bs;
697: MKReal(t,rp);
698: *b = (Num)rp;
699: }
700: }
701:
1.10 ! kondoh 702: void absintvald(IntervalDouble a, IntervalDouble *b)
! 703: {
! 704: double ai,as,inf,sup;
! 705:
! 706: ai = INF(a);
! 707: as = SUP(a);
! 708: if ( as < 0 ) {
! 709: sup = -ai;
! 710: inf = -as;
! 711: } else if (ai < 0) {
! 712: inf = 0.0;
! 713: sup = MAX(as, -ai);
! 714: } else {
! 715: inf = ai;
! 716: sup = as;
! 717: }
! 718: MKIntervalDouble(inf,sup,*b);
! 719: }
! 720:
1.8 noro 721: void distanceitvd(IntervalDouble a, IntervalDouble b, Num *c)
722: {
723: double ai,as,bi,bs;
724: double di,ds;
725: double t;
726: Real rp;
727:
728: ai = INF(a);
729: as = SUP(a);
730: bi = INF(b);
731: bs = SUP(b);
732: di = bi - ai;
733: ds = bs - as;
734:
735: if (di < 0) di = -di;
736: if (ds < 0) ds = -ds;
737: if (di > ds) t = di; else t = ds;
738: MKReal(t,rp);
739: *c = (Num)rp;
1.1 saito 740: }
741:
742: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to d-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