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