Annotation of OpenXM_contrib2/asir2000/engine/t-itv.c, Revision 1.6
1.1 saito 1: /*
1.6 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/t-itv.c,v 1.5 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: void itvtois(Itv a, Num *inf, Num *sup)
12: {
1.6 ! noro 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: }
1.1 saito 21: }
22:
23: void istoitv(Num inf, Num sup, Itv *rp)
24: {
1.6 ! noro 25: Itv c;
1.1 saito 26:
1.6 ! noro 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: }
1.1 saito 38: }
39:
40: void additvp(Itv a, Itv b, Itv *c)
41: {
1.6 ! noro 42: Num ai,as,bi,bs;
! 43: Num inf,sup;
1.1 saito 44:
1.6 ! noro 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,(IntervalDouble *)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: }
1.1 saito 61: }
62:
63: void subitvp(Itv a, Itv b, Itv *c)
64: {
1.6 ! noro 65: Num ai,as,bi,bs;
! 66: Num inf,sup;
1.1 saito 67:
1.6 ! noro 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,(IntervalDouble *)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: }
1.1 saito 84: }
85:
86: void mulitvp(Itv a, Itv b, Itv *c)
87: {
1.6 ! noro 88: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t;
! 89: int ba,bb;
1.1 saito 90:
1.6 ! noro 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,(IntervalDouble *)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: }
1.1 saito 144: }
145:
146: int initvp(Num a, Itv b)
147: {
1.6 ! noro 148: Num inf, sup;
1.1 saito 149:
1.6 ! noro 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;
1.1 saito 154: }
155:
156: int itvinitvp(Itv a, Itv b)
157: {
1.6 ! noro 158: Num ai,as,bi,bs;
1.1 saito 159:
1.6 ! noro 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;
1.1 saito 164: }
165:
166: void divitvp(Itv a, Itv b, Itv *c)
167: {
1.6 ! noro 168: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t;
! 169: int ba,bb;
1.1 saito 170:
1.6 ! noro 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,(IntervalDouble *)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: }
1.1 saito 216: }
217:
218: void pwritvp(Itv a, Num e, Itv *c)
219: {
1.6 ! noro 220: int ei;
! 221: Itv t;
1.1 saito 222:
1.6 ! noro 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) ) {
1.3 ohara 230: #if defined(PARI) && 0
1.6 ! noro 231: gpui_ri((Obj)a,(Obj)c,(Obj *)c);
1.1 saito 232: #else
1.6 ! noro 233: error("pwritv : can't calculate a fractional power");
1.1 saito 234: #endif
1.6 ! noro 235: } else {
! 236: ei = QTOS((Q)e);
! 237: pwritv0p(a,ei,&t);
! 238: if ( SGN((Q)e) < 0 )
! 239: divnum(0,(Num)ONE,(Num)t,c);
! 240: else
! 241: *c = t;
! 242: }
1.1 saito 243: }
244:
245: void pwritv0p(Itv a, int e, Itv *c)
246: {
1.6 ! noro 247: Num inf, sup;
! 248: Num ai,Xmin,Xmax;
! 249: Q ne;
! 250:
! 251: if ( e == 1 )
! 252: *c = a;
! 253: else {
! 254: STOQ(e,ne);
! 255: if ( !(e%2) ) {
! 256: if ( initvp(0,a) ) {
! 257: Xmin = 0;
! 258: chsgnnum(INF(a),&ai);
! 259: if ( compnum(0,ai,SUP(a)) > 0 ) {
! 260: Xmax = ai;
! 261: } else {
! 262: Xmax = SUP(a);
! 263: }
! 264: } else {
! 265: if ( compnum(0,INF(a),0) > 0 ) {
! 266: Xmin = INF(a);
! 267: Xmax = SUP(a);
! 268: } else {
! 269: Xmin = SUP(a);
! 270: Xmax = INF(a);
! 271: }
! 272: }
! 273: } else {
! 274: Xmin = INF(a);
! 275: Xmax = SUP(a);
! 276: }
! 277: if ( ! Xmin ) inf = 0;
! 278: else pwrnum(0,Xmin,ne,&inf);
! 279: pwrnum(0,Xmax,ne,&sup);
! 280: istoitv(inf,sup,c);
! 281: }
1.1 saito 282: }
283:
284: void chsgnitvp(Itv a, Itv *c)
285: {
1.6 ! noro 286: Num inf,sup;
1.1 saito 287:
1.6 ! noro 288: if ( !a )
! 289: *c = 0;
! 290: else if ( NID(a) <= N_B )
! 291: chsgnnum(a,c);
! 292: else {
! 293: chsgnnum(INF((Itv)a),&inf);
! 294: chsgnnum(SUP((Itv)a),&sup);
! 295: istoitv(inf,sup,c);
! 296: }
1.1 saito 297: }
298:
299: int cmpitvp(Itv a, Itv b)
300: {
1.6 ! noro 301: Num ai,as,bi,bs;
! 302: int s,t;
1.1 saito 303:
1.6 ! noro 304: if ( !a ) {
! 305: if ( !b || (NID(b)<=N_B) )
! 306: return compnum(0,a,b);
! 307: else
! 308: return -1;
! 309: } else if ( !b ) {
! 310: if ( !a || (NID(a)<=N_B) )
! 311: return compnum(0,a,b);
! 312: else
! 313: return initvp((Num)b,a);
! 314: } else {
! 315: itvtois(a,&ai,&as);
! 316: itvtois(b,&bi,&bs);
! 317: s = compnum(0,ai,bi) ;
! 318: t = compnum(0,as,bs) ;
! 319: if ( !s && !t ) return 0;
! 320: else return -1;
! 321: }
1.1 saito 322: }
323:
324: void miditvp(Itv a, Num *b)
325: {
1.6 ! noro 326: Num ai,as;
! 327: Num t;
1.1 saito 328:
1.6 ! noro 329: if ( ! a ) *b = 0;
! 330: else if ( (NID(a) <= N_B) )
! 331: *b = (Num)a;
! 332: else {
! 333: STOQ(2,TWO);
! 334: itvtois(a,&ai,&as);
! 335: addnum(0,ai,as,&t);
! 336: divnum(0,t,TWO,b);
! 337: }
1.1 saito 338: }
339:
340: void cupitvp(Itv a, Itv b, Itv *c)
341: {
1.6 ! noro 342: Num ai,as,bi,bs;
! 343: Num inf,sup;
1.1 saito 344:
1.6 ! noro 345: itvtois(a,&ai,&as);
! 346: itvtois(b,&bi,&bs);
! 347: if ( compnum(0,ai,bi) > 0 ) inf = bi;
! 348: else inf = ai;
! 349: if ( compnum(0,as,bs) > 0 ) sup = as;
! 350: else sup = bs;
! 351: istoitv(inf,sup,c);
1.1 saito 352: }
353:
354: void capitvp(Itv a, Itv b, Itv *c)
355: {
1.6 ! noro 356: Num ai,as,bi,bs;
! 357: Num inf,sup;
1.1 saito 358:
1.6 ! noro 359: itvtois(a,&ai,&as);
! 360: itvtois(b,&bi,&bs);
! 361: if ( compnum(0,ai,bi) > 0 ) inf = ai;
! 362: else inf = bi;
! 363: if ( compnum(0,as,bs) > 0 ) sup = bs;
! 364: else sup = as;
! 365: if ( compnum(0,inf,sup) > 0 ) *c = 0;
! 366: else istoitv(inf,sup,c);
1.1 saito 367: }
368:
369:
370: void widthitvp(Itv a, Num *b)
371: {
1.6 ! noro 372: Num ai,as;
1.1 saito 373:
1.6 ! noro 374: if ( ! a ) *b = 0;
! 375: else if ( (NID(a) <= N_B) )
! 376: *b = (Num)a;
! 377: else {
! 378: itvtois(a,&ai,&as);
! 379: subnum(0,as,ai,b);
! 380: }
1.1 saito 381: }
382:
383: void absitvp(Itv a, Num *b)
384: {
1.6 ! noro 385: Num ai,as,bi,bs;
1.1 saito 386:
1.6 ! noro 387: if ( ! a ) *b = 0;
! 388: else if ( (NID(a) <= N_B) ) {
! 389: if ( compnum(0,a,0) < 0 ) chsgnnum(a,b);
! 390: else *b = (Num)a;
! 391: } else {
! 392: itvtois(a,&ai,&as);
! 393: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
! 394: else bi = ai;
! 395: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
! 396: else bs = as;
! 397: if ( compnum(0,bi,bs) > 0 ) *b = bi;
! 398: else *b = bs;
! 399: }
1.1 saito 400: }
401:
402: void distanceitvp(Itv a, Itv b, Num *c)
403: {
1.6 ! noro 404: Num ai,as,bi,bs;
! 405: Num di,ds;
! 406: Itv d;
! 407:
! 408: itvtois(a,&ai,&as);
! 409: itvtois(b,&bi,&bs);
! 410: subnum(0,ai,bi,&di);
! 411: subnum(0,as,bs,&ds);
! 412: istoitv(di,ds,&d);
! 413: absitvp(d,c);
1.1 saito 414: }
415:
416: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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