Reduce diffs against vendor source (Sun fdlibm 5.3).

author: das <das@FreeBSD.org> 2005-02-04 18:26:06 +0000
committer: das <das@FreeBSD.org> 2005-02-04 18:26:06 +0000
commit: a47af911ba88a733157f219bcea374aafc0e6a34 (patch)
tree: 358c37c6dfb32659eea19d2a72397142931bc562 /lib/msun/src
parent: e8403b23681c8ac8f2c86499943b43c359a4e301 (diff)
download: FreeBSD-src-a47af911ba88a733157f219bcea374aafc0e6a34.zip
FreeBSD-src-a47af911ba88a733157f219bcea374aafc0e6a34.tar.gz
28 files changed, 476 insertions, 428 deletions
diff --git a/lib/msun/src/e_acos.c b/lib/msun/src/e_acos.c
index ff6269b..b021a1e 100644
--- a/lib/msun/src/e_acos.c
+++ b/lib/msun/src/e_acos.c
@@ -1,11 +1,12 @@
-/* @(#)e_acos.c 5.1 93/09/24 */
+
+/* @(#)e_acos.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -15,14 +16,14 @@ static char rcsid[] = "$FreeBSD$";
 #endif
 
 /* __ieee754_acos(x)
- * Method :
+ * Method :                  
  *	acos(x)  = pi/2 - asin(x)
  *	acos(-x) = pi/2 + asin(x)
  * For |x|<=0.5
  *	acos(x) = pi/2 - (x + x*x^2*R(x^2))	(see asin.c)
  * For x>0.5
  * 	acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
- *		= 2asin(sqrt((1-x)/2))
+ *		= 2asin(sqrt((1-x)/2))  
  *		= 2s + 2s*z*R(z) 	...z=(1-x)/2, s=sqrt(z)
  *		= 2f + (2c + 2s*z*R(z))
  *     where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
@@ -35,7 +36,7 @@ static char rcsid[] = "$FreeBSD$";
  *	if x is NaN, return x itself;
  *	if |x|>1, return NaN with invalid signal.
  *
- * Function needed: __ieee754_sqrt
+ * Function needed: sqrt
  */
 
 #include "math.h"
@@ -84,13 +85,13 @@ __ieee754_acos(double x)
 	    z = (one+x)*0.5;
 	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
 	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
-	    s = __ieee754_sqrt(z);
+	    s = sqrt(z);
 	    r = p/q;
 	    w = r*s-pio2_lo;
 	    return pi - 2.0*(s+w);
 	} else {			/* x > 0.5 */
 	    z = (one-x)*0.5;
-	    s = __ieee754_sqrt(z);
+	    s = sqrt(z);
 	    df = s;
 	    SET_LOW_WORD(df,0);
 	    c  = (z-df*df)/(s+df);
diff --git a/lib/msun/src/e_acosh.c b/lib/msun/src/e_acosh.c
index 98938f9..25e11c3 100644
--- a/lib/msun/src/e_acosh.c
+++ b/lib/msun/src/e_acosh.c
@@ -1,13 +1,15 @@
-/* @(#)e_acosh.c 5.1 93/09/24 */
+
+/* @(#)e_acosh.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
+ *
  */
 
 #ifndef lint
@@ -16,7 +18,7 @@ static char rcsid[] = "$FreeBSD$";
 
 /* __ieee754_acosh(x)
  * Method :
- *	Based on
+ *	Based on 
  *		acosh(x) = log [ x + sqrt(x*x-1) ]
  *	we have
  *		acosh(x) := log(x)+ln2,	if x is large; else
@@ -47,15 +49,15 @@ __ieee754_acosh(double x)
 	} else if(hx >=0x41b00000) {	/* x > 2**28 */
 	    if(hx >=0x7ff00000) {	/* x is inf of NaN */
 	        return x+x;
-	    } else
+	    } else 
 		return __ieee754_log(x)+ln2;	/* acosh(huge)=log(2x) */
 	} else if(((hx-0x3ff00000)|lx)==0) {
 	    return 0.0;			/* acosh(1) = 0 */
 	} else if (hx > 0x40000000) {	/* 2**28 > x > 2 */
 	    t=x*x;
-	    return __ieee754_log(2.0*x-one/(x+__ieee754_sqrt(t-one)));
+	    return __ieee754_log(2.0*x-one/(x+sqrt(t-one)));
 	} else {			/* 1<x<2 */
 	    t = x-one;
-	    return log1p(t+__ieee754_sqrt(2.0*t+t*t));
+	    return log1p(t+sqrt(2.0*t+t*t));
 	}
 }
diff --git a/lib/msun/src/e_asin.c b/lib/msun/src/e_asin.c
index f4c4f40..044c1be 100644
--- a/lib/msun/src/e_asin.c
+++ b/lib/msun/src/e_asin.c
@@ -1,11 +1,12 @@
-/* @(#)e_asin.c 5.1 93/09/24 */
+
+/* @(#)e_asin.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -15,12 +16,12 @@ static char rcsid[] = "$FreeBSD$";
 #endif
 
 /* __ieee754_asin(x)
- * Method :
+ * Method :                  
  *	Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
  *	we approximate asin(x) on [0,0.5] by
  *		asin(x) = x + x*x^2*R(x^2)
  *	where
- *		R(x^2) is a rational approximation of (asin(x)-x)/x^3
+ *		R(x^2) is a rational approximation of (asin(x)-x)/x^3 
  *	and its remez error is bounded by
  *		|(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
  *
@@ -78,12 +79,12 @@ __ieee754_asin(double x)
 	    GET_LOW_WORD(lx,x);
 	    if(((ix-0x3ff00000)|lx)==0)
 		    /* asin(1)=+-pi/2 with inexact */
-		return x*pio2_hi+x*pio2_lo;
-	    return (x-x)/(x-x);		/* asin(|x|>1) is NaN */
+		return x*pio2_hi+x*pio2_lo;	
+	    return (x-x)/(x-x);		/* asin(|x|>1) is NaN */   
 	} else if (ix<0x3fe00000) {	/* |x|<0.5 */
 	    if(ix<0x3e400000) {		/* if |x| < 2**-27 */
 		if(huge+x>one) return x;/* return x with inexact if x!=0*/
-	    } else
+	    } else 
 		t = x*x;
 		p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
 		q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
@@ -95,7 +96,7 @@ __ieee754_asin(double x)
 	t = w*0.5;
 	p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
 	q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
-	s = __ieee754_sqrt(t);
+	s = sqrt(t);
 	if(ix>=0x3FEF3333) { 	/* if |x| > 0.975 */
 	    w = p/q;
 	    t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
@@ -107,6 +108,6 @@ __ieee754_asin(double x)
 	    p  = 2.0*s*r-(pio2_lo-2.0*c);
 	    q  = pio4_hi-2.0*w;
 	    t  = pio4_hi-(p-q);
-	}
-	if(hx>0) return t; else return -t;
+	}    
+	if(hx>0) return t; else return -t;    
 }
diff --git a/lib/msun/src/e_atan2.c b/lib/msun/src/e_atan2.c
index e9c9858..a05a840 100644
--- a/lib/msun/src/e_atan2.c
+++ b/lib/msun/src/e_atan2.c
@@ -1,13 +1,15 @@
-/* @(#)e_atan2.c 5.1 93/09/24 */
+
+/* @(#)e_atan2.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
+ *
  */
 
 #ifndef lint
@@ -17,7 +19,7 @@ static char rcsid[] = "$FreeBSD$";
 /* __ieee754_atan2(y,x)
  * Method :
  *	1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
- *	2. Reduce x to positive by (if x and y are unexceptional):
+ *	2. Reduce x to positive by (if x and y are unexceptional): 
  *		ARG (x+iy) = arctan(y/x)   	   ... if x > 0,
  *		ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
  *
@@ -35,9 +37,9 @@ static char rcsid[] = "$FreeBSD$";
  *	ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
  *
  * Constants:
- * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough
+ * The hexadecimal values are the intended ones for the following 
+ * constants. The decimal values may be used, provided that the 
+ * compiler will convert from decimal to binary accurately enough 
  * to produce the hexadecimal values shown.
  */
 
@@ -66,13 +68,13 @@ __ieee754_atan2(double y, double x)
 	if(((ix|((lx|-lx)>>31))>0x7ff00000)||
 	   ((iy|((ly|-ly)>>31))>0x7ff00000))	/* x or y is NaN */
 	   return x+y;
