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