![[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 / asir2018 / engine|
Return to p-itv.c CVS log | Up to [local] / OpenXM_contrib2 / asir2018 / engine |
1.1 noro 1: /*
1.2 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2018/engine/p-itv.c,v 1.1 2018/09/19 05:45:07 noro Exp $
1.1 noro 3: */
4: #if defined(INTERVAL)
5: #include "ca.h"
6: #include "base.h"
7: #if defined(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 ) {
54: // NEWIntervalBigFloat((IntervalBigFloat)c);
55: c=MALLOC(sizeof(struct oIntervalBigFloat));
56: OID(c)=O_N;
57: NID(c)=N_IntervalBigFloat;
58: } else
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) )
88: additvd((Num)a,(Num)b,(IntervalDouble *)c);
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) )
111: subitvd((Num)a,(Num)b,(IntervalDouble *)c);
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) )
132: mulitvd((Num)a,(Num)b,(IntervalDouble *)c);
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) )
214: divitvd((Num)a,(Num)b,(IntervalDouble *)c);
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) ) {
265: #if defined(PARI) && 0
266: gpui_ri((Obj)a,(Obj)c,(Obj *)c);
267: #else
268: error("pwritv : can't calculate a fractional power");
269: #endif
270: } else {
1.2 ! noro 271: ei = ZTOS((Q)e);
1.1 noro 272: pwritv0p(a,ei,&t);
273: if ( SGN((Q)e) < 0 )
274: divnum(0,(Num)ONE,(Num)t,(Num *)c);
275: else
276: *c = t;
277: }
278: }
279:
280: void pwritv0p(Itv a, int e, Itv *c)
281: {
282: Num inf, sup;
283: Num ai,Xmin,Xmax;
284: Q ne;
285:
286: if ( e == 1 )
287: *c = a;
288: else {
1.2 ! noro 289: STOZ(e,ne);
1.1 noro 290: if ( !(e%2) ) {
291: if ( initvp(0,a) ) {
292: Xmin = 0;
293: chsgnnum(INF(a),&ai);
294: if ( compnum(0,ai,SUP(a)) > 0 ) {
295: Xmax = ai;
296: } else {
297: Xmax = SUP(a);
298: }
299: } else {
300: if ( compnum(0,INF(a),0) > 0 ) {
301: Xmin = INF(a);
302: Xmax = SUP(a);
303: } else {
304: Xmin = SUP(a);
305: Xmax = INF(a);
306: }
307: }
308: } else {
309: Xmin = INF(a);
310: Xmax = SUP(a);
311: }
312: if ( ! Xmin ) inf = 0;
313: else pwrnum(0,Xmin,(Num)ne,&inf);
314: pwrnum(0,Xmax,(Num)ne,&sup);
315: istoitv(inf,sup,c);
316: }
317: }
318:
319: void chsgnitvp(Itv a, Itv *c)
320: {
321: Num inf,sup;
322:
323: if ( !a )
324: *c = 0;
325: else if ( NID(a) <= N_B )
326: chsgnnum((Num)a,(Num *)c);
327: else {
328: chsgnnum(INF((Itv)a),&inf);
329: chsgnnum(SUP((Itv)a),&sup);
330: istoitv(inf,sup,c);
331: }
332: }
333:
334: int cmpitvp(Itv a, Itv b)
335: {
336: Num ai,as,bi,bs;
337: int s,t;
338:
339: if ( !a ) {
340: if ( !b || (NID(b)<=N_B) ) {
341: return compnum(0,(Num)a,(Num)b);
342: } else {
343: itvtois(b,&bi,&bs);
344: if ( compnum(0,(Num)a,bs) > 0 ) return 1;
345: else if ( compnum(0,bi,(Num)a) > 0 ) return -1;
346: else return 0;
347: }
348: } else if ( !b ) {
349: if ( !a || (NID(a)<=N_B) ) {
350: return compnum(0,(Num)a,(Num)b);
351: } else {
352: itvtois(a,&ai,&as);
353: if ( compnum(0,ai,(Num)b) > 0 ) return 1;
354: else if ( compnum(0,(Num)b,as) > 0 ) return -1;
355: else return 0;
356: }
357: } else {
358: itvtois(a,&ai,&as);
359: itvtois(b,&bi,&bs);
360: s = compnum(0,ai,bs) ;
361: t = compnum(0,bi,as) ;
362: if ( s > 0 ) return 1;
363: else if ( t > 0 ) return -1;
364: else return 0;
365: }
366: }
367:
368: void miditvp(Itv a, Num *b)
369: {
370: Num ai,as;
371: Num t;
372:
373: if ( ! a ) *b = 0;
374: else if ( (NID(a) <= N_B) )
375: *b = (Num)a;
376: else {
1.2 ! noro 377: STOZ(2,TWO);
1.1 noro 378: itvtois(a,&ai,&as);
379: addnum(0,ai,as,&t);
380: divnum(0,t,(Num)TWO,b);
381: }
382: }
383:
384: void cupitvp(Itv a, Itv b, Itv *c)
385: {
386: Num ai,as,bi,bs;
387: Num inf,sup;
388:
389: itvtois(a,&ai,&as);
390: itvtois(b,&bi,&bs);
391: if ( compnum(0,ai,bi) > 0 ) inf = bi;
392: else inf = ai;
393: if ( compnum(0,as,bs) > 0 ) sup = as;
394: else sup = bs;
395: istoitv(inf,sup,c);
396: }
397:
398: void capitvp(Itv a, Itv b, Itv *c)
399: {
400: Num ai,as,bi,bs;
401: Num inf,sup;
402:
403: itvtois(a,&ai,&as);
404: itvtois(b,&bi,&bs);
405: if ( compnum(0,ai,bi) > 0 ) inf = ai;
406: else inf = bi;
407: if ( compnum(0,as,bs) > 0 ) sup = bs;
408: else sup = as;
409: if ( compnum(0,inf,sup) > 0 ) *c = 0;
410: else istoitv(inf,sup,c);
411: }
412:
413:
414: void widthitvp(Itv a, Num *b)
415: {
416: Num ai,as;
417:
418: if ( ! a ) *b = 0;
419: else if ( (NID(a) <= N_B) )
420: *b = (Num)a;
421: else {
422: itvtois(a,&ai,&as);
423: subnum(0,as,ai,b);
424: }
425: }
426:
427: void absitvp(Itv a, Num *b)
428: {
429: Num ai,as,bi,bs;
430:
431: if ( ! a ) *b = 0;
432: else if ( (NID(a) <= N_B) ) {
433: if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
434: else *b = (Num)a;
435: } else {
436: itvtois(a,&ai,&as);
437: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
438: else bi = ai;
439: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
440: else bs = as;
441: if ( compnum(0,bi,bs) > 0 ) *b = bi;
442: else *b = bs;
443: }
444: }
445:
446: void distanceitvp(Itv a, Itv b, Num *c)
447: {
448: Num ai,as,bi,bs;
449: Num di,ds;
450: Itv d;
451:
452: itvtois(a,&ai,&as);
453: itvtois(b,&bi,&bs);
454: subnum(0,ai,bi,&di);
455: subnum(0,as,bs,&ds);
456: istoitv(di,ds,&d);
457: absitvp(d,c);
458: }
459:
460: #endif