-	if(((hx-0x3ff00000)|lx)==0) return atan(y);   /* x=1.0 */
+	if((hx-0x3ff00000|lx)==0) return atan(y);   /* x=1.0 */
 	m = ((hy>>31)&1)|((hx>>30)&2);	/* 2*sign(x)+sign(y) */
 
     /* when y = 0 */
 	if((iy|ly)==0) {
 	    switch(m) {
-		case 0:
+		case 0: 
 		case 1: return y; 	/* atan(+-0,+anything)=+-0 */
 		case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
 		case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
@@ -80,7 +82,7 @@ __ieee754_atan2(double y, double x)
 	}
     /* when x = 0 */
 	if((ix|lx)==0) return (hy<0)?  -pi_o_2-tiny: pi_o_2+tiny;
-
+	    
     /* when x is INF */
 	if(ix==0x7ff00000) {
 	    if(iy==0x7ff00000) {
diff --git a/lib/msun/src/e_atanh.c b/lib/msun/src/e_atanh.c
index d838c13..9a89669 100644
--- a/lib/msun/src/e_atanh.c
+++ b/lib/msun/src/e_atanh.c
@@ -1,13 +1,15 @@
-/* @(#)e_atanh.c 5.1 93/09/24 */
+
+/* @(#)e_atanh.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
+ *
  */
 
 #ifndef lint
@@ -21,7 +23,7 @@ static char rcsid[] = "$FreeBSD$";
  *                  1              2x                          x
  *	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
  *                  2             1 - x                      1 - x
- *
+ *	
  * 	For x<0.5
  *	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
  *
@@ -48,14 +50,14 @@ __ieee754_atanh(double x)
 	ix = hx&0x7fffffff;
 	if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
 	    return (x-x)/(x-x);
-	if(ix==0x3ff00000)
+	if(ix==0x3ff00000) 
 	    return x/zero;
 	if(ix<0x3e300000&&(huge+x)>zero) return x;	/* x<2**-28 */
 	SET_HIGH_WORD(x,ix);
 	if(ix<0x3fe00000) {		/* x < 0.5 */
 	    t = x+x;
 	    t = 0.5*log1p(t+t*x/(one-x));
-	} else
+	} else 
 	    t = 0.5*log1p((x+x)/(one-x));
 	if(hx>=0) return t; else return -t;
 }
diff --git a/lib/msun/src/e_cosh.c b/lib/msun/src/e_cosh.c
index 1a67f9a..f8a2b9b 100644
--- a/lib/msun/src/e_cosh.c
+++ b/lib/msun/src/e_cosh.c
@@ -1,11 +1,12 @@
-/* @(#)e_cosh.c 5.1 93/09/24 */
+
+/* @(#)e_cosh.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -15,18 +16,18 @@ static char rcsid[] = "$FreeBSD$";
 #endif
 
 /* __ieee754_cosh(x)
- * Method :
+ * Method : 
  * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
- *	1. Replace x by |x| (cosh(x) = cosh(-x)).
- *	2.
- *		                                        [ exp(x) - 1 ]^2
+ *	1. Replace x by |x| (cosh(x) = cosh(-x)). 
+ *	2. 
+ *		                                        [ exp(x) - 1 ]^2 
  *	    0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
  *			       			           2*exp(x)
  *
  *		                                  exp(x) +  1/exp(x)
  *	    ln2/2    <= x <= 22     :  cosh(x) := -------------------
  *			       			          2
- *	    22       <= x <= lnovft :  cosh(x) := exp(x)/2
+ *	    22       <= x <= lnovft :  cosh(x) := exp(x)/2 
  *	    lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
  *	    ln2ovft  <  x	    :  cosh(x) := huge*huge (overflow)
  *
@@ -52,7 +53,7 @@ __ieee754_cosh(double x)
 	ix &= 0x7fffffff;
 
     /* x is INF or NaN */
-	if(ix>=0x7ff00000) return x*x;
+	if(ix>=0x7ff00000) return x*x;	
 
     /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
 	if(ix<0x3fd62e43) {
diff --git a/lib/msun/src/e_exp.c b/lib/msun/src/e_exp.c
index 7ce259c..60b8b2a 100644
--- a/lib/msun/src/e_exp.c
+++ b/lib/msun/src/e_exp.c
@@ -1,11 +1,11 @@
-/* @(#)e_exp.c 5.1 93/09/24 */
+
+/* @(#)e_exp.c 1.6 04/04/22 */
 /*
  * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -22,36 +22,36 @@ static char rcsid[] = "$FreeBSD$";
  *      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
  *	Given x, find r and integer k such that
  *
- *               x = k*ln2 + r,  |r| <= 0.5*ln2.
+ *               x = k*ln2 + r,  |r| <= 0.5*ln2.  
  *
- *      Here r will be represented as r = hi-lo for better
+ *      Here r will be represented as r = hi-lo for better 
  *	accuracy.
  *
  *   2. Approximation of exp(r) by a special rational function on
  *	the interval [0,0.34658]:
  *	Write
  *	    R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
- *      We use a special Reme algorithm on [0,0.34658] to generate
- * 	a polynomial of degree 5 to approximate R. The maximum error
+ *      We use a special Remes algorithm on [0,0.34658] to generate 
+ * 	a polynomial of degree 5 to approximate R. The maximum error 
  *	of this polynomial approximation is bounded by 2**-59. In
  *	other words,
  *	    R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
  *  	(where z=r*r, and the values of P1 to P5 are listed below)
  *	and
  *	    |                  5          |     -59
- *	    | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
+ *	    | 2.0+P1*z+...+P5*z   -  R(z) | <= 2 
  *	    |                             |
  *	The computation of exp(r) thus becomes
  *                             2*r
  *		exp(r) = 1 + -------
  *		              R - r
- *                                 r*R1(r)
+ *                                 r*R1(r)	
  *		       = 1 + r + ----------- (for better accuracy)
  *		                  2 - R1(r)
  *	where
  *			         2       4             10
  *		R1(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
- *
+ *	
  *   3. Scale back to obtain exp(x):
  *	From step 1, we have
  *	   exp(x) = 2^k * exp(r)
@@ -66,13 +66,13 @@ static char rcsid[] = "$FreeBSD$";
  *	1 ulp (unit in the last place).
  *
  * Misc. info.
- *	For IEEE double
+ *	For IEEE double 
  *	    if x >  7.09782712893383973096e+02 then exp(x) overflow
  *	    if x < -7.45133219101941108420e+02 then exp(x) underflow
  *
  * Constants:
- * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the
+ * The hexadecimal values are the intended ones for the following 
+ * constants. The decimal values may be used, provided that the 
  * compiler will convert from decimal to binary accurately enough
  * to produce the hexadecimal values shown.
  */
@@ -124,17 +124,17 @@ __ieee754_exp(double x)	/* default IEEE double exp */
 	}
 
     /* argument reduction */
-	if(hx > 0x3fd62e42) {		/* if  |x| > 0.5 ln2 */
+	if(hx > 0x3fd62e42) {		/* if  |x| > 0.5 ln2 */ 
 	    if(hx < 0x3FF0A2B2) {	/* and |x| < 1.5 ln2 */
 		hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
 	    } else {
-		k  = invln2*x+halF[xsb];
+		k  = (int)(invln2*x+halF[xsb]);
 		t  = k;
 		hi = x - t*ln2HI[0];	/* t*ln2HI is exact here */
 		lo = t*ln2LO[0];
 	    }
 	    x  = hi - lo;
-	}
+	} 
 	else if(hx < 0x3e300000)  {	/* when |x|<2**-28 */
 	    if(huge+x>one) return one+x;/* trigger inexact */
 	}
@@ -143,7 +143,7 @@ __ieee754_exp(double x)	/* default IEEE double exp */
     /* x is now in primary range */
 	t  = x*x;
 	c  = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
-	if(k==0) 	return one-((x*c)/(c-2.0)-x);
+	if(k==0) 	return one-((x*c)/(c-2.0)-x); 
 	else 		y = one-((lo-(x*c)/(2.0-c))-hi);
 	if(k >= -1021) {
 	    u_int32_t hy;
diff --git a/lib/msun/src/e_fmod.c b/lib/msun/src/e_fmod.c
index 06e7c15..b731ce1 100644
--- a/lib/msun/src/e_fmod.c
+++ b/lib/msun/src/e_fmod.c
@@ -1,11 +1,12 @@
-/* @(#)e_fmod.c 5.1 93/09/24 */
+
+/* @(#)e_fmod.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -14,7 +15,7 @@
 static char rcsid[] = "$FreeBSD$";
 #endif
 
-/*
+/* 
  * __ieee754_fmod(x,y)
  * Return x mod y in exact arithmetic
  * Method: shift and subtract
@@ -43,7 +44,7 @@ __ieee754_fmod(double x, double y)
 	    return (x*y)/(x*y);
 	if(hx<=hy) {
 	    if((hx<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
-	    if(lx==ly)
+	    if(lx==ly) 
 		return Zero[(u_int32_t)sx>>31];	/* |x|=|y| return x*0*/
 	}
 
@@ -66,7 +67,7 @@ __ieee754_fmod(double x, double y)
 	} else iy = (hy>>20)-1023;
 
     /* set up {hx,lx}, {hy,ly} and align y to x */
-	if(ix >= -1022)
+	if(ix >= -1022) 
 	    hx = 0x00100000|(0x000fffff&hx);
 	else {		/* subnormal x, shift x to normal */
 	    n = -1022-ix;
@@ -78,7 +79,7 @@ __ieee754_fmod(double x, double y)
 		lx = 0;
 	    }
 	}
-	if(iy >= -1022)
+	if(iy >= -1022) 
 	    hy = 0x00100000|(0x000fffff&hy);
 	else {		/* subnormal y, shift y to normal */
 	    n = -1022-iy;
diff --git a/lib/msun/src/e_gamma.c b/lib/msun/src/e_gamma.c
index 6cfbe17..e3ebd7c 100644
--- a/lib/msun/src/e_gamma.c
+++ b/lib/msun/src/e_gamma.c
@@ -1,11 +1,12 @@
-/* @(#)e_gamma.c 5.1 93/09/24 */
+
+/* @(#)e_gamma.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  *
diff --git a/lib/msun/src/e_gamma_r.c b/lib/msun/src/e_gamma_r.c
index 79ae4ea..48b981a 100644
--- a/lib/msun/src/e_gamma_r.c
+++ b/lib/msun/src/e_gamma_r.c
@@ -1,13 +1,15 @@
-/* @(#)er_gamma.c 5.1 93/09/24 */
+
+/* @(#)e_gamma_r.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
+ *
  */
 
 #ifndef lint
@@ -15,8 +17,8 @@ static char rcsid[] = "$FreeBSD$";
 #endif
 
 /* __ieee754_gamma_r(x, signgamp)
- * Reentrant version of the logarithm of the Gamma function
- * with user provide pointer for the sign of Gamma(x).
+ * Reentrant version of the logarithm of the Gamma function 
+ * with user provide pointer for the sign of Gamma(x). 
  *
  * Method: See __ieee754_lgamma_r
  */
diff --git a/lib/msun/src/e_hypot.c b/lib/msun/src/e_hypot.c
index aa50051..e2f82ac 100644
--- a/lib/msun/src/e_hypot.c
+++ b/lib/msun/src/e_hypot.c
@@ -1,11 +1,12 @@
-/* @(#)e_hypot.c 5.1 93/09/24 */
+
+/* @(#)e_hypot.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -16,12 +17,12 @@ static char rcsid[] = "$FreeBSD$";
 
 /* __ieee754_hypot(x,y)
  *
- * Method :
- *	If (assume round-to-nearest) z=x*x+y*y
- *	has error less than sqrt(2)/2 ulp, than
+ * Method :                  
+ *	If (assume round-to-nearest) z=x*x+y*y 
+ *	has error less than sqrt(2)/2 ulp, than 
  *	sqrt(z) has error less than 1 ulp (exercise).
  *
- *	So, compute sqrt(x*x+y*y) with some care as
+ *	So, compute sqrt(x*x+y*y) with some care as 
  *	follows to get the error below 1 ulp:
  *
  *	Assume x>y>0;
@@ -31,10 +32,10 @@ static char rcsid[] = "$FreeBSD$";
  *	where x1 = x with lower 32 bits cleared, x2 = x-x1; else
  *	2. if x <= 2y use
  *		t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
- *	where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
+ *	where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, 
  *	y1= y with lower 32 bits chopped, y2 = y-y1.
- *
- *	NOTE: scaling may be necessary if some argument is too
+ *		
+ *	NOTE: scaling may be necessary if some argument is too 
  *	      large or too tiny
  *
  * Special cases:
@@ -42,8 +43,8 @@ static char rcsid[] = "$FreeBSD$";
  *	hypot(x,y) is NAN if x or y is NAN.
  *
  * Accuracy:
- * 	hypot(x,y) returns sqrt(x^2+y^2) with error less
- * 	than 1 ulps (units in the last place)
+ * 	hypot(x,y) returns sqrt(x^2+y^2) with error less 
+ * 	than 1 ulps (units in the last place) 
  */
 
 #include "math.h"
@@ -103,7 +104,7 @@ __ieee754_hypot(double x, double y)
 	    t1 = 0;
 	    SET_HIGH_WORD(t1,ha);
 	    t2 = a-t1;
-	    w  = __ieee754_sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
+	    w  = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
 	} else {
 	    a  = a+a;
 	    y1 = 0;
@@ -112,7 +113,7 @@ __ieee754_hypot(double x, double y)
 	    t1 = 0;
 	    SET_HIGH_WORD(t1,ha+0x00100000);
 	    t2 = a - t1;
-	    w  = __ieee754_sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
+	    w  = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
 	}
 	if(k!=0) {
 	    u_int32_t high;
diff --git a/lib/msun/src/e_j0.c b/lib/msun/src/e_j0.c
index ca08293..493f0ef 100644
--- a/lib/msun/src/e_j0.c
+++ b/lib/msun/src/e_j0.c
@@ -1,11 +1,12 @@
-/* @(#)e_j0.c 5.1 93/09/24 */
+
+/* @(#)e_j0.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -33,20 +34,20 @@ static char rcsid[] = "$FreeBSD$";
  * 	   (To avoid cancellation, use
  *		sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  * 	    to compute the worse one.)
- *
+ *	   
  *	3 Special cases
  *		j0(nan)= nan
  *		j0(0) = 1
  *		j0(inf) = 0
- *
+ *		
  * Method -- y0(x):
  *	1. For x<2.
- *	   Since
+ *	   Since 
  *		y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
  *	   therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
  *	   We use the following function to approximate y0,
  *		y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
- *	   where
+ *	   where 
  *		U(z) = u00 + u01*z + ... + u06*z^6
  *		V(z) = 1  + v01*z + ... + v04*z^4
  *	   with absolute approximation error bounded by 2**-72.
@@ -151,7 +152,7 @@ __ieee754_y0(double x)
 	EXTRACT_WORDS(hx,lx,x);
         ix = 0x7fffffff&hx;
     /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
-	if(ix>=0x7ff00000) return  one/(x+x*x);
+	if(ix>=0x7ff00000) return  one/(x+x*x); 
         if((ix|lx)==0) return -one/zero;
         if(hx<0) return zero/zero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
@@ -284,7 +285,7 @@ static const double pS2[5] = {
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
 	return one+ r/s;
 }
-
+		
 
 /* For x >= 8, the asymptotic expansions of qzero is
  *	-1/8 s + 75/1024 s^3 - ..., where s = 1/x.
diff --git a/lib/msun/src/e_j1.c b/lib/msun/src/e_j1.c
index 97fdcf8..a01b812 100644
--- a/lib/msun/src/e_j1.c
+++ b/lib/msun/src/e_j1.c
@@ -1,11 +1,12 @@
-/* @(#)e_j1.c 5.1 93/09/24 */
+
+/* @(#)e_j1.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -34,16 +35,16 @@ static char rcsid[] = "$FreeBSD$";
  * 	   (To avoid cancellation, use
  *		sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  * 	    to compute the worse one.)
- *
+ *	   
  *	3 Special cases
  *		j1(nan)= nan
  *		j1(0) = 0
  *		j1(inf) = 0
- *
+ *		
  * Method -- y1(x):
- *	1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
+ *	1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN 
  *	2. For x<2.
- *	   Since
+ *	   Since 
  *		y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
  *	   therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
  *	   We use the following function to approximate y1,
@@ -148,7 +149,7 @@ __ieee754_y1(double x)
 	EXTRACT_WORDS(hx,lx,x);
         ix = 0x7fffffff&hx;
     /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
-	if(ix>=0x7ff00000) return  one/(x+x*x);
+	if(ix>=0x7ff00000) return  one/(x+x*x); 
         if((ix|lx)==0) return -one/zero;
         if(hx<0) return zero/zero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
@@ -178,10 +179,10 @@ __ieee754_y1(double x)
                     z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
                 }
                 return z;
-        }
+        } 
         if(ix<=0x3c900000) {    /* x < 2**-54 */
             return(-tpi/x);
-        }
+        } 
         z = x*x;
         u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
         v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
