![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to f-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 f-itv.c CVS log | Up to [local] / OpenXM_contrib2 / asir2018 / engine |
1.1 noro 1: /*
1.2 ! noro 2: * $OpenXM: OpenXM_contrib2/asir2018/engine/f-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: #include "itv-pari.h"
10: long get_pariprec();
11: #endif
12:
13: void ToBf(Num a, BF *rp)
14: {
15: GEN pa, pb, pc;
16: BF bn, bd;
17: BF c;
18: long ltop, lbot;
19:
20: if ( ! a ) {
21: *rp = 0;
22: }
23: else switch ( NID(a) ) {
24: case N_Q:
25: ltop = avma;
26: ritopa_i(NM((Q)a), SGN((Q)a), &pa);
27: pb = cgetr(get_pariprec());
28: mpaff(pa, pb);
29: if ( INT((Q)a) ) {
30: lbot = avma;
31: pb = gerepile(ltop, lbot, pb);
32: patori(pb, &c);
33: cgiv(pb);
34: *rp = c;
35: } else {
36: patori(pb, &bn);
37: ritopa_i(DN((Q)a), 1, &pa);
38: pb = cgetr(get_pariprec());
39: mpaff(pa, pb);
40: lbot = avma;
41: pb = gerepile(ltop, lbot, pb);
42: patori(pb, &bd);
43: cgiv(pb);
44: divbf((Num)bn,(Num)bd, (Num *)&c);
45: *rp = c;
46: }
47: break;
48: case N_R:
49: pb = dbltor(BDY((Real)a));
50: patori(pb, &c);
51: cgiv(pb);
52: *rp = c;
53: break;
54: case N_B:
55: *rp = (BF)a;
56: break;
57: default:
58: error("ToBf : not implemented");
59: break;
60: }
61: }
62:
63: void double2bf(double d, BF *rp)
64: {
65: GEN p;
66:
67: p = dbltor(d);
68: patori(p, rp);
69: cgiv(p);
70: }
71:
72: double bf2double(BF a)
73: {
74: GEN p;
75:
76: ritopa(a, &p);
77: return rtodbl(p);
78: }
79:
80: void getulp(BF a, BF *au)
81: {
82: GEN b, c;
83: long lb, i;
84:
85: ritopa(a, &b);
86: lb = lg(b);
87: c = cgetr(lb);
88: c[1] = b[1];
89: for (i=2;i<lb-1;i++) c[i] = 0;
90: c[i] = 1;
91: setsigne(c, 1);
92: patori(c,au);
93: cgiv(c);
94: cgiv(b);
95: }
96:
97: void addulp(IntervalBigFloat a, IntervalBigFloat *c)
98: {
99: Num ai, as, aiu, asu, inf, sup;
100:
101: itvtois((Itv)a,&ai,&as);
102: if ( ai ) {
103: getulp((BF)ai, (BF *)&aiu);
104: subbf(ai,aiu,&inf);
105: } else inf = ai;
106: if ( as ) {
107: getulp((BF)as, (BF *)&asu);
108: addbf(as,asu,&sup);
109: } else sup = as;
110: istoitv(inf,sup, (Itv *)c);
111: }
112:
113: void additvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
114: {
115: Num ai,as,bi,bs,mas,mbs,tmp;
116: Num inf,sup;
117: GEN pa, pb, z;
118: long ltop, lbot;
119:
120: if ( !a )
121: *c = b;
122: else if ( !b )
123: *c = a;
124: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
125: addnum(0,(Num)a,(Num)b,(Num *)c);
126: else {
127: itvtois((Itv)a,&inf,&sup);
128: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
129: itvtois((Itv)b,&inf,&sup);
130: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
131: #if 0
132: printexpr(CO, ai);
133: printexpr(CO, as);
134: printexpr(CO, bi);
135: printexpr(CO, bs);
136: #endif
137:
138: addnum(0,ai,bi,&inf);
139: addnum(0,as,bs,&sup);
140: istoitv(inf,sup,(Itv *)&tmp);
141: addulp((IntervalBigFloat)tmp, c);
142: return;
143: }
144: }
145:
146: void subitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
147: {
148: Num ai,as,bi,bs,mas, mbs;
149: Num inf,sup,tmp;
150: GEN pa, pb, z;
151: long ltop, lbot;
152:
153: if ( !a )
154: chsgnnum((Num)b,(Num *)c);
155: else if ( !b )
156: *c = a;
157: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
158: subnum(0,(Num)a,(Num)b,(Num *)c);
159: else {
160: itvtois((Itv)a,&inf,&sup);
161: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
162: itvtois((Itv)b,&inf,&sup);
163: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
164: subnum(0,ai,bs,&inf);
165: subnum(0,as,bi,&sup);
166: istoitv(inf,sup,(Itv *)&tmp);
167: addulp((IntervalBigFloat)tmp, c);
168: }
169: }
170:
171: void mulitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
172: {
173: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t,tmp;
174: Num inf, sup;
175: int ba,bb;
176:
177: if ( !a || !b )
178: *c = 0;
179: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
180: mulnum(0,(Num)a,(Num)b,(Num *)c);
181: else {
182: itvtois((Itv)a,&inf,&sup);
183: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
184: itvtois((Itv)b,&inf,&sup);
185: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
186:
187: /* MUST check if ai, as, bi, bs are bigfloat. */
188: chsgnnum(as,&a2);
189: chsgnnum(bs,&b2);
190:
191:
192: if ( ba = (compnum(0,0,a2) > 0) ) {
193: a1 = ai;
194: } else {
195: a1 = a2;
196: a2 = ai;
197: }
198: if ( bb = (compnum(0,0,b2) > 0) ) {
199: b1 = bi;
200: } else {
201: b1 = b2;
202: b2 = bi;
203: }
204: mulnum(0,a2,b2,&t);
205: subnum(0,0,t,&c2);
206: if ( compnum(0,0,b1) > 0 ) {
207: mulnum(0,a2,b1,&t);
208: subnum(0,0,t,&c1);
209: if ( compnum(0,0,a1) > 0 ) {
210: mulnum(0,a1,b2,&t);
211: subnum(0,0,t,&p);
212: if ( compnum(0,c1,p) > 0 ) c1 = p;
213: mulnum(0,a1,b1,&t);
214: subnum(0,0,t,&p);
215: if ( compnum(0,c2,p) > 0 ) c2 = p;
216: }
217: } else {
218: if ( compnum(0,0,a1) > 0 ) {
219: mulnum(0,a1,b2,&t);
220: subnum(0,0,t,&c1);
221: } else {
222: mulnum(0,a1,b1,&c1);
223: }
224: }
225: if ( ba == bb ) {
226: subnum(0,0,c2,&t);
227: istoitv(c1,t,(Itv *)&tmp);
228: } else {
229: subnum(0,0,c1,&t);
230: istoitv(c2,t,(Itv *)&tmp);
231: }
232: addulp((IntervalBigFloat)tmp, c);
233: }
234: }
235:
236: int initvf(Num a, Itv b)
237: {
238: Num inf, sup;
239:
240: itvtois(b, &inf, &sup);
241: if ( compnum(0,inf,a) > 0 ) return 0;
242: else if ( compnum(0,a,sup) > 0 ) return 0;
243: else return 1;
244: }
245:
246: int itvinitvf(Itv a, Itv b)
247: {
248: Num ai,as,bi,bs;
249:
250: itvtois(a, &ai, &as);
251: itvtois(b, &bi, &bs);
252: if ( compnum(0,ai,bi) >= 0 && compnum(0,bs,as) >= 0 ) return 1;
253: else return 0;
254: }
255:
256: void divitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
257: {
258: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t,tmp;
259: Num inf, sup;
260: int ba,bb;
261:
262: if ( !b )
263: error("divitv : division by 0");
264: else if ( !a )
265: *c = 0;
266: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
267: divnum(0,(Num)a,(Num)b,(Num *)c);
268: else {
269: itvtois((Itv)a,&inf,&sup);
270: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
271: itvtois((Itv)b,&inf,&sup);
272: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
273: /* MUST check if ai, as, bi, bs are bigfloat. */
274: chsgnnum(as,&a2);
275: chsgnnum(bs,&b2);
276: if ( ba = (compnum(0,0,a2) > 0) ) {
277: a1 = ai;
278: } else {
279: a1 = a2;
280: a2 = ai;
281: }
282: if ( bb = (compnum(0,0,b2) > 0) ) {
283: b1 = bi;
284: } else {
285: b1 = b2;
286: b2 = bi;
287: }
288: if ( compnum(0,0,b1) >= 0 ) {
289: error("divitv : division by interval including 0.");
290: } else {
291: divnum(0,a2,b1,&c2);
292: if ( compnum(0,0,a1) > 0 ) {
293: divnum(0,a1,b1,&c1);
294: } else {
295: divnum(0,a1,b2,&t);
296: subnum(0,0,t,&c1);
297: }
298: }
299: if ( ba == bb ) {
300: subnum(0,0,c2,&t);
301: istoitv(c1,t,(Itv *)&tmp);
302: } else {
303: subnum(0,0,c1,&t);
304: istoitv(c2,t,(Itv *)&tmp);
305: }
306: addulp((IntervalBigFloat)tmp, c);
307: }
308: }
309:
310: void pwritvf(Itv a, Num e, Itv *c)
311: {
312: int ei;
313: Itv t;
314:
315: if ( !e )
316: *c = (Itv)ONE;
317: else if ( !a )
318: *c = 0;
319: else if ( NID(a) <= N_B )
320: pwrnum(0,(Num)a,e,(Num *)c);
321: else if ( !INT(e) ) {
322: #if defined(PARI) && 0
323: gpui_ri((Obj)a,(Obj)c,(Obj *)c);
324: #else
325: error("pwritv : can't calculate a fractional power");
326: #endif
327: } else {
1.2 ! noro 328: ei = ZTOS((Q)e);
1.1 noro 329: pwritv0f(a,ei,&t);
330: if ( SGN((Q)e) < 0 )
331: divbf((Num)ONE,(Num)t,(Num *)c);
332: else
333: *c = t;
334: }
335: }
336:
337: void pwritv0f(Itv a, int e, Itv *c)
338: {
339: Num inf, sup;
340: Num ai,Xmin,Xmax;
341: IntervalBigFloat tmp;
342: Q ne;
343:
344: if ( e == 1 )
345: *c = a;
346: else {
1.2 ! noro 347: STOZ(e,ne);
1.1 noro 348: if ( !(e%2) ) {
349: if ( initvp(0,a) ) {
350: Xmin = 0;
351: chsgnnum(INF(a),&ai);
352: if ( compnum(0,ai,SUP(a)) > 0 ) {
353: Xmax = ai;
354: } else {
355: Xmax = SUP(a);
356: }
357: } else {
358: if ( compnum(0,INF(a),0) > 0 ) {
359: Xmin = INF(a);
360: Xmax = SUP(a);
361: } else {
362: Xmin = SUP(a);
363: Xmax = INF(a);
364: }
365: }
366: } else {
367: Xmin = INF(a);
368: Xmax = SUP(a);
369: }
370: if ( ! Xmin ) inf = 0;
371: else pwrbf(Xmin,(Num)ne,&inf);
372: if ( ! Xmax ) sup = 0;
373: else pwrbf(Xmax,(Num)ne,&sup);
374: istoitv(inf,sup,(Itv *)&tmp);
375: addulp(tmp, (IntervalBigFloat *)c);
376: }
377: }
378:
379: void chsgnitvf(Itv a, Itv *c)
380: {
381: Num inf,sup;
382:
383: if ( !a )
384: *c = 0;
385: else if ( NID(a) <= N_B )
386: chsgnnum((Num)a,(Num *)c);
387: else {
388: chsgnnum(INF((Itv)a),&inf);
389: chsgnnum(SUP((Itv)a),&sup);
390: istoitv(inf,sup,c);
391: }
392: }
393:
394: int cmpitvf(Itv a, Itv b)
395: {
396: Num ai,as,bi,bs;
397: int s,t;
398:
399: if ( !a ) {
400: if ( !b || (NID(b)<=N_B) ) {
401: return compnum(0,(Num)a,(Num)b);
402: } else {
403: itvtois(b,&bi,&bs);
404: if ( compnum(0,(Num)a,bs) > 0 ) return 1;
405: else if ( compnum(0,bi,(Num)a) > 0 ) return -1;
406: else return 0;
407: }
408: } else if ( !b ) {
409: if ( !a || (NID(a)<=N_B) ) {
410: return compnum(0,(Num)a,(Num)b);
411: } else {
412: itvtois(a,&ai,&as);
413: if ( compnum(0,ai,(Num)b) > 0 ) return 1;
414: else if ( compnum(0,(Num)b,as) > 0 ) return -1;
415: else return 0;
416: }
417: } else {
418: itvtois(a,&ai,&as);
419: itvtois(b,&bi,&bs);
420: s = compnum(0,ai,bs) ;
421: t = compnum(0,bi,as) ;
422: if ( s > 0 ) return 1;
423: else if ( t > 0 ) return -1;
424: else return 0;
425: }
426: }
427:
428: void miditvf(Itv a, Num *b)
429: {
430: Num ai,as;
431: Num t;
432:
433: if ( ! a ) *b = 0;
434: else if ( (NID(a) <= N_B) )
435: *b = (Num)a;
436: else {
1.2 ! noro 437: STOZ(2,TWO);
1.1 noro 438: itvtois(a,&ai,&as);
439: addbf(ai,as,&t);
440: divbf(t,(Num)TWO,b);
441: }
442: }
443:
444: void cupitvf(Itv a, Itv b, Itv *c)
445: {
446: Num ai,as,bi,bs;
447: Num inf,sup;
448:
449: itvtois(a,&ai,&as);
450: itvtois(b,&bi,&bs);
451: if ( compnum(0,ai,bi) > 0 ) inf = bi;
452: else inf = ai;
453: if ( compnum(0,as,bs) > 0 ) sup = as;
454: else sup = bs;
455: istoitv(inf,sup,c);
456: }
457:
458: void capitvf(Itv a, Itv b, Itv *c)
459: {
460: Num ai,as,bi,bs;
461: Num inf,sup;
462:
463: itvtois(a,&ai,&as);
464: itvtois(b,&bi,&bs);
465: if ( compnum(0,ai,bi) > 0 ) inf = ai;
466: else inf = bi;
467: if ( compnum(0,as,bs) > 0 ) sup = bs;
468: else sup = as;
469: if ( compnum(0,inf,sup) > 0 ) *c = 0;
470: else istoitv(inf,sup,c);
471: }
472:
473:
474: void widthitvf(Itv a, Num *b)
475: {
476: Num ai,as;
477:
478: if ( ! a ) *b = 0;
479: else if ( (NID(a) <= N_B) )
480: *b = (Num)a;
481: else {
482: itvtois(a,&ai,&as);
483: subnum(0,as,ai,b);
484: }
485: }
486:
487: void absitvf(Itv a, Num *b)
488: {
489: Num ai,as,bi,bs;
490:
491: if ( ! a ) *b = 0;
492: else if ( (NID(a) <= N_B) ) {
493: if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
494: else *b = (Num)a;
495: } else {
496: itvtois(a,&ai,&as);
497: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
498: else bi = ai;
499: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
500: else bs = as;
501: if ( compnum(0,bi,bs) > 0 ) *b = bi;
502: else *b = bs;
503: }
504: }
505:
506: void distanceitvf(Itv a, Itv b, Num *c)
507: {
508: Num ai,as,bi,bs;
509: Num di,ds;
510: Itv d;
511:
512: itvtois(a,&ai,&as);
513: itvtois(b,&bi,&bs);
514: subnum(0,ai,bi,&di);
515: subnum(0,as,bs,&ds);
516: istoitv(di,ds,&d);
517: absitvf(d,c);
518: }
519:
520: #endif