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