@@ -278,7 +279,7 @@ static const double ps2[5] = {
         s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
         return one+ r/s;
 }
-
+		
 
 /* For x >= 8, the asymptotic expansions of qone is
  *	3/8 s - 105/1024 s^3 - ..., where s = 1/x.
diff --git a/lib/msun/src/e_jn.c b/lib/msun/src/e_jn.c
index d0162cd..ca26230 100644
--- a/lib/msun/src/e_jn.c
+++ b/lib/msun/src/e_jn.c
@@ -1,11 +1,12 @@
-/* @(#)e_jn.c 5.1 93/09/24 */
+
+/* @(#)e_jn.c 1.4 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -18,7 +19,7 @@ static char rcsid[] = "$FreeBSD$";
  * __ieee754_jn(n, x), __ieee754_yn(n, x)
  * floating point Bessel's function of the 1st and 2nd kind
  * of order n
- *
+ *          
  * Special cases:
  *	y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  *	y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
@@ -37,7 +38,7 @@ static char rcsid[] = "$FreeBSD$";
  *	yn(n,x) is similar in all respects, except
  *	that forward recursion is used for all
  *	values of n>1.
- *
+ *	
  */
 
 #include "math.h"
@@ -64,7 +65,7 @@ __ieee754_jn(int n, double x)
 	ix = 0x7fffffff&hx;
     /* if J(n,NaN) is NaN */
 	if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
