![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to d-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 d-itv.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 saito 1: /*
1.3 ! ohara 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/d-itv.c,v 1.2 2002/01/08 04:14:37 kondoh Exp $
1.1 saito 3: */
4: #if defined(INTERVAL)
5: #include <float.h>
6: #include "ca.h"
7: #include "base.h"
1.3 ! ohara 8: #if defined(PARI)
1.1 saito 9: #include "genpari.h"
10: #endif
11:
12: #if defined(ITVDEBUG)
13: void printbinint(int d)
14: {
15: int i, j, mask;
16: union {
17: int x;
18: char c[sizeof(int)];
19: } a;
20:
21: a.x = d;
22: for(i=sizeof(int)-1;i>=0;i--) {
23: mask = 0x80;
24: for(j=0;j<8;j++) {
25: if (a.c[i] & mask) fprintf(stderr,"1");
26: else fprintf(stderr,"0");
27: mask >>= 1;
28: }
29: }
30: fprintf(stderr,"\n");
31: }
32: #endif
33:
34: double NatToRealUp(N a, int *expo)
35: {
36: int *p;
37: int l,all,i,j,s,tb,top,tail;
38: unsigned int t,m[2];
39: N b,c;
40: int kk, carry, rem;
41:
42: p = BD(a); l = PL(a) - 1;
43: for ( top = 0, t = p[l]; t; t >>= 1, top++ );
44: all = top + BSH*l; tail = (53-top)%BSH; j = l-(53-top)/BSH-1;
45:
46:
47: if ( j >= 0 ) {
48:
49: #if defined(ITVDEBUG)
50: fprintf(stderr," %d : tail = %d\n", j, tail);
51: printbinint(p[j]);
52: #endif
53: kk = (1<< (BSH - tail)) - 1;
54: rem = p[j] & kk;
55: if ( !rem )
56: for (i=0;i<j;i++)
57: if ( p[j] ) {
58: carry = 1;
59: break;
60: }
61: else carry = 1;
62: if ( carry ) {
63: b = NALLOC(j+1+1);
64: PL(b) = j+1;
65: p = BD(b);
66: for(i=0;i<j;i++) p[i] = 0;
67: p[j]=(1<< (BSH - tail));
68:
69: addn(a,b,&c);
70:
71: p = BD(c); l = PL(c) - 1;
72: for ( top = 0, t = p[l]; t; t >>= 1, top++ );
73: all = top + BSH*l; tail = (53-top)%BSH; i = l-(53-top)/BSH-1;
74: } else i = j;
75: } else i = j;
76:
77:
78: m[1] = i < 0 ? 0 : p[i]>>(BSH-tail);
79: for ( j = 1, i++, tb = tail; i <= l; i++ ) {
80: s = 32-tb; t = i < 0 ? 0 : p[i];
81: if ( BSH > s ) {
82: m[j] |= ((t&((1<<s)-1))<<tb);
83: if ( !j )
84: break;
85: else {
86: j--; m[j] = t>>s; tb = BSH-s;
87: }
88: } else {
89: m[j] |= (t<<tb); tb += BSH;
90: }
91: }
92: s = (all-1)+1023;
93: m[0] = (m[0]&((1<<20)-1))|(MIN(2046,s)<<20); *expo = MAX(s-2046,0);
94: #ifdef vax
95: t = m[0]; m[0] = m[1]; m[1] = t; itod(m);
96: #endif
97: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__)
98: t = m[0]; m[0] = m[1]; m[1] = t;
99: #endif
100: return *((double *)m);
101: }
102:
103: static double Q2doubleDown(Q a)
104: {
105: double nm,dn,t,snm,rd;
106: int enm,edn,e;
107: unsigned int *p,s;
108:
109: nm = (SGN(a)>0)? NatToReal(NM(a),&enm): NatToRealUp(NM(a),&enm);
110:
111: if ( INT(a) )
112: if ( enm >= 1 )
113: error("Q2doubleDown : Overflow");
114: else
115: return SGN(a)>0 ? nm : -nm;
116: else {
117: dn = (SGN(a)>0)? NatToReal(DN(a),&edn): NatToRealUp(DN(a),&edn);
118:
119: if ( ((e = enm - edn) >= 1024) || (e <= -1023) )
120: error("Q2doubleDown : Overflow"); /* XXX */
121: else {
122: t = 0.0; p = (unsigned int *)&t; /* Machine */
123: *p = ((e+1023)<<20);
124: #ifdef vax
125: s = p[0]; p[0] = p[1]; p[1] = s; itod(p);
126: #endif
127: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__)
128: s = p[0]; p[0] = p[1]; p[1] = s;
129: #endif
130: FPMINUSINF
131: snm = (SGN(a)>0) ? nm : -nm;
132: rd = snm / dn * t;
133: FPNEAREST
134: return rd;
135: }
136: }
137: }
138:
139:
140: static double Q2doubleUp(Q a)
141: {
142: double nm,dn,t,snm,rd;
143: int enm,edn,e;
144: unsigned int *p,s;
145:
146: nm = (SGN(a)>0)? NatToRealUp(NM(a),&enm): NatToReal(NM(a),&enm);
147:
148: if ( INT(a) )
149: if ( enm >= 1 )
150: error("Q2doubleUp : Overflow");
151: else
152: return SGN(a)>0 ? nm : -nm;
153: else {
154: dn = (SGN(a)>0)? NatToRealUp(DN(a),&edn): NatToReal(DN(a),&edn);
155:
156: if ( ((e = enm - edn) >= 1024) || (e <= -1023) )
157: error("Q2doubleUp : Overflow"); /* XXX */
158: else {
159: t = 0.0; p = (unsigned int *)&t; /* Machine */
160: *p = ((e+1023)<<20);
161: #ifdef vax
162: s = p[0]; p[0] = p[1]; p[1] = s; itod(p);
163: #endif
164: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__)
165: s = p[0]; p[0] = p[1]; p[1] = s;
166: #endif
167: #if 0
168: FPPLUSINF
169: snm = (SGN(a)>0) ? nm : -nm;
170: rd = snm / dn * t;
171: #endif
172:
173: FPMINUSINF
174: snm = (SGN(a)>0) ? -nm : nm;
175: rd = snm / dn * t;
176: FPNEAREST
177: return (-rd);
178: }
179: }
180: }
181:
182: static double PARI2doubleDown(BF a)
183: {
184: GEN p;
185: double d;
186:
187: ritopa(a, &p);
188: d = rtodbl(p);
189: cgiv(p);
190: return d;
191: }
192:
193: static double PARI2doubleUp(BF a)
194: {
195: return PARI2doubleDown(a);
196: }
197:
198: double subulpd(double d)
199: {
200: double inf;
201:
202: FPMINUSINF
203: inf = d - DBL_MIN;
204: FPNEAREST
205: return inf;
206: }
207:
208: double addulpd(double d)
209: {
210: double md, sup;
211:
212: md = -d;
213: FPMINUSINF
214: sup = md - DBL_MIN;
215: FPNEAREST
216: return (-sup);
217: }
218:
219: double ToRealDown(Num a)
220: {
221: double inf;
222:
223: if ( ! a ) return 0.0;
224: switch ( NID(a) ) {
225: case N_Q:
226: inf = Q2doubleDown((Q)a); break;
227: case N_R:
228: inf = subulpd(BDY((Real)a)); break;
229: case N_B:
230: inf = PARI2doubleDown((BF)a); break;
231: case N_IP:
232: inf = ToRealDown(INF((Itv)a));
233: break;
1.2 kondoh 234: case N_IntervalDouble:
235: inf = INF((IntervalDouble)a); break;
1.1 saito 236: case N_A:
237: default:
238: inf = 0.0;
239: error("ToRealDown: not supported operands.");
240: break;
241: }
242: return inf;
243: }
244:
245: double ToRealUp(Num a)
246: {
247: double sup;
248:
249: if ( ! a ) return 0.0;
250: switch ( NID(a) ) {
251: case N_Q:
252: sup = Q2doubleUp((Q)a); break;
253: case N_R:
254: sup = addulpd(BDY((Real)a)); break;
255: case N_B:
256: sup = PARI2doubleUp((BF)a); break;
257: case N_IP:
258: sup = ToRealUp(SUP((Itv)a)); break;
1.2 kondoh 259: case N_IntervalDouble:
260: sup = SUP((IntervalDouble)a); break;
1.1 saito 261: case N_A:
262: default:
263: sup = 0.0;
264: error("ToRealUp: not supported operands.");
265: break;
266: }
267: return sup;
268: }
269:
270:
271: void Num2double(Num a, double *inf, double *sup)
272: {
273: switch ( NID(a) ) {
274: case N_Q:
275: *inf = Q2doubleDown((Q)a);
276: *sup = Q2doubleUp((Q)a);
277: break;
278: case N_R:
279: *inf = BDY((Real)a);
280: *sup = BDY((Real)a);
281: break;
282: case N_B:
283: *inf = PARI2doubleDown((BF)a);
284: *sup = PARI2doubleUp((BF)a);
285: break;
286: case N_IP:
287: *inf = ToRealDown(INF((Itv)a));
288: *sup = ToRealUp(SUP((Itv)a));
289: break;
1.2 kondoh 290: case N_IntervalDouble:
291: *inf = INF((IntervalDouble)a);
292: *sup = SUP((IntervalDouble)a);
1.1 saito 293: break;
294: case N_A:
295: default:
296: *inf = 0.0;
297: *sup = 0.0;
298: error("Num2double: not supported operands.");
299: break;
300: }
301: }
302:
303:
1.2 kondoh 304: void additvd(Num a, Num b, IntervalDouble *c)
1.1 saito 305: {
306: double ai,as,mas, bi,bs;
307: double inf,sup;
308:
309: if ( !a ) {
1.2 kondoh 310: *c = (IntervalDouble)b;
1.1 saito 311: } else if ( !b ) {
1.2 kondoh 312: *c = (IntervalDouble)a;
1.1 saito 313: #if 0
314: } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
315: addnum(0,a,b,c);
316: #endif
317: } else {
318:
319: Num2double((Num)a,&ai,&as);
320: Num2double((Num)b,&bi,&bs);
321: mas = -as;
322: FPMINUSINF
323: inf = ai + bi;
324: sup = mas - bs; /* as + bs = ( - ( (-as) - bs ) ) */
325: FPNEAREST
1.2 kondoh 326: MKIntervalDouble(inf,(-sup),*c);
1.1 saito 327: }
328: }
329:
1.2 kondoh 330: void subitvd(Num a, Num b, IntervalDouble *c)
1.1 saito 331: {
332: double ai,as,mas, bi,bs;
333: double inf,sup;
334:
335: if ( !a ) {
1.2 kondoh 336: *c = (IntervalDouble)b;
1.1 saito 337: } else if ( !b ) {
1.2 kondoh 338: *c = (IntervalDouble)a;
1.1 saito 339: #if 0
340: } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
341: addnum(0,a,b,c);
342: #endif
343: } else {
344: Num2double(a,&ai,&as);
345: Num2double(b,&bi,&bs);
346: FPMINUSINF
347: inf = ai - bs;
348: sup = ( bi - as ); /* as - bi = ( - ( bi - as ) ) */
349: FPNEAREST
1.2 kondoh 350: MKIntervalDouble(inf,(-sup),*c);
1.1 saito 351: }
352: }
353:
1.2 kondoh 354: void mulitvd(Num a, Num b, IntervalDouble *c)
1.1 saito 355: {
356: double ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p;
357: double inf, sup;
358: int ba,bb;
359:
360: if ( !a || !b ) {
361: *c = 0;
362: #if 0
363: } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
364: mulnum(0,a,b,c);
365: #endif
366: } else {
367: Num2double(a,&ai,&as);
368: Num2double(b,&bi,&bs);
369:
370: FPMINUSINF
371:
372: a2 = -as;
373: b2 = -bs;
374:
375: if ( ba = ( a2 < 0.0 ) ) {
376: a1 = ai;
377: } else {
378: a1 = a2;
379: a2 = ai;
380: }
381: if ( bb = ( b2 < 0.0 ) ) {
382: b1 = bi;
383: } else {
384: b1 = b2;
385: b2 = bi;
386: }
387:
388: c2 = - a2 * b2;
389: if ( b1 < 0.0 ) {
390: c1 = - a2 * b1;
391: if ( a1 < 0.0 ) {
392: p = - a1 * b2;
393: if ( p < c1 ) c1 = p;
394: p = - a1 * b1;
395: if ( p < c2 ) c2 = p;
396: }
397: } else {
398: c1 = (a1<0.0)?( - a1 * b2 ):( a1 * b1 );
399: }
400: if ( ba == bb ) {
401: inf = c1;
402: sup = - c2;
403: } else {
404: inf = c2;
405: sup = - c1;
406: }
407: FPNEAREST
1.2 kondoh 408: MKIntervalDouble(inf,sup,*c);
1.1 saito 409: }
410: }
411:
1.2 kondoh 412: int initvd(Num a, IntervalDouble b)
1.1 saito 413: {
414: Real inf, sup;
415:
1.2 kondoh 416: if ( NID(b) == N_IntervalDouble ) {
1.1 saito 417: MKReal(INF(b), inf);
418: MKReal(SUP(b), sup);
419: } else return 0;
420: if ( compnum(0,(Num)inf,a) > 0 ) return 0;
421: else if ( compnum(0,a,(Num)sup) > 0 ) return 0;
422: else return 1;
423: }
424:
1.2 kondoh 425: void divitvd(Num a, Num b, IntervalDouble *c)
1.1 saito 426: {
427: double ai,as,bi,bs,a1,a2,b1,b2,c1,c2;
428: double inf, sup;
429: int ba,bb;
430:
431: if ( !b ) {
432: *c = 0;
433: error("divitvd : division by 0");
434: } else if ( !a ) {
435: *c = 0;
436: #if 0
437: } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) {
438: divnum(0,a,b,c);
439: #endif
440: } else {
441: Num2double(a,&ai,&as);
442: Num2double(b,&bi,&bs);
443:
444: FPMINUSINF
445:
446: a2 = -as;
447: b2 = -bs;
448:
449: if ( ba = ( a2 < 0.0 ) ) {
450: a1 = ai;
451: } else {
452: a1 = a2;
453: a2 = ai;
454: }
455: if ( bb = ( b2 < 0.0 ) ) {
456: b1 = bi;
457: } else {
458: b1 = b2;
459: b2 = bi;
460: }
461:
462: if ( b1 <= 0.0 ) {
463: *c = 0;
464: error("divitvd : division by 0 interval");
465: } else {
466: c2 = a2 / b1;
467: c1 = (a1<0.0)?( a1 / b1 ):( - a1 / b2 );
468: }
469: if ( ba == bb ) {
470: inf = c1;
471: sup = - c2;
472: } else {
473: inf = c2;
474: sup = - c1;
475: }
476: FPNEAREST
1.2 kondoh 477: MKIntervalDouble(inf,sup,*c);
1.1 saito 478: }
479: }
480:
1.2 kondoh 481: void pwritvd(Num a, Num e, IntervalDouble *c)
1.1 saito 482: {
483: int ei;
1.2 kondoh 484: IntervalDouble t;
1.1 saito 485:
486: if ( !e ) {
1.2 kondoh 487: *c = (IntervalDouble)ONE;
1.1 saito 488: } else if ( !a ) {
489: *c = 0;
490: #if 0
491: } else if ( NID(a) <= N_IP ) {
492: pwrnum(0,a,e,c);
493: #endif
494: } else if ( !INT(e) ) {
1.3 ! ohara 495: #if defined(PARI) && 0
1.1 saito 496: GEN pa,pe,z;
497: int ltop,lbot;
498:
499: ltop = avma; ritopa(a,&pa); ritopa(e,&pe); lbot = avma;
500: z = gerepile(ltop,lbot,gpui(pa,pe,prec));
501: patori(z,c); cgiv(z);
502: #else
503: error("pwritvd : can't calculate a fractional power");
504: #endif
505: } else {
506: ei = QTOS((Q)e);
1.2 kondoh 507: pwritv0d((IntervalDouble)a,ei,&t);
1.1 saito 508: if ( SGN((Q)e) < 0 )
509: divnum(0,(Num)ONE,(Num)t,(Num *)c);
510: else
511: *c = t;
512: }
513: }
514:
1.2 kondoh 515: void pwritv0d(IntervalDouble a, int e, IntervalDouble *c)
1.1 saito 516: {
517: double inf, sup;
518: double t, Xmin, Xmax;
519: int i;
520:
521: if ( e == 1 )
522: *c = a;
523: else {
524: if ( !(e%2) ) {
525: if ( initvd(0,a) ) {
526: Xmin = 0.0;
527: t = - INF(a);
528: Xmax = SUP(a);
529: if ( t > Xmax ) Xmax = t;
530: } else {
531: if ( INF(a) > 0.0 ) {
532: Xmin = INF(a);
533: Xmax = SUP(a);
534: } else {
535: Xmin = - SUP(a);
536: Xmax = - INF(a);
537: }
538: }
539: } else {
540: Xmin = INF(a);
541: Xmax = SUP(a);
542: }
543: FPPLUSINF
544: sup = (Xmax!=0.0)?pwrreal0(Xmax,e):0.0;
545: FPMINUSINF
546: inf = (Xmin!=0.0)?pwrreal0(Xmin,e):0.0;
547: FPNEAREST
1.2 kondoh 548: MKIntervalDouble(inf,sup,*c);
1.1 saito 549: }
550: }
551:
1.2 kondoh 552: void chsgnitvd(IntervalDouble a, IntervalDouble *c)
1.1 saito 553: {
554: double inf,sup;
555:
556: if ( !a )
557: *c = 0;
558: #if 0
559: else if ( NID(a) <= N_IP )
560: chsgnnum(a,c);
561: #endif
562: else {
563: inf = - SUP(a);
564: sup = - INF(a);
1.2 kondoh 565: MKIntervalDouble(inf,sup,*c);
1.1 saito 566: }
567: }
568:
1.2 kondoh 569: int cmpitvd(IntervalDouble a, IntervalDouble b)
570: /* a > b: 1 */
571: /* a = b: 0 */
572: /* a < b: -1 */
1.1 saito 573: {
574: double ai,as,bi,bs;
575:
576: if ( !a ) {
577: if ( !b || (NID(b)<=N_IP) ) {
578: return compnum(0,(Num)a,(Num)b);
579: } else {
580: bi = INF(b);
581: bs = SUP(b);
582: if ( bi > 0.0 ) return -1;
583: else if ( bs < 0.0 ) return 1;
584: else return 0;
585: }
586: } else if ( !b ) {
587: if ( !a || (NID(a)<=N_IP) ) {
588: return compnum(0,(Num)a,(Num)b);
589: } else {
590: ai = INF(a);
591: as = SUP(a);
592: if ( ai > 0.0 ) return 1;
593: else if ( as < 0.0 ) return -1;
594: else return 0;
595: }
596: } else {
597: bi = INF(b);
598: bs = SUP(b);
599: ai = INF(a);
600: as = SUP(a);
601: if ( ai > bs ) return 1;
602: else if ( bi > as ) return -1;
603: else return 0;
604: }
605: }
606:
1.2 kondoh 607: void miditvd(IntervalDouble a, Num *b)
1.1 saito 608: {
609: double t;
610: Real rp;
611:
612: if ( ! a ) *b = 0;
613: else if ( (NID(a) <= N_IP) )
614: *b = (Num)a;
615: else {
616: t = (INF(a) + SUP(a))/2.0;
617: MKReal(t,rp);
618: *b = (Num)rp;
619: }
620: }
621:
1.2 kondoh 622: void cupitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
1.1 saito 623: {
624: double ai,as,bi,bs;
625: double inf,sup;
626:
627: ai = INF(a);
628: as = SUP(a);
629: bi = INF(b);
630: bs = SUP(b);
631: inf = MIN(ai,bi);
632: sup = MAX(as,bs);
1.2 kondoh 633: MKIntervalDouble(inf,sup,*c);
1.1 saito 634: }
635:
1.2 kondoh 636: void capitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
1.1 saito 637: {
638: double ai,as,bi,bs;
639: double inf,sup;
640:
641: ai = INF(a);
642: as = SUP(a);
643: bi = INF(b);
644: bs = SUP(b);
645: inf = MAX(ai,bi);
646: sup = MIN(as,bs);
647: if ( inf > sup ) *c = 0;
1.2 kondoh 648: else MKIntervalDouble(inf,sup,*c);
1.1 saito 649: }
650:
651:
1.2 kondoh 652: void widthitvd(IntervalDouble a, Num *b)
1.1 saito 653: {
654: double t;
655: Real rp;
656:
657: if ( ! a ) *b = 0;
658: else {
659: t = SUP(a) - INF(a);
660: MKReal(t,rp);
661: *b = (Num)rp;
662: }
663: }
664:
1.2 kondoh 665: void absitvd(IntervalDouble a, Num *b)
1.1 saito 666: {
667: double ai,as,bi,bs;
668: double t;
669: Real rp;
670:
671: if ( ! a ) *b = 0;
672: else {
673: ai = INF(a);
674: as = SUP(a);
675: if (ai < 0) bi = -ai; else bi = ai;
676: if (as < 0) bs = -as; else bs = as;
677: if ( bi > bs ) t = bi; else t = bs;
678: MKReal(t,rp);
679: *b = (Num)rp;
680: }
681: }
682:
1.2 kondoh 683: void distanceitvd(IntervalDouble a, IntervalDouble b, Num *c)
1.1 saito 684: {
685: double ai,as,bi,bs;
686: double di,ds;
687: double t;
688: Real rp;
689:
690: ai = INF(a);
691: as = SUP(a);
692: bi = INF(b);
693: bs = SUP(b);
694: di = bi - ai;
695: ds = bs - as;
696:
697: if (di < 0) di = -di;
698: if (ds < 0) ds = -ds;
699: if (di > ds) t = di; else t = ds;
700: MKReal(t,rp);
701: *c = (Num)rp;
702: }
703:
704: #endif