![[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|
Return to f-itv.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! saito 1: /*
! 2: * $OpenXM: $
! 3: */
! 4: #if defined(INTERVAL)
! 5: #include "ca.h"
! 6: #include "base.h"
! 7: #if PARI
! 8: #include "genpari.h"
! 9: #include "itv-pari.h"
! 10: extern long prec;
! 11: #endif
! 12:
! 13: void ToBf(Num a, BF *rp)
! 14: {
! 15: GEN pa, pb, pc;
! 16: BF bn, bd;
! 17: BF c;
! 18: long ltop, lbot;
! 19:
! 20: if ( ! a ) {
! 21: *rp = 0;
! 22: }
! 23: else switch ( NID(a) ) {
! 24: case N_Q:
! 25: ltop = avma;
! 26: ritopa_i(NM((Q)a), SGN((Q)a), &pa);
! 27: pb = cgetr(prec);
! 28: mpaff(pa, pb);
! 29: if ( INT((Q)a) ) {
! 30: lbot = avma;
! 31: pb = gerepile(ltop, lbot, pb);
! 32: patori(pb, &c);
! 33: cgiv(pb);
! 34: *rp = c;
! 35: } else {
! 36: patori(pb, &bn);
! 37: ritopa_i(DN((Q)a), 1, &pa);
! 38: pb = cgetr(prec);
! 39: mpaff(pa, pb);
! 40: lbot = avma;
! 41: pb = gerepile(ltop, lbot, pb);
! 42: patori(pb, &bd);
! 43: cgiv(pb);
! 44: divbf((Num)bn,(Num)bd, (Num *)&c);
! 45: *rp = c;
! 46: }
! 47: break;
! 48: case N_R:
! 49: pb = dbltor(BDY((Real)a));
! 50: patori(pb, &c);
! 51: cgiv(pb);
! 52: *rp = c;
! 53: break;
! 54: case N_B:
! 55: *rp = (BF)a;
! 56: break;
! 57: default:
! 58: error("ToBf : not implemented");
! 59: break;
! 60: }
! 61: }
! 62:
! 63: void double2bf(double d, BF *rp)
! 64: {
! 65: GEN p;
! 66:
! 67: p = dbltor(d);
! 68: patori(p, rp);
! 69: cgiv(p);
! 70: }
! 71:
! 72: double bf2double(BF a)
! 73: {
! 74: GEN p;
! 75:
! 76: ritopa(a, &p);
! 77: return rtodbl(p);
! 78: }
! 79:
! 80: void getulp(BF a, BF *au)
! 81: {
! 82: GEN b, c;
! 83: long lb, i;
! 84:
! 85: ritopa(a, &b);
! 86: lb = lg(b);
! 87: c = cgetr(lb);
! 88: c[1] = b[1];
! 89: for (i=2;i<lb-1;i++) c[i] = 0;
! 90: c[i] = 1;
! 91: setsigne(c, 1);
! 92: patori(c,au);
! 93: cgiv(c);
! 94: cgiv(b);
! 95: }
! 96:
! 97: void addulp(ItvF a, ItvF *c)
! 98: {
! 99: Num ai, as, aiu, asu, inf, sup;
! 100:
! 101: itvtois((Itv)a,&ai,&as);
! 102: if ( ai ) {
! 103: getulp((BF)ai, (BF *)&aiu);
! 104: subbf(ai,aiu,&inf);
! 105: } else inf = ai;
! 106: if ( as ) {
! 107: getulp((BF)as, (BF *)&asu);
! 108: addbf(as,asu,&sup);
! 109: } else sup = as;
! 110: istoitv(inf,sup, (Itv *)c);
! 111: }
! 112:
! 113: void additvf(ItvF a, ItvF b, ItvF *c)
! 114: {
! 115: Num ai,as,bi,bs,mas,mbs,tmp;
! 116: Num inf,sup;
! 117: GEN pa, pb, z;
! 118: long ltop, lbot;
! 119:
! 120: if ( !a )
! 121: *c = b;
! 122: else if ( !b )
! 123: *c = a;
! 124: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
! 125: addnum(0,(Num)a,(Num)b,(Num *)c);
! 126: else {
! 127: itvtois((Itv)a,&inf,&sup);
! 128: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
! 129: itvtois((Itv)b,&inf,&sup);
! 130: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
! 131: #if 0
! 132: printexpr(CO, ai);
! 133: printexpr(CO, as);
! 134: printexpr(CO, bi);
! 135: printexpr(CO, bs);
! 136: #endif
! 137:
! 138: #if 1
! 139: addnum(0,ai,bi,&inf);
! 140: addnum(0,as,bs,&sup);
! 141: istoitv(inf,sup,(Itv *)&tmp);
! 142: addulp((ItvF)tmp, c);
! 143: return;
! 144: #else
! 145: ltop = avma;
! 146: ritopa(ai,&pa);
! 147: ritopa(bi,&pb);
! 148: lbot = avma;
! 149: z = gerepile(ltop,lbot,PariAddDown(pa,pb));
! 150: patori(z,&inf); cgiv(z);
! 151:
! 152: /* MUST check if ai, as, bi, bs are bigfloat. */
! 153:
! 154: /* as + bs = ( - ( (-as) + (-bs) ) ) */
! 155: chsgnbf(as,&mas);
! 156: chsgnbf(bs,&mbs);
! 157: ltop = avma;
! 158: ritopa(mas,&pa);
! 159: ritopa(mbs,&pb);
! 160: lbot = avma;
! 161: z = gerepile(ltop,lbot,PariAddDown(pa,pb));
! 162: patori(z,&tmp); cgiv(z);
! 163:
! 164: chsgnbf(tmp,&sup);
! 165: istoitv(inf,sup,c);
! 166: #endif
! 167: }
! 168: }
! 169:
! 170: void subitvf(ItvF a, ItvF b, ItvF *c)
! 171: {
! 172: Num ai,as,bi,bs,mas, mbs;
! 173: Num inf,sup,tmp;
! 174: GEN pa, pb, z;
! 175: long ltop, lbot;
! 176:
! 177: if ( !a )
! 178: chsgnnum((Num)b,(Num *)c);
! 179: else if ( !b )
! 180: *c = a;
! 181: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
! 182: subnum(0,(Num)a,(Num)b,(Num *)c);
! 183: else {
! 184: itvtois((Itv)a,&inf,&sup);
! 185: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
! 186: itvtois((Itv)b,&inf,&sup);
! 187: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
! 188: #if 1
! 189: subnum(0,ai,bs,&inf);
! 190: subnum(0,as,bi,&sup);
! 191: istoitv(inf,sup,(Itv *)&tmp);
! 192: addulp((ItvF)tmp, c);
! 193: #else
! 194:
! 195: /* MUST check if ai, as, bi, bs are bigfloat. */
! 196: /* ai - bs = ai + (-bs) */
! 197: chsgnbf(bs,&mbs);
! 198: ltop = avma;
! 199: ritopa(ai,&pa);
! 200: ritopa(mbs,&pb);
! 201: lbot = avma;
! 202: z = gerepile(ltop,lbot,PariAddDown(pa,pb));
! 203: patori(z,&inf); cgiv(z);
! 204:
! 205: /* as - bi = ( - ( bi + (-as) ) ) */
! 206: chsgnbf(as,&mas);
! 207: ltop = avma;
! 208: ritopa(mas,&pa);
! 209: ritopa(bi,&pb);
! 210: lbot = avma;
! 211: z = gerepile(ltop,lbot,PariAddDown(pa,pb));
! 212: patori(z,&tmp); cgiv(z);
! 213:
! 214: chsgnbf(tmp,&sup);
! 215: istoitv(inf,sup,c);
! 216: #endif
! 217: }
! 218: }
! 219:
! 220: void mulitvf(ItvF a, ItvF b, ItvF *c)
! 221: {
! 222: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t,tmp;
! 223: Num inf, sup;
! 224: int ba,bb;
! 225:
! 226: if ( !a || !b )
! 227: *c = 0;
! 228: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
! 229: mulnum(0,(Num)a,(Num)b,(Num *)c);
! 230: else {
! 231: itvtois((Itv)a,&inf,&sup);
! 232: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
! 233: itvtois((Itv)b,&inf,&sup);
! 234: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
! 235:
! 236: /* MUST check if ai, as, bi, bs are bigfloat. */
! 237: chsgnnum(as,&a2);
! 238: chsgnnum(bs,&b2);
! 239:
! 240:
! 241: if ( ba = (compnum(0,0,a2) > 0) ) {
! 242: a1 = ai;
! 243: } else {
! 244: a1 = a2;
! 245: a2 = ai;
! 246: }
! 247: if ( bb = (compnum(0,0,b2) > 0) ) {
! 248: b1 = bi;
! 249: } else {
! 250: b1 = b2;
! 251: b2 = bi;
! 252: }
! 253: mulnum(0,a2,b2,&t);
! 254: subnum(0,0,t,&c2);
! 255: if ( compnum(0,0,b1) > 0 ) {
! 256: mulnum(0,a2,b1,&t);
! 257: subnum(0,0,t,&c1);
! 258: if ( compnum(0,0,a1) > 0 ) {
! 259: mulnum(0,a1,b2,&t);
! 260: subnum(0,0,t,&p);
! 261: if ( compnum(0,c1,p) > 0 ) c1 = p;
! 262: mulnum(0,a1,b1,&t);
! 263: subnum(0,0,t,&p);
! 264: if ( compnum(0,c2,p) > 0 ) c2 = p;
! 265: }
! 266: } else {
! 267: if ( compnum(0,0,a1) > 0 ) {
! 268: mulnum(0,a1,b2,&t);
! 269: subnum(0,0,t,&c1);
! 270: } else {
! 271: mulnum(0,a1,b1,&c1);
! 272: }
! 273: }
! 274: if ( ba == bb ) {
! 275: subnum(0,0,c2,&t);
! 276: istoitv(c1,t,(Itv *)&tmp);
! 277: } else {
! 278: subnum(0,0,c1,&t);
! 279: istoitv(c2,t,(Itv *)&tmp);
! 280: }
! 281: addulp((ItvF)tmp, c);
! 282: }
! 283: }
! 284:
! 285: int initvf(Num a, Itv b)
! 286: {
! 287: Num inf, sup;
! 288:
! 289: itvtois(b, &inf, &sup);
! 290: if ( compnum(0,inf,a) > 0 ) return 0;
! 291: else if ( compnum(0,a,sup) > 0 ) return 0;
! 292: else return 1;
! 293: }
! 294:
! 295: int itvinitvf(Itv a, Itv b)
! 296: {
! 297: Num ai,as,bi,bs;
! 298:
! 299: itvtois(a, &ai, &as);
! 300: itvtois(b, &bi, &bs);
! 301: if ( compnum(0,ai,bi) >= 0 && compnum(0,bs,as) >= 0 ) return 1;
! 302: else return 0;
! 303: }
! 304:
! 305: void divitvf(ItvF a, ItvF b, ItvF *c)
! 306: {
! 307: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t,tmp;
! 308: Num inf, sup;
! 309: int ba,bb;
! 310:
! 311: if ( !b )
! 312: error("divitv : division by 0");
! 313: else if ( !a )
! 314: *c = 0;
! 315: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
! 316: divnum(0,(Num)a,(Num)b,(Num *)c);
! 317: else {
! 318: itvtois((Itv)a,&inf,&sup);
! 319: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
! 320: itvtois((Itv)b,&inf,&sup);
! 321: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
! 322: /* MUST check if ai, as, bi, bs are bigfloat. */
! 323: chsgnnum(as,&a2);
! 324: chsgnnum(bs,&b2);
! 325: if ( ba = (compnum(0,0,a2) > 0) ) {
! 326: a1 = ai;
! 327: } else {
! 328: a1 = a2;
! 329: a2 = ai;
! 330: }
! 331: if ( bb = (compnum(0,0,b2) > 0) ) {
! 332: b1 = bi;
! 333: } else {
! 334: b1 = b2;
! 335: b2 = bi;
! 336: }
! 337: if ( compnum(0,0,b1) >= 0 ) {
! 338: error("divitv : division by interval including 0.");
! 339: } else {
! 340: divnum(0,a2,b1,&c2);
! 341: if ( compnum(0,0,a1) > 0 ) {
! 342: divnum(0,a1,b1,&c1);
! 343: } else {
! 344: divnum(0,a1,b2,&t);
! 345: subnum(0,0,t,&c1);
! 346: }
! 347: }
! 348: if ( ba == bb ) {
! 349: subnum(0,0,c2,&t);
! 350: istoitv(c1,t,(Itv *)&tmp);
! 351: } else {
! 352: subnum(0,0,c1,&t);
! 353: istoitv(c2,t,(Itv *)&tmp);
! 354: }
! 355: addulp((ItvF)tmp, c);
! 356: }
! 357: }
! 358:
! 359: void pwritvf(Itv a, Num e, Itv *c)
! 360: {
! 361: int ei;
! 362: Itv t;
! 363:
! 364: if ( !e )
! 365: *c = (Itv)ONE;
! 366: else if ( !a )
! 367: *c = 0;
! 368: else if ( NID(a) <= N_B )
! 369: pwrnum(0,(Num)a,e,(Num *)c);
! 370: else if ( !INT(e) ) {
! 371: #if PARI && 0
! 372: GEN pa,pe,z;
! 373: int ltop,lbot;
! 374:
! 375: ltop = avma; ritopa(a,&pa); ritopa(e,&pe); lbot = avma;
! 376: z = gerepile(ltop,lbot,gpui(pa,pe,prec));
! 377: patori(z,c); cgiv(z);
! 378: #else
! 379: error("pwritv : can't calculate a fractional power");
! 380: #endif
! 381: } else {
! 382: ei = QTOS((Q)e);
! 383: pwritv0f(a,ei,&t);
! 384: if ( SGN((Q)e) < 0 )
! 385: divbf((Num)ONE,(Num)t,(Num *)c);
! 386: else
! 387: *c = t;
! 388: }
! 389: }
! 390:
! 391: void pwritv0f(Itv a, int e, Itv *c)
! 392: {
! 393: Num inf, sup;
! 394: Num ai,Xmin,Xmax;
! 395: ItvF tmp;
! 396: Q ne;
! 397:
! 398: if ( e == 1 )
! 399: *c = a;
! 400: else {
! 401: STOQ(e,ne);
! 402: if ( !(e%2) ) {
! 403: if ( initvp(0,a) ) {
! 404: Xmin = 0;
! 405: chsgnnum(INF(a),&ai);
! 406: if ( compnum(0,ai,SUP(a)) > 0 ) {
! 407: Xmax = ai;
! 408: } else {
! 409: Xmax = SUP(a);
! 410: }
! 411: } else {
! 412: if ( compnum(0,INF(a),0) > 0 ) {
! 413: Xmin = INF(a);
! 414: Xmax = SUP(a);
! 415: } else {
! 416: Xmin = SUP(a);
! 417: Xmax = INF(a);
! 418: }
! 419: }
! 420: } else {
! 421: Xmin = INF(a);
! 422: Xmax = SUP(a);
! 423: }
! 424: if ( ! Xmin ) inf = 0;
! 425: else pwrbf(Xmin,(Num)ne,&inf);
! 426: if ( ! Xmax ) sup = 0;
! 427: else pwrbf(Xmax,(Num)ne,&sup);
! 428: istoitv(inf,sup,(Itv *)&tmp);
! 429: addulp(tmp, (ItvF *)c);
! 430: }
! 431: }
! 432:
! 433: void chsgnitvf(Itv a, Itv *c)
! 434: {
! 435: Num inf,sup;
! 436:
! 437: if ( !a )
! 438: *c = 0;
! 439: else if ( NID(a) <= N_B )
! 440: chsgnnum((Num)a,(Num *)c);
! 441: else {
! 442: chsgnnum(INF((Itv)a),&inf);
! 443: chsgnnum(SUP((Itv)a),&sup);
! 444: istoitv(inf,sup,c);
! 445: }
! 446: }
! 447:
! 448: int cmpitvf(Itv a, Itv b)
! 449: {
! 450: Num ai,as,bi,bs;
! 451: int s,t;
! 452:
! 453: if ( !a ) {
! 454: if ( !b || (NID(b)<=N_B) )
! 455: return compnum(0,(Num)a,(Num)b);
! 456: else
! 457: return -1;
! 458: } else if ( !b ) {
! 459: if ( !a || (NID(a)<=N_B) )
! 460: return compnum(0,(Num)a,(Num)b);
! 461: else
! 462: return initvp((Num)b,a);
! 463: } else {
! 464: itvtois(a,&ai,&as);
! 465: itvtois(b,&bi,&bs);
! 466: s = compnum(0,ai,bi);
! 467: t = compnum(0,as,bs);
! 468: if ( !s && !t ) return 0;
! 469: else return -1;
! 470: }
! 471: }
! 472:
! 473: void miditvf(Itv a, Num *b)
! 474: {
! 475: Num ai,as;
! 476: Num t;
! 477: Q TWO;
! 478:
! 479: if ( ! a ) *b = 0;
! 480: else if ( (NID(a) <= N_B) )
! 481: *b = (Num)a;
! 482: else {
! 483: STOQ(2,TWO);
! 484: itvtois(a,&ai,&as);
! 485: addbf(ai,as,&t);
! 486: divbf(t,(Num)TWO,b);
! 487: }
! 488: }
! 489:
! 490: void cupitvf(Itv a, Itv b, Itv *c)
! 491: {
! 492: Num ai,as,bi,bs;
! 493: Num inf,sup;
! 494:
! 495: itvtois(a,&ai,&as);
! 496: itvtois(b,&bi,&bs);
! 497: if ( compnum(0,ai,bi) > 0 ) inf = bi;
! 498: else inf = ai;
! 499: if ( compnum(0,as,bs) > 0 ) sup = as;
! 500: else sup = bs;
! 501: istoitv(inf,sup,c);
! 502: }
! 503:
! 504: void capitvf(Itv a, Itv b, Itv *c)
! 505: {
! 506: Num ai,as,bi,bs;
! 507: Num inf,sup;
! 508:
! 509: itvtois(a,&ai,&as);
! 510: itvtois(b,&bi,&bs);
! 511: if ( compnum(0,ai,bi) > 0 ) inf = ai;
! 512: else inf = bi;
! 513: if ( compnum(0,as,bs) > 0 ) sup = bs;
! 514: else sup = as;
! 515: if ( compnum(0,inf,sup) > 0 ) *c = 0;
! 516: else istoitv(inf,sup,c);
! 517: }
! 518:
! 519:
! 520: void widthitvf(Itv a, Num *b)
! 521: {
! 522: Num ai,as;
! 523:
! 524: if ( ! a ) *b = 0;
! 525: else if ( (NID(a) <= N_B) )
! 526: *b = (Num)a;
! 527: else {
! 528: itvtois(a,&ai,&as);
! 529: subnum(0,as,ai,b);
! 530: }
! 531: }
! 532:
! 533: void absitvf(Itv a, Num *b)
! 534: {
! 535: Num ai,as,bi,bs;
! 536:
! 537: if ( ! a ) *b = 0;
! 538: else if ( (NID(a) <= N_B) ) {
! 539: if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
! 540: else *b = (Num)a;
! 541: } else {
! 542: itvtois(a,&ai,&as);
! 543: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
! 544: else bi = ai;
! 545: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
! 546: else bs = as;
! 547: if ( compnum(0,bi,bs) > 0 ) *b = bi;
! 548: else *b = bs;
! 549: }
! 550: }
! 551:
! 552: void distanceitvf(Itv a, Itv b, Num *c)
! 553: {
! 554: Num ai,as,bi,bs;
! 555: Num di,ds;
! 556: Itv d;
! 557:
! 558: itvtois(a,&ai,&as);
! 559: itvtois(b,&bi,&bs);
! 560: subnum(0,ai,bi,&di);
! 561: subnum(0,as,bs,&ds);
! 562: istoitv(di,ds,&d);
! 563: absitvf(d,c);
! 564: }
! 565:
! 566: #endif