diff options
Diffstat (limited to 'sysdeps/ieee754/support.c')
| -rw-r--r-- | sysdeps/ieee754/support.c | 110 |
1 files changed, 53 insertions, 57 deletions
diff --git a/sysdeps/ieee754/support.c b/sysdeps/ieee754/support.c index e976839421..2651ea8be4 100644 --- a/sysdeps/ieee754/support.c +++ b/sysdeps/ieee754/support.c @@ -10,10 +10,6 @@ * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. * 4. Neither the name of the University nor the names of its contributors * may be used to endorse or promote products derived from this software * without specific prior written permission. @@ -35,37 +31,37 @@ static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 6/4/93"; #endif /* not lint */ -/* - * Some IEEE standard 754 recommended functions and remainder and sqrt for +/* + * Some IEEE standard 754 recommended functions and remainder and sqrt for * supporting the C elementary functions. ****************************************************************************** * WARNING: * These codes are developed (in double) to support the C elementary * functions temporarily. They are not universal, and some of them are very - * slow (in particular, drem and sqrt is extremely inefficient). Each - * computer system should have its implementation of these functions using + * slow (in particular, drem and sqrt is extremely inefficient). Each + * computer system should have its implementation of these functions using * its own assembler. ****************************************************************************** * * IEEE 754 required operations: - * drem(x,p) + * drem(x,p) * returns x REM y = x - [x/y]*y , where [x/y] is the integer * nearest x/y; in half way case, choose the even one. - * sqrt(x) - * returns the square root of x correctly rounded according to + * sqrt(x) + * returns the square root of x correctly rounded according to * the rounding mod. * * IEEE 754 recommended functions: - * (a) copysign(x,y) - * returns x with the sign of y. - * (b) scalb(x,N) + * (a) copysign(x,y) + * returns x with the sign of y. + * (b) scalb(x,N) * returns x * (2**N), for integer values N. - * (c) logb(x) - * returns the unbiased exponent of x, a signed integer in - * double precision, except that logb(0) is -INF, logb(INF) + * (c) logb(x) + * returns the unbiased exponent of x, a signed integer in + * double precision, except that logb(0) is -INF, logb(INF) * is +INF, and logb(NAN) is that NAN. - * (d) finite(x) - * returns the value TRUE if -INF < x < +INF and returns + * (d) finite(x) + * returns the value TRUE if -INF < x < +INF and returns * FALSE otherwise. * * @@ -78,7 +74,7 @@ static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 6/4/93"; #if defined(vax)||defined(tahoe) /* VAX D format */ #include <errno.h> static const unsigned short msign=0x7fff , mexp =0x7f80 ; - static const short prep1=57, gap=7, bias=129 ; + static const short prep1=57, gap=7, bias=129 ; static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ; #else /* defined(vax)||defined(tahoe) */ static const unsigned short msign=0x7fff, mexp =0x7ff0 ; @@ -97,7 +93,7 @@ double x; int N; unsigned short *px=(unsigned short *) &x; #endif /* national */ - if( x == zero ) return(x); + if( x == zero ) return(x); #if defined(vax)||defined(tahoe) if( (k= *px & mexp ) != ~msign ) { @@ -117,7 +113,7 @@ double x; int N; if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap); else x=novf+novf; /* overflow */ else - if( k > -prep1 ) + if( k > -prep1 ) /* gradual underflow */ {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);} else @@ -147,7 +143,7 @@ double x,y; } double logb(x) -double x; +double x; { #ifdef national @@ -164,8 +160,8 @@ double x; return ( (k>>gap) - bias ); else if( x != zero) return ( -1022.0 ); - else - return(-(1.0/zero)); + else + return(-(1.0/zero)); else if(x != x) return(x); else @@ -174,7 +170,7 @@ double x; } finite(x) -double x; +double x; { #if defined(vax)||defined(tahoe) return(1); @@ -192,16 +188,16 @@ double x,p; { short sign; double hp,dp,tmp; - unsigned short k; + unsigned short k; #ifdef national unsigned short - *px=(unsigned short *) &x +3, + *px=(unsigned short *) &x +3, *pp=(unsigned short *) &p +3, *pd=(unsigned short *) &dp +3, *pt=(unsigned short *) &tmp+3; #else /* national */ unsigned short - *px=(unsigned short *) &x , + *px=(unsigned short *) &x , *pp=(unsigned short *) &p , *pd=(unsigned short *) &dp , *pt=(unsigned short *) &tmp; @@ -230,13 +226,13 @@ double x,p; #endif /* defined(vax)||defined(tahoe) */ { if (p != p) return p; else return x;} - else if ( ((*pp & mexp)>>gap) <= 1 ) + else if ( ((*pp & mexp)>>gap) <= 1 ) /* subnormal p, or almost subnormal p */ { double b; b=scalb(1.0,(int)prep1); p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);} else if ( p >= novf/2) { p /= 2 ; x /= 2; return(drem(x,p)*2);} - else + else { dp=p+p; hp=p/2; sign= *px & ~msign ; @@ -294,13 +290,13 @@ double x; } /* sqrt(INF) is INF */ - if(!finite(x)) return(x); + if(!finite(x)) return(x); /* scale x to [1,4) */ n=logb(x); x=scalb(x,-n); if((m=logb(x))!=0) x=scalb(x,-m); /* subnormal number */ - m += n; + m += n; n = m/2; if((n+n)!=m) {x *= 2; m -=1; n=m/2;} @@ -313,23 +309,23 @@ double x; else s *= 2; } - + /* generate the last bit and determine the final rounding */ - r/=2; x *= 4; + r/=2; x *= 4; if(x==zero) goto end; 100+r; /* trigger inexact flag */ if(s<x) { q+=r; x -=s; s += 2; s *= 2; x *= 4; - t = (x-s)-5; + t = (x-s)-5; b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */ b=1.0+r/4; if(b>1.0) t=1; /* b>1 : Round-to-(+INF) */ if(t>=0) q+=r; } /* else: Round-to-nearest */ - else { - s *= 2; x *= 4; - t = (x-s)-1; + else { + s *= 2; x *= 4; + t = (x-s)-1; b=1.0+3*r/4; if(b==1.0) goto end; b=1.0+r/4; if(b>1.0) t=1; if(t>=0) q+=r; } - + end: return(scalb(q,n)); } @@ -342,13 +338,13 @@ end: return(scalb(q,n)); * * Warning: this code should not get compiled in unless ALL of * the following machine-dependent routines are supplied. - * + * * Required machine dependent functions (not on a VAX): * swapINX(i): save inexact flag and reset it to "i" * swapENI(e): save inexact enable and reset it to "e" */ -double drem(x,y) +double drem(x,y) double x,y; { @@ -363,8 +359,8 @@ double x,y; double hy,y1,t,t1; short k; long n; - int i,e; - unsigned short xexp,yexp, *px =(unsigned short *) &x , + int i,e; + unsigned short xexp,yexp, *px =(unsigned short *) &x , nx,nf, *py =(unsigned short *) &y , sign, *pt =(unsigned short *) &t , *pt1 =(unsigned short *) &t1 ; @@ -381,7 +377,7 @@ double x,y; /* save the inexact flag and inexact enable in i and e respectively * and reset them to zero */ - i=swapINX(0); e=swapENI(0); + i=swapINX(0); e=swapENI(0); /* subnormal number */ nx=0; @@ -391,7 +387,7 @@ double x,y; if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;} nf=nx; - py[n0] &= 0x7fff; + py[n0] &= 0x7fff; px[n0] &= 0x7fff; /* mask off the least significant 27 bits of y */ @@ -408,7 +404,7 @@ loop: if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */ {pt[n0]+=k;pt1[n0]+=k;} n=x/t; x=(x-n*t1)-n*(t-t1); - } + } /* end while (x > y) */ if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;} @@ -416,14 +412,14 @@ loop: /* final adjustment */ hy=y/2.0; - if(x>hy||((x==hy)&&n%2==1)) x-=y; + if(x>hy||((x==hy)&&n%2==1)) x-=y; px[n0] ^= sign; if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;} /* restore inexact flag and inexact enable */ - swapINX(i); swapENI(e); + swapINX(i); swapENI(e); - return(x); + return(x); } #endif @@ -435,7 +431,7 @@ loop: * * Warning: this code should not get compiled in unless ALL of * the following machine-dependent routines are supplied. - * + * * Required machine dependent functions: * swapINX(i) ...return the status of INEXACT flag and reset it to "i" * swapRM(r) ...return the current Rounding Mode and reset it to "r" @@ -456,16 +452,16 @@ double x; double const b54=134217728.*134217728.; /* b54=2**54 */ long mx,scalx; long const mexp=0x7ff00000; - int i,j,r,e,swapINX(),swapRM(),swapENI(); + int i,j,r,e,swapINX(),swapRM(),swapENI(); unsigned long *py=(unsigned long *) &y , *pt=(unsigned long *) &t , *px=(unsigned long *) &x ; #ifdef national /* ordering of word in a floating point number */ - const int n0=1, n1=0; + const int n0=1, n1=0; #else - const int n0=0, n1=1; + const int n0=0, n1=1; #endif -/* Rounding Mode: RN ...round-to-nearest +/* Rounding Mode: RN ...round-to-nearest * RZ ...round-towards 0 * RP ...round-towards +INF * RM ...round-towards -INF @@ -502,10 +498,10 @@ double x; t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */ - t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; + t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; -/* twiddle last bit to force y correctly rounded */ +/* twiddle last bit to force y correctly rounded */ swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */ swapINX(0); /* ...clear INEXACT flag */ swapENI(e); /* ...restore inexact enable status */ |