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