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