diff options
author | rgrimes <rgrimes@FreeBSD.org> | 1995-05-30 05:51:47 +0000 |
---|---|---|
committer | rgrimes <rgrimes@FreeBSD.org> | 1995-05-30 05:51:47 +0000 |
commit | f05428e4cd63dde97bac14b84dd146a5c00455e3 (patch) | |
tree | e1331adb5d216f2b3fa6baa6491752348d2e5f10 /lib/libm/ieee | |
parent | 6de57e42c294763c78d77b0a9a7c5a08008a378a (diff) | |
download | FreeBSD-src-f05428e4cd63dde97bac14b84dd146a5c00455e3.zip FreeBSD-src-f05428e4cd63dde97bac14b84dd146a5c00455e3.tar.gz |
Remove trailing whitespace.
Diffstat (limited to 'lib/libm/ieee')
-rw-r--r-- | lib/libm/ieee/cabs.c | 26 | ||||
-rw-r--r-- | lib/libm/ieee/cbrt.c | 12 | ||||
-rw-r--r-- | lib/libm/ieee/support.c | 106 |
3 files changed, 72 insertions, 72 deletions
diff --git a/lib/libm/ieee/cabs.c b/lib/libm/ieee/cabs.c index c4a6fb0..5f1fc62 100644 --- a/lib/libm/ieee/cabs.c +++ b/lib/libm/ieee/cabs.c @@ -38,7 +38,7 @@ static char sccsid[] = "@(#)cabs.c 8.1 (Berkeley) 6/4/93"; /* HYPOT(X,Y) * RETURN THE SQUARE ROOT OF X^2 + Y^2 WHERE Z=X+iY * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 11/28/84; + * CODED IN C BY K.C. NG, 11/28/84; * REVISED BY K.C. NG, 7/12/85. * * Required system supported functions : @@ -52,16 +52,16 @@ static char sccsid[] = "@(#)cabs.c 8.1 (Berkeley) 6/4/93"; * y if y > x (hence x is never smaller than y). * 2. Hypot(x,y) is computed by: * Case I, x/y > 2 - * + * * y * hypot = x + ----------------------------- * 2 * sqrt ( 1 + [x/y] ) + x/y * - * Case II, x/y <= 2 + * Case II, x/y <= 2 * y * hypot = x + -------------------------------------------------- - * 2 + * 2 * [x/y] - 2 * (sqrt(2)+1) + (x-y)/y + ----------------------------- * 2 @@ -107,7 +107,7 @@ double hypot(x,y) double x, y; { - static const double zero=0, one=1, + static const double zero=0, one=1, small=1.0E-18; /* fl(1+small)==1 */ static const ibig=30; /* fl(1+2**(2*ibig))==1 */ double t,r; @@ -115,15 +115,15 @@ double x, y; if(finite(x)) if(finite(y)) - { + { x=copysign(x,one); y=copysign(y,one); - if(y > x) + if(y > x) { t=x; x=y; y=t; } if(x == zero) return(zero); if(y == zero) return(x); exp= logb(x); - if(exp-(int)logb(y) > ibig ) + if(exp-(int)logb(y) > ibig ) /* raise inexact flag and return |x| */ { one+small; return(x); } @@ -144,7 +144,7 @@ double x, y; else if(y==y) /* y is +-INF */ return(copysign(y,one)); - else + else return(y); /* y is NaN and x is finite */ else if(x==x) /* x is +-INF */ @@ -199,16 +199,16 @@ double x, y; if(finite(x)) if(finite(y)) - { + { x=copysign(x,one); y=copysign(y,one); - if(y > x) + if(y > x) { temp=x; x=y; y=temp; } if(x == zero) return(zero); if(y == zero) return(x); exp= logb(x); x=scalb(x,-exp); - if(exp-(int)logb(y) > ibig ) + if(exp-(int)logb(y) > ibig ) /* raise inexact flag and return |x| */ { one+small; return(scalb(x,exp)); } else y=scalb(y,-exp); @@ -217,7 +217,7 @@ double x, y; else if(y==y) /* y is +-INF */ return(copysign(y,one)); - else + else return(y); /* y is NaN and x is finite */ else if(x==x) /* x is +-INF */ diff --git a/lib/libm/ieee/cbrt.c b/lib/libm/ieee/cbrt.c index fe5fb95..6d73779 100644 --- a/lib/libm/ieee/cbrt.c +++ b/lib/libm/ieee/cbrt.c @@ -50,8 +50,8 @@ static char sccsid[] = "@(#)cbrt.c 8.1 (Berkeley) 6/4/93"; * long interger at the address of a floating point number will be the * leading 32 bits of that floating point number (i.e., sign, exponent, * and the 20 most significant bits). - * On a National machine, it has different ordering; therefore, this code - * must be compiled with flag -DNATIONAL. + * On a National machine, it has different ordering; therefore, this code + * must be compiled with flag -DNATIONAL. */ #if !defined(vax)&&!defined(tahoe) @@ -65,7 +65,7 @@ static const double F= 45./28., G= 5./14.; -double cbrt(x) +double cbrt(x) double x; { double r,s,t=0.0,w; @@ -93,15 +93,15 @@ double x; t*=x; pt[n0]=pt[n0]/3+B2; } else - pt[n0]=px[n0]/3+B1; + pt[n0]=px[n0]/3+B1; /* new cbrt to 23 bits, may be implemented in single precision */ r=t*t/x; s=C+r*t; - t*=G+F/(s+E+D/s); + t*=G+F/(s+E+D/s); - /* chopped to 20 bits and make it larger than cbrt(x) */ + /* chopped to 20 bits and make it larger than cbrt(x) */ pt[n1]=0; pt[n0]+=0x00000001; diff --git a/lib/libm/ieee/support.c b/lib/libm/ieee/support.c index e976839..0ae10bb 100644 --- a/lib/libm/ieee/support.c +++ b/lib/libm/ieee/support.c @@ -35,37 +35,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 +78,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 +97,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 +117,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 +147,7 @@ double x,y; } double logb(x) -double x; +double x; { #ifdef national @@ -164,8 +164,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 +174,7 @@ double x; } finite(x) -double x; +double x; { #if defined(vax)||defined(tahoe) return(1); @@ -192,16 +192,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 +230,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 +294,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 +313,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 +342,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 +363,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 +381,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 +391,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 +408,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 +416,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 +435,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 +456,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 +502,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 */ |