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