Annotation of OpenXM_contrib2/asir2018/engine/f-itv.c, Revision 1.4
1.1 noro 1: /*
1.4 ! kondoh 2: * $OpenXM: OpenXM_contrib2/asir2018/engine/f-itv.c,v 1.3 2019/06/04 07:11:23 kondoh 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: {
1.4 ! kondoh 390: long ei;
1.1 noro 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.4 ! kondoh 406: //ei = QTOS((Q)e);
! 407: ei = mpz_get_si(BDY((Q)e));
1.3 kondoh 408: if (ei<0) ei = -ei;
1.1 noro 409: pwritv0f(a,ei,&t);
1.4 ! kondoh 410: // if ( SGN((Q)e) < 0 )
! 411: if ( sgnq((Q)e) < 0 )
1.3 kondoh 412: divitvf((IntervalBigFloat)ONE,(IntervalBigFloat)t,(IntervalBigFloat *)c); /* bug ?? */
1.1 noro 413: else
414: *c = t;
415: }
416: }
417:
1.4 ! kondoh 418: void pwritv0f(Itv a, long e, Itv *c)
1.1 noro 419: {
420: Num inf, sup;
421: Num ai,Xmin,Xmax;
422: IntervalBigFloat tmp;
423: Q ne;
1.3 kondoh 424: int current_roundmode;
1.1 noro 425:
426: if ( e == 1 )
427: *c = a;
428: else {
1.4 ! kondoh 429: STOZ(e,ne);
1.1 noro 430: if ( !(e%2) ) {
431: if ( initvp(0,a) ) {
432: Xmin = 0;
433: chsgnnum(INF(a),&ai);
434: if ( compnum(0,ai,SUP(a)) > 0 ) {
435: Xmax = ai;
436: } else {
437: Xmax = SUP(a);
438: }
439: } else {
440: if ( compnum(0,INF(a),0) > 0 ) {
441: Xmin = INF(a);
442: Xmax = SUP(a);
443: } else {
444: Xmin = SUP(a);
445: Xmax = INF(a);
446: }
447: }
448: } else {
449: Xmin = INF(a);
450: Xmax = SUP(a);
451: }
1.3 kondoh 452:
453: current_roundmode = mpfr_roundmode;
1.1 noro 454: if ( ! Xmin ) inf = 0;
1.3 kondoh 455: else {
456: mpfr_roundmode = MPFR_RNDD;
457: pwrbf(Xmin,(Num)ne,&inf);
458: }
1.1 noro 459: if ( ! Xmax ) sup = 0;
1.3 kondoh 460: else {
461: mpfr_roundmode = MPFR_RNDU;
462: pwrbf(Xmax,(Num)ne,&sup);
463: }
1.1 noro 464: istoitv(inf,sup,(Itv *)&tmp);
1.3 kondoh 465: mpfr_roundmode = current_roundmode;
466: *c = (Itv)tmp;
467: // addulp(tmp, (IntervalBigFloat *)c);
1.1 noro 468: }
469: }
470:
1.3 kondoh 471: void chsgnitvf(IntervalBigFloat a, IntervalBigFloat *rp)
1.1 noro 472: {
473: Num inf,sup;
1.3 kondoh 474: IntervalBigFloat c;
475: int current_roundmode;
1.1 noro 476:
477: if ( !a )
1.3 kondoh 478: *rp = 0;
479: else if ( NID(a) < N_IntervalBigFloat )
480: chsgnnum((Num)a,(Num *)rp);
1.1 noro 481: else {
1.3 kondoh 482: current_roundmode = mpfr_roundmode;
483:
484: mpfr_roundmode = MPFR_RNDU;
485: chsgnnum((Num)INF(a),&sup);
486: mpfr_roundmode = MPFR_RNDD;
487: chsgnnum((Num)SUP(a),&inf);
488: //MKIntervalBigFloat((BF)inf,(BF)sup,c);
489: istoitv(inf,sup,(Itv *)&c);
490: *rp = c;
491: mpfr_roundmode = current_roundmode;
1.1 noro 492: }
493: }
494:
495: int cmpitvf(Itv a, Itv b)
496: {
497: Num ai,as,bi,bs;
498: int s,t;
499:
500: if ( !a ) {
501: if ( !b || (NID(b)<=N_B) ) {
502: return compnum(0,(Num)a,(Num)b);
503: } else {
504: itvtois(b,&bi,&bs);
505: if ( compnum(0,(Num)a,bs) > 0 ) return 1;
506: else if ( compnum(0,bi,(Num)a) > 0 ) return -1;
507: else return 0;
508: }
509: } else if ( !b ) {
510: if ( !a || (NID(a)<=N_B) ) {
511: return compnum(0,(Num)a,(Num)b);
512: } else {
513: itvtois(a,&ai,&as);
514: if ( compnum(0,ai,(Num)b) > 0 ) return 1;
515: else if ( compnum(0,(Num)b,as) > 0 ) return -1;
516: else return 0;
517: }
518: } else {
519: itvtois(a,&ai,&as);
520: itvtois(b,&bi,&bs);
521: s = compnum(0,ai,bs) ;
522: t = compnum(0,bi,as) ;
523: if ( s > 0 ) return 1;
524: else if ( t > 0 ) return -1;
525: else return 0;
526: }
527: }
528:
529: void miditvf(Itv a, Num *b)
530: {
531: Num ai,as;
532: Num t;
533:
534: if ( ! a ) *b = 0;
535: else if ( (NID(a) <= N_B) )
536: *b = (Num)a;
537: else {
1.4 ! kondoh 538: //STOQ(2,TWO);
1.1 noro 539: itvtois(a,&ai,&as);
540: addbf(ai,as,&t);
541: divbf(t,(Num)TWO,b);
542: }
543: }
544:
545: void cupitvf(Itv a, Itv b, Itv *c)
546: {
547: Num ai,as,bi,bs;
548: Num inf,sup;
549:
550: itvtois(a,&ai,&as);
551: itvtois(b,&bi,&bs);
552: if ( compnum(0,ai,bi) > 0 ) inf = bi;
553: else inf = ai;
554: if ( compnum(0,as,bs) > 0 ) sup = as;
555: else sup = bs;
556: istoitv(inf,sup,c);
557: }
558:
559: void capitvf(Itv a, Itv b, Itv *c)
560: {
561: Num ai,as,bi,bs;
562: Num inf,sup;
563:
564: itvtois(a,&ai,&as);
565: itvtois(b,&bi,&bs);
566: if ( compnum(0,ai,bi) > 0 ) inf = ai;
567: else inf = bi;
568: if ( compnum(0,as,bs) > 0 ) sup = bs;
569: else sup = as;
570: if ( compnum(0,inf,sup) > 0 ) *c = 0;
571: else istoitv(inf,sup,c);
572: }
573:
574:
575: void widthitvf(Itv a, Num *b)
576: {
577: Num ai,as;
578:
579: if ( ! a ) *b = 0;
580: else if ( (NID(a) <= N_B) )
581: *b = (Num)a;
582: else {
583: itvtois(a,&ai,&as);
584: subnum(0,as,ai,b);
585: }
586: }
587:
588: void absitvf(Itv a, Num *b)
589: {
590: Num ai,as,bi,bs;
591:
592: if ( ! a ) *b = 0;
593: else if ( (NID(a) <= N_B) ) {
594: if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
595: else *b = (Num)a;
596: } else {
597: itvtois(a,&ai,&as);
598: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
599: else bi = ai;
600: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
601: else bs = as;
602: if ( compnum(0,bi,bs) > 0 ) *b = bi;
603: else *b = bs;
604: }
605: }
606:
607: void distanceitvf(Itv a, Itv b, Num *c)
608: {
609: Num ai,as,bi,bs;
610: Num di,ds;
611: Itv d;
612:
613: itvtois(a,&ai,&as);
614: itvtois(b,&bi,&bs);
615: subnum(0,ai,bi,&di);
616: subnum(0,as,bs,&ds);
617: istoitv(di,ds,&d);
618: absitvf(d,c);
619: }
620:
621: #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