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