![[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.4 ! kondoh 2: * $OpenXM: OpenXM_contrib2/asir2018/engine/p-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: #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: {
1.4 ! kondoh 275: long ei;
1.1 noro 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.4 ! kondoh 291: //ei = QTOS((Q)e);
! 292: ei = mpz_get_si(BDY((Q)e));
1.1 noro 293: pwritv0p(a,ei,&t);
1.4 ! kondoh 294: // if ( SGN((Q)e) < 0 )
! 295: if ( sgnq((Q)e) < 0 )
1.1 noro 296: divnum(0,(Num)ONE,(Num)t,(Num *)c);
297: else
298: *c = t;
299: }
300: }
301:
1.4 ! kondoh 302: void pwritv0p(Itv a, long e, Itv *c)
1.1 noro 303: {
304: Num inf, sup;
305: Num ai,Xmin,Xmax;
306: Q ne;
307:
308: if ( e == 1 )
309: *c = a;
310: else {
1.4 ! kondoh 311: STOZ(e,ne);
1.1 noro 312: if ( !(e%2) ) {
313: if ( initvp(0,a) ) {
314: Xmin = 0;
315: chsgnnum(INF(a),&ai);
316: if ( compnum(0,ai,SUP(a)) > 0 ) {
317: Xmax = ai;
318: } else {
319: Xmax = SUP(a);
320: }
321: } else {
322: if ( compnum(0,INF(a),0) > 0 ) {
323: Xmin = INF(a);
324: Xmax = SUP(a);
325: } else {
326: Xmin = SUP(a);
327: Xmax = INF(a);
328: }
329: }
330: } else {
331: Xmin = INF(a);
332: Xmax = SUP(a);
333: }
334: if ( ! Xmin ) inf = 0;
335: else pwrnum(0,Xmin,(Num)ne,&inf);
336: pwrnum(0,Xmax,(Num)ne,&sup);
337: istoitv(inf,sup,c);
338: }
339: }
340:
341: void chsgnitvp(Itv a, Itv *c)
342: {
343: Num inf,sup;
344:
345: if ( !a )
346: *c = 0;
347: else if ( NID(a) <= N_B )
348: chsgnnum((Num)a,(Num *)c);
349: else {
350: chsgnnum(INF((Itv)a),&inf);
351: chsgnnum(SUP((Itv)a),&sup);
352: istoitv(inf,sup,c);
353: }
354: }
355:
356: int cmpitvp(Itv a, Itv b)
357: {
358: Num ai,as,bi,bs;
359: int s,t;
360:
361: if ( !a ) {
362: if ( !b || (NID(b)<=N_B) ) {
363: return compnum(0,(Num)a,(Num)b);
364: } else {
365: itvtois(b,&bi,&bs);
366: if ( compnum(0,(Num)a,bs) > 0 ) return 1;
367: else if ( compnum(0,bi,(Num)a) > 0 ) return -1;
368: else return 0;
369: }
370: } else if ( !b ) {
371: if ( !a || (NID(a)<=N_B) ) {
372: return compnum(0,(Num)a,(Num)b);
373: } else {
374: itvtois(a,&ai,&as);
375: if ( compnum(0,ai,(Num)b) > 0 ) return 1;
376: else if ( compnum(0,(Num)b,as) > 0 ) return -1;
377: else return 0;
378: }
379: } else {
380: itvtois(a,&ai,&as);
381: itvtois(b,&bi,&bs);
382: s = compnum(0,ai,bs) ;
383: t = compnum(0,bi,as) ;
384: if ( s > 0 ) return 1;
385: else if ( t > 0 ) return -1;
386: else return 0;
387: }
388: }
389:
390: void miditvp(Itv a, Num *b)
391: {
392: Num ai,as;
393: Num t;
394:
395: if ( ! a ) *b = 0;
396: else if ( (NID(a) <= N_B) )
397: *b = (Num)a;
398: else {
1.4 ! kondoh 399: //STOZ(2,TWO);
1.1 noro 400: itvtois(a,&ai,&as);
401: addnum(0,ai,as,&t);
402: divnum(0,t,(Num)TWO,b);
403: }
404: }
405:
406: void cupitvp(Itv a, Itv b, Itv *c)
407: {
408: Num ai,as,bi,bs;
409: Num inf,sup;
410:
411: itvtois(a,&ai,&as);
412: itvtois(b,&bi,&bs);
413: if ( compnum(0,ai,bi) > 0 ) inf = bi;
414: else inf = ai;
415: if ( compnum(0,as,bs) > 0 ) sup = as;
416: else sup = bs;
417: istoitv(inf,sup,c);
418: }
419:
420: void capitvp(Itv a, Itv b, Itv *c)
421: {
422: Num ai,as,bi,bs;
423: Num inf,sup;
424:
425: itvtois(a,&ai,&as);
426: itvtois(b,&bi,&bs);
427: if ( compnum(0,ai,bi) > 0 ) inf = ai;
428: else inf = bi;
429: if ( compnum(0,as,bs) > 0 ) sup = bs;
430: else sup = as;
431: if ( compnum(0,inf,sup) > 0 ) *c = 0;
432: else istoitv(inf,sup,c);
433: }
434:
435:
436: void widthitvp(Itv a, Num *b)
437: {
438: Num ai,as;
439:
440: if ( ! a ) *b = 0;
441: else if ( (NID(a) <= N_B) )
442: *b = (Num)a;
443: else {
444: itvtois(a,&ai,&as);
445: subnum(0,as,ai,b);
446: }
447: }
448:
449: void absitvp(Itv a, Num *b)
450: {
451: Num ai,as,bi,bs;
452:
453: if ( ! a ) *b = 0;
454: else if ( (NID(a) <= N_B) ) {
455: if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
456: else *b = (Num)a;
457: } else {
458: itvtois(a,&ai,&as);
459: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
460: else bi = ai;
461: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
462: else bs = as;
463: if ( compnum(0,bi,bs) > 0 ) *b = bi;
464: else *b = bs;
465: }
466: }
467:
468: void distanceitvp(Itv a, Itv b, Num *c)
469: {
470: Num ai,as,bi,bs;
471: Num di,ds;
472: Itv d;
473:
474: itvtois(a,&ai,&as);
475: itvtois(b,&bi,&bs);
476: subnum(0,ai,bi,&di);
477: subnum(0,as,bs,&ds);
478: istoitv(di,ds,&d);
479: absitvp(d,c);
480: }
481:
482: #endif