Annotation of OpenXM_contrib2/asir2000/engine/f-itv.c, Revision 1.8
1.1 saito 1: /*
1.8 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/f-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: #include "itv-pari.h"
1.7 ohara 10: long get_pariprec();
1.1 saito 11: #endif
12:
13: void ToBf(Num a, BF *rp)
14: {
1.8 ! noro 15: GEN pa, pb, pc;
! 16: BF bn, bd;
! 17: BF c;
! 18: long ltop, lbot;
! 19:
! 20: if ( ! a ) {
! 21: *rp = 0;
! 22: }
! 23: else switch ( NID(a) ) {
! 24: case N_Q:
! 25: ltop = avma;
! 26: ritopa_i(NM((Q)a), SGN((Q)a), &pa);
! 27: pb = cgetr(get_pariprec());
! 28: mpaff(pa, pb);
! 29: if ( INT((Q)a) ) {
! 30: lbot = avma;
! 31: pb = gerepile(ltop, lbot, pb);
! 32: patori(pb, &c);
! 33: cgiv(pb);
! 34: *rp = c;
! 35: } else {
! 36: patori(pb, &bn);
! 37: ritopa_i(DN((Q)a), 1, &pa);
! 38: pb = cgetr(get_pariprec());
! 39: mpaff(pa, pb);
! 40: lbot = avma;
! 41: pb = gerepile(ltop, lbot, pb);
! 42: patori(pb, &bd);
! 43: cgiv(pb);
! 44: divbf((Num)bn,(Num)bd, (Num *)&c);
! 45: *rp = c;
! 46: }
! 47: break;
! 48: case N_R:
! 49: pb = dbltor(BDY((Real)a));
! 50: patori(pb, &c);
! 51: cgiv(pb);
! 52: *rp = c;
! 53: break;
! 54: case N_B:
! 55: *rp = (BF)a;
! 56: break;
! 57: default:
! 58: error("ToBf : not implemented");
! 59: break;
! 60: }
1.1 saito 61: }
62:
63: void double2bf(double d, BF *rp)
64: {
1.8 ! noro 65: GEN p;
1.1 saito 66:
1.8 ! noro 67: p = dbltor(d);
! 68: patori(p, rp);
! 69: cgiv(p);
1.1 saito 70: }
71:
1.8 ! noro 72: double bf2double(BF a)
1.1 saito 73: {
1.8 ! noro 74: GEN p;
1.1 saito 75:
1.8 ! noro 76: ritopa(a, &p);
! 77: return rtodbl(p);
1.1 saito 78: }
79:
80: void getulp(BF a, BF *au)
81: {
1.8 ! noro 82: GEN b, c;
! 83: long lb, i;
1.1 saito 84:
1.8 ! noro 85: ritopa(a, &b);
! 86: lb = lg(b);
! 87: c = cgetr(lb);
! 88: c[1] = b[1];
! 89: for (i=2;i<lb-1;i++) c[i] = 0;
! 90: c[i] = 1;
! 91: setsigne(c, 1);
! 92: patori(c,au);
! 93: cgiv(c);
! 94: cgiv(b);
1.1 saito 95: }
96:
1.2 kondoh 97: void addulp(IntervalBigFloat a, IntervalBigFloat *c)
1.1 saito 98: {
1.8 ! noro 99: Num ai, as, aiu, asu, inf, sup;
1.1 saito 100:
1.8 ! noro 101: itvtois((Itv)a,&ai,&as);
! 102: if ( ai ) {
! 103: getulp((BF)ai, (BF *)&aiu);
! 104: subbf(ai,aiu,&inf);
! 105: } else inf = ai;
! 106: if ( as ) {
! 107: getulp((BF)as, (BF *)&asu);
! 108: addbf(as,asu,&sup);
! 109: } else sup = as;
! 110: istoitv(inf,sup, (Itv *)c);
1.1 saito 111: }
112:
1.2 kondoh 113: void additvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 114: {
1.8 ! noro 115: Num ai,as,bi,bs,mas,mbs,tmp;
! 116: Num inf,sup;
! 117: GEN pa, pb, z;
! 118: long ltop, lbot;
! 119:
! 120: if ( !a )
! 121: *c = b;
! 122: else if ( !b )
! 123: *c = a;
! 124: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
! 125: addnum(0,(Num)a,(Num)b,(Num *)c);
! 126: else {
! 127: itvtois((Itv)a,&inf,&sup);
! 128: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
! 129: itvtois((Itv)b,&inf,&sup);
! 130: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
1.1 saito 131: #if 0
132: printexpr(CO, ai);
133: printexpr(CO, as);
134: printexpr(CO, bi);
135: printexpr(CO, bs);
136: #endif
137:
1.8 ! noro 138: addnum(0,ai,bi,&inf);
! 139: addnum(0,as,bs,&sup);
! 140: istoitv(inf,sup,(Itv *)&tmp);
! 141: addulp((IntervalBigFloat)tmp, c);
! 142: return;
! 143: }
1.1 saito 144: }
145:
1.2 kondoh 146: void subitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 147: {
1.8 ! noro 148: Num ai,as,bi,bs,mas, mbs;
! 149: Num inf,sup,tmp;
! 150: GEN pa, pb, z;
! 151: long ltop, lbot;
! 152:
! 153: if ( !a )
! 154: chsgnnum((Num)b,(Num *)c);
! 155: else if ( !b )
! 156: *c = a;
! 157: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
! 158: subnum(0,(Num)a,(Num)b,(Num *)c);
! 159: else {
! 160: itvtois((Itv)a,&inf,&sup);
! 161: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
! 162: itvtois((Itv)b,&inf,&sup);
! 163: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
! 164: subnum(0,ai,bs,&inf);
! 165: subnum(0,as,bi,&sup);
! 166: istoitv(inf,sup,(Itv *)&tmp);
! 167: addulp((IntervalBigFloat)tmp, c);
! 168: }
1.1 saito 169: }
170:
1.2 kondoh 171: void mulitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 172: {
1.8 ! noro 173: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t,tmp;
! 174: Num inf, sup;
! 175: int ba,bb;
! 176:
! 177: if ( !a || !b )
! 178: *c = 0;
! 179: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
! 180: mulnum(0,(Num)a,(Num)b,(Num *)c);
! 181: else {
! 182: itvtois((Itv)a,&inf,&sup);
! 183: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
! 184: itvtois((Itv)b,&inf,&sup);
! 185: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
! 186:
! 187: /* MUST check if ai, as, bi, bs are bigfloat. */
! 188: chsgnnum(as,&a2);
! 189: chsgnnum(bs,&b2);
! 190:
! 191:
! 192: if ( ba = (compnum(0,0,a2) > 0) ) {
! 193: a1 = ai;
! 194: } else {
! 195: a1 = a2;
! 196: a2 = ai;
! 197: }
! 198: if ( bb = (compnum(0,0,b2) > 0) ) {
! 199: b1 = bi;
! 200: } else {
! 201: b1 = b2;
! 202: b2 = bi;
! 203: }
! 204: mulnum(0,a2,b2,&t);
! 205: subnum(0,0,t,&c2);
! 206: if ( compnum(0,0,b1) > 0 ) {
! 207: mulnum(0,a2,b1,&t);
! 208: subnum(0,0,t,&c1);
! 209: if ( compnum(0,0,a1) > 0 ) {
! 210: mulnum(0,a1,b2,&t);
! 211: subnum(0,0,t,&p);
! 212: if ( compnum(0,c1,p) > 0 ) c1 = p;
! 213: mulnum(0,a1,b1,&t);
! 214: subnum(0,0,t,&p);
! 215: if ( compnum(0,c2,p) > 0 ) c2 = p;
! 216: }
! 217: } else {
! 218: if ( compnum(0,0,a1) > 0 ) {
! 219: mulnum(0,a1,b2,&t);
! 220: subnum(0,0,t,&c1);
! 221: } else {
! 222: mulnum(0,a1,b1,&c1);
! 223: }
! 224: }
! 225: if ( ba == bb ) {
! 226: subnum(0,0,c2,&t);
! 227: istoitv(c1,t,(Itv *)&tmp);
! 228: } else {
! 229: subnum(0,0,c1,&t);
! 230: istoitv(c2,t,(Itv *)&tmp);
! 231: }
! 232: addulp((IntervalBigFloat)tmp, c);
! 233: }
1.1 saito 234: }
235:
236: int initvf(Num a, Itv b)
237: {
1.8 ! noro 238: Num inf, sup;
1.1 saito 239:
1.8 ! noro 240: itvtois(b, &inf, &sup);
! 241: if ( compnum(0,inf,a) > 0 ) return 0;
! 242: else if ( compnum(0,a,sup) > 0 ) return 0;
! 243: else return 1;
1.1 saito 244: }
245:
246: int itvinitvf(Itv a, Itv b)
247: {
1.8 ! noro 248: Num ai,as,bi,bs;
1.1 saito 249:
1.8 ! noro 250: itvtois(a, &ai, &as);
! 251: itvtois(b, &bi, &bs);
! 252: if ( compnum(0,ai,bi) >= 0 && compnum(0,bs,as) >= 0 ) return 1;
! 253: else return 0;
1.1 saito 254: }
255:
1.2 kondoh 256: void divitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 257: {
1.8 ! noro 258: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t,tmp;
! 259: Num inf, sup;
! 260: int ba,bb;
! 261:
! 262: if ( !b )
! 263: error("divitv : division by 0");
! 264: else if ( !a )
! 265: *c = 0;
! 266: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
! 267: divnum(0,(Num)a,(Num)b,(Num *)c);
! 268: else {
! 269: itvtois((Itv)a,&inf,&sup);
! 270: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
! 271: itvtois((Itv)b,&inf,&sup);
! 272: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
1.1 saito 273: /* MUST check if ai, as, bi, bs are bigfloat. */
1.8 ! noro 274: chsgnnum(as,&a2);
! 275: chsgnnum(bs,&b2);
! 276: if ( ba = (compnum(0,0,a2) > 0) ) {
! 277: a1 = ai;
! 278: } else {
! 279: a1 = a2;
! 280: a2 = ai;
! 281: }
! 282: if ( bb = (compnum(0,0,b2) > 0) ) {
! 283: b1 = bi;
! 284: } else {
! 285: b1 = b2;
! 286: b2 = bi;
! 287: }
! 288: if ( compnum(0,0,b1) >= 0 ) {
! 289: error("divitv : division by interval including 0.");
! 290: } else {
! 291: divnum(0,a2,b1,&c2);
! 292: if ( compnum(0,0,a1) > 0 ) {
! 293: divnum(0,a1,b1,&c1);
! 294: } else {
! 295: divnum(0,a1,b2,&t);
! 296: subnum(0,0,t,&c1);
! 297: }
! 298: }
! 299: if ( ba == bb ) {
! 300: subnum(0,0,c2,&t);
! 301: istoitv(c1,t,(Itv *)&tmp);
! 302: } else {
! 303: subnum(0,0,c1,&t);
! 304: istoitv(c2,t,(Itv *)&tmp);
! 305: }
! 306: addulp((IntervalBigFloat)tmp, c);
! 307: }
1.1 saito 308: }
309:
310: void pwritvf(Itv a, Num e, Itv *c)
311: {
1.8 ! noro 312: int ei;
! 313: Itv t;
1.1 saito 314:
1.8 ! noro 315: if ( !e )
! 316: *c = (Itv)ONE;
! 317: else if ( !a )
! 318: *c = 0;
! 319: else if ( NID(a) <= N_B )
! 320: pwrnum(0,(Num)a,e,(Num *)c);
! 321: else if ( !INT(e) ) {
1.3 ohara 322: #if defined(PARI) && 0
1.8 ! noro 323: gpui_ri((Obj)a,(Obj)c,(Obj *)c);
1.1 saito 324: #else
1.8 ! noro 325: error("pwritv : can't calculate a fractional power");
1.1 saito 326: #endif
1.8 ! noro 327: } else {
! 328: ei = QTOS((Q)e);
! 329: pwritv0f(a,ei,&t);
! 330: if ( SGN((Q)e) < 0 )
! 331: divbf((Num)ONE,(Num)t,(Num *)c);
! 332: else
! 333: *c = t;
! 334: }
1.1 saito 335: }
336:
337: void pwritv0f(Itv a, int e, Itv *c)
338: {
1.8 ! noro 339: Num inf, sup;
! 340: Num ai,Xmin,Xmax;
! 341: IntervalBigFloat tmp;
! 342: Q ne;
! 343:
! 344: if ( e == 1 )
! 345: *c = a;
! 346: else {
! 347: STOQ(e,ne);
! 348: if ( !(e%2) ) {
! 349: if ( initvp(0,a) ) {
! 350: Xmin = 0;
! 351: chsgnnum(INF(a),&ai);
! 352: if ( compnum(0,ai,SUP(a)) > 0 ) {
! 353: Xmax = ai;
! 354: } else {
! 355: Xmax = SUP(a);
! 356: }
! 357: } else {
! 358: if ( compnum(0,INF(a),0) > 0 ) {
! 359: Xmin = INF(a);
! 360: Xmax = SUP(a);
! 361: } else {
! 362: Xmin = SUP(a);
! 363: Xmax = INF(a);
! 364: }
! 365: }
! 366: } else {
! 367: Xmin = INF(a);
! 368: Xmax = SUP(a);
! 369: }
! 370: if ( ! Xmin ) inf = 0;
! 371: else pwrbf(Xmin,(Num)ne,&inf);
! 372: if ( ! Xmax ) sup = 0;
! 373: else pwrbf(Xmax,(Num)ne,&sup);
! 374: istoitv(inf,sup,(Itv *)&tmp);
! 375: addulp(tmp, (IntervalBigFloat *)c);
! 376: }
1.1 saito 377: }
378:
379: void chsgnitvf(Itv a, Itv *c)
380: {
1.8 ! noro 381: Num inf,sup;
1.1 saito 382:
1.8 ! noro 383: if ( !a )
! 384: *c = 0;
! 385: else if ( NID(a) <= N_B )
! 386: chsgnnum((Num)a,(Num *)c);
! 387: else {
! 388: chsgnnum(INF((Itv)a),&inf);
! 389: chsgnnum(SUP((Itv)a),&sup);
! 390: istoitv(inf,sup,c);
! 391: }
1.1 saito 392: }
393:
394: int cmpitvf(Itv a, Itv b)
395: {
1.8 ! noro 396: Num ai,as,bi,bs;
! 397: int s,t;
1.1 saito 398:
1.8 ! noro 399: if ( !a ) {
! 400: if ( !b || (NID(b)<=N_B) ) {
! 401: return compnum(0,(Num)a,(Num)b);
! 402: } else {
! 403: itvtois(b,&bi,&bs);
! 404: if ( compnum(0,(Num)a,bs) > 0 ) return 1;
! 405: else if ( compnum(0,bi,(Num)a) > 0 ) return -1;
! 406: else return 0;
! 407: }
! 408: } else if ( !b ) {
! 409: if ( !a || (NID(a)<=N_B) ) {
! 410: return compnum(0,(Num)a,(Num)b);
! 411: } else {
! 412: itvtois(a,&ai,&as);
! 413: if ( compnum(0,ai,(Num)b) > 0 ) return 1;
! 414: else if ( compnum(0,(Num)b,as) > 0 ) return -1;
! 415: else return 0;
! 416: }
! 417: } else {
! 418: itvtois(a,&ai,&as);
! 419: itvtois(b,&bi,&bs);
! 420: s = compnum(0,ai,bs) ;
! 421: t = compnum(0,bi,as) ;
! 422: if ( s > 0 ) return 1;
! 423: else if ( t > 0 ) return -1;
! 424: else return 0;
! 425: }
1.1 saito 426: }
427:
428: void miditvf(Itv a, Num *b)
429: {
1.8 ! noro 430: Num ai,as;
! 431: Num t;
1.1 saito 432:
1.8 ! noro 433: if ( ! a ) *b = 0;
! 434: else if ( (NID(a) <= N_B) )
! 435: *b = (Num)a;
! 436: else {
! 437: STOQ(2,TWO);
! 438: itvtois(a,&ai,&as);
! 439: addbf(ai,as,&t);
! 440: divbf(t,(Num)TWO,b);
! 441: }
1.1 saito 442: }
443:
444: void cupitvf(Itv a, Itv b, Itv *c)
445: {
1.8 ! noro 446: Num ai,as,bi,bs;
! 447: Num inf,sup;
1.1 saito 448:
1.8 ! noro 449: itvtois(a,&ai,&as);
! 450: itvtois(b,&bi,&bs);
! 451: if ( compnum(0,ai,bi) > 0 ) inf = bi;
! 452: else inf = ai;
! 453: if ( compnum(0,as,bs) > 0 ) sup = as;
! 454: else sup = bs;
! 455: istoitv(inf,sup,c);
1.1 saito 456: }
457:
458: void capitvf(Itv a, Itv b, Itv *c)
459: {
1.8 ! noro 460: Num ai,as,bi,bs;
! 461: Num inf,sup;
1.1 saito 462:
1.8 ! noro 463: itvtois(a,&ai,&as);
! 464: itvtois(b,&bi,&bs);
! 465: if ( compnum(0,ai,bi) > 0 ) inf = ai;
! 466: else inf = bi;
! 467: if ( compnum(0,as,bs) > 0 ) sup = bs;
! 468: else sup = as;
! 469: if ( compnum(0,inf,sup) > 0 ) *c = 0;
! 470: else istoitv(inf,sup,c);
1.1 saito 471: }
472:
473:
474: void widthitvf(Itv a, Num *b)
475: {
1.8 ! noro 476: Num ai,as;
1.1 saito 477:
1.8 ! noro 478: if ( ! a ) *b = 0;
! 479: else if ( (NID(a) <= N_B) )
! 480: *b = (Num)a;
! 481: else {
! 482: itvtois(a,&ai,&as);
! 483: subnum(0,as,ai,b);
! 484: }
1.1 saito 485: }
486:
487: void absitvf(Itv a, Num *b)
488: {
1.8 ! noro 489: Num ai,as,bi,bs;
1.1 saito 490:
1.8 ! noro 491: if ( ! a ) *b = 0;
! 492: else if ( (NID(a) <= N_B) ) {
! 493: if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
! 494: else *b = (Num)a;
! 495: } else {
! 496: itvtois(a,&ai,&as);
! 497: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
! 498: else bi = ai;
! 499: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
! 500: else bs = as;
! 501: if ( compnum(0,bi,bs) > 0 ) *b = bi;
! 502: else *b = bs;
! 503: }
1.1 saito 504: }
505:
506: void distanceitvf(Itv a, Itv b, Num *c)
507: {
1.8 ! noro 508: Num ai,as,bi,bs;
! 509: Num di,ds;
! 510: Itv d;
! 511:
! 512: itvtois(a,&ai,&as);
! 513: itvtois(b,&bi,&bs);
! 514: subnum(0,ai,bi,&di);
! 515: subnum(0,as,bs,&ds);
! 516: istoitv(di,ds,&d);
! 517: absitvf(d,c);
1.1 saito 518: }
519:
520: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to f-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