![[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.7 ! ohara 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/d-itv.c,v 1.6 2015/08/14 13:51:54 fujimoto 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
1.7 ! ohara 97: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID)
1.1 saito 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
1.7 ! ohara 127: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) | defined(ANDROID)
1.1 saito 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
1.7 ! ohara 164: #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID)
1.1 saito 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.4 ohara 496: gpui_ri((Obj)a,(Obj)c,(Obj *)c);
1.1 saito 497: #else
498: error("pwritvd : can't calculate a fractional power");
499: #endif
500: } else {
501: ei = QTOS((Q)e);
1.2 kondoh 502: pwritv0d((IntervalDouble)a,ei,&t);
1.1 saito 503: if ( SGN((Q)e) < 0 )
504: divnum(0,(Num)ONE,(Num)t,(Num *)c);
505: else
506: *c = t;
507: }
508: }
509:
1.2 kondoh 510: void pwritv0d(IntervalDouble a, int e, IntervalDouble *c)
1.1 saito 511: {
512: double inf, sup;
513: double t, Xmin, Xmax;
514: int i;
515:
516: if ( e == 1 )
517: *c = a;
518: else {
519: if ( !(e%2) ) {
520: if ( initvd(0,a) ) {
521: Xmin = 0.0;
522: t = - INF(a);
523: Xmax = SUP(a);
524: if ( t > Xmax ) Xmax = t;
525: } else {
526: if ( INF(a) > 0.0 ) {
527: Xmin = INF(a);
528: Xmax = SUP(a);
529: } else {
530: Xmin = - SUP(a);
531: Xmax = - INF(a);
532: }
533: }
534: } else {
535: Xmin = INF(a);
536: Xmax = SUP(a);
537: }
538: FPPLUSINF
539: sup = (Xmax!=0.0)?pwrreal0(Xmax,e):0.0;
540: FPMINUSINF
541: inf = (Xmin!=0.0)?pwrreal0(Xmin,e):0.0;
542: FPNEAREST
1.2 kondoh 543: MKIntervalDouble(inf,sup,*c);
1.1 saito 544: }
545: }
546:
1.2 kondoh 547: void chsgnitvd(IntervalDouble a, IntervalDouble *c)
1.1 saito 548: {
549: double inf,sup;
550:
551: if ( !a )
552: *c = 0;
553: #if 0
554: else if ( NID(a) <= N_IP )
555: chsgnnum(a,c);
556: #endif
557: else {
558: inf = - SUP(a);
559: sup = - INF(a);
1.2 kondoh 560: MKIntervalDouble(inf,sup,*c);
1.1 saito 561: }
562: }
563:
1.2 kondoh 564: int cmpitvd(IntervalDouble a, IntervalDouble b)
565: /* a > b: 1 */
566: /* a = b: 0 */
567: /* a < b: -1 */
1.1 saito 568: {
569: double ai,as,bi,bs;
570:
571: if ( !a ) {
572: if ( !b || (NID(b)<=N_IP) ) {
573: return compnum(0,(Num)a,(Num)b);
574: } else {
575: bi = INF(b);
576: bs = SUP(b);
577: if ( bi > 0.0 ) return -1;
578: else if ( bs < 0.0 ) return 1;
579: else return 0;
580: }
581: } else if ( !b ) {
582: if ( !a || (NID(a)<=N_IP) ) {
583: return compnum(0,(Num)a,(Num)b);
584: } else {
585: ai = INF(a);
586: as = SUP(a);
587: if ( ai > 0.0 ) return 1;
588: else if ( as < 0.0 ) return -1;
589: else return 0;
590: }
591: } else {
592: bi = INF(b);
593: bs = SUP(b);
594: ai = INF(a);
595: as = SUP(a);
596: if ( ai > bs ) return 1;
597: else if ( bi > as ) return -1;
598: else return 0;
599: }
600: }
601:
1.2 kondoh 602: void miditvd(IntervalDouble a, Num *b)
1.1 saito 603: {
604: double t;
605: Real rp;
606:
607: if ( ! a ) *b = 0;
608: else if ( (NID(a) <= N_IP) )
609: *b = (Num)a;
610: else {
611: t = (INF(a) + SUP(a))/2.0;
612: MKReal(t,rp);
613: *b = (Num)rp;
614: }
615: }
616:
1.2 kondoh 617: void cupitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
1.1 saito 618: {
619: double ai,as,bi,bs;
620: double inf,sup;
621:
622: ai = INF(a);
623: as = SUP(a);
624: bi = INF(b);
625: bs = SUP(b);
626: inf = MIN(ai,bi);
627: sup = MAX(as,bs);
1.2 kondoh 628: MKIntervalDouble(inf,sup,*c);
1.1 saito 629: }
630:
1.2 kondoh 631: void capitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c)
1.1 saito 632: {
633: double ai,as,bi,bs;
634: double inf,sup;
635:
636: ai = INF(a);
637: as = SUP(a);
638: bi = INF(b);
639: bs = SUP(b);
640: inf = MAX(ai,bi);
641: sup = MIN(as,bs);
642: if ( inf > sup ) *c = 0;
1.2 kondoh 643: else MKIntervalDouble(inf,sup,*c);
1.1 saito 644: }
645:
646:
1.2 kondoh 647: void widthitvd(IntervalDouble a, Num *b)
1.1 saito 648: {
649: double t;
650: Real rp;
651:
652: if ( ! a ) *b = 0;
653: else {
654: t = SUP(a) - INF(a);
655: MKReal(t,rp);
656: *b = (Num)rp;
657: }
658: }
659:
1.2 kondoh 660: void absitvd(IntervalDouble a, Num *b)
1.1 saito 661: {
662: double ai,as,bi,bs;
663: double t;
664: Real rp;
665:
666: if ( ! a ) *b = 0;
667: else {
668: ai = INF(a);
669: as = SUP(a);
670: if (ai < 0) bi = -ai; else bi = ai;
671: if (as < 0) bs = -as; else bs = as;
672: if ( bi > bs ) t = bi; else t = bs;
673: MKReal(t,rp);
674: *b = (Num)rp;
675: }
676: }
677:
1.2 kondoh 678: void distanceitvd(IntervalDouble a, IntervalDouble b, Num *c)
1.1 saito 679: {
680: double ai,as,bi,bs;
681: double di,ds;
682: double t;
683: Real rp;
684:
685: ai = INF(a);
686: as = SUP(a);
687: bi = INF(b);
688: bs = SUP(b);
689: di = bi - ai;
690: ds = bs - as;
691:
692: if (di < 0) di = -di;
693: if (ds < 0) ds = -ds;
694: if (di > ds) t = di; else t = ds;
695: MKReal(t,rp);
696: *c = (Num)rp;
697: }
698:
699: #endif