=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/d-itv.c,v retrieving revision 1.7 retrieving revision 1.8 diff -u -p -r1.7 -r1.8 --- OpenXM_contrib2/asir2000/engine/d-itv.c 2016/06/29 08:16:11 1.7 +++ OpenXM_contrib2/asir2000/engine/d-itv.c 2018/03/29 01:32:51 1.8 @@ -1,5 +1,5 @@ /* - * $OpenXM: OpenXM_contrib2/asir2000/engine/d-itv.c,v 1.6 2015/08/14 13:51:54 fujimoto Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/engine/d-itv.c,v 1.7 2016/06/29 08:16:11 ohara Exp $ */ #if defined(INTERVAL) #include @@ -12,688 +12,688 @@ #if defined(ITVDEBUG) void printbinint(int d) { - int i, j, mask; - union { - int x; - char c[sizeof(int)]; - } a; + int i, j, mask; + union { + int x; + char c[sizeof(int)]; + } a; - a.x = d; - for(i=sizeof(int)-1;i>=0;i--) { - mask = 0x80; - for(j=0;j<8;j++) { - if (a.c[i] & mask) fprintf(stderr,"1"); - else fprintf(stderr,"0"); - mask >>= 1; - } - } - fprintf(stderr,"\n"); + a.x = d; + for(i=sizeof(int)-1;i>=0;i--) { + mask = 0x80; + for(j=0;j<8;j++) { + if (a.c[i] & mask) fprintf(stderr,"1"); + else fprintf(stderr,"0"); + mask >>= 1; + } + } + fprintf(stderr,"\n"); } #endif double NatToRealUp(N a, int *expo) { - int *p; - int l,all,i,j,s,tb,top,tail; - unsigned int t,m[2]; - N b,c; - int kk, carry, rem; - - p = BD(a); l = PL(a) - 1; - for ( top = 0, t = p[l]; t; t >>= 1, top++ ); - all = top + BSH*l; tail = (53-top)%BSH; j = l-(53-top)/BSH-1; + int *p; + int l,all,i,j,s,tb,top,tail; + unsigned int t,m[2]; + N b,c; + int kk, carry, rem; + + p = BD(a); l = PL(a) - 1; + for ( top = 0, t = p[l]; t; t >>= 1, top++ ); + all = top + BSH*l; tail = (53-top)%BSH; j = l-(53-top)/BSH-1; - if ( j >= 0 ) { + if ( j >= 0 ) { #if defined(ITVDEBUG) - fprintf(stderr," %d : tail = %d\n", j, tail); - printbinint(p[j]); + fprintf(stderr," %d : tail = %d\n", j, tail); + printbinint(p[j]); #endif - kk = (1<< (BSH - tail)) - 1; - rem = p[j] & kk; - if ( !rem ) - for (i=0;i>= 1, top++ ); - all = top + BSH*l; tail = (53-top)%BSH; i = l-(53-top)/BSH-1; - } else i = j; - } else i = j; + p = BD(c); l = PL(c) - 1; + for ( top = 0, t = p[l]; t; t >>= 1, top++ ); + all = top + BSH*l; tail = (53-top)%BSH; i = l-(53-top)/BSH-1; + } else i = j; + } else i = j; - m[1] = i < 0 ? 0 : p[i]>>(BSH-tail); - for ( j = 1, i++, tb = tail; i <= l; i++ ) { - s = 32-tb; t = i < 0 ? 0 : p[i]; - if ( BSH > s ) { - m[j] |= ((t&((1<>s; tb = BSH-s; - } - } else { - m[j] |= (t<>(BSH-tail); + for ( j = 1, i++, tb = tail; i <= l; i++ ) { + s = 32-tb; t = i < 0 ? 0 : p[i]; + if ( BSH > s ) { + m[j] |= ((t&((1<>s; tb = BSH-s; + } + } else { + m[j] |= (t<0)? NatToReal(NM(a),&enm): NatToRealUp(NM(a),&enm); + nm = (SGN(a)>0)? NatToReal(NM(a),&enm): NatToRealUp(NM(a),&enm); - if ( INT(a) ) - if ( enm >= 1 ) - error("Q2doubleDown : Overflow"); - else - return SGN(a)>0 ? nm : -nm; - else { - dn = (SGN(a)>0)? NatToReal(DN(a),&edn): NatToRealUp(DN(a),&edn); + if ( INT(a) ) + if ( enm >= 1 ) + error("Q2doubleDown : Overflow"); + else + return SGN(a)>0 ? nm : -nm; + else { + dn = (SGN(a)>0)? NatToReal(DN(a),&edn): NatToRealUp(DN(a),&edn); - if ( ((e = enm - edn) >= 1024) || (e <= -1023) ) - error("Q2doubleDown : Overflow"); /* XXX */ - else { - t = 0.0; p = (unsigned int *)&t; /* Machine */ - *p = ((e+1023)<<20); + if ( ((e = enm - edn) >= 1024) || (e <= -1023) ) + error("Q2doubleDown : Overflow"); /* XXX */ + else { + t = 0.0; p = (unsigned int *)&t; /* Machine */ + *p = ((e+1023)<<20); #ifdef vax - s = p[0]; p[0] = p[1]; p[1] = s; itod(p); + s = p[0]; p[0] = p[1]; p[1] = s; itod(p); #endif #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) | defined(ANDROID) - s = p[0]; p[0] = p[1]; p[1] = s; + s = p[0]; p[0] = p[1]; p[1] = s; #endif - FPMINUSINF - snm = (SGN(a)>0) ? nm : -nm; - rd = snm / dn * t; - FPNEAREST - return rd; - } - } + FPMINUSINF + snm = (SGN(a)>0) ? nm : -nm; + rd = snm / dn * t; + FPNEAREST + return rd; + } + } } -static double Q2doubleUp(Q a) +static double Q2doubleUp(Q a) { - double nm,dn,t,snm,rd; - int enm,edn,e; - unsigned int *p,s; + double nm,dn,t,snm,rd; + int enm,edn,e; + unsigned int *p,s; - nm = (SGN(a)>0)? NatToRealUp(NM(a),&enm): NatToReal(NM(a),&enm); + nm = (SGN(a)>0)? NatToRealUp(NM(a),&enm): NatToReal(NM(a),&enm); - if ( INT(a) ) - if ( enm >= 1 ) - error("Q2doubleUp : Overflow"); - else - return SGN(a)>0 ? nm : -nm; - else { - dn = (SGN(a)>0)? NatToRealUp(DN(a),&edn): NatToReal(DN(a),&edn); + if ( INT(a) ) + if ( enm >= 1 ) + error("Q2doubleUp : Overflow"); + else + return SGN(a)>0 ? nm : -nm; + else { + dn = (SGN(a)>0)? NatToRealUp(DN(a),&edn): NatToReal(DN(a),&edn); - if ( ((e = enm - edn) >= 1024) || (e <= -1023) ) - error("Q2doubleUp : Overflow"); /* XXX */ - else { - t = 0.0; p = (unsigned int *)&t; /* Machine */ - *p = ((e+1023)<<20); + if ( ((e = enm - edn) >= 1024) || (e <= -1023) ) + error("Q2doubleUp : Overflow"); /* XXX */ + else { + t = 0.0; p = (unsigned int *)&t; /* Machine */ + *p = ((e+1023)<<20); #ifdef vax - s = p[0]; p[0] = p[1]; p[1] = s; itod(p); + s = p[0]; p[0] = p[1]; p[1] = s; itod(p); #endif #if defined(MIPSEL) || defined(TOWNS) || defined(linux) || defined(__alpha) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(ANDROID) - s = p[0]; p[0] = p[1]; p[1] = s; + s = p[0]; p[0] = p[1]; p[1] = s; #endif #if 0 - FPPLUSINF - snm = (SGN(a)>0) ? nm : -nm; - rd = snm / dn * t; + FPPLUSINF + snm = (SGN(a)>0) ? nm : -nm; + rd = snm / dn * t; #endif - FPMINUSINF - snm = (SGN(a)>0) ? -nm : nm; - rd = snm / dn * t; - FPNEAREST - return (-rd); - } - } + FPMINUSINF + snm = (SGN(a)>0) ? -nm : nm; + rd = snm / dn * t; + FPNEAREST + return (-rd); + } + } } -static double PARI2doubleDown(BF a) +static double PARI2doubleDown(BF a) { GEN p; - double d; + double d; ritopa(a, &p); d = rtodbl(p); - cgiv(p); - return d; + cgiv(p); + return d; } -static double PARI2doubleUp(BF a) +static double PARI2doubleUp(BF a) { - return PARI2doubleDown(a); + return PARI2doubleDown(a); } -double subulpd(double d) +double subulpd(double d) { - double inf; + double inf; - FPMINUSINF - inf = d - DBL_MIN; - FPNEAREST - return inf; + FPMINUSINF + inf = d - DBL_MIN; + FPNEAREST + return inf; } -double addulpd(double d) +double addulpd(double d) { - double md, sup; + double md, sup; - md = -d; - FPMINUSINF - sup = md - DBL_MIN; - FPNEAREST - return (-sup); + md = -d; + FPMINUSINF + sup = md - DBL_MIN; + FPNEAREST + return (-sup); } -double ToRealDown(Num a) +double ToRealDown(Num a) { - double inf; + double inf; - if ( ! a ) return 0.0; - switch ( NID(a) ) { - case N_Q: - inf = Q2doubleDown((Q)a); break; - case N_R: - inf = subulpd(BDY((Real)a)); break; - case N_B: - inf = PARI2doubleDown((BF)a); break; - case N_IP: - inf = ToRealDown(INF((Itv)a)); - break; - case N_IntervalDouble: - inf = INF((IntervalDouble)a); break; - case N_A: - default: - inf = 0.0; - error("ToRealDown: not supported operands."); - break; - } - return inf; + if ( ! a ) return 0.0; + switch ( NID(a) ) { + case N_Q: + inf = Q2doubleDown((Q)a); break; + case N_R: + inf = subulpd(BDY((Real)a)); break; + case N_B: + inf = PARI2doubleDown((BF)a); break; + case N_IP: + inf = ToRealDown(INF((Itv)a)); + break; + case N_IntervalDouble: + inf = INF((IntervalDouble)a); break; + case N_A: + default: + inf = 0.0; + error("ToRealDown: not supported operands."); + break; + } + return inf; } -double ToRealUp(Num a) +double ToRealUp(Num a) { - double sup; + double sup; - if ( ! a ) return 0.0; - switch ( NID(a) ) { - case N_Q: - sup = Q2doubleUp((Q)a); break; - case N_R: - sup = addulpd(BDY((Real)a)); break; - case N_B: - sup = PARI2doubleUp((BF)a); break; - case N_IP: - sup = ToRealUp(SUP((Itv)a)); break; - case N_IntervalDouble: - sup = SUP((IntervalDouble)a); break; - case N_A: - default: - sup = 0.0; - error("ToRealUp: not supported operands."); - break; - } - return sup; + if ( ! a ) return 0.0; + switch ( NID(a) ) { + case N_Q: + sup = Q2doubleUp((Q)a); break; + case N_R: + sup = addulpd(BDY((Real)a)); break; + case N_B: + sup = PARI2doubleUp((BF)a); break; + case N_IP: + sup = ToRealUp(SUP((Itv)a)); break; + case N_IntervalDouble: + sup = SUP((IntervalDouble)a); break; + case N_A: + default: + sup = 0.0; + error("ToRealUp: not supported operands."); + break; + } + return sup; } -void Num2double(Num a, double *inf, double *sup) +void Num2double(Num a, double *inf, double *sup) { - switch ( NID(a) ) { - case N_Q: - *inf = Q2doubleDown((Q)a); - *sup = Q2doubleUp((Q)a); - break; - case N_R: - *inf = BDY((Real)a); - *sup = BDY((Real)a); - break; - case N_B: - *inf = PARI2doubleDown((BF)a); - *sup = PARI2doubleUp((BF)a); - break; - case N_IP: - *inf = ToRealDown(INF((Itv)a)); - *sup = ToRealUp(SUP((Itv)a)); - break; - case N_IntervalDouble: - *inf = INF((IntervalDouble)a); - *sup = SUP((IntervalDouble)a); - break; - case N_A: - default: - *inf = 0.0; - *sup = 0.0; - error("Num2double: not supported operands."); - break; - } + switch ( NID(a) ) { + case N_Q: + *inf = Q2doubleDown((Q)a); + *sup = Q2doubleUp((Q)a); + break; + case N_R: + *inf = BDY((Real)a); + *sup = BDY((Real)a); + break; + case N_B: + *inf = PARI2doubleDown((BF)a); + *sup = PARI2doubleUp((BF)a); + break; + case N_IP: + *inf = ToRealDown(INF((Itv)a)); + *sup = ToRealUp(SUP((Itv)a)); + break; + case N_IntervalDouble: + *inf = INF((IntervalDouble)a); + *sup = SUP((IntervalDouble)a); + break; + case N_A: + default: + *inf = 0.0; + *sup = 0.0; + error("Num2double: not supported operands."); + break; + } } void additvd(Num a, Num b, IntervalDouble *c) { - double ai,as,mas, bi,bs; - double inf,sup; + double ai,as,mas, bi,bs; + double inf,sup; - if ( !a ) { - *c = (IntervalDouble)b; - } else if ( !b ) { - *c = (IntervalDouble)a; -#if 0 - } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) { - addnum(0,a,b,c); + if ( !a ) { + *c = (IntervalDouble)b; + } else if ( !b ) { + *c = (IntervalDouble)a; +#if 0 + } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) { + addnum(0,a,b,c); #endif - } else { + } else { - Num2double((Num)a,&ai,&as); - Num2double((Num)b,&bi,&bs); - mas = -as; - FPMINUSINF - inf = ai + bi; - sup = mas - bs; /* as + bs = ( - ( (-as) - bs ) ) */ - FPNEAREST - MKIntervalDouble(inf,(-sup),*c); - } + Num2double((Num)a,&ai,&as); + Num2double((Num)b,&bi,&bs); + mas = -as; + FPMINUSINF + inf = ai + bi; + sup = mas - bs; /* as + bs = ( - ( (-as) - bs ) ) */ + FPNEAREST + MKIntervalDouble(inf,(-sup),*c); + } } void subitvd(Num a, Num b, IntervalDouble *c) { - double ai,as,mas, bi,bs; - double inf,sup; + double ai,as,mas, bi,bs; + double inf,sup; - if ( !a ) { - *c = (IntervalDouble)b; - } else if ( !b ) { - *c = (IntervalDouble)a; -#if 0 - } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) { - addnum(0,a,b,c); + if ( !a ) { + *c = (IntervalDouble)b; + } else if ( !b ) { + *c = (IntervalDouble)a; +#if 0 + } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) { + addnum(0,a,b,c); #endif - } else { - Num2double(a,&ai,&as); - Num2double(b,&bi,&bs); - FPMINUSINF - inf = ai - bs; - sup = ( bi - as ); /* as - bi = ( - ( bi - as ) ) */ - FPNEAREST - MKIntervalDouble(inf,(-sup),*c); - } + } else { + Num2double(a,&ai,&as); + Num2double(b,&bi,&bs); + FPMINUSINF + inf = ai - bs; + sup = ( bi - as ); /* as - bi = ( - ( bi - as ) ) */ + FPNEAREST + MKIntervalDouble(inf,(-sup),*c); + } } -void mulitvd(Num a, Num b, IntervalDouble *c) +void mulitvd(Num a, Num b, IntervalDouble *c) { - double ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p; - double inf, sup; - int ba,bb; + double ai,as,bi,bs,a1,a2,b1,b2,c1,c2,p; + double inf, sup; + int ba,bb; - if ( !a || !b ) { - *c = 0; -#if 0 - } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) { - mulnum(0,a,b,c); + if ( !a || !b ) { + *c = 0; +#if 0 + } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) { + mulnum(0,a,b,c); #endif - } else { - Num2double(a,&ai,&as); - Num2double(b,&bi,&bs); + } else { + Num2double(a,&ai,&as); + Num2double(b,&bi,&bs); - FPMINUSINF + FPMINUSINF - a2 = -as; - b2 = -bs; + a2 = -as; + b2 = -bs; - if ( ba = ( a2 < 0.0 ) ) { - a1 = ai; - } else { - a1 = a2; - a2 = ai; - } - if ( bb = ( b2 < 0.0 ) ) { - b1 = bi; - } else { - b1 = b2; - b2 = bi; - } + if ( ba = ( a2 < 0.0 ) ) { + a1 = ai; + } else { + a1 = a2; + a2 = ai; + } + if ( bb = ( b2 < 0.0 ) ) { + b1 = bi; + } else { + b1 = b2; + b2 = bi; + } - c2 = - a2 * b2; - if ( b1 < 0.0 ) { - c1 = - a2 * b1; - if ( a1 < 0.0 ) { - p = - a1 * b2; - if ( p < c1 ) c1 = p; - p = - a1 * b1; - if ( p < c2 ) c2 = p; - } - } else { - c1 = (a1<0.0)?( - a1 * b2 ):( a1 * b1 ); - } - if ( ba == bb ) { - inf = c1; - sup = - c2; - } else { - inf = c2; - sup = - c1; - } - FPNEAREST - MKIntervalDouble(inf,sup,*c); - } + c2 = - a2 * b2; + if ( b1 < 0.0 ) { + c1 = - a2 * b1; + if ( a1 < 0.0 ) { + p = - a1 * b2; + if ( p < c1 ) c1 = p; + p = - a1 * b1; + if ( p < c2 ) c2 = p; + } + } else { + c1 = (a1<0.0)?( - a1 * b2 ):( a1 * b1 ); + } + if ( ba == bb ) { + inf = c1; + sup = - c2; + } else { + inf = c2; + sup = - c1; + } + FPNEAREST + MKIntervalDouble(inf,sup,*c); + } } int initvd(Num a, IntervalDouble b) { - Real inf, sup; + Real inf, sup; - if ( NID(b) == N_IntervalDouble ) { - MKReal(INF(b), inf); - MKReal(SUP(b), sup); - } else return 0; - if ( compnum(0,(Num)inf,a) > 0 ) return 0; - else if ( compnum(0,a,(Num)sup) > 0 ) return 0; - else return 1; + if ( NID(b) == N_IntervalDouble ) { + MKReal(INF(b), inf); + MKReal(SUP(b), sup); + } else return 0; + if ( compnum(0,(Num)inf,a) > 0 ) return 0; + else if ( compnum(0,a,(Num)sup) > 0 ) return 0; + else return 1; } -void divitvd(Num a, Num b, IntervalDouble *c) +void divitvd(Num a, Num b, IntervalDouble *c) { - double ai,as,bi,bs,a1,a2,b1,b2,c1,c2; - double inf, sup; - int ba,bb; + double ai,as,bi,bs,a1,a2,b1,b2,c1,c2; + double inf, sup; + int ba,bb; - if ( !b ) { - *c = 0; - error("divitvd : division by 0"); - } else if ( !a ) { - *c = 0; -#if 0 - } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) { - divnum(0,a,b,c); + if ( !b ) { + *c = 0; + error("divitvd : division by 0"); + } else if ( !a ) { + *c = 0; +#if 0 + } else if ( (NID(a) <= N_IP) && (NID(b) <= N_IP ) ) { + divnum(0,a,b,c); #endif - } else { - Num2double(a,&ai,&as); - Num2double(b,&bi,&bs); + } else { + Num2double(a,&ai,&as); + Num2double(b,&bi,&bs); - FPMINUSINF + FPMINUSINF - a2 = -as; - b2 = -bs; + a2 = -as; + b2 = -bs; - if ( ba = ( a2 < 0.0 ) ) { - a1 = ai; - } else { - a1 = a2; - a2 = ai; - } - if ( bb = ( b2 < 0.0 ) ) { - b1 = bi; - } else { - b1 = b2; - b2 = bi; - } + if ( ba = ( a2 < 0.0 ) ) { + a1 = ai; + } else { + a1 = a2; + a2 = ai; + } + if ( bb = ( b2 < 0.0 ) ) { + b1 = bi; + } else { + b1 = b2; + b2 = bi; + } - if ( b1 <= 0.0 ) { - *c = 0; - error("divitvd : division by 0 interval"); - } else { - c2 = a2 / b1; - c1 = (a1<0.0)?( a1 / b1 ):( - a1 / b2 ); - } - if ( ba == bb ) { - inf = c1; - sup = - c2; - } else { - inf = c2; - sup = - c1; - } - FPNEAREST - MKIntervalDouble(inf,sup,*c); - } + if ( b1 <= 0.0 ) { + *c = 0; + error("divitvd : division by 0 interval"); + } else { + c2 = a2 / b1; + c1 = (a1<0.0)?( a1 / b1 ):( - a1 / b2 ); + } + if ( ba == bb ) { + inf = c1; + sup = - c2; + } else { + inf = c2; + sup = - c1; + } + FPNEAREST + MKIntervalDouble(inf,sup,*c); + } } -void pwritvd(Num a, Num e, IntervalDouble *c) +void pwritvd(Num a, Num e, IntervalDouble *c) { - int ei; - IntervalDouble t; + int ei; + IntervalDouble t; - if ( !e ) { - *c = (IntervalDouble)ONE; - } else if ( !a ) { - *c = 0; -#if 0 - } else if ( NID(a) <= N_IP ) { - pwrnum(0,a,e,c); + if ( !e ) { + *c = (IntervalDouble)ONE; + } else if ( !a ) { + *c = 0; +#if 0 + } else if ( NID(a) <= N_IP ) { + pwrnum(0,a,e,c); #endif - } else if ( !INT(e) ) { + } else if ( !INT(e) ) { #if defined(PARI) && 0 - gpui_ri((Obj)a,(Obj)c,(Obj *)c); + gpui_ri((Obj)a,(Obj)c,(Obj *)c); #else - error("pwritvd : can't calculate a fractional power"); + error("pwritvd : can't calculate a fractional power"); #endif - } else { - ei = QTOS((Q)e); - pwritv0d((IntervalDouble)a,ei,&t); - if ( SGN((Q)e) < 0 ) - divnum(0,(Num)ONE,(Num)t,(Num *)c); - else - *c = t; - } + } else { + ei = QTOS((Q)e); + pwritv0d((IntervalDouble)a,ei,&t); + if ( SGN((Q)e) < 0 ) + divnum(0,(Num)ONE,(Num)t,(Num *)c); + else + *c = t; + } } -void pwritv0d(IntervalDouble a, int e, IntervalDouble *c) +void pwritv0d(IntervalDouble a, int e, IntervalDouble *c) { - double inf, sup; - double t, Xmin, Xmax; - int i; + double inf, sup; + double t, Xmin, Xmax; + int i; - if ( e == 1 ) - *c = a; - else { - if ( !(e%2) ) { - if ( initvd(0,a) ) { - Xmin = 0.0; - t = - INF(a); - Xmax = SUP(a); - if ( t > Xmax ) Xmax = t; - } else { - if ( INF(a) > 0.0 ) { - Xmin = INF(a); - Xmax = SUP(a); - } else { - Xmin = - SUP(a); - Xmax = - INF(a); - } - } - } else { - Xmin = INF(a); - Xmax = SUP(a); - } - FPPLUSINF - sup = (Xmax!=0.0)?pwrreal0(Xmax,e):0.0; - FPMINUSINF - inf = (Xmin!=0.0)?pwrreal0(Xmin,e):0.0; - FPNEAREST - MKIntervalDouble(inf,sup,*c); - } + if ( e == 1 ) + *c = a; + else { + if ( !(e%2) ) { + if ( initvd(0,a) ) { + Xmin = 0.0; + t = - INF(a); + Xmax = SUP(a); + if ( t > Xmax ) Xmax = t; + } else { + if ( INF(a) > 0.0 ) { + Xmin = INF(a); + Xmax = SUP(a); + } else { + Xmin = - SUP(a); + Xmax = - INF(a); + } + } + } else { + Xmin = INF(a); + Xmax = SUP(a); + } + FPPLUSINF + sup = (Xmax!=0.0)?pwrreal0(Xmax,e):0.0; + FPMINUSINF + inf = (Xmin!=0.0)?pwrreal0(Xmin,e):0.0; + FPNEAREST + MKIntervalDouble(inf,sup,*c); + } } -void chsgnitvd(IntervalDouble a, IntervalDouble *c) +void chsgnitvd(IntervalDouble a, IntervalDouble *c) { - double inf,sup; + double inf,sup; - if ( !a ) - *c = 0; -#if 0 - else if ( NID(a) <= N_IP ) - chsgnnum(a,c); + if ( !a ) + *c = 0; +#if 0 + else if ( NID(a) <= N_IP ) + chsgnnum(a,c); #endif - else { - inf = - SUP(a); - sup = - INF(a); - MKIntervalDouble(inf,sup,*c); - } + else { + inf = - SUP(a); + sup = - INF(a); + MKIntervalDouble(inf,sup,*c); + } } -int cmpitvd(IntervalDouble a, IntervalDouble b) +int cmpitvd(IntervalDouble a, IntervalDouble b) /* a > b: 1 */ /* a = b: 0 */ /* a < b: -1 */ { - double ai,as,bi,bs; + double ai,as,bi,bs; - if ( !a ) { - if ( !b || (NID(b)<=N_IP) ) { - return compnum(0,(Num)a,(Num)b); - } else { - bi = INF(b); - bs = SUP(b); - if ( bi > 0.0 ) return -1; - else if ( bs < 0.0 ) return 1; - else return 0; - } - } else if ( !b ) { - if ( !a || (NID(a)<=N_IP) ) { - return compnum(0,(Num)a,(Num)b); - } else { - ai = INF(a); - as = SUP(a); - if ( ai > 0.0 ) return 1; - else if ( as < 0.0 ) return -1; - else return 0; - } - } else { - bi = INF(b); - bs = SUP(b); - ai = INF(a); - as = SUP(a); - if ( ai > bs ) return 1; - else if ( bi > as ) return -1; - else return 0; - } + if ( !a ) { + if ( !b || (NID(b)<=N_IP) ) { + return compnum(0,(Num)a,(Num)b); + } else { + bi = INF(b); + bs = SUP(b); + if ( bi > 0.0 ) return -1; + else if ( bs < 0.0 ) return 1; + else return 0; + } + } else if ( !b ) { + if ( !a || (NID(a)<=N_IP) ) { + return compnum(0,(Num)a,(Num)b); + } else { + ai = INF(a); + as = SUP(a); + if ( ai > 0.0 ) return 1; + else if ( as < 0.0 ) return -1; + else return 0; + } + } else { + bi = INF(b); + bs = SUP(b); + ai = INF(a); + as = SUP(a); + if ( ai > bs ) return 1; + else if ( bi > as ) return -1; + else return 0; + } } -void miditvd(IntervalDouble a, Num *b) +void miditvd(IntervalDouble a, Num *b) { - double t; - Real rp; + double t; + Real rp; - if ( ! a ) *b = 0; - else if ( (NID(a) <= N_IP) ) - *b = (Num)a; - else { - t = (INF(a) + SUP(a))/2.0; - MKReal(t,rp); - *b = (Num)rp; - } + if ( ! a ) *b = 0; + else if ( (NID(a) <= N_IP) ) + *b = (Num)a; + else { + t = (INF(a) + SUP(a))/2.0; + MKReal(t,rp); + *b = (Num)rp; + } } -void cupitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c) +void cupitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c) { - double ai,as,bi,bs; - double inf,sup; + double ai,as,bi,bs; + double inf,sup; - ai = INF(a); - as = SUP(a); - bi = INF(b); - bs = SUP(b); - inf = MIN(ai,bi); - sup = MAX(as,bs); - MKIntervalDouble(inf,sup,*c); + ai = INF(a); + as = SUP(a); + bi = INF(b); + bs = SUP(b); + inf = MIN(ai,bi); + sup = MAX(as,bs); + MKIntervalDouble(inf,sup,*c); } -void capitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c) +void capitvd(IntervalDouble a, IntervalDouble b, IntervalDouble *c) { - double ai,as,bi,bs; - double inf,sup; + double ai,as,bi,bs; + double inf,sup; - ai = INF(a); - as = SUP(a); - bi = INF(b); - bs = SUP(b); - inf = MAX(ai,bi); - sup = MIN(as,bs); - if ( inf > sup ) *c = 0; - else MKIntervalDouble(inf,sup,*c); + ai = INF(a); + as = SUP(a); + bi = INF(b); + bs = SUP(b); + inf = MAX(ai,bi); + sup = MIN(as,bs); + if ( inf > sup ) *c = 0; + else MKIntervalDouble(inf,sup,*c); } -void widthitvd(IntervalDouble a, Num *b) +void widthitvd(IntervalDouble a, Num *b) { - double t; - Real rp; + double t; + Real rp; - if ( ! a ) *b = 0; - else { - t = SUP(a) - INF(a); - MKReal(t,rp); - *b = (Num)rp; - } + if ( ! a ) *b = 0; + else { + t = SUP(a) - INF(a); + MKReal(t,rp); + *b = (Num)rp; + } } -void absitvd(IntervalDouble a, Num *b) +void absitvd(IntervalDouble a, Num *b) { - double ai,as,bi,bs; - double t; - Real rp; + double ai,as,bi,bs; + double t; + Real rp; - if ( ! a ) *b = 0; - else { - ai = INF(a); - as = SUP(a); - if (ai < 0) bi = -ai; else bi = ai; - if (as < 0) bs = -as; else bs = as; - if ( bi > bs ) t = bi; else t = bs; - MKReal(t,rp); - *b = (Num)rp; - } + if ( ! a ) *b = 0; + else { + ai = INF(a); + as = SUP(a); + if (ai < 0) bi = -ai; else bi = ai; + if (as < 0) bs = -as; else bs = as; + if ( bi > bs ) t = bi; else t = bs; + MKReal(t,rp); + *b = (Num)rp; + } } -void distanceitvd(IntervalDouble a, IntervalDouble b, Num *c) +void distanceitvd(IntervalDouble a, IntervalDouble b, Num *c) { - double ai,as,bi,bs; - double di,ds; - double t; - Real rp; + double ai,as,bi,bs; + double di,ds; + double t; + Real rp; - ai = INF(a); - as = SUP(a); - bi = INF(b); - bs = SUP(b); - di = bi - ai; - ds = bs - as; + ai = INF(a); + as = SUP(a); + bi = INF(b); + bs = SUP(b); + di = bi - ai; + ds = bs - as; - if (di < 0) di = -di; - if (ds < 0) ds = -ds; - if (di > ds) t = di; else t = ds; - MKReal(t,rp); - *c = (Num)rp; + if (di < 0) di = -di; + if (ds < 0) ds = -ds; + if (di > ds) t = di; else t = ds; + MKReal(t,rp); + *c = (Num)rp; } #endif