![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to p-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 p-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: #endif
! 10:
! 11: extern int zerorewrite;
! 12:
! 13: void itvtois(Itv a, Num *inf, Num *sup)
! 14: {
! 15: if ( !a )
! 16: *inf = *sup = (Num)0;
! 17: else if ( NID(a) <= N_B ) {
! 18: *inf = (Num)a; *sup = (Num)a;
! 19: } else {
! 20: *inf = INF(a);
! 21: *sup = SUP(a);
! 22: }
! 23: }
! 24:
! 25: void istoitv(Num inf, Num sup, Itv *rp)
! 26: {
! 27: Num i, s;
! 28: Itv c;
! 29: int type=0;
! 30:
! 31: if ( !inf && !sup ) {
! 32: *rp = 0;
! 33: return;
! 34: }
! 35: if ( !sup ) {
! 36: s = 0;
! 37: }
! 38: else if ( NID( sup )==N_B ) {
! 39: type = 1;
! 40: ToBf(sup, (BF *)&s);
! 41: } else {
! 42: s = sup;
! 43: }
! 44: if ( !inf ) {
! 45: i = 0;
! 46: }
! 47: else if ( NID( inf )==N_B ) {
! 48: type = 1;
! 49: ToBf(inf, (BF *)&i);
! 50: } else {
! 51: i = inf;
! 52: }
! 53: if ( type )
! 54: NEWItvF((ItvF)c);
! 55: else
! 56: NEWItvP(c);
! 57:
! 58: if ( compnum(0,i,s) >= 0 ) {
! 59: INF(c) = s; SUP(c) = i;
! 60: } else {
! 61: INF(c) = i; SUP(c) = s;
! 62: }
! 63:
! 64: if (zerorewrite)
! 65: if ( initvp(0,c) ) {
! 66: *rp = 0;
! 67: return;
! 68: }
! 69: *rp = c;
! 70: }
! 71:
! 72: void additvp(Itv a, Itv b, Itv *c)
! 73: {
! 74: Num ai,as,bi,bs;
! 75: Num inf,sup;
! 76:
! 77: if ( !a )
! 78: *c = b;
! 79: else if ( !b )
! 80: *c = a;
! 81: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
! 82: addnum(0,(Num)a,(Num)b,(Num *)c);
! 83: else if ( (NID(a) == N_IP) && (NID(b) == N_R )
! 84: ||(NID(a) == N_R ) && (NID(b) == N_IP) )
! 85: additvd((Num)a,(Num)b,(ItvD *)c);
! 86: else {
! 87: itvtois(a,&ai,&as);
! 88: itvtois(b,&bi,&bs);
! 89: addnum(0,ai,bi,&inf);
! 90: addnum(0,as,bs,&sup);
! 91: istoitv(inf,sup,c);
! 92: }
! 93: }
! 94:
! 95: void subitvp(Itv a, Itv b, Itv *c)
! 96: {
! 97: Num ai,as,bi,bs;
! 98: Num inf,sup;
! 99:
! 100: if ( !a )
! 101: chsgnnum((Num)b,(Num *)c);
! 102: else if ( !b )
! 103: *c = a;
! 104: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
! 105: subnum(0,(Num)a,(Num)b,(Num *)c);
! 106: else if ( (NID(a) == N_IP) && (NID(b) == N_R )
! 107: ||(NID(a) == N_R ) && (NID(b) == N_IP) )
! 108: subitvd((Num)a,(Num)b,(ItvD *)c);
! 109: else {
! 110: itvtois(a,&ai,&as);
! 111: itvtois(b,&bi,&bs);
! 112: subnum(0,ai,bs,&inf);
! 113: subnum(0,as,bi,&sup);
! 114: istoitv(inf,sup,c);
! 115: }
! 116: }
! 117:
! 118: void mulitvp(Itv a, Itv b, Itv *c)
! 119: {
! 120: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t;
! 121: int ba,bb;
! 122:
! 123: if ( !a || !b )
! 124: *c = 0;
! 125: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
! 126: mulnum(0,(Num)a,(Num)b,(Num *)c);
! 127: else if ( (NID(a) == N_IP) && (NID(b) == N_R )
! 128: ||(NID(a) == N_R ) && (NID(b) == N_IP) )
! 129: mulitvd((Num)a,(Num)b,(ItvD *)c);
! 130: else {
! 131: itvtois(a,&ai,&as);
! 132: itvtois(b,&bi,&bs);
! 133: chsgnnum(as,&a2);
! 134: chsgnnum(bs,&b2);
! 135: if ( ba = (compnum(0,0,a2) > 0) ) {
! 136: a1 = ai;
! 137: } else {
! 138: a1 = a2;
! 139: a2 = ai;
! 140: }
! 141: if ( bb = (compnum(0,0,b2) > 0) ) {
! 142: b1 = bi;
! 143: } else {
! 144: b1 = b2;
! 145: b2 = bi;
! 146: }
! 147: mulnum(0,a2,b2,&t);
! 148: subnum(0,0,t,&c2);
! 149: if ( compnum(0,0,b1) > 0 ) {
! 150: mulnum(0,a2,b1,&t);
! 151: subnum(0,0,t,&c1);
! 152: if ( compnum(0,0,a1) > 0 ) {
! 153: mulnum(0,a1,b2,&t);
! 154: subnum(0,0,t,&p);
! 155: if ( compnum(0,c1,p) > 0 ) c1 = p;
! 156: mulnum(0,a1,b1,&t);
! 157: subnum(0,0,t,&p);
! 158: if ( compnum(0,c2,p) > 0 ) c2 = p;
! 159: }
! 160: } else {
! 161: if ( compnum(0,0,a1) > 0 ) {
! 162: mulnum(0,a1,b2,&t);
! 163: subnum(0,0,t,&c1);
! 164: } else {
! 165: mulnum(0,a1,b1,&c1);
! 166: }
! 167: }
! 168: if ( ba == bb ) {
! 169: subnum(0,0,c2,&t);
! 170: istoitv(c1,t,c);
! 171: } else {
! 172: subnum(0,0,c1,&t);
! 173: istoitv(c2,t,c);
! 174: }
! 175: }
! 176: }
! 177:
! 178: int initvp(Num a, Itv b)
! 179: {
! 180: Num inf, sup;
! 181:
! 182: itvtois(b, &inf, &sup);
! 183: if ( compnum(0,inf,a) > 0 ) return 0;
! 184: else if ( compnum(0,a,sup) > 0 ) return 0;
! 185: else return 1;
! 186: }
! 187:
! 188: int itvinitvp(Itv a, Itv b)
! 189: {
! 190: Num ai,as,bi,bs;
! 191:
! 192: itvtois(a, &ai, &as);
! 193: itvtois(b, &bi, &bs);
! 194: if ( compnum(0,ai,bi) >= 0 && compnum(0,bs,as) >= 0 ) return 1;
! 195: else return 0;
! 196: }
! 197:
! 198: void divitvp(Itv a, Itv b, Itv *c)
! 199: {
! 200: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t;
! 201: int ba,bb;
! 202:
! 203: if ( !b )
! 204: error("divitv : division by 0");
! 205: else if ( !a )
! 206: *c = 0;
! 207: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
! 208: divnum(0,(Num)a,(Num)b,(Num *)c);
! 209: else if ( (NID(a) == N_IP) && (NID(b) == N_R )
! 210: ||(NID(a) == N_R ) && (NID(b) == N_IP) )
! 211: divitvd((Num)a,(Num)b,(ItvD *)c);
! 212: else {
! 213: itvtois(a,&ai,&as);
! 214: itvtois(b,&bi,&bs);
! 215: chsgnnum(as,&a2);
! 216: chsgnnum(bs,&b2);
! 217: if ( ba = (compnum(0,0,a2) > 0) ) {
! 218: a1 = ai;
! 219: } else {
! 220: a1 = a2;
! 221: a2 = ai;
! 222: }
! 223: if ( bb = (compnum(0,0,b2) > 0) ) {
! 224: b1 = bi;
! 225: } else {
! 226: b1 = b2;
! 227: b2 = bi;
! 228: }
! 229: if ( compnum(0,0,b1) >= 0 ) {
! 230: error("divitv : division by interval including 0.");
! 231: } else {
! 232: divnum(0,a2,b1,&c2);
! 233: if ( compnum(0,0,a1) > 0 ) {
! 234: divnum(0,a1,b1,&c1);
! 235: } else {
! 236: divnum(0,a1,b2,&t);
! 237: subnum(0,0,t,&c1);
! 238: }
! 239: }
! 240: if ( ba == bb ) {
! 241: subnum(0,0,c2,&t);
! 242: istoitv(c1,t,c);
! 243: } else {
! 244: subnum(0,0,c1,&t);
! 245: istoitv(c2,t,c);
! 246: }
! 247: }
! 248: }
! 249:
! 250: void pwritvp(Itv a, Num e, Itv *c)
! 251: {
! 252: int ei;
! 253: Itv t;
! 254:
! 255: if ( !e )
! 256: *c = (Itv)ONE;
! 257: else if ( !a )
! 258: *c = 0;
! 259: else if ( NID(a) <= N_B )
! 260: pwrnum(0,(Num)a,(Num)e,(Num *)c);
! 261: else if ( !INT(e) ) {
! 262: #if PARI && 0
! 263: GEN pa,pe,z;
! 264: int ltop,lbot;
! 265:
! 266: ltop = avma; ritopa(a,&pa); ritopa(e,&pe); lbot = avma;
! 267: z = gerepile(ltop,lbot,gpui(pa,pe,prec));
! 268: patori(z,c); cgiv(z);
! 269: #else
! 270: error("pwritv : can't calculate a fractional power");
! 271: #endif
! 272: } else {
! 273: ei = QTOS((Q)e);
! 274: pwritv0p(a,ei,&t);
! 275: if ( SGN((Q)e) < 0 )
! 276: divnum(0,(Num)ONE,(Num)t,(Num *)c);
! 277: else
! 278: *c = t;
! 279: }
! 280: }
! 281:
! 282: void pwritv0p(Itv a, int e, Itv *c)
! 283: {
! 284: Num inf, sup;
! 285: Num ai,Xmin,Xmax;
! 286: Q ne;
! 287:
! 288: if ( e == 1 )
! 289: *c = a;
! 290: else {
! 291: STOQ(e,ne);
! 292: if ( !(e%2) ) {
! 293: if ( initvp(0,a) ) {
! 294: Xmin = 0;
! 295: chsgnnum(INF(a),&ai);
! 296: if ( compnum(0,ai,SUP(a)) > 0 ) {
! 297: Xmax = ai;
! 298: } else {
! 299: Xmax = SUP(a);
! 300: }
! 301: } else {
! 302: if ( compnum(0,INF(a),0) > 0 ) {
! 303: Xmin = INF(a);
! 304: Xmax = SUP(a);
! 305: } else {
! 306: Xmin = SUP(a);
! 307: Xmax = INF(a);
! 308: }
! 309: }
! 310: } else {
! 311: Xmin = INF(a);
! 312: Xmax = SUP(a);
! 313: }
! 314: if ( ! Xmin ) inf = 0;
! 315: else pwrnum(0,Xmin,(Num)ne,&inf);
! 316: pwrnum(0,Xmax,(Num)ne,&sup);
! 317: istoitv(inf,sup,c);
! 318: }
! 319: }
! 320:
! 321: void chsgnitvp(Itv a, Itv *c)
! 322: {
! 323: Num inf,sup;
! 324:
! 325: if ( !a )
! 326: *c = 0;
! 327: else if ( NID(a) <= N_B )
! 328: chsgnnum((Num)a,(Num *)c);
! 329: else {
! 330: chsgnnum(INF((Itv)a),&inf);
! 331: chsgnnum(SUP((Itv)a),&sup);
! 332: istoitv(inf,sup,c);
! 333: }
! 334: }
! 335:
! 336: int cmpitvp(Itv a, Itv b)
! 337: {
! 338: Num ai,as,bi,bs;
! 339: int s,t;
! 340:
! 341: if ( !a ) {
! 342: if ( !b || (NID(b)<=N_B) )
! 343: return compnum(0,(Num)a,(Num)b);
! 344: else
! 345: return -1;
! 346: } else if ( !b ) {
! 347: if ( !a || (NID(a)<=N_B) )
! 348: return compnum(0,(Num)a,(Num)b);
! 349: else
! 350: return initvp((Num)b,a);
! 351: } else {
! 352: itvtois(a,&ai,&as);
! 353: itvtois(b,&bi,&bs);
! 354: s = compnum(0,ai,bi) ;
! 355: t = compnum(0,as,bs) ;
! 356: if ( !s && !t ) return 0;
! 357: else return -1;
! 358: }
! 359: }
! 360:
! 361: void miditvp(Itv a, Num *b)
! 362: {
! 363: Num ai,as;
! 364: Num t;
! 365: Q TWO;
! 366:
! 367: if ( ! a ) *b = 0;
! 368: else if ( (NID(a) <= N_B) )
! 369: *b = (Num)a;
! 370: else {
! 371: STOQ(2,TWO);
! 372: itvtois(a,&ai,&as);
! 373: addnum(0,ai,as,&t);
! 374: divnum(0,t,(Num)TWO,b);
! 375: }
! 376: }
! 377:
! 378: void cupitvp(Itv a, Itv b, Itv *c)
! 379: {
! 380: Num ai,as,bi,bs;
! 381: Num inf,sup;
! 382:
! 383: itvtois(a,&ai,&as);
! 384: itvtois(b,&bi,&bs);
! 385: if ( compnum(0,ai,bi) > 0 ) inf = bi;
! 386: else inf = ai;
! 387: if ( compnum(0,as,bs) > 0 ) sup = as;
! 388: else sup = bs;
! 389: istoitv(inf,sup,c);
! 390: }
! 391:
! 392: void capitvp(Itv a, Itv b, Itv *c)
! 393: {
! 394: Num ai,as,bi,bs;
! 395: Num inf,sup;
! 396:
! 397: itvtois(a,&ai,&as);
! 398: itvtois(b,&bi,&bs);
! 399: if ( compnum(0,ai,bi) > 0 ) inf = ai;
! 400: else inf = bi;
! 401: if ( compnum(0,as,bs) > 0 ) sup = bs;
! 402: else sup = as;
! 403: if ( compnum(0,inf,sup) > 0 ) *c = 0;
! 404: else istoitv(inf,sup,c);
! 405: }
! 406:
! 407:
! 408: void widthitvp(Itv a, Num *b)
! 409: {
! 410: Num ai,as;
! 411:
! 412: if ( ! a ) *b = 0;
! 413: else if ( (NID(a) <= N_B) )
! 414: *b = (Num)a;
! 415: else {
! 416: itvtois(a,&ai,&as);
! 417: subnum(0,as,ai,b);
! 418: }
! 419: }
! 420:
! 421: void absitvp(Itv a, Num *b)
! 422: {
! 423: Num ai,as,bi,bs;
! 424:
! 425: if ( ! a ) *b = 0;
! 426: else if ( (NID(a) <= N_B) ) {
! 427: if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
! 428: else *b = (Num)a;
! 429: } else {
! 430: itvtois(a,&ai,&as);
! 431: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
! 432: else bi = ai;
! 433: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
! 434: else bs = as;
! 435: if ( compnum(0,bi,bs) > 0 ) *b = bi;
! 436: else *b = bs;
! 437: }
! 438: }
! 439:
! 440: void distanceitvp(Itv a, Itv b, Num *c)
! 441: {
! 442: Num ai,as,bi,bs;
! 443: Num di,ds;
! 444: Itv d;
! 445:
! 446: itvtois(a,&ai,&as);
! 447: itvtois(b,&bi,&bs);
! 448: subnum(0,ai,bi,&di);
! 449: subnum(0,as,bs,&ds);
! 450: istoitv(di,ds,&d);
! 451: absitvp(d,c);
! 452: }
! 453:
! 454: #endif