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