![[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 / asir2000 / engine|
Return to f-itv.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 saito 1: /*
1.2 ! kondoh 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/f-itv.c,v 1.1 2000/12/22 10:03:28 saito Exp $
1.1 saito 3: */
4: #if defined(INTERVAL)
5: #include "ca.h"
6: #include "base.h"
7: #if PARI
8: #include "genpari.h"
9: #include "itv-pari.h"
10: extern long prec;
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(prec);
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(prec);
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:
1.2 ! kondoh 97: void addulp(IntervalBigFloat a, IntervalBigFloat *c)
1.1 saito 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:
1.2 ! kondoh 113: void additvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 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: #if 1
139: addnum(0,ai,bi,&inf);
140: addnum(0,as,bs,&sup);
141: istoitv(inf,sup,(Itv *)&tmp);
1.2 ! kondoh 142: addulp((IntervalBigFloat)tmp, c);
1.1 saito 143: return;
144: #else
145: ltop = avma;
146: ritopa(ai,&pa);
147: ritopa(bi,&pb);
148: lbot = avma;
149: z = gerepile(ltop,lbot,PariAddDown(pa,pb));
150: patori(z,&inf); cgiv(z);
151:
152: /* MUST check if ai, as, bi, bs are bigfloat. */
153:
154: /* as + bs = ( - ( (-as) + (-bs) ) ) */
155: chsgnbf(as,&mas);
156: chsgnbf(bs,&mbs);
157: ltop = avma;
158: ritopa(mas,&pa);
159: ritopa(mbs,&pb);
160: lbot = avma;
161: z = gerepile(ltop,lbot,PariAddDown(pa,pb));
162: patori(z,&tmp); cgiv(z);
163:
164: chsgnbf(tmp,&sup);
165: istoitv(inf,sup,c);
166: #endif
167: }
168: }
169:
1.2 ! kondoh 170: void subitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 171: {
172: Num ai,as,bi,bs,mas, mbs;
173: Num inf,sup,tmp;
174: GEN pa, pb, z;
175: long ltop, lbot;
176:
177: if ( !a )
178: chsgnnum((Num)b,(Num *)c);
179: else if ( !b )
180: *c = a;
181: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
182: subnum(0,(Num)a,(Num)b,(Num *)c);
183: else {
184: itvtois((Itv)a,&inf,&sup);
185: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
186: itvtois((Itv)b,&inf,&sup);
187: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
188: #if 1
189: subnum(0,ai,bs,&inf);
190: subnum(0,as,bi,&sup);
191: istoitv(inf,sup,(Itv *)&tmp);
1.2 ! kondoh 192: addulp((IntervalBigFloat)tmp, c);
1.1 saito 193: #else
194:
195: /* MUST check if ai, as, bi, bs are bigfloat. */
196: /* ai - bs = ai + (-bs) */
197: chsgnbf(bs,&mbs);
198: ltop = avma;
199: ritopa(ai,&pa);
200: ritopa(mbs,&pb);
201: lbot = avma;
202: z = gerepile(ltop,lbot,PariAddDown(pa,pb));
203: patori(z,&inf); cgiv(z);
204:
205: /* as - bi = ( - ( bi + (-as) ) ) */
206: chsgnbf(as,&mas);
207: ltop = avma;
208: ritopa(mas,&pa);
209: ritopa(bi,&pb);
210: lbot = avma;
211: z = gerepile(ltop,lbot,PariAddDown(pa,pb));
212: patori(z,&tmp); cgiv(z);
213:
214: chsgnbf(tmp,&sup);
215: istoitv(inf,sup,c);
216: #endif
217: }
218: }
219:
1.2 ! kondoh 220: void mulitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 221: {
222: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p,t,tmp;
223: Num inf, sup;
224: int ba,bb;
225:
226: if ( !a || !b )
227: *c = 0;
228: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
229: mulnum(0,(Num)a,(Num)b,(Num *)c);
230: else {
231: itvtois((Itv)a,&inf,&sup);
232: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
233: itvtois((Itv)b,&inf,&sup);
234: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
235:
236: /* MUST check if ai, as, bi, bs are bigfloat. */
237: chsgnnum(as,&a2);
238: chsgnnum(bs,&b2);
239:
240:
241: if ( ba = (compnum(0,0,a2) > 0) ) {
242: a1 = ai;
243: } else {
244: a1 = a2;
245: a2 = ai;
246: }
247: if ( bb = (compnum(0,0,b2) > 0) ) {
248: b1 = bi;
249: } else {
250: b1 = b2;
251: b2 = bi;
252: }
253: mulnum(0,a2,b2,&t);
254: subnum(0,0,t,&c2);
255: if ( compnum(0,0,b1) > 0 ) {
256: mulnum(0,a2,b1,&t);
257: subnum(0,0,t,&c1);
258: if ( compnum(0,0,a1) > 0 ) {
259: mulnum(0,a1,b2,&t);
260: subnum(0,0,t,&p);
261: if ( compnum(0,c1,p) > 0 ) c1 = p;
262: mulnum(0,a1,b1,&t);
263: subnum(0,0,t,&p);
264: if ( compnum(0,c2,p) > 0 ) c2 = p;
265: }
266: } else {
267: if ( compnum(0,0,a1) > 0 ) {
268: mulnum(0,a1,b2,&t);
269: subnum(0,0,t,&c1);
270: } else {
271: mulnum(0,a1,b1,&c1);
272: }
273: }
274: if ( ba == bb ) {
275: subnum(0,0,c2,&t);
276: istoitv(c1,t,(Itv *)&tmp);
277: } else {
278: subnum(0,0,c1,&t);
279: istoitv(c2,t,(Itv *)&tmp);
280: }
1.2 ! kondoh 281: addulp((IntervalBigFloat)tmp, c);
1.1 saito 282: }
283: }
284:
285: int initvf(Num a, Itv b)
286: {
287: Num inf, sup;
288:
289: itvtois(b, &inf, &sup);
290: if ( compnum(0,inf,a) > 0 ) return 0;
291: else if ( compnum(0,a,sup) > 0 ) return 0;
292: else return 1;
293: }
294:
295: int itvinitvf(Itv a, Itv b)
296: {
297: Num ai,as,bi,bs;
298:
299: itvtois(a, &ai, &as);
300: itvtois(b, &bi, &bs);
301: if ( compnum(0,ai,bi) >= 0 && compnum(0,bs,as) >= 0 ) return 1;
302: else return 0;
303: }
304:
1.2 ! kondoh 305: void divitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 306: {
307: Num ai,as,bi,bs,a1,a2,b1,b2,c1,c2,t,tmp;
308: Num inf, sup;
309: int ba,bb;
310:
311: if ( !b )
312: error("divitv : division by 0");
313: else if ( !a )
314: *c = 0;
315: else if ( (NID(a) <= N_B) && (NID(b) <= N_B ) )
316: divnum(0,(Num)a,(Num)b,(Num *)c);
317: else {
318: itvtois((Itv)a,&inf,&sup);
319: ToBf(inf, (BF *)&ai); ToBf(sup, (BF *)&as);
320: itvtois((Itv)b,&inf,&sup);
321: ToBf(inf, (BF *)&bi); ToBf(sup, (BF *)&bs);
322: /* MUST check if ai, as, bi, bs are bigfloat. */
323: chsgnnum(as,&a2);
324: chsgnnum(bs,&b2);
325: if ( ba = (compnum(0,0,a2) > 0) ) {
326: a1 = ai;
327: } else {
328: a1 = a2;
329: a2 = ai;
330: }
331: if ( bb = (compnum(0,0,b2) > 0) ) {
332: b1 = bi;
333: } else {
334: b1 = b2;
335: b2 = bi;
336: }
337: if ( compnum(0,0,b1) >= 0 ) {
338: error("divitv : division by interval including 0.");
339: } else {
340: divnum(0,a2,b1,&c2);
341: if ( compnum(0,0,a1) > 0 ) {
342: divnum(0,a1,b1,&c1);
343: } else {
344: divnum(0,a1,b2,&t);
345: subnum(0,0,t,&c1);
346: }
347: }
348: if ( ba == bb ) {
349: subnum(0,0,c2,&t);
350: istoitv(c1,t,(Itv *)&tmp);
351: } else {
352: subnum(0,0,c1,&t);
353: istoitv(c2,t,(Itv *)&tmp);
354: }
1.2 ! kondoh 355: addulp((IntervalBigFloat)tmp, c);
1.1 saito 356: }
357: }
358:
359: void pwritvf(Itv a, Num e, Itv *c)
360: {
361: int ei;
362: Itv t;
363:
364: if ( !e )
365: *c = (Itv)ONE;
366: else if ( !a )
367: *c = 0;
368: else if ( NID(a) <= N_B )
369: pwrnum(0,(Num)a,e,(Num *)c);
370: else if ( !INT(e) ) {
371: #if PARI && 0
372: GEN pa,pe,z;
373: int ltop,lbot;
374:
375: ltop = avma; ritopa(a,&pa); ritopa(e,&pe); lbot = avma;
376: z = gerepile(ltop,lbot,gpui(pa,pe,prec));
377: patori(z,c); cgiv(z);
378: #else
379: error("pwritv : can't calculate a fractional power");
380: #endif
381: } else {
382: ei = QTOS((Q)e);
383: pwritv0f(a,ei,&t);
384: if ( SGN((Q)e) < 0 )
385: divbf((Num)ONE,(Num)t,(Num *)c);
386: else
387: *c = t;
388: }
389: }
390:
391: void pwritv0f(Itv a, int e, Itv *c)
392: {
393: Num inf, sup;
394: Num ai,Xmin,Xmax;
1.2 ! kondoh 395: IntervalBigFloat tmp;
1.1 saito 396: Q ne;
397:
398: if ( e == 1 )
399: *c = a;
400: else {
401: STOQ(e,ne);
402: if ( !(e%2) ) {
403: if ( initvp(0,a) ) {
404: Xmin = 0;
405: chsgnnum(INF(a),&ai);
406: if ( compnum(0,ai,SUP(a)) > 0 ) {
407: Xmax = ai;
408: } else {
409: Xmax = SUP(a);
410: }
411: } else {
412: if ( compnum(0,INF(a),0) > 0 ) {
413: Xmin = INF(a);
414: Xmax = SUP(a);
415: } else {
416: Xmin = SUP(a);
417: Xmax = INF(a);
418: }
419: }
420: } else {
421: Xmin = INF(a);
422: Xmax = SUP(a);
423: }
424: if ( ! Xmin ) inf = 0;
425: else pwrbf(Xmin,(Num)ne,&inf);
426: if ( ! Xmax ) sup = 0;
427: else pwrbf(Xmax,(Num)ne,&sup);
428: istoitv(inf,sup,(Itv *)&tmp);
1.2 ! kondoh 429: addulp(tmp, (IntervalBigFloat *)c);
1.1 saito 430: }
431: }
432:
433: void chsgnitvf(Itv a, Itv *c)
434: {
435: Num inf,sup;
436:
437: if ( !a )
438: *c = 0;
439: else if ( NID(a) <= N_B )
440: chsgnnum((Num)a,(Num *)c);
441: else {
442: chsgnnum(INF((Itv)a),&inf);
443: chsgnnum(SUP((Itv)a),&sup);
444: istoitv(inf,sup,c);
445: }
446: }
447:
448: int cmpitvf(Itv a, Itv b)
449: {
450: Num ai,as,bi,bs;
451: int s,t;
452:
453: if ( !a ) {
454: if ( !b || (NID(b)<=N_B) )
455: return compnum(0,(Num)a,(Num)b);
456: else
457: return -1;
458: } else if ( !b ) {
459: if ( !a || (NID(a)<=N_B) )
460: return compnum(0,(Num)a,(Num)b);
461: else
462: return initvp((Num)b,a);
463: } else {
464: itvtois(a,&ai,&as);
465: itvtois(b,&bi,&bs);
1.2 ! kondoh 466: s = compnum(0,ai,bs) ;
! 467: t = compnum(0,bi,as) ;
! 468: if ( s > 0 ) return 1;
! 469: else if ( t > 0 ) return -1;
! 470: else return 0;
1.1 saito 471: }
472: }
473:
474: void miditvf(Itv a, Num *b)
475: {
476: Num ai,as;
477: Num t;
478: Q TWO;
479:
480: if ( ! a ) *b = 0;
481: else if ( (NID(a) <= N_B) )
482: *b = (Num)a;
483: else {
484: STOQ(2,TWO);
485: itvtois(a,&ai,&as);
486: addbf(ai,as,&t);
487: divbf(t,(Num)TWO,b);
488: }
489: }
490:
491: void cupitvf(Itv a, Itv b, Itv *c)
492: {
493: Num ai,as,bi,bs;
494: Num inf,sup;
495:
496: itvtois(a,&ai,&as);
497: itvtois(b,&bi,&bs);
498: if ( compnum(0,ai,bi) > 0 ) inf = bi;
499: else inf = ai;
500: if ( compnum(0,as,bs) > 0 ) sup = as;
501: else sup = bs;
502: istoitv(inf,sup,c);
503: }
504:
505: void capitvf(Itv a, Itv b, Itv *c)
506: {
507: Num ai,as,bi,bs;
508: Num inf,sup;
509:
510: itvtois(a,&ai,&as);
511: itvtois(b,&bi,&bs);
512: if ( compnum(0,ai,bi) > 0 ) inf = ai;
513: else inf = bi;
514: if ( compnum(0,as,bs) > 0 ) sup = bs;
515: else sup = as;
516: if ( compnum(0,inf,sup) > 0 ) *c = 0;
517: else istoitv(inf,sup,c);
518: }
519:
520:
521: void widthitvf(Itv a, Num *b)
522: {
523: Num ai,as;
524:
525: if ( ! a ) *b = 0;
526: else if ( (NID(a) <= N_B) )
527: *b = (Num)a;
528: else {
529: itvtois(a,&ai,&as);
530: subnum(0,as,ai,b);
531: }
532: }
533:
534: void absitvf(Itv a, Num *b)
535: {
536: Num ai,as,bi,bs;
537:
538: if ( ! a ) *b = 0;
539: else if ( (NID(a) <= N_B) ) {
540: if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
541: else *b = (Num)a;
542: } else {
543: itvtois(a,&ai,&as);
544: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
545: else bi = ai;
546: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
547: else bs = as;
548: if ( compnum(0,bi,bs) > 0 ) *b = bi;
549: else *b = bs;
550: }
551: }
552:
553: void distanceitvf(Itv a, Itv b, Num *c)
554: {
555: Num ai,as,bi,bs;
556: Num di,ds;
557: Itv d;
558:
559: itvtois(a,&ai,&as);
560: itvtois(b,&bi,&bs);
561: subnum(0,ai,bi,&di);
562: subnum(0,as,bs,&ds);
563: istoitv(di,ds,&d);
564: absitvf(d,c);
565: }
566:
567: #endif