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