-	if(n<0){
+	if(n<0){		
 		n = -n;
 		x = -x;
 		hx ^= 0x80000000;
@@ -75,13 +76,13 @@ __ieee754_jn(int n, double x)
 	x = fabs(x);
 	if((ix|lx)==0||ix>=0x7ff00000) 	/* if x is 0 or inf */
 	    b = zero;
-	else if((double)n<=x) {
+	else if((double)n<=x) {   
 		/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
 	    if(ix>=0x52D00000) { /* x > 2**302 */
-    /* (x >> n**2)
+    /* (x >> n**2) 
      *	    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
      *	    Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-     *	    Let s=sin(x), c=cos(x),
+     *	    Let s=sin(x), c=cos(x), 
      *		xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
      *
      *		   n	sin(xn)*sqt2	cos(xn)*sqt2
@@ -98,7 +99,7 @@ __ieee754_jn(int n, double x)
 		    case 3: temp =  cos(x)-sin(x); break;
 		}
 		b = invsqrtpi*temp/sqrt(x);
-	    } else {
+	    } else {	
 	        a = __ieee754_j0(x);
 	        b = __ieee754_j1(x);
 	        for(i=1;i<n;i++){
@@ -109,7 +110,7 @@ __ieee754_jn(int n, double x)
 	    }
 	} else {
 	    if(ix<0x3e100000) {	/* x < 2**-29 */
-    /* x is tiny, return the first Taylor expansion of J(n,x)
+    /* x is tiny, return the first Taylor expansion of J(n,x) 
      * J(n,x) = 1/n!*(x/2)^n  - ...
      */
 		if(n>33)	/* underflow */
@@ -124,14 +125,14 @@ __ieee754_jn(int n, double x)
 		}
 	    } else {
 		/* use backward recurrence */
-		/* 			x      x^2      x^2
+		/* 			x      x^2      x^2       
 		 *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 		 *			2n  - 2(n+1) - 2(n+2)
 		 *
-		 * 			1      1        1
+		 * 			1      1        1       
 		 *  (for large x)   =  ----  ------   ------   .....
 		 *			2n   2(n+1)   2(n+2)
-		 *			-- - ------ - ------ -
+		 *			-- - ------ - ------ - 
 		 *			 x     x         x
 		 *
 		 * Let w = 2n/x and h=2/x, then the above quotient
@@ -147,9 +148,9 @@ __ieee754_jn(int n, double x)
 		 * To determine how many terms needed, let
 		 * Q(0) = w, Q(1) = w(w+h) - 1,
 		 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
-		 * When Q(k) > 1e4	good for single
-		 * When Q(k) > 1e9	good for double
-		 * When Q(k) > 1e17	good for quadruple
+		 * When Q(k) > 1e4	good for single 
+		 * When Q(k) > 1e9	good for double 
+		 * When Q(k) > 1e17	good for quadruple 
 		 */
 	    /* determine k */
 		double t,v;
@@ -228,10 +229,10 @@ __ieee754_yn(int n, double x)
 	if(n==1) return(sign*__ieee754_y1(x));
 	if(ix==0x7ff00000) return zero;
 	if(ix>=0x52D00000) { /* x > 2**302 */
-    /* (x >> n**2)
+    /* (x >> n**2) 
      *	    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
      *	    Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-     *	    Let s=sin(x), c=cos(x),
+     *	    Let s=sin(x), c=cos(x), 
      *		xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
      *
      *		   n	sin(xn)*sqt2	cos(xn)*sqt2
diff --git a/lib/msun/src/e_lgamma.c b/lib/msun/src/e_lgamma.c
index 8158cbe..39d3f85 100644
--- a/lib/msun/src/e_lgamma.c
+++ b/lib/msun/src/e_lgamma.c
@@ -1,13 +1,15 @@
-/* @(#)e_lgamma.c 5.1 93/09/24 */
+
+/* @(#)e_lgamma.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
+ *
  */
 
 #ifndef lint
diff --git a/lib/msun/src/e_lgamma_r.c b/lib/msun/src/e_lgamma_r.c
index ea2d97e..91989d5 100644
--- a/lib/msun/src/e_lgamma_r.c
+++ b/lib/msun/src/e_lgamma_r.c
@@ -1,13 +1,15 @@
-/* @(#)er_lgamma.c 5.1 93/09/24 */
+
+/* @(#)e_lgamma_r.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
+ *
  */
 
 #ifndef lint
@@ -15,12 +17,12 @@ static char rcsid[] = "$FreeBSD$";
 #endif
 
 /* __ieee754_lgamma_r(x, signgamp)
- * Reentrant version of the logarithm of the Gamma function
- * with user provide pointer for the sign of Gamma(x).
+ * Reentrant version of the logarithm of the Gamma function 
+ * with user provide pointer for the sign of Gamma(x). 
  *
  * Method:
  *   1. Argument Reduction for 0 < x <= 8
- * 	Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
+ * 	Since gamma(1+s)=s*gamma(s), for x in [0,8], we may 
  * 	reduce x to a number in [1.5,2.5] by
  * 		lgamma(1+s) = log(s) + lgamma(s)
  *	for example,
@@ -58,27 +60,27 @@ static char rcsid[] = "$FreeBSD$";
  *	by
  *	  			    3       5             11
  *		w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
- *	where
+ *	where 
  *		|w - f(z)| < 2**-58.74
- *
+ *		
  *   4. For negative x, since (G is gamma function)
  *		-x*G(-x)*G(x) = pi/sin(pi*x),
  * 	we have
  * 		G(x) = pi/(sin(pi*x)*(-x)*G(-x))
  *	since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
- *	Hence, for x<0, signgam = sign(sin(pi*x)) and
+ *	Hence, for x<0, signgam = sign(sin(pi*x)) and 
  *		lgamma(x) = log(|Gamma(x)|)
  *			  = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
- *	Note: one should avoid compute pi*(-x) directly in the
+ *	Note: one should avoid compute pi*(-x) directly in the 
  *	      computation of sin(pi*(-x)).
- *
+ *		
  *   5. Special Cases
  *		lgamma(2+s) ~ s*(1-Euler) for tiny s
  *		lgamma(1)=lgamma(2)=0
  *		lgamma(x) ~ -log(x) for tiny x
  *		lgamma(0) = lgamma(inf) = inf
  *	 	lgamma(-integer) = +-inf
- *
+ *	
  */
 
 #include "math.h"
@@ -187,9 +189,9 @@ static const double zero=  0.00000000000000000000e+00;
         }
 	switch (n) {
 	    case 0:   y =  __kernel_sin(pi*y,zero,0); break;
-	    case 1:
+	    case 1:   
 	    case 2:   y =  __kernel_cos(pi*(0.5-y),zero); break;
-	    case 3:
+	    case 3:  
 	    case 4:   y =  __kernel_sin(pi*(one-y),zero,0); break;
 	    case 5:
 	    case 6:   y = -__kernel_cos(pi*(y-1.5),zero); break;
@@ -258,7 +260,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
 		p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
 		p  = z*p1-(tt-w*(p2+y*p3));
 		r += (tf + p); break;
-	      case 2:
+	      case 2:	
 		p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
 		p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
 		r += (-0.5*y + p1/p2);
@@ -287,7 +289,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
 	    y = z*z;
 	    w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
 	    r = (x-half)*(t-one)+w;
-	} else
+	} else 
     /* 2**58 <= x <= inf */
 	    r =  x*(__ieee754_log(x)-one);
 	if(hx<0) r = nadj - r;
diff --git a/lib/msun/src/e_log.c b/lib/msun/src/e_log.c
index aa1a70b..eee7cf9 100644
--- a/lib/msun/src/e_log.c
+++ b/lib/msun/src/e_log.c
@@ -1,11 +1,12 @@
-/* @(#)e_log.c 5.1 93/09/24 */
+
+/* @(#)e_log.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -17,17 +18,17 @@ static char rcsid[] = "$FreeBSD$";
 /* __ieee754_log(x)
  * Return the logrithm of x
  *
- * Method :
- *   1. Argument Reduction: find k and f such that
- *			x = 2^k * (1+f),
+ * Method :                  
+ *   1. Argument Reduction: find k and f such that 
+ *			x = 2^k * (1+f), 
  *	   where  sqrt(2)/2 < 1+f < sqrt(2) .
  *
  *   2. Approximation of log(1+f).
  *	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
  *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
  *	     	 = 2s + s*R
- *      We use a special Reme algorithm on [0,0.1716] to generate
- * 	a polynomial of degree 14 to approximate R The maximum error
+ *      We use a special Reme algorithm on [0,0.1716] to generate 
+ * 	a polynomial of degree 14 to approximate R The maximum error 
  *	of this polynomial approximation is bounded by 2**-58.45. In
  *	other words,
  *		        2      4      6      8      10      12      14
@@ -35,22 +36,22 @@ static char rcsid[] = "$FreeBSD$";
  *  	(the values of Lg1 to Lg7 are listed in the program)
  *	and
  *	    |      2          14          |     -58.45
- *	    | Lg1*s +...+Lg7*s    -  R(z) | <= 2
+ *	    | Lg1*s +...+Lg7*s    -  R(z) | <= 2 
  *	    |                             |
  *	Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
  *	In order to guarantee error in log below 1ulp, we compute log
  *	by
  *		log(1+f) = f - s*(f - R)	(if f is not too large)
  *		log(1+f) = f - (hfsq - s*(hfsq+R)).	(better accuracy)
- *
- *	3. Finally,  log(x) = k*ln2 + log(1+f).
+ *	
+ *	3. Finally,  log(x) = k*ln2 + log(1+f).  
  *			    = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
- *	   Here ln2 is split into two floating point number:
+ *	   Here ln2 is split into two floating point number: 
  *			ln2_hi + ln2_lo,
  *	   where n*ln2_hi is always exact for |n| < 2000.
  *
  * Special cases:
- *	log(x) is NaN with signal if x < 0 (including -INF) ;
+ *	log(x) is NaN with signal if x < 0 (including -INF) ; 
  *	log(+INF) is +INF; log(0) is -INF with signal;
  *	log(NaN) is that NaN with no signal.
  *
@@ -59,9 +60,9 @@ static char rcsid[] = "$FreeBSD$";
  *	1 ulp (unit in the last place).
  *
  * Constants:
- * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough
+ * The hexadecimal values are the intended ones for the following 
+ * constants. The decimal values may be used, provided that the 
+ * compiler will convert from decimal to binary accurately enough 
  * to produce the hexadecimal values shown.
  */
 
@@ -93,12 +94,12 @@ __ieee754_log(double x)
 
 	k=0;
 	if (hx < 0x00100000) {			/* x < 2**-1022  */
-	    if (((hx&0x7fffffff)|lx)==0)
+	    if (((hx&0x7fffffff)|lx)==0) 
 		return -two54/zero;		/* log(+-0)=-inf */
 	    if (hx<0) return (x-x)/zero;	/* log(-#) = NaN */
 	    k -= 54; x *= two54; /* subnormal number, scale up x */
 	    GET_HIGH_WORD(hx,x);
-	}
+	} 
 	if (hx >= 0x7ff00000) return x+x;
 	k += (hx>>20)-1023;
 	hx &= 0x000fffff;
@@ -113,14 +114,14 @@ __ieee754_log(double x)
 	    if(k==0) return f-R; else {dk=(double)k;
 	    	     return dk*ln2_hi-((R-dk*ln2_lo)-f);}
 	}
- 	s = f/(2.0+f);
+ 	s = f/(2.0+f); 
 	dk = (double)k;
 	z = s*s;
 	i = hx-0x6147a;
 	w = z*z;
 	j = 0x6b851-hx;
-	t1= w*(Lg2+w*(Lg4+w*Lg6));
-	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
+	t1= w*(Lg2+w*(Lg4+w*Lg6)); 
+	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 
 	i |= j;
 	R = t2+t1;
 	if(i>0) {
diff --git a/lib/msun/src/e_log10.c b/lib/msun/src/e_log10.c
index 7a76def..871a44a 100644
--- a/lib/msun/src/e_log10.c
+++ b/lib/msun/src/e_log10.c
@@ -1,11 +1,12 @@
-/* @(#)e_log10.c 5.1 93/09/24 */
+
+/* @(#)e_log10.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -16,26 +17,26 @@ static char rcsid[] = "$FreeBSD$";
 
 /* __ieee754_log10(x)
  * Return the base 10 logarithm of x
- *
+ * 
  * Method :
  *	Let log10_2hi = leading 40 bits of log10(2) and
  *	    log10_2lo = log10(2) - log10_2hi,
  *	    ivln10   = 1/log(10) rounded.
  *	Then
- *		n = ilogb(x),
+ *		n = ilogb(x), 
  *		if(n<0)  n = n+1;
  *		x = scalbn(x,-n);
  *		log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
  *
  * Note 1:
- *	To guarantee log10(10**n)=n, where 10**n is normal, the rounding
+ *	To guarantee log10(10**n)=n, where 10**n is normal, the rounding 
  *	mode must set to Round-to-Nearest.
  * Note 2:
  *	[1/log(10)] rounded to 53 bits has error  .198   ulps;
  *	log10 is monotonic at all binary break points.
  *
  * Special cases:
- *	log10(x) is NaN with signal if x < 0;
+ *	log10(x) is NaN with signal if x < 0; 
  *	log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
  *	log10(NaN) is that NaN with no signal;
  *	log10(10**N) = N  for N=0,1,...,22.
diff --git a/lib/msun/src/e_pow.c b/lib/msun/src/e_pow.c
index 04a2da5..95d0855 100644
--- a/lib/msun/src/e_pow.c
+++ b/lib/msun/src/e_pow.c
@@ -1,11 +1,10 @@
-/* @(#)e_pow.c 5.1 93/09/24 */
+/* @(#)e_pow.c 1.5 04/04/22 SMI */
 /*
  * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -21,7 +20,7 @@ static char rcsid[] = "$FreeBSD$";
  *	1. Compute and return log2(x) in two pieces:
  *		log2(x) = w1 + w2,
  *	   where w1 has 53-24 = 29 bit trailing zeros.
- *	2. Perform y*log2(x) = n+y' by simulating muti-precision
+ *	2. Perform y*log2(x) = n+y' by simulating muti-precision 
  *	   arithmetic, where |y'|<=0.5.
  *	3. Return x**y = 2**n*exp(y'*log2)
  *
@@ -49,13 +48,13 @@ static char rcsid[] = "$FreeBSD$";
  * Accuracy:
  *	pow(x,y) returns x**y nearly rounded. In particular
  *			pow(integer,integer)
- *	always returns the correct integer provided it is
+ *	always returns the correct integer provided it is 
  *	representable.
  *
  * Constants :
- * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough
+ * The hexadecimal values are the intended ones for the following 
+ * constants. The decimal values may be used, provided that the 
+ * compiler will convert from decimal to binary accurately enough 
  * to produce the hexadecimal values shown.
  */
 
@@ -99,7 +98,7 @@ double
 __ieee754_pow(double x, double y)
 {
 	double z,ax,z_h,z_l,p_h,p_l;
-	double y1,t1,t2,r,s,sn,t,u,v,w;
+	double y1,t1,t2,r,s,t,u,v,w;
 	int32_t i,j,k,yisint,n;
 	int32_t hx,hy,ix,iy;
 	u_int32_t lx,ly;
@@ -109,12 +108,12 @@ __ieee754_pow(double x, double y)
 	ix = hx&0x7fffffff;  iy = hy&0x7fffffff;
 
     /* y==zero: x**0 = 1 */
-	if((iy|ly)==0) return one;
+	if((iy|ly)==0) return one; 	
 
     /* +-NaN return x+y */
 	if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
-	   iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
-		return x+y;
+	   iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) 
+		return x+y;	
 
     /* determine if y is an odd int when x < 0
      * yisint = 0	... y is not an integer
@@ -122,7 +121,7 @@ __ieee754_pow(double x, double y)
      * yisint = 2	... y is an even int
      */
 	yisint  = 0;
-	if(hx<0) {
+	if(hx<0) {	
 	    if(iy>=0x43400000) yisint = 2; /* even integer y */
 	    else if(iy>=0x3ff00000) {
 		k = (iy>>20)-0x3ff;	   /* exponent */
@@ -133,11 +132,11 @@ __ieee754_pow(double x, double y)
 		    j = iy>>(20-k);
 		    if((j<<(20-k))==iy) yisint = 2-(j&1);
 		}
-	    }
-	}
+	    }		
+	} 
 
     /* special value of y */
-	if(ly==0) {
+	if(ly==0) { 	
 	    if (iy==0x7ff00000) {	/* y is +-inf */
 	        if(((ix-0x3ff00000)|lx)==0)
 		    return  y - y;	/* inf**+-1 is NaN */
@@ -145,14 +144,14 @@ __ieee754_pow(double x, double y)
 		    return (hy>=0)? y: zero;
 	        else			/* (|x|<1)**-,+inf = inf,0 */
 		    return (hy<0)?-y: zero;
-	    }
+	    } 
 	    if(iy==0x3ff00000) {	/* y is  +-1 */
 		if(hy<0) return one/x; else return x;
 	    }
 	    if(hy==0x40000000) return x*x; /* y is  2 */
 	    if(hy==0x3fe00000) {	/* y is  0.5 */
 		if(hx>=0)	/* x >= +0 */
-		return __ieee754_sqrt(x);
+		return sqrt(x);	
 	    }
 	}
 
