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