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