Annotation of OpenXM_contrib2/asir2000/engine/p-itv.c, Revision 1.2
1.1 saito 1: /*
1.2 ! kondoh 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/p-itv.c,v 1.1 2000/12/22 10:03:28 saito Exp $
1.1 saito 3: */
4: #if defined(INTERVAL)
5: #include "ca.h"
6: #include "base.h"
7: #if PARI
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) ) {
262: #if PARI && 0
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 ) {
342: if ( !b || (NID(b)<=N_B) )
343: return compnum(0,(Num)a,(Num)b);
344: else
345: return -1;
346: } else if ( !b ) {
347: if ( !a || (NID(a)<=N_B) )
348: return compnum(0,(Num)a,(Num)b);
349: else
350: return initvp((Num)b,a);
351: } else {
352: itvtois(a,&ai,&as);
353: itvtois(b,&bi,&bs);
1.2 ! kondoh 354: s = compnum(0,ai,bs) ;
! 355: t = compnum(0,bi,as) ;
! 356: if ( s > 0 ) return 1;
! 357: else if ( t > 0 ) return -1;
! 358: else return 0;
1.1 saito 359: }
360: }
361:
362: void miditvp(Itv a, Num *b)
363: {
364: Num ai,as;
365: Num t;
366: Q TWO;
367:
368: if ( ! a ) *b = 0;
369: else if ( (NID(a) <= N_B) )
370: *b = (Num)a;
371: else {
372: STOQ(2,TWO);
373: itvtois(a,&ai,&as);
374: addnum(0,ai,as,&t);
375: divnum(0,t,(Num)TWO,b);
376: }
377: }
378:
379: void cupitvp(Itv a, Itv b, Itv *c)
380: {
381: Num ai,as,bi,bs;
382: Num inf,sup;
383:
384: itvtois(a,&ai,&as);
385: itvtois(b,&bi,&bs);
386: if ( compnum(0,ai,bi) > 0 ) inf = bi;
387: else inf = ai;
388: if ( compnum(0,as,bs) > 0 ) sup = as;
389: else sup = bs;
390: istoitv(inf,sup,c);
391: }
392:
393: void capitvp(Itv a, Itv b, Itv *c)
394: {
395: Num ai,as,bi,bs;
396: Num inf,sup;
397:
398: itvtois(a,&ai,&as);
399: itvtois(b,&bi,&bs);
400: if ( compnum(0,ai,bi) > 0 ) inf = ai;
401: else inf = bi;
402: if ( compnum(0,as,bs) > 0 ) sup = bs;
403: else sup = as;
404: if ( compnum(0,inf,sup) > 0 ) *c = 0;
405: else istoitv(inf,sup,c);
406: }
407:
408:
409: void widthitvp(Itv a, Num *b)
410: {
411: Num ai,as;
412:
413: if ( ! a ) *b = 0;
414: else if ( (NID(a) <= N_B) )
415: *b = (Num)a;
416: else {
417: itvtois(a,&ai,&as);
418: subnum(0,as,ai,b);
419: }
420: }
421:
422: void absitvp(Itv a, Num *b)
423: {
424: Num ai,as,bi,bs;
425:
426: if ( ! a ) *b = 0;
427: else if ( (NID(a) <= N_B) ) {
428: if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
429: else *b = (Num)a;
430: } else {
431: itvtois(a,&ai,&as);
432: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
433: else bi = ai;
434: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
435: else bs = as;
436: if ( compnum(0,bi,bs) > 0 ) *b = bi;
437: else *b = bs;
438: }
439: }
440:
441: void distanceitvp(Itv a, Itv b, Num *c)
442: {
443: Num ai,as,bi,bs;
444: Num di,ds;
445: Itv d;
446:
447: itvtois(a,&ai,&as);
448: itvtois(b,&bi,&bs);
449: subnum(0,ai,bi,&di);
450: subnum(0,as,bs,&ds);
451: istoitv(di,ds,&d);
452: absitvp(d,c);
453: }
454:
455: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>