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