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