@@ -165,13 +164,13 @@ __ieee754_pow(double x, double y)
 		if(hx<0) {
 		    if(((ix-0x3ff00000)|yisint)==0) {
 			z = (z-z)/(z-z); /* (-1)**non-int is NaN */
-		    } else if(yisint==1)
+		    } else if(yisint==1) 
 			z = -z;		/* (x<0)**odd = -(|x|**odd) */
 		}
 		return z;
 	    }
 	}
-
+    
     /* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be
 	n = (hx>>31)+1;
        but ANSI C says a right shift of a signed negative quantity is
@@ -181,8 +180,8 @@ __ieee754_pow(double x, double y)
     /* (x<0)**(non-int) is NaN */
 	if((n|yisint)==0) return (x-x)/(x-x);
 
-	sn = one; /* s (sign of result -ve**odd) = -1 else = 1 */
-	if((n|(yisint-1))==0) sn = -one;/* (-ve)**(odd int) */
+	s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
+	if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */
 
     /* |y| is huge */
 	if(iy>0x41e00000) { /* if |y| > 2**31 */
@@ -191,11 +190,11 @@ __ieee754_pow(double x, double y)
 		if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
 	    }
 	/* over/underflow if x is not close to one */
-	    if(ix<0x3fefffff) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
-	    if(ix>0x3ff00000) return (hy>0)? sn*huge*huge:sn*tiny*tiny;
-	/* now |1-x| is tiny <= 2**-20, suffice to compute
+	    if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny;
+	    if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*tiny*tiny;
+	/* now |1-x| is tiny <= 2**-20, suffice to compute 
 	   log(x) by x-x^2/2+x^3/3-x^4/4 */
-	    t = ax-1;		/* t has 20 trailing zeros */
+	    t = ax-one;		/* t has 20 trailing zeros */
 	    w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
 	    u = ivln2_h*t;	/* ivln2_h has 21 sig. bits */
 	    v = t*ivln2_l-w*ivln2;
@@ -203,7 +202,7 @@ __ieee754_pow(double x, double y)
 	    SET_LOW_WORD(t1,0);
 	    t2 = v-(t1-u);
 	} else {
-	    double s2,s_h,s_l,t_h,t_l;
+	    double ss,s2,s_h,s_l,t_h,t_l;
 	    n = 0;
 	/* take care subnormal number */
 	    if(ix<0x00100000)
@@ -217,11 +216,11 @@ __ieee754_pow(double x, double y)
 	    else {k=0;n+=1;ix -= 0x00100000;}
 	    SET_HIGH_WORD(ax,ix);
 
-	/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
+	/* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
 	    u = ax-bp[k];		/* bp[0]=1.0, bp[1]=1.5 */
 	    v = one/(ax+bp[k]);
-	    s = u*v;
-	    s_h = s;
+	    ss = u*v;
+	    s_h = ss;
 	    SET_LOW_WORD(s_h,0);
 	/* t_h=ax+bp[k] High */
 	    t_h = zero;
@@ -229,23 +228,23 @@ __ieee754_pow(double x, double y)
 	    t_l = ax - (t_h-bp[k]);
 	    s_l = v*((u-s_h*t_h)-s_h*t_l);
 	/* compute log(ax) */
-	    s2 = s*s;
+	    s2 = ss*ss;
 	    r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
-	    r += s_l*(s_h+s);
+	    r += s_l*(s_h+ss);
 	    s2  = s_h*s_h;
 	    t_h = 3.0+s2+r;
 	    SET_LOW_WORD(t_h,0);
 	    t_l = r-((t_h-3.0)-s2);
-	/* u+v = s*(1+...) */
+	/* u+v = ss*(1+...) */
 	    u = s_h*t_h;
-	    v = s_l*t_h+t_l*s;
-	/* 2/(3log2)*(s+...) */
+	    v = s_l*t_h+t_l*ss;
+	/* 2/(3log2)*(ss+...) */
 	    p_h = u+v;
 	    SET_LOW_WORD(p_h,0);
 	    p_l = v-(p_h-u);
 	    z_h = cp_h*p_h;		/* cp_h+cp_l = 2/(3*log2) */
 	    z_l = cp_l*p_h+p_l*cp+dp_l[k];
-	/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
+	/* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */
 	    t = (double)n;
 	    t1 = (((z_h+z_l)+dp_h[k])+t);
 	    SET_LOW_WORD(t1,0);
@@ -261,15 +260,15 @@ __ieee754_pow(double x, double y)
 	EXTRACT_WORDS(j,i,z);
 	if (j>=0x40900000) {				/* z >= 1024 */
 	    if(((j-0x40900000)|i)!=0)			/* if z > 1024 */
-		return sn*huge*huge;			/* overflow */
+		return s*huge*huge;			/* overflow */
 	    else {
-		if(p_l+ovt>z-p_h) return sn*huge*huge;	/* overflow */
+		if(p_l+ovt>z-p_h) return s*huge*huge;	/* overflow */
 	    }
 	} else if((j&0x7fffffff)>=0x4090cc00 ) {	/* z <= -1075 */
-	    if(((j-0xc090cc00)|i)!=0) 			/* z < -1075 */
-		return sn*tiny*tiny;			/* underflow */
+	    if(((j-0xc090cc00)|i)!=0) 		/* z < -1075 */
+		return s*tiny*tiny;		/* underflow */
 	    else {
-		if(p_l<=z-p_h) return sn*tiny*tiny;	/* underflow */
+		if(p_l<=z-p_h) return s*tiny*tiny;	/* underflow */
 	    }
 	}
     /*
@@ -286,7 +285,7 @@ __ieee754_pow(double x, double y)
 	    n = ((n&0x000fffff)|0x00100000)>>(20-k);
 	    if(j<0) n = -n;
 	    p_h -= t;
-	}
+	} 
 	t = p_l+p_h;
 	SET_LOW_WORD(t,0);
 	u = t*lg2_h;
@@ -301,5 +300,5 @@ __ieee754_pow(double x, double y)
 	j += (n<<20);
 	if((j>>20)<=0) z = scalbn(z,n);	/* subnormal output */
 	else SET_HIGH_WORD(z,j);
-	return sn*z;
+	return s*z;
 }
diff --git a/lib/msun/src/e_rem_pio2.c b/lib/msun/src/e_rem_pio2.c
index 668c1ea..5b4a08d 100644
--- a/lib/msun/src/e_rem_pio2.c
+++ b/lib/msun/src/e_rem_pio2.c
@@ -1,13 +1,15 @@
-/* @(#)e_rem_pio2.c 5.1 93/09/24 */
+
+/* @(#)e_rem_pio2.c 1.4 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
+ *
  */
 
 #ifndef lint
@@ -15,8 +17,8 @@ static char rcsid[] = "$FreeBSD$";
 #endif
 
 /* __ieee754_rem_pio2(x,y)
- *
- * return the remainder of x rem pi/2 in y[0]+y[1]
+ * 
+ * return the remainder of x rem pi/2 in y[0]+y[1] 
  * use __kernel_rem_pio2()
  */
 
