Annotation of OpenXM_contrib2/asir2018/engine/t-itv.c, Revision 1.1
1.1 ! noro 1: /*
! 2: * $OpenXM$
! 3: */
! 4: #if defined(INTERVAL)
! 5: #include "ca.h"
! 6: #include "base.h"
! 7: #if defined(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,(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: }
! 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,(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: }
! 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,(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: }
! 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,(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: }
! 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 defined(PARI) && 0
! 231: gpui_ri((Obj)a,(Obj)c,(Obj *)c);
! 232: #else
! 233: error("pwritv : can't calculate a fractional power");
! 234: #endif
! 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: }
! 243: }
! 244:
! 245: void pwritv0p(Itv a, int e, Itv *c)
! 246: {
! 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: }
! 282: }
! 283:
! 284: void chsgnitvp(Itv a, Itv *c)
! 285: {
! 286: Num inf,sup;
! 287:
! 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: }
! 297: }
! 298:
! 299: int cmpitvp(Itv a, Itv b)
! 300: {
! 301: Num ai,as,bi,bs;
! 302: int s,t;
! 303:
! 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: }
! 322: }
! 323:
! 324: void miditvp(Itv a, Num *b)
! 325: {
! 326: Num ai,as;
! 327: Num t;
! 328:
! 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: }
! 338: }
! 339:
! 340: void cupitvp(Itv a, Itv b, Itv *c)
! 341: {
! 342: Num ai,as,bi,bs;
! 343: Num inf,sup;
! 344:
! 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);
! 352: }
! 353:
! 354: void capitvp(Itv a, Itv b, Itv *c)
! 355: {
! 356: Num ai,as,bi,bs;
! 357: Num inf,sup;
! 358:
! 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);
! 367: }
! 368:
! 369:
! 370: void widthitvp(Itv a, Num *b)
! 371: {
! 372: Num ai,as;
! 373:
! 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: }
! 381: }
! 382:
! 383: void absitvp(Itv a, Num *b)
! 384: {
! 385: Num ai,as,bi,bs;
! 386:
! 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: }
! 400: }
! 401:
! 402: void distanceitvp(Itv a, Itv b, Num *c)
! 403: {
! 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);
! 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 / asir2018 / engine