![[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 / asir2000 / engine|
Return to t-itv.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 saito 1: /*
1.5 ! ohara 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/t-itv.c,v 1.4 2003/10/20 07:18:42 saito Exp $
1.1 saito 3: */
4: #if defined(INTERVAL)
5: #include "ca.h"
6: #include "base.h"
1.3 ohara 7: #if defined(PARI)
1.1 saito 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) )
1.2 kondoh 53: additvd((Num)a,(Num)b,(IntervalDouble *)c);
1.1 saito 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) )
1.2 kondoh 76: subitvd((Num)a,(Num)b,(IntervalDouble *)c);
1.1 saito 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) )
1.2 kondoh 97: mulitvd((Num)a,(Num)b,(IntervalDouble *)c);
1.1 saito 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) )
1.2 kondoh 179: divitvd((Num)a,(Num)b,(IntervalDouble *)c);
1.1 saito 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: {
220: int ei;
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) ) {
1.3 ohara 230: #if defined(PARI) && 0
1.5 ! ohara 231: gpui_ri((Obj)a,(Obj)c,(Obj *)c);
1.1 saito 232: #else
233: error("pwritv : can't calculate a fractional power");
234: #endif
235: } else {
236: ei = QTOS((Q)e);
237: pwritv0p(a,ei,&t);
238: if ( SGN((Q)e) < 0 )
239: divnum(0,(Num)ONE,(Num)t,c);
240: else
241: *c = t;
242: }
243: }
244:
245: void pwritv0p(Itv a, int e, Itv *c)
246: {
247: Num inf, sup;
248: Num ai,Xmin,Xmax;
249: Q ne;
250:
251: if ( e == 1 )
252: *c = a;
253: else {
254: STOQ(e,ne);
255: if ( !(e%2) ) {
256: if ( initvp(0,a) ) {
257: Xmin = 0;
258: chsgnnum(INF(a),&ai);
259: if ( compnum(0,ai,SUP(a)) > 0 ) {
260: Xmax = ai;
261: } else {
262: Xmax = SUP(a);
263: }
264: } else {
265: if ( compnum(0,INF(a),0) > 0 ) {
266: Xmin = INF(a);
267: Xmax = SUP(a);
268: } else {
269: Xmin = SUP(a);
270: Xmax = INF(a);
271: }
272: }
273: } else {
274: Xmin = INF(a);
275: Xmax = SUP(a);
276: }
277: if ( ! Xmin ) inf = 0;
278: else pwrnum(0,Xmin,ne,&inf);
279: pwrnum(0,Xmax,ne,&sup);
280: istoitv(inf,sup,c);
281: }
282: }
283:
284: void chsgnitvp(Itv a, Itv *c)
285: {
286: Num inf,sup;
287:
288: if ( !a )
289: *c = 0;
290: else if ( NID(a) <= N_B )
291: chsgnnum(a,c);
292: else {
293: chsgnnum(INF((Itv)a),&inf);
294: chsgnnum(SUP((Itv)a),&sup);
295: istoitv(inf,sup,c);
296: }
297: }
298:
299: int cmpitvp(Itv a, Itv b)
300: {
301: Num ai,as,bi,bs;
302: int s,t;
303:
304: if ( !a ) {
305: if ( !b || (NID(b)<=N_B) )
306: return compnum(0,a,b);
307: else
308: return -1;
309: } else if ( !b ) {
310: if ( !a || (NID(a)<=N_B) )
311: return compnum(0,a,b);
312: else
313: return initvp((Num)b,a);
314: } else {
315: itvtois(a,&ai,&as);
316: itvtois(b,&bi,&bs);
317: s = compnum(0,ai,bi) ;
318: t = compnum(0,as,bs) ;
319: if ( !s && !t ) return 0;
320: else return -1;
321: }
322: }
323:
324: void miditvp(Itv a, Num *b)
325: {
326: Num ai,as;
327: Num t;
328:
329: if ( ! a ) *b = 0;
330: else if ( (NID(a) <= N_B) )
331: *b = (Num)a;
332: else {
333: STOQ(2,TWO);
334: itvtois(a,&ai,&as);
335: addnum(0,ai,as,&t);
336: divnum(0,t,TWO,b);
337: }
338: }
339:
340: void cupitvp(Itv a, Itv b, Itv *c)
341: {
342: Num ai,as,bi,bs;
343: Num inf,sup;
344:
345: itvtois(a,&ai,&as);
346: itvtois(b,&bi,&bs);
347: if ( compnum(0,ai,bi) > 0 ) inf = bi;
348: else inf = ai;
349: if ( compnum(0,as,bs) > 0 ) sup = as;
350: else sup = bs;
351: istoitv(inf,sup,c);
352: }
353:
354: void capitvp(Itv a, Itv b, Itv *c)
355: {
356: Num ai,as,bi,bs;
357: Num inf,sup;
358:
359: itvtois(a,&ai,&as);
360: itvtois(b,&bi,&bs);
361: if ( compnum(0,ai,bi) > 0 ) inf = ai;
362: else inf = bi;
363: if ( compnum(0,as,bs) > 0 ) sup = bs;
364: else sup = as;
365: if ( compnum(0,inf,sup) > 0 ) *c = 0;
366: else istoitv(inf,sup,c);
367: }
368:
369:
370: void widthitvp(Itv a, Num *b)
371: {
372: Num ai,as;
373:
374: if ( ! a ) *b = 0;
375: else if ( (NID(a) <= N_B) )
376: *b = (Num)a;
377: else {
378: itvtois(a,&ai,&as);
379: subnum(0,as,ai,b);
380: }
381: }
382:
383: void absitvp(Itv a, Num *b)
384: {
385: Num ai,as,bi,bs;
386:
387: if ( ! a ) *b = 0;
388: else if ( (NID(a) <= N_B) ) {
389: if ( compnum(0,a,0) < 0 ) chsgnnum(a,b);
390: else *b = (Num)a;
391: } else {
392: itvtois(a,&ai,&as);
393: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
394: else bi = ai;
395: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
396: else bs = as;
397: if ( compnum(0,bi,bs) > 0 ) *b = bi;
398: else *b = bs;
399: }
400: }
401:
402: void distanceitvp(Itv a, Itv b, Num *c)
403: {
404: Num ai,as,bi,bs;
405: Num di,ds;
406: Itv d;
407:
408: itvtois(a,&ai,&as);
409: itvtois(b,&bi,&bs);
410: subnum(0,ai,bi,&di);
411: subnum(0,as,bs,&ds);
412: istoitv(di,ds,&d);
413: absitvp(d,c);
414: }
415:
416: #endif