![[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|
Return to d-itv.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! saito 1: /*
! 2: * $OpenXM: $
! 3: */
! 4: #if defined(INTERVAL)
! 5: #include <float.h>
! 6: #include "ca.h"
! 7: #include "base.h"
! 8: #if PARI
! 9: #include "genpari.h"
! 10: #endif
! 11:
! 12: #if defined(ITVDEBUG)
! 13: void printbinint(int d)
! 14: {
! 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");
! 31: }
! 32: #endif
! 33:
! 34: double NatToRealUp(N a, int *expo)
! 35: {
! 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;
! 45:
! 46:
! 47: if ( j >= 0 ) {
! 48:
! 49: #if defined(ITVDEBUG)
! 50: fprintf(stderr," %d : tail = %d\n", j, tail);
! 51: printbinint(p[j]);
! 52: #endif
! 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);
! 94: #ifdef vax
! 95: t = m[0]; m[0] = m[1]; m[1] = t; itod(m);
! 96: #endif
! 97: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__)
! 98: t = m[0]; m[0] = m[1]; m[1] = t;
! 99: #endif
! 100: return *((double *)m);
! 101: }
! 102:
! 103: static double Q2doubleDown(Q a)
! 104: {
! 105: double nm,dn,t,snm,rd;
! 106: int enm,edn,e;
! 107: unsigned int *p,s;
! 108:
! 109: nm = (SGN(a)>0)? NatToReal(NM(a),&enm): NatToRealUp(NM(a),&enm);
! 110:
! 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);
! 118:
! 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);
! 124: #ifdef vax
! 125: s = p[0]; p[0] = p[1]; p[1] = s; itod(p);
! 126: #endif
! 127: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__)
! 128: s = p[0]; p[0] = p[1]; p[1] = s;
! 129: #endif
! 130: FPMINUSINF
! 131: snm = (SGN(a)>0) ? nm : -nm;
! 132: rd = snm / dn * t;
! 133: FPNEAREST
! 134: return rd;
! 135: }
! 136: }
! 137: }
! 138:
! 139:
! 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);
! 161: #ifdef vax
! 162: s = p[0]; p[0] = p[1]; p[1] = s; itod(p);
! 163: #endif
! 164: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__)
! 165: s = p[0]; p[0] = p[1]; p[1] = s;
! 166: #endif
! 167: #if 0
! 168: FPPLUSINF
! 169: snm = (SGN(a)>0) ? nm : -nm;
! 170: rd = snm / dn * t;
! 171: #endif
! 172:
! 173: FPMINUSINF
! 174: snm = (SGN(a)>0) ? -nm : nm;
! 175: rd = snm / dn * t;
! 176: FPNEAREST
! 177: return (-rd);
! 178: }
! 179: }
! 180: }
! 181:
! 182: static double PARI2doubleDown(BF a)
! 183: {
! 184: GEN p;
! 185: double d;
! 186:
! 187: ritopa(a, &p);
! 188: d = rtodbl(p);
! 189: cgiv(p);
! 190: return d;
! 191: }
! 192:
! 193: static double PARI2doubleUp(BF a)
! 194: {
! 195: return PARI2doubleDown(a);
! 196: }
! 197:
! 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_ID:
! 235: inf = INF((ItvD)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_ID:
! 260: sup = SUP((ItvD)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_ID:
! 291: *inf = INF((ItvD)a);
! 292: *sup = SUP((ItvD)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: }
! 301: }
! 302:
! 303:
! 304: void additvd(Num a, Num b, ItvD *c)
! 305: {
! 306: double ai,as,mas, bi,bs;
! 307: double inf,sup;
! 308:
! 309: if ( !a ) {
! 310: *c = (ItvD)b;
! 311: } else if ( !b ) {
! 312: *c = (ItvD)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: MKItvD(inf,(-sup),*c);
! 327: }
! 328: }
! 329:
! 330: void subitvd(Num a, Num b, ItvD *c)
! 331: {
! 332: double ai,as,mas, bi,bs;
! 333: double inf,sup;
! 334:
! 335: if ( !a ) {
! 336: *c = (ItvD)b;
! 337: } else if ( !b ) {
! 338: *c = (ItvD)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: MKItvD(inf,(-sup),*c);
! 351: }
! 352: }
! 353:
! 354: void mulitvd(Num a, Num b, ItvD *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: MKItvD(inf,sup,*c);
! 409: }
! 410: }
! 411:
! 412: int initvd(Num a, ItvD b)
! 413: {
! 414: Real inf, sup;
! 415:
! 416: if ( NID(b) == N_ID ) {
! 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, ItvD *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: MKItvD(inf,sup,*c);
! 478: }
! 479: }
! 480:
! 481: void pwritvd(Num a, Num e, ItvD *c)
! 482: {
! 483: int ei;
! 484: ItvD t;
! 485:
! 486: if ( !e ) {
! 487: *c = (ItvD)ONE;
! 488: } else if ( !a ) {
! 489: *c = 0;
! 490: #if 0
! 491: } else if ( NID(a) <= N_IP ) {
! 492: pwrnum(0,a,e,c);
! 493: #endif
! 494: } else if ( !INT(e) ) {
! 495: #if PARI && 0
! 496: GEN pa,pe,z;
! 497: int ltop,lbot;
! 498:
! 499: ltop = avma; ritopa(a,&pa); ritopa(e,&pe); lbot = avma;
! 500: z = gerepile(ltop,lbot,gpui(pa,pe,prec));
! 501: patori(z,c); cgiv(z);
! 502: #else
! 503: error("pwritvd : can't calculate a fractional power");
! 504: #endif
! 505: } else {
! 506: ei = QTOS((Q)e);
! 507: pwritv0d((ItvD)a,ei,&t);
! 508: if ( SGN((Q)e) < 0 )
! 509: divnum(0,(Num)ONE,(Num)t,(Num *)c);
! 510: else
! 511: *c = t;
! 512: }
! 513: }
! 514:
! 515: void pwritv0d(ItvD a, int e, ItvD *c)
! 516: {
! 517: double inf, sup;
! 518: double t, Xmin, Xmax;
! 519: int i;
! 520:
! 521: if ( e == 1 )
! 522: *c = a;
! 523: else {
! 524: if ( !(e%2) ) {
! 525: if ( initvd(0,a) ) {
! 526: Xmin = 0.0;
! 527: t = - INF(a);
! 528: Xmax = SUP(a);
! 529: if ( t > Xmax ) Xmax = t;
! 530: } else {
! 531: if ( INF(a) > 0.0 ) {
! 532: Xmin = INF(a);
! 533: Xmax = SUP(a);
! 534: } else {
! 535: Xmin = - SUP(a);
! 536: Xmax = - INF(a);
! 537: }
! 538: }
! 539: } else {
! 540: Xmin = INF(a);
! 541: Xmax = SUP(a);
! 542: }
! 543: FPPLUSINF
! 544: sup = (Xmax!=0.0)?pwrreal0(Xmax,e):0.0;
! 545: FPMINUSINF
! 546: inf = (Xmin!=0.0)?pwrreal0(Xmin,e):0.0;
! 547: FPNEAREST
! 548: MKItvD(inf,sup,*c);
! 549: }
! 550: }
! 551:
! 552: void chsgnitvd(ItvD a, ItvD *c)
! 553: {
! 554: double inf,sup;
! 555:
! 556: if ( !a )
! 557: *c = 0;
! 558: #if 0
! 559: else if ( NID(a) <= N_IP )
! 560: chsgnnum(a,c);
! 561: #endif
! 562: else {
! 563: inf = - SUP(a);
! 564: sup = - INF(a);
! 565: MKItvD(inf,sup,*c);
! 566: }
! 567: }
! 568:
! 569: int cmpitvd(ItvD a, ItvD b)
! 570: {
! 571: double ai,as,bi,bs;
! 572:
! 573: if ( !a ) {
! 574: if ( !b || (NID(b)<=N_IP) ) {
! 575: return compnum(0,(Num)a,(Num)b);
! 576: } else {
! 577: bi = INF(b);
! 578: bs = SUP(b);
! 579: if ( bi > 0.0 ) return -1;
! 580: else if ( bs < 0.0 ) return 1;
! 581: else return 0;
! 582: }
! 583: } else if ( !b ) {
! 584: if ( !a || (NID(a)<=N_IP) ) {
! 585: return compnum(0,(Num)a,(Num)b);
! 586: } else {
! 587: ai = INF(a);
! 588: as = SUP(a);
! 589: if ( ai > 0.0 ) return 1;
! 590: else if ( as < 0.0 ) return -1;
! 591: else return 0;
! 592: }
! 593: } else {
! 594: bi = INF(b);
! 595: bs = SUP(b);
! 596: ai = INF(a);
! 597: as = SUP(a);
! 598: if ( ai > bs ) return 1;
! 599: else if ( bi > as ) return -1;
! 600: else return 0;
! 601: }
! 602: }
! 603:
! 604: void miditvd(ItvD a, Num *b)
! 605: {
! 606: double t;
! 607: Real rp;
! 608:
! 609: if ( ! a ) *b = 0;
! 610: else if ( (NID(a) <= N_IP) )
! 611: *b = (Num)a;
! 612: else {
! 613: t = (INF(a) + SUP(a))/2.0;
! 614: MKReal(t,rp);
! 615: *b = (Num)rp;
! 616: }
! 617: }
! 618:
! 619: void cupitvd(ItvD a, ItvD b, ItvD *c)
! 620: {
! 621: double ai,as,bi,bs;
! 622: double inf,sup;
! 623:
! 624: ai = INF(a);
! 625: as = SUP(a);
! 626: bi = INF(b);
! 627: bs = SUP(b);
! 628: inf = MIN(ai,bi);
! 629: sup = MAX(as,bs);
! 630: MKItvD(inf,sup,*c);
! 631: }
! 632:
! 633: void capitvd(ItvD a, ItvD b, ItvD *c)
! 634: {
! 635: double ai,as,bi,bs;
! 636: double inf,sup;
! 637:
! 638: ai = INF(a);
! 639: as = SUP(a);
! 640: bi = INF(b);
! 641: bs = SUP(b);
! 642: inf = MAX(ai,bi);
! 643: sup = MIN(as,bs);
! 644: if ( inf > sup ) *c = 0;
! 645: else MKItvD(inf,sup,*c);
! 646: }
! 647:
! 648:
! 649: void widthitvd(ItvD a, Num *b)
! 650: {
! 651: double t;
! 652: Real rp;
! 653:
! 654: if ( ! a ) *b = 0;
! 655: else {
! 656: t = SUP(a) - INF(a);
! 657: MKReal(t,rp);
! 658: *b = (Num)rp;
! 659: }
! 660: }
! 661:
! 662: void absitvd(ItvD a, Num *b)
! 663: {
! 664: double ai,as,bi,bs;
! 665: double t;
! 666: Real rp;
! 667:
! 668: if ( ! a ) *b = 0;
! 669: else {
! 670: ai = INF(a);
! 671: as = SUP(a);
! 672: if (ai < 0) bi = -ai; else bi = ai;
! 673: if (as < 0) bs = -as; else bs = as;
! 674: if ( bi > bs ) t = bi; else t = bs;
! 675: MKReal(t,rp);
! 676: *b = (Num)rp;
! 677: }
! 678: }
! 679:
! 680: void distanceitvd(ItvD a, ItvD b, Num *c)
! 681: {
! 682: double ai,as,bi,bs;
! 683: double di,ds;
! 684: double t;
! 685: Real rp;
! 686:
! 687: ai = INF(a);
! 688: as = SUP(a);
! 689: bi = INF(b);
! 690: bs = SUP(b);
! 691: di = bi - ai;
! 692: ds = bs - as;
! 693:
! 694: if (di < 0) di = -di;
! 695: if (ds < 0) ds = -ds;
! 696: if (di > ds) t = di; else t = ds;
! 697: MKReal(t,rp);
! 698: *c = (Num)rp;
! 699: }
! 700:
! 701: #endif