![[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.6 ! saito 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/f-itv.c,v 1.5 2003/10/20 07:18:42 saito Exp $
1.1 saito 3: */
4: #if defined(INTERVAL)
5: #include "ca.h"
6: #include "base.h"
1.3 ohara 7: #if defined(PARI)
1.1 saito 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: addnum(0,ai,bi,&inf);
139: addnum(0,as,bs,&sup);
140: istoitv(inf,sup,(Itv *)&tmp);
1.2 kondoh 141: addulp((IntervalBigFloat)tmp, c);
1.1 saito 142: return;
143: }
144: }
145:
1.2 kondoh 146: void subitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 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);
1.2 kondoh 167: addulp((IntervalBigFloat)tmp, c);
1.1 saito 168: }
169: }
170:
1.2 kondoh 171: void mulitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 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: }
1.2 kondoh 232: addulp((IntervalBigFloat)tmp, c);
1.1 saito 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:
1.2 kondoh 256: void divitvf(IntervalBigFloat a, IntervalBigFloat b, IntervalBigFloat *c)
1.1 saito 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: }
1.2 kondoh 306: addulp((IntervalBigFloat)tmp, c);
1.1 saito 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) ) {
1.3 ohara 322: #if defined(PARI) && 0
1.1 saito 323: GEN pa,pe,z;
324: int ltop,lbot;
325:
326: ltop = avma; ritopa(a,&pa); ritopa(e,&pe); lbot = avma;
327: z = gerepile(ltop,lbot,gpui(pa,pe,prec));
328: patori(z,c); cgiv(z);
329: #else
330: error("pwritv : can't calculate a fractional power");
331: #endif
332: } else {
333: ei = QTOS((Q)e);
334: pwritv0f(a,ei,&t);
335: if ( SGN((Q)e) < 0 )
336: divbf((Num)ONE,(Num)t,(Num *)c);
337: else
338: *c = t;
339: }
340: }
341:
342: void pwritv0f(Itv a, int e, Itv *c)
343: {
344: Num inf, sup;
345: Num ai,Xmin,Xmax;
1.2 kondoh 346: IntervalBigFloat tmp;
1.1 saito 347: Q ne;
348:
349: if ( e == 1 )
350: *c = a;
351: else {
352: STOQ(e,ne);
353: if ( !(e%2) ) {
354: if ( initvp(0,a) ) {
355: Xmin = 0;
356: chsgnnum(INF(a),&ai);
357: if ( compnum(0,ai,SUP(a)) > 0 ) {
358: Xmax = ai;
359: } else {
360: Xmax = SUP(a);
361: }
362: } else {
363: if ( compnum(0,INF(a),0) > 0 ) {
364: Xmin = INF(a);
365: Xmax = SUP(a);
366: } else {
367: Xmin = SUP(a);
368: Xmax = INF(a);
369: }
370: }
371: } else {
372: Xmin = INF(a);
373: Xmax = SUP(a);
374: }
375: if ( ! Xmin ) inf = 0;
376: else pwrbf(Xmin,(Num)ne,&inf);
377: if ( ! Xmax ) sup = 0;
378: else pwrbf(Xmax,(Num)ne,&sup);
379: istoitv(inf,sup,(Itv *)&tmp);
1.2 kondoh 380: addulp(tmp, (IntervalBigFloat *)c);
1.1 saito 381: }
382: }
383:
384: void chsgnitvf(Itv a, Itv *c)
385: {
386: Num inf,sup;
387:
388: if ( !a )
389: *c = 0;
390: else if ( NID(a) <= N_B )
391: chsgnnum((Num)a,(Num *)c);
392: else {
393: chsgnnum(INF((Itv)a),&inf);
394: chsgnnum(SUP((Itv)a),&sup);
395: istoitv(inf,sup,c);
396: }
397: }
398:
399: int cmpitvf(Itv a, Itv b)
400: {
401: Num ai,as,bi,bs;
402: int s,t;
403:
404: if ( !a ) {
1.4 kondoh 405: if ( !b || (NID(b)<=N_B) ) {
1.1 saito 406: return compnum(0,(Num)a,(Num)b);
1.4 kondoh 407: } else {
408: itvtois(b,&bi,&bs);
409: if ( compnum(0,(Num)a,bs) > 0 ) return 1;
410: else if ( compnum(0,bi,(Num)a) > 0 ) return -1;
411: else return 0;
412: }
1.1 saito 413: } else if ( !b ) {
1.4 kondoh 414: if ( !a || (NID(a)<=N_B) ) {
1.1 saito 415: return compnum(0,(Num)a,(Num)b);
1.4 kondoh 416: } else {
417: itvtois(a,&ai,&as);
418: if ( compnum(0,ai,(Num)b) > 0 ) return 1;
419: else if ( compnum(0,(Num)b,as) > 0 ) return -1;
420: else return 0;
421: }
1.1 saito 422: } else {
423: itvtois(a,&ai,&as);
424: itvtois(b,&bi,&bs);
1.2 kondoh 425: s = compnum(0,ai,bs) ;
426: t = compnum(0,bi,as) ;
427: if ( s > 0 ) return 1;
428: else if ( t > 0 ) return -1;
429: else return 0;
1.1 saito 430: }
431: }
432:
433: void miditvf(Itv a, Num *b)
434: {
435: Num ai,as;
436: Num t;
437:
438: if ( ! a ) *b = 0;
439: else if ( (NID(a) <= N_B) )
440: *b = (Num)a;
441: else {
442: STOQ(2,TWO);
443: itvtois(a,&ai,&as);
444: addbf(ai,as,&t);
445: divbf(t,(Num)TWO,b);
446: }
447: }
448:
449: void cupitvf(Itv a, Itv b, Itv *c)
450: {
451: Num ai,as,bi,bs;
452: Num inf,sup;
453:
454: itvtois(a,&ai,&as);
455: itvtois(b,&bi,&bs);
456: if ( compnum(0,ai,bi) > 0 ) inf = bi;
457: else inf = ai;
458: if ( compnum(0,as,bs) > 0 ) sup = as;
459: else sup = bs;
460: istoitv(inf,sup,c);
461: }
462:
463: void capitvf(Itv a, Itv b, Itv *c)
464: {
465: Num ai,as,bi,bs;
466: Num inf,sup;
467:
468: itvtois(a,&ai,&as);
469: itvtois(b,&bi,&bs);
470: if ( compnum(0,ai,bi) > 0 ) inf = ai;
471: else inf = bi;
472: if ( compnum(0,as,bs) > 0 ) sup = bs;
473: else sup = as;
474: if ( compnum(0,inf,sup) > 0 ) *c = 0;
475: else istoitv(inf,sup,c);
476: }
477:
478:
479: void widthitvf(Itv a, Num *b)
480: {
481: Num ai,as;
482:
483: if ( ! a ) *b = 0;
484: else if ( (NID(a) <= N_B) )
485: *b = (Num)a;
486: else {
487: itvtois(a,&ai,&as);
488: subnum(0,as,ai,b);
489: }
490: }
491:
492: void absitvf(Itv a, Num *b)
493: {
494: Num ai,as,bi,bs;
495:
496: if ( ! a ) *b = 0;
497: else if ( (NID(a) <= N_B) ) {
498: if ( compnum(0,(Num)a,0) < 0 ) chsgnnum((Num)a,b);
499: else *b = (Num)a;
500: } else {
501: itvtois(a,&ai,&as);
502: if ( compnum(0,ai,0) < 0 ) chsgnnum(ai,&bi);
503: else bi = ai;
504: if ( compnum(0,as,0) < 0 ) chsgnnum(as,&bs);
505: else bs = as;
506: if ( compnum(0,bi,bs) > 0 ) *b = bi;
507: else *b = bs;
508: }
509: }
510:
511: void distanceitvf(Itv a, Itv b, Num *c)
512: {
513: Num ai,as,bi,bs;
514: Num di,ds;
515: Itv d;
516:
517: itvtois(a,&ai,&as);
518: itvtois(b,&bi,&bs);
519: subnum(0,ai,bi,&di);
520: subnum(0,as,bs,&ds);
521: istoitv(di,ds,&d);
522: absitvf(d,c);
523: }
524:
525: #endif