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