@@ -24,20 +26,20 @@ static char rcsid[] = "$FreeBSD$";
 #include "math_private.h"
 
 /*
- * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
+ * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 
  */
 static const int32_t two_over_pi[] = {
-0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
-0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
-0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
-0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
-0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
-0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
-0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
-0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
-0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
-0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
-0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
+0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 
+0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 
+0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 
+0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 
+0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 
+0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 
+0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 
+0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 
+0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 
+0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 
+0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 
 };
 
 static const int32_t npio2_hw[] = {
@@ -83,7 +85,7 @@ pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 	if(ix<=0x3fe921fb)   /* |x| ~<= pi/4 , no need for reduction */
 	    {y[0] = x; y[1] = 0; return 0;}
 	if(ix<0x4002d97c) {  /* |x| < 3pi/4, special case with n=+-1 */
-	    if(hx>0) {
+	    if(hx>0) { 
 		z = x - pio2_1;
 		if(ix!=0x3ff921fb) { 	/* 33+53 bit pi is good enough */
 		    y[0] = z - pio2_1t;
@@ -113,27 +115,27 @@ pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 	    fn = (double)n;
 	    r  = t-fn*pio2_1;
 	    w  = fn*pio2_1t;	/* 1st round good to 85 bit */
-	    if(n<32&&ix!=npio2_hw[n-1]) {
+	    if(n<32&&ix!=npio2_hw[n-1]) {	
 		y[0] = r-w;	/* quick check no cancellation */
 	    } else {
 	        u_int32_t high;
 	        j  = ix>>20;
-	        y[0] = r-w;
+	        y[0] = r-w; 
 		GET_HIGH_WORD(high,y[0]);
 	        i = j-((high>>20)&0x7ff);
 	        if(i>16) {  /* 2nd iteration needed, good to 118 */
 		    t  = r;
-		    w  = fn*pio2_2;
+		    w  = fn*pio2_2;	
 		    r  = t-w;
-		    w  = fn*pio2_2t-((t-r)-w);
+		    w  = fn*pio2_2t-((t-r)-w);	
 		    y[0] = r-w;
 		    GET_HIGH_WORD(high,y[0]);
 		    i = j-((high>>20)&0x7ff);
 		    if(i>49)  {	/* 3rd iteration need, 151 bits acc */
 		    	t  = r;	/* will cover all possible cases */
-		    	w  = fn*pio2_3;
+		    	w  = fn*pio2_3;	
 		    	r  = t-w;
-		    	w  = fn*pio2_3t-((t-r)-w);
+		    	w  = fn*pio2_3t-((t-r)-w);	
 		    	y[0] = r-w;
 		    }
 		}
@@ -142,7 +144,7 @@ pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 	    if(hx<0) 	{y[0] = -y[0]; y[1] = -y[1]; return -n;}
 	    else	 return n;
 	}
-    /*
+    /* 
      * all other (large) arguments
      */
 	if(ix>=0x7ff00000) {		/* x is inf or NaN */
diff --git a/lib/msun/src/e_remainder.c b/lib/msun/src/e_remainder.c
index ba226b7..1cea6c1 100644
--- a/lib/msun/src/e_remainder.c
+++ b/lib/msun/src/e_remainder.c
@@ -1,11 +1,12 @@
-/* @(#)e_remainder.c 5.1 93/09/24 */
+
+/* @(#)e_remainder.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -15,11 +16,11 @@ static char rcsid[] = "$FreeBSD$";
 #endif
 
 /* __ieee754_remainder(x,p)
- * Return :
- * 	returns  x REM p  =  x - [x/p]*p as if in infinite
- * 	precise arithmetic, where [x/p] is the (infinite bit)
+ * Return :                  
+ * 	returns  x REM p  =  x - [x/p]*p as if in infinite 
+ * 	precise arithmetic, where [x/p] is the (infinite bit) 
  *	integer nearest x/p (in half way case choose the even one).
- * Method :
+ * Method : 
  *	Based on fmod() return x-[x/p]chopped*p exactlp.
  */
 
diff --git a/lib/msun/src/e_scalb.c b/lib/msun/src/e_scalb.c
index 6553d50..f7d46d0 100644
--- a/lib/msun/src/e_scalb.c
+++ b/lib/msun/src/e_scalb.c
@@ -1,11 +1,12 @@
-/* @(#)e_scalb.c 5.1 93/09/24 */
+
+/* @(#)e_scalb.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -16,7 +17,7 @@ static char rcsid[] = "$FreeBSD$";
 
 /*
  * __ieee754_scalb(x, fn) is provide for
- * passing various standard test suite. One
+ * passing various standard test suite. One 
  * should use scalbn() instead.
  */
 
@@ -34,7 +35,7 @@ __ieee754_scalb(double x, double fn)
 #ifdef _SCALB_INT
 	return scalbn(x,fn);
 #else
-	if ((isnan)(x)||(isnan)(fn)) return x*fn;
+	if (isnan(x)||isnan(fn)) return x*fn;
 	if (!finite(fn)) {
 	    if(fn>0.0) return x*fn;
 	    else       return x/(-fn);
diff --git a/lib/msun/src/e_sinh.c b/lib/msun/src/e_sinh.c
index 76c4576..09b0597 100644
--- a/lib/msun/src/e_sinh.c
+++ b/lib/msun/src/e_sinh.c
@@ -1,11 +1,12 @@
-/* @(#)e_sinh.c 5.1 93/09/24 */
+
+/* @(#)e_sinh.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -15,15 +16,15 @@ static char rcsid[] = "$FreeBSD$";
 #endif
 
 /* __ieee754_sinh(x)
- * Method :
+ * Method : 
  * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
- *	1. Replace x by |x| (sinh(-x) = -sinh(x)).
- *	2.
+ *	1. Replace x by |x| (sinh(-x) = -sinh(x)). 
+ *	2. 
  *		                                    E + E/(E+1)
  *	    0        <= x <= 22     :  sinh(x) := --------------, E=expm1(x)
  *			       			        2
  *
- *	    22       <= x <= lnovft :  sinh(x) := exp(x)/2
+ *	    22       <= x <= lnovft :  sinh(x) := exp(x)/2 
  *	    lnovft   <= x <= ln2ovft:  sinh(x) := exp(x/2)/2 * exp(x/2)
  *	    ln2ovft  <  x	    :  sinh(x) := x*shuge (overflow)
  *
@@ -49,7 +50,7 @@ __ieee754_sinh(double x)
 	ix = jx&0x7fffffff;
 
     /* x is INF or NaN */
-	if(ix>=0x7ff00000) return x+x;
+	if(ix>=0x7ff00000) return x+x;	
 
 	h = 0.5;
 	if (jx<0) h = -h;
diff --git a/lib/msun/src/e_sqrt.c b/lib/msun/src/e_sqrt.c
index 5aeadee..3b3b2a7 100644
--- a/lib/msun/src/e_sqrt.c
+++ b/lib/msun/src/e_sqrt.c
@@ -1,11 +1,12 @@
-/* @(#)e_sqrt.c 5.1 93/09/24 */
+
+/* @(#)e_sqrt.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -19,10 +20,10 @@ static char rcsid[] = "$FreeBSD$";
  *           ------------------------------------------
  *	     |  Use the hardware sqrt if you have one |
  *           ------------------------------------------
- * Method:
- *   Bit by bit method using integer arithmetic. (Slow, but portable)
+ * Method: 
+ *   Bit by bit method using integer arithmetic. (Slow, but portable) 
  *   1. Normalization
- *	Scale x to y in [1,4) with even powers of 2:
+ *	Scale x to y in [1,4) with even powers of 2: 
  *	find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
  *		sqrt(x) = 2^k * sqrt(y)
  *   2. Bit by bit computation
@@ -31,9 +32,9 @@ static char rcsid[] = "$FreeBSD$";
  *                                     i+1         2
  *	    s  = 2*q , and	y  =  2   * ( y - q  ).		(1)
  *	     i      i            i                 i
- *
- *	To compute q    from q , one checks whether
- *		    i+1       i
+ *                                                        
+ *	To compute q    from q , one checks whether 
+ *		    i+1       i                       
  *
  *			      -(i+1) 2
  *			(q + 2      ) <= y.			(2)
@@ -43,12 +44,12 @@ static char rcsid[] = "$FreeBSD$";
  *		 	       i+1   i             i+1   i
  *
  *	With some algebric manipulation, it is not difficult to see
- *	that (2) is equivalent to
+ *	that (2) is equivalent to 
  *                             -(i+1)
  *			s  +  2       <= y			(3)
  *			 i                i
  *
- *	The advantage of (3) is that s  and y  can be computed by
+ *	The advantage of (3) is that s  and y  can be computed by 
  *				      i      i
  *	the following recurrence formula:
  *	    if (3) is false
@@ -60,10 +61,10 @@ static char rcsid[] = "$FreeBSD$";
  *                         -i                     -(i+1)
  *	    s	  =  s  + 2  ,  y    = y  -  s  - 2  		(5)
  *           i+1      i          i+1    i     i
- *
- *	One may easily use induction to prove (4) and (5).
+ *				
+ *	One may easily use induction to prove (4) and (5). 
  *	Note. Since the left hand side of (3) contain only i+2 bits,
- *	      it does not necessary to do a full (53-bit) comparison
+ *	      it does not necessary to do a full (53-bit) comparison 
  *	      in (3).
  *   3. Final rounding
  *	After generating the 53 bits result, we compute one more bit.
@@ -73,7 +74,7 @@ static char rcsid[] = "$FreeBSD$";
  *	The rounding mode can be detected by checking whether
  *	huge + tiny is equal to huge, and whether huge - tiny is
  *	equal to huge for some floating point number "huge" and "tiny".
- *
+ *		
  * Special cases:
  *	sqrt(+-0) = +-0 	... exact
  *	sqrt(inf) = inf
@@ -100,10 +101,10 @@ __ieee754_sqrt(double x)
 	EXTRACT_WORDS(ix0,ix1,x);
 
     /* take care of Inf and NaN */
-	if((ix0&0x7ff00000)==0x7ff00000) {
+	if((ix0&0x7ff00000)==0x7ff00000) {			
 	    return x*x+x;		/* sqrt(NaN)=NaN, sqrt(+inf)=+inf
 					   sqrt(-inf)=sNaN */
-	}
+	} 
     /* take care of zero */
 	if(ix0<=0) {
 	    if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
@@ -137,12 +138,12 @@ __ieee754_sqrt(double x)
 	r = 0x00200000;		/* r = moving bit from right to left */
 
 	while(r!=0) {
-	    t = s0+r;
-	    if(t<=ix0) {
-		s0   = t+r;
-		ix0 -= t;
-		q   += r;
-	    }
+	    t = s0+r; 
+	    if(t<=ix0) { 
+		s0   = t+r; 
+		ix0 -= t; 
+		q   += r; 
+	    } 
 	    ix0 += ix0 + ((ix1&sign)>>31);
 	    ix1 += ix1;
 	    r>>=1;
@@ -150,9 +151,9 @@ __ieee754_sqrt(double x)
 
 	r = sign;
 	while(r!=0) {
-	    t1 = s1+r;
+	    t1 = s1+r; 
 	    t  = s0;
-	    if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
+	    if((t<ix0)||((t==ix0)&&(t1<=ix1))) { 
 		s1  = t1+r;
 		if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
 		ix0 -= t;
@@ -173,7 +174,7 @@ __ieee754_sqrt(double x)
 	        if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
 		else if (z>one) {
 		    if (q1==(u_int32_t)0xfffffffe) q+=1;
-		    q1+=2;
+		    q1+=2; 
 		} else
 	            q1 += (q1&1);
 	    }
@@ -189,18 +190,18 @@ __ieee754_sqrt(double x)
 /*
 Other methods  (use floating-point arithmetic)
 -------------
-(This is a copy of a drafted paper by Prof W. Kahan
+(This is a copy of a drafted paper by Prof W. Kahan 
 and K.C. Ng, written in May, 1986)
 
-	Two algorithms are given here to implement sqrt(x)
+	Two algorithms are given here to implement sqrt(x) 
 	(IEEE double precision arithmetic) in software.
 	Both supply sqrt(x) correctly rounded. The first algorithm (in
 	Section A) uses newton iterations and involves four divisions.
 	The second one uses reciproot iterations to avoid division, but
 	requires more multiplications. Both algorithms need the ability
-	to chop results of arithmetic operations instead of round them,
+	to chop results of arithmetic operations instead of round them, 
 	and the INEXACT flag to indicate when an arithmetic operation
-	is executed exactly with no roundoff error, all part of the
+	is executed exactly with no roundoff error, all part of the 
 	standard (IEEE 754-1985). The ability to perform shift, add,
 	subtract and logical AND operations upon 32-bit words is needed
 	too, though not part of the standard.
@@ -210,7 +211,7 @@ A.  sqrt(x) by Newton Iteration
    (1)	Initial approximation
 
 	Let x0 and x1 be the leading and the trailing 32-bit words of
-	a floating point number x (in IEEE double format) respectively
+	a floating point number x (in IEEE double format) respectively 
 
 	    1    11		     52				  ...widths
 	   ------------------------------------------------------
@@ -218,7 +219,7 @@ A.  sqrt(x) by Newton Iteration
 	   ------------------------------------------------------
 	      msb    lsb  msb				      lsb ...order
 
-
+ 
 	     ------------------------  	     ------------------------
 	x0:  |s|   e    |    f1     |	 x1: |          f2           |
 	     ------------------------  	     ------------------------
@@ -243,7 +244,7 @@ A.  sqrt(x) by Newton Iteration
 
     (2)	Iterative refinement
 
-	Apply Heron's rule three times to y, we have y approximates
+	Apply Heron's rule three times to y, we have y approximates 
 	sqrt(x) to within 1 ulp (Unit in the Last Place):
 
 		y := (y+x/y)/2		... almost 17 sig. bits
@@ -268,12 +269,12 @@ A.  sqrt(x) by Newton Iteration
 	it requires more multiplications and additions. Also x must be
 	scaled in advance to avoid spurious overflow in evaluating the
 	expression 3y*y+x. Hence it is not recommended uless division
-	is slow. If division is very slow, then one should use the
+	is slow. If division is very slow, then one should use the 
 	reciproot algorithm given in section B.
 
     (3) Final adjustment
 
-	By twiddling y's last bit it is possible to force y to be
+	By twiddling y's last bit it is possible to force y to be 
 	correctly rounded according to the prevailing rounding mode
 	as follows. Let r and i be copies of the rounding mode and
 	inexact flag before entering the square root program. Also we
@@ -304,7 +305,7 @@ A.  sqrt(x) by Newton Iteration
 	        I := i;	 		... restore inexact flag
 	        R := r;  		... restore rounded mode
 	        return sqrt(x):=y.
-
+		    
     (4)	Special cases
 
 	Square root of +inf, +-0, or NaN is itself;
@@ -323,7 +324,7 @@ B.  sqrt(x) by Reciproot Iteration
 	    k := 0x5fe80000 - (x0>>1);
 	    y0:= k - T2[63&(k>>14)].	... y ~ 1/sqrt(x) to 7.8 bits
 
-	Here k is a 32-bit integer and T2[] is an integer array
+	Here k is a 32-bit integer and T2[] is an integer array 
 	containing correction terms. Now magically the floating
 	value of y (y's leading 32-bit word is y0, the value of
 	its trailing word y1 is set to zero) approximates 1/sqrt(x)
@@ -344,9 +345,9 @@ B.  sqrt(x) by Reciproot Iteration
 
 	Apply Reciproot iteration three times to y and multiply the
 	result by x to get an approximation z that matches sqrt(x)
-	to about 1 ulp. To be exact, we will have
+	to about 1 ulp. To be exact, we will have 
 		-1ulp < sqrt(x)-z<1.0625ulp.
-
+	
 	... set rounding mode to Round-to-nearest
 	   y := y*(1.5-0.5*x*y*y)	... almost 15 sig. bits to 1/sqrt(x)
 	   y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
@@ -355,14 +356,14 @@ B.  sqrt(x) by Reciproot Iteration
 	   z := z + 0.5*z*(1-z*y)	... about 1 ulp to sqrt(x)
 
 	Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
-	(a) the term z*y in the final iteration is always less than 1;
+	(a) the term z*y in the final iteration is always less than 1; 
 	(b) the error in the final result is biased upward so that
 		-1 ulp < sqrt(x) - z < 1.0625 ulp
 	    instead of |sqrt(x)-z|<1.03125ulp.
 
     (3)	Final adjustment
 
-	By twiddling y's last bit it is possible to force y to be
+	By twiddling y's last bit it is possible to force y to be 
 	correctly rounded according to the prevailing rounding mode
 	as follows. Let r and i be copies of the rounding mode and
 	inexact flag before entering the square root program. Also we
@@ -402,27 +403,27 @@ B.  sqrt(x) by Reciproot Iteration
 	    I := 1;		... Raise Inexact flag: z is not exact
 	else {
 	    j := 1 - [(x0>>20)&1]	... j = logb(x) mod 2
-	    k := z1 >> 26;		... get z's 25-th and 26-th
+	    k := z1 >> 26;		... get z's 25-th and 26-th 
 					    fraction bits
 	    I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
 	}
 	R:= r		... restore rounded mode
 	return sqrt(x):=z.
 
-	If multiplication is cheaper then the foregoing red tape, the
+	If multiplication is cheaper then the foregoing red tape, the 
 	Inexact flag can be evaluated by
 
 	    I := i;
 	    I := (z*z!=x) or I.
 
-	Note that z*z can overwrite I; this value must be sensed if it is
+	Note that z*z can overwrite I; this value must be sensed if it is 
 	True.
 
 	Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
 	zero.
 
 		    --------------------
-		z1: |        f2        |
+		z1: |        f2        | 
 		    --------------------
 		bit 31		   bit 0
 
@@ -439,7 +440,7 @@ B.  sqrt(x) by Reciproot Iteration
 	11			01		even
 	-------------------------------------------------
 
-    (4)	Special cases (see (4) of Section A).
-
+    (4)	Special cases (see (4) of Section A).	
+ 
  */
-
+ 
diff --git a/lib/msun/src/k_cos.c b/lib/msun/src/k_cos.c
index f8fe2fd..ad78f8e 100644
--- a/lib/msun/src/k_cos.c
+++ b/lib/msun/src/k_cos.c
@@ -1,11 +1,12 @@
-/* @(#)k_cos.c 5.1 93/09/24 */
+
+/* @(#)k_cos.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -18,7 +19,7 @@ static char rcsid[] = "$FreeBSD$";
  * __kernel_cos( x,  y )
  * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
  * Input x is assumed to be bounded by ~pi/4 in magnitude.
- * Input y is the tail of x.
+ * Input y is the tail of x. 
  *
  * Algorithm
  *	1. Since cos(-x) = cos(x), we need only to consider positive x.
@@ -28,15 +29,15 @@ static char rcsid[] = "$FreeBSD$";
  *		  	                 4            14
  *	   	cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
  *	   where the remez error is
- *
+ *	
  * 	|              2     4     6     8     10    12     14 |     -58
  * 	|cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  )| <= 2
- * 	|    					               |
- *
- * 	               4     6     8     10    12     14
+ * 	|    					               | 
+ * 
+ * 	               4     6     8     10    12     14 
  *	4. let r = C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  , then
  *	       cos(x) = 1 - x*x/2 + r
- *	   since cos(x+y) ~ cos(x) - sin(x)*y
+ *	   since cos(x+y) ~ cos(x) - sin(x)*y 
  *			  ~ cos(x) - x*y,
  *	   a correction term is necessary in cos(x) and hence
  *		cos(x+y) = 1 - (x*x/2 - (r - x*y))
@@ -73,7 +74,7 @@ __kernel_cos(double x, double y)
 	}
 	z  = x*x;
 	r  = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
-	if(ix < 0x3FD33333) 			/* if |x| < 0.3 */
+	if(ix < 0x3FD33333) 			/* if |x| < 0.3 */ 
 	    return one - (0.5*z - (z*r - x*y));
 	else {
 	    if(ix > 0x3fe90000) {		/* x > 0.78125 */
diff --git a/lib/msun/src/k_rem_pio2.c b/lib/msun/src/k_rem_pio2.c
index 3613c38..61f7dce 100644
--- a/lib/msun/src/k_rem_pio2.c
+++ b/lib/msun/src/k_rem_pio2.c
@@ -1,11 +1,12 @@
-/* @(#)k_rem_pio2.c 5.1 93/09/24 */
+
+/* @(#)k_rem_pio2.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -17,12 +18,12 @@ static char rcsid[] = "$FreeBSD$";
 /*
  * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
  * double x[],y[]; int e0,nx,prec; int ipio2[];
- *
- * __kernel_rem_pio2 return the last three digits of N with
+ * 
+ * __kernel_rem_pio2 return the last three digits of N with 
  *		y = x - N*pi/2
  * so that |y| < pi/2.
  *
- * The method is to compute the integer (mod 8) and fraction parts of
+ * The method is to compute the integer (mod 8) and fraction parts of 
  * (2/pi)*x without doing the full multiplication. In general we
  * skip the part of the product that are known to be a huge integer (
  * more accurately, = 0 mod 8 ). Thus the number of operations are
@@ -31,10 +32,10 @@ static char rcsid[] = "$FreeBSD$";
  * (2/pi) is represented by an array of 24-bit integers in ipio2[].
  *
  * Input parameters:
- * 	x[]	The input value (must be positive) is broken into nx
+ * 	x[]	The input value (must be positive) is broken into nx 
  *		pieces of 24-bit integers in double precision format.
- *		x[i] will be the i-th 24 bit of x. The scaled exponent
- *		of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
+ *		x[i] will be the i-th 24 bit of x. The scaled exponent 
+ *		of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 
  *		match x's up to 24 bits.
  *
  *		Example of breaking a double positive z into x[0]+x[1]+x[2]:
@@ -71,8 +72,8 @@ static char rcsid[] = "$FreeBSD$";
  *			3	113 bits (quad)
  *
  *	ipio2[]
- *		integer array, contains the (24*i)-th to (24*i+23)-th
- *		bit of 2/pi after binary point. The corresponding
+ *		integer array, contains the (24*i)-th to (24*i+23)-th 
+ *		bit of 2/pi after binary point. The corresponding 
  *		floating value is
  *
  *			ipio2[i] * 2^(-24(i+1)).
@@ -87,8 +88,8 @@ static char rcsid[] = "$FreeBSD$";
  *		in the computation. The recommended value is 2,3,4,
  *		6 for single, double, extended,and quad.
  *
- * 	jz	local integer variable indicating the number of
- *		terms of ipio2[] used.
+ * 	jz	local integer variable indicating the number of 
+ *		terms of ipio2[] used. 
  *
  *	jx	nx - 1
  *
@@ -108,9 +109,9 @@ static char rcsid[] = "$FreeBSD$";
  *		exponent for q[i] would be q0-24*i.
  *
  *	PIo2[]	double precision array, obtained by cutting pi/2
- *		into 24 bits chunks.
+ *		into 24 bits chunks. 
  *
- *	f[]	ipio2[] in floating point
+ *	f[]	ipio2[] in floating point 
  *
  *	iq[]	integer array by breaking up q[] in 24-bits chunk.
  *
@@ -124,9 +125,9 @@ static char rcsid[] = "$FreeBSD$";
 
 /*
  * Constants:
- * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough
+ * The hexadecimal values are the intended ones for the following 
+ * constants. The decimal values may be used, provided that the 
+ * compiler will convert from decimal to binary accurately enough 
  * to produce the hexadecimal values shown.
  */
 
@@ -146,7 +147,7 @@ static const double PIo2[] = {
   2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
 };
 
-static const double
+static const double			
 zero   = 0.0,
 one    = 1.0,
 two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
@@ -194,7 +195,7 @@ recompute:
 	    i  = (iq[jz-1]>>(24-q0)); n += i;
 	    iq[jz-1] -= i<<(24-q0);
 	    ih = iq[jz-1]>>(23-q0);
-	}
+	} 
 	else if(q0==0) ih = iq[jz-1]>>23;
 	else if(z>=0.5) ih=2;
 
@@ -245,7 +246,7 @@ recompute:
 	    while(iq[jz]==0) { jz--; q0-=24;}
 	} else { /* break z into 24-bit if necessary */
 	    z = scalbn(z,-q0);
-	    if(z>=two24) {
+	    if(z>=two24) { 
 		fw = (double)((int32_t)(twon24*z));
 		iq[jz] = (int32_t)(z-two24*fw);
 		jz += 1; q0 += 24;
@@ -270,29 +271,29 @@ recompute:
 	    case 0:
 		fw = 0.0;
 		for (i=jz;i>=0;i--) fw += fq[i];
-		y[0] = (ih==0)? fw: -fw;
+		y[0] = (ih==0)? fw: -fw; 
 		break;
 	    case 1:
 	    case 2:
 		fw = 0.0;
-		for (i=jz;i>=0;i--) fw += fq[i];
-		y[0] = (ih==0)? fw: -fw;
+		for (i=jz;i>=0;i--) fw += fq[i]; 
+		y[0] = (ih==0)? fw: -fw; 
 		fw = fq[0]-fw;
 		for (i=1;i<=jz;i++) fw += fq[i];
-		y[1] = (ih==0)? fw: -fw;
+		y[1] = (ih==0)? fw: -fw; 
 		break;
 	    case 3:	/* painful */
 		for (i=jz;i>0;i--) {
-		    fw      = fq[i-1]+fq[i];
+		    fw      = fq[i-1]+fq[i]; 
 		    fq[i]  += fq[i-1]-fw;
 		    fq[i-1] = fw;
 		}
 		for (i=jz;i>1;i--) {
-		    fw      = fq[i-1]+fq[i];
+		    fw      = fq[i-1]+fq[i]; 
 		    fq[i]  += fq[i-1]-fw;
 		    fq[i-1] = fw;
 		}
-		for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
+		for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; 
 		if(ih==0) {
 		    y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
 		} else {
diff --git a/lib/msun/src/k_sin.c b/lib/msun/src/k_sin.c
index 52d3faa..55e2bdc 100644
--- a/lib/msun/src/k_sin.c
+++ b/lib/msun/src/k_sin.c
@@ -1,11 +1,12 @@
-/* @(#)k_sin.c 5.1 93/09/24 */
+
+/* @(#)k_sin.c 1.3 95/01/18 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
  */
@@ -18,24 +19,24 @@ static char rcsid[] = "$FreeBSD$";
  * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
  * Input x is assumed to be bounded by ~pi/4 in magnitude.
  * Input y is the tail of x.
- * Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
+ * Input iy indicates whether y is 0. (if iy=0, y assume to be 0). 
  *
  * Algorithm
- *	1. Since sin(-x) = -sin(x), we need only to consider positive x.
+ *	1. Since sin(-x) = -sin(x), we need only to consider positive x. 
  *	2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
  *	3. sin(x) is approximated by a polynomial of degree 13 on
  *	   [0,pi/4]
  *		  	         3            13
  *	   	sin(x) ~ x + S1*x + ... + S6*x
  *	   where
- *
+ *	
  * 	|sin(x)         2     4     6     8     10     12  |     -58
  * 	|----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
- * 	|  x 					           |
- *
+ * 	|  x 					           | 
+ * 
  *	4. sin(x+y) = sin(x) + sin'(x')*y
  *		    ~ sin(x) + (1-x*x/2)*y
- *	   For better accuracy, let
+ *	   For better accuracy, let 
  *		     3      2      2      2      2
  *		r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
  *	   then                   3    2
diff --git a/lib/msun/src/k_tan.c b/lib/msun/src/k_tan.c
index b7cdf35..16b2dd0 100644
--- a/lib/msun/src/k_tan.c
+++ b/lib/msun/src/k_tan.c
@@ -1,4 +1,5 @@
-/* @(#)k_tan.c 5.1 93/09/24 */
+/* @(#)k_tan.c 1.5 04/04/22 SMI */
+
 /*
  * ====================================================
  * Copyright 2004 Sun Microsystems, Inc.  All Rights Reserved.
@@ -9,6 +10,7 @@
  * ====================================================
  */
 
+/* INDENT OFF */
 #ifndef lint
 static char rcsid[] = "$FreeBSD$";
 #endif
@@ -17,13 +19,12 @@ static char rcsid[] = "$FreeBSD$";
  * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
  * Input x is assumed to be bounded by ~pi/4 in magnitude.
  * Input y is the tail of x.
- * Input k indicates whether tan (if k=1) or
- * -1/tan (if k= -1) is returned.
+ * Input k indicates whether tan (if k = 1) or -1/tan (if k = -1) is returned.
  *
  * Algorithm
  *	1. Since tan(-x) = -tan(x), we need only to consider positive x.
  *	2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
- *	3. tan(x) is approximated by an odd polynomial of degree 27 on
+ *	3. tan(x) is approximated by a odd polynomial of degree 27 on
  *	   [0,0.67434]
  *		  	         3             27
  *	   	tan(x) ~ x + T1*x + ... + T13*x
@@ -49,34 +50,38 @@ static char rcsid[] = "$FreeBSD$";
 
 #include "math.h"
 #include "math_private.h"
-static const double
-one   =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
-pio4  =  7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
-pio4lo=  3.06161699786838301793e-17, /* 0x3C81A626, 0x33145C07 */
-T[] =  {
-  3.33333333333334091986e-01, /* 0x3FD55555, 0x55555563 */
-  1.33333333333201242699e-01, /* 0x3FC11111, 0x1110FE7A */
-  5.39682539762260521377e-02, /* 0x3FABA1BA, 0x1BB341FE */
-  2.18694882948595424599e-02, /* 0x3F9664F4, 0x8406D637 */
-  8.86323982359930005737e-03, /* 0x3F8226E3, 0xE96E8493 */
-  3.59207910759131235356e-03, /* 0x3F6D6D22, 0xC9560328 */
-  1.45620945432529025516e-03, /* 0x3F57DBC8, 0xFEE08315 */
-  5.88041240820264096874e-04, /* 0x3F4344D8, 0xF2F26501 */
-  2.46463134818469906812e-04, /* 0x3F3026F7, 0x1A8D1068 */
-  7.81794442939557092300e-05, /* 0x3F147E88, 0xA03792A6 */
-  7.14072491382608190305e-05, /* 0x3F12B80F, 0x32F0A7E9 */
- -1.85586374855275456654e-05, /* 0xBEF375CB, 0xDB605373 */
-  2.59073051863633712884e-05, /* 0x3EFB2A70, 0x74BF7AD4 */
+static const double xxx[] = {
+		 3.33333333333334091986e-01,	/* 3FD55555, 55555563 */
+		 1.33333333333201242699e-01,	/* 3FC11111, 1110FE7A */
+		 5.39682539762260521377e-02,	/* 3FABA1BA, 1BB341FE */
+		 2.18694882948595424599e-02,	/* 3F9664F4, 8406D637 */
+		 8.86323982359930005737e-03,	/* 3F8226E3, E96E8493 */
+		 3.59207910759131235356e-03,	/* 3F6D6D22, C9560328 */
+		 1.45620945432529025516e-03,	/* 3F57DBC8, FEE08315 */
+		 5.88041240820264096874e-04,	/* 3F4344D8, F2F26501 */
+		 2.46463134818469906812e-04,	/* 3F3026F7, 1A8D1068 */
+		 7.81794442939557092300e-05,	/* 3F147E88, A03792A6 */
+		 7.14072491382608190305e-05,	/* 3F12B80F, 32F0A7E9 */
+		-1.85586374855275456654e-05,	/* BEF375CB, DB605373 */
+		 2.59073051863633712884e-05,	/* 3EFB2A70, 74BF7AD4 */
+/* one */	 1.00000000000000000000e+00,	/* 3FF00000, 00000000 */
+/* pio4 */	 7.85398163397448278999e-01,	/* 3FE921FB, 54442D18 */
+/* pio4lo */	 3.06161699786838301793e-17	/* 3C81A626, 33145C07 */
 };
+#define	one	xxx[13]
+#define	pio4	xxx[14]
+#define	pio4lo	xxx[15]
+#define	T	xxx
+/* INDENT ON */
 
 double
-__kernel_tan(double x, double y, int iy)
-{
-	double z,r,v,w,s;
-	int32_t ix,hx;
+__kernel_tan(double x, double y, int iy) {
+	double z, r, v, w, s;
+	int32_t ix, hx;
+
 	GET_HIGH_WORD(hx,x);
-	ix = hx&0x7fffffff;	/* high word of |x| */
-	if(ix<0x3e300000) {			/* x < 2**-28 */
+	ix = hx & 0x7fffffff;			/* high word of |x| */
+	if (ix < 0x3e300000) {			/* x < 2**-28 */
 		if ((int) x == 0) {		/* generate inexact */
 			u_int32_t low;
 			GET_LOW_WORD(low,x);
@@ -99,39 +104,51 @@ __kernel_tan(double x, double y, int iy)
 			}
 		}
 	}
-	if(ix>=0x3FE59428) { 			/* |x|>=0.6744 */
-	    if(hx<0) {x = -x; y = -y;}
-	    z = pio4-x;
-	    w = pio4lo-y;
-	    x = z+w; y = 0.0;
+	if (ix >= 0x3FE59428) {	/* |x| >= 0.6744 */
+		if (hx < 0) {
+			x = -x;
+			y = -y;
+		}
+		z = pio4 - x;
+		w = pio4lo - y;
+		x = z + w;
+		y = 0.0;
 	}
-	z	=  x*x;
-	w 	=  z*z;
-    /* Break x^5*(T[1]+x^2*T[2]+...) into
-     *	  x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
-     *	  x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
-     */
-	r = T[1]+w*(T[3]+w*(T[5]+w*(T[7]+w*(T[9]+w*T[11]))));
-	v = z*(T[2]+w*(T[4]+w*(T[6]+w*(T[8]+w*(T[10]+w*T[12])))));
-	s = z*x;
-	r = y + z*(s*(r+v)+y);
-	r += T[0]*s;
-	w = x+r;
-	if(ix>=0x3FE59428) {
-	    v = (double)iy;
-	    return (double)(1-((hx>>30)&2))*(v-2.0*(x-(w*w/(w+v)-r)));
+	z = x * x;
+	w = z * z;
+	/*
+	 * Break x^5*(T[1]+x^2*T[2]+...) into
+	 * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
+	 * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
+	 */
+	r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] +
+		w * T[11]))));
+	v = z * (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] +
+		w * T[12])))));
+	s = z * x;
+	r = y + z * (s * (r + v) + y);
+	r += T[0] * s;
+	w = x + r;
+	if (ix >= 0x3FE59428) {
+		v = (double) iy;
+		return (double) (1 - ((hx >> 30) & 2)) *
+			(v - 2.0 * (x - (w * w / (w + v) - r)));
 	}
-	if(iy==1) return w;
-	else {		/* if allow error up to 2 ulp,
-			   simply return -1.0/(x+r) here */
-     /*  compute -1.0/(x+r) accurately */
-	    double a,t;
-	    z  = w;
-	    SET_LOW_WORD(z,0);
-	    v  = r-(z - x); 	/* z+v = r+x */
-	    t = a  = -1.0/w;	/* a = -1.0/w */
-	    SET_LOW_WORD(t,0);
-	    s  = 1.0+t*z;
-	    return t+a*(s+t*v);
+	if (iy == 1)
+		return w;
+	else {
+		/*
+		 * if allow error up to 2 ulp, simply return
+		 * -1.0 / (x+r) here
+		 */
+		/* compute -1.0 / (x+r) accurately */
+		double a, t;
+		z = w;
+		SET_LOW_WORD(z,0);
+		v = r - (z - x);	/* z+v = r+x */
+		t = a = -1.0 / w;	/* a = -1.0/w */
+		SET_LOW_WORD(t,0);
+		s = 1.0 + t * z;
+		return t + a * (s + t * v);
 	}
 }
author	das <das@FreeBSD.org>	2005-02-04 18:26:06 +0000
committer	das <das@FreeBSD.org>	2005-02-04 18:26:06 +0000
commit	a47af911ba88a733157f219bcea374aafc0e6a34 (patch)
tree	358c37c6dfb32659eea19d2a72397142931bc562 /lib/msun/src
parent	e8403b23681c8ac8f2c86499943b43c359a4e301 (diff)
download	FreeBSD-src-a47af911ba88a733157f219bcea374aafc0e6a34.zip FreeBSD-src-a47af911ba88a733157f219bcea374aafc0e6a34.tar.gz