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