Specialize fdlibm to jerry

Keep IEEE code paths only: * removed SVID, XOPEN, POSIX code paths from everywhere. * deleted s_lib_version.c, as version is only useful if multiple standards are supported. * deleted k_standard.c, as it handles non-IEEE exception cases only. * renamed the e_{acos,asin,atan2,exp,fmod,log,pow,sqrt}.c sources as s_.*, dropped the __ieee754_ prefix from the names of the appropriate functions therein, and deleted the w_{acos,asin,atan2,exp,fmod,log,pow,sqrt}.c wrapper code. Keep C99 declaration variants only: * removed old C-style function declaration variants. * removed data declaration variants where const qualifier was not used. Clean unused sources/functions: * removed s_{rint,significand,tanh}.c and the appropriate functions defined therein. JerryScript-DCO-1.0-Signed-off-by: Akos Kiss akiss@inf.u-szeged.hu
2016-03-13 22:30:47 +01:00
parent d674b92f26
commit 397dff81ee
38 changed files with 266 additions and 1905 deletions
@@ -13,39 +13,31 @@
 */
 /* __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()
 */
 #include "fdlibm.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
 */
 #ifdef __STDC__
 static const int two_over_pi[] = {
-#else
+0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
-static int two_over_pi[] = {
+0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
-#endif
+0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
-0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 
+0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
-0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 
+0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
-0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 
+0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
-0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 
+0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
-0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 
+0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
-0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 
+0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
-0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 
+0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
-0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 
+0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 
 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 
 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 
 };
 #ifdef __STDC__
 static const int npio2_hw[] = {
 #else
 static int npio2_hw[] = {
 #endif
 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
@@ -64,11 +56,7 @@ static int npio2_hw[] = {
 * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
 */
-#ifdef __STDC__
+static const double
 static const double 
 #else
 static double 
 #endif
 zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
 half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
 two24 =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
@@ -80,12 +68,7 @@ pio2_2t =  2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
 pio2_3  =  2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
 pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
-#ifdef __STDC__
+int __ieee754_rem_pio2(double x, double *y)
 	int __ieee754_rem_pio2(double x, double *y)
 #else
 	int __ieee754_rem_pio2(x,y)
 	double x,y[];
 #endif
 {
 	double z,w,t,r,fn;
 	double tx[3];
@@ -126,7 +109,7 @@ 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 {
 	        j  = ix>>20;
@@ -134,16 +117,16 @@ pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 	        i = j-(((__HI(y[0]))>>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;
 		    i = j-(((__HI(y[0]))>>20)&0x7ff);
 		    if(i>49)  {	/* 3rd iteration need, 151 bits acc */
 		    	t  = r;	/* will cover all possible cases */
 		    	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;
 		    }
 		}
@@ -152,7 +135,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 */
@@ -24,56 +24,17 @@
 #ifdef __LITTLE_ENDIAN
 #define __HI(x) *(1+(int*)&x)
 #define __LO(x) *(int*)&x
 #define __HIp(x) *(1+(int*)x)
 #define __LOp(x) *(int*)x
 #else
 #define __HI(x) *(int*)&x
 #define __LO(x) *(1+(int*)&x)
 #define __HIp(x) *(int*)x
 #define __LOp(x) *(1+(int*)x)
 #endif
 #ifdef __STDC__
 #define	__P(p)	p
 #else
 #define	__P(p)	()
 #endif
 /*
 * ANSI/POSIX
 */
 extern int signgam;
 #define	MAXFLOAT	((float)3.40282346638528860e+38)
 enum fdversion {fdlibm_ieee = -1, fdlibm_svid, fdlibm_xopen, fdlibm_posix};
 #define _LIB_VERSION_TYPE enum fdversion
 #define _LIB_VERSION _fdlib_version
 /* if global variable _LIB_VERSION is not desirable, one may
 * change the following to be a constant by:
 *	#define _LIB_VERSION_TYPE const enum version
 * In that case, after one initializes the value _LIB_VERSION (see
 * s_lib_version.c) during compile time, it cannot be modified
 * in the middle of a program
 */
 extern  _LIB_VERSION_TYPE  _LIB_VERSION;
 #define _IEEE_  fdlibm_ieee
 #define _SVID_  fdlibm_svid
 #define _XOPEN_ fdlibm_xopen
 #define _POSIX_ fdlibm_posix
 struct exception {
 	int type;
 	char *name;
 	double arg1;
 	double arg2;
 	double retval;
 };
 #define	HUGE		MAXFLOAT
 /*
@@ -93,126 +54,39 @@ struct exception {
 /*
 * ANSI/POSIX
 */
-extern double acos __P((double));
+extern double acos (double);
-extern double asin __P((double));
+extern double asin (double);
-extern double atan __P((double));
+extern double atan (double);
-extern double atan2 __P((double, double));
+extern double atan2 (double, double);
-extern double cos __P((double));
+extern double cos (double);
-extern double sin __P((double));
+extern double sin (double);
-extern double tan __P((double));
+extern double tan (double);
-extern double cosh __P((double));
+extern double exp (double);
-extern double sinh __P((double));
+extern double log (double);
 extern double tanh __P((double));
-extern double exp __P((double));
+extern double pow (double, double);
-extern double frexp __P((double, int *));
+extern double sqrt (double);
 extern double ldexp __P((double, int));
 extern double log __P((double));
 extern double log10 __P((double));
 extern double modf __P((double, double *));
-extern double pow __P((double, double));
+extern double ceil (double);
-extern double sqrt __P((double));
+extern double fabs (double);
 extern double floor (double);
 extern double fmod (double, double);
-extern double ceil __P((double));
+extern int isnan (double);
-extern double fabs __P((double));
+extern int finite (double);
 extern double floor __P((double));
 extern double fmod __P((double, double));
 extern double erf __P((double));
 extern double erfc __P((double));
 extern double gamma __P((double));
 extern double hypot __P((double, double));
 extern int isnan __P((double));
 extern int finite __P((double));
 extern double j0 __P((double));
 extern double j1 __P((double));
 extern double jn __P((int, double));
 extern double lgamma __P((double));
 extern double y0 __P((double));
 extern double y1 __P((double));
 extern double yn __P((int, double));
 extern double acosh __P((double));
 extern double asinh __P((double));
 extern double atanh __P((double));
 extern double cbrt __P((double));
 extern double logb __P((double));
 extern double nextafter __P((double, double));
 extern double remainder __P((double, double));
 #ifdef _SCALB_INT
 extern double scalb __P((double, int));
 #else
 extern double scalb __P((double, double));
 #endif
 extern int matherr __P((struct exception *));
 /*
 * IEEE Test Vector
 */
 extern double significand __P((double));
 /*
 * Functions callable from C, intended to support IEEE arithmetic.
 */
-extern double copysign __P((double, double));
+extern double copysign (double, double);
-extern int ilogb __P((double));
+extern double scalbn (double, int);
 extern double rint __P((double));
 extern double scalbn __P((double, int));
 /*
 * BSD math library entry points
 */
 extern double expm1 __P((double));
 extern double log1p __P((double));
 /*
 * Reentrant version of gamma & lgamma; passes signgam back by reference
 * as the second argument; user must allocate space for signgam.
 */
 #ifdef _REENTRANT
 extern double gamma_r __P((double, int *));
 extern double lgamma_r __P((double, int *));
 #endif	/* _REENTRANT */
 /* ieee style elementary functions */
-extern double __ieee754_sqrt __P((double));
+extern int    __ieee754_rem_pio2 (double,double*);
 extern double __ieee754_acos __P((double));
 extern double __ieee754_acosh __P((double));
 extern double __ieee754_log __P((double));
 extern double __ieee754_atanh __P((double));
 extern double __ieee754_asin __P((double));
 extern double __ieee754_atan2 __P((double,double));
 extern double __ieee754_exp __P((double));
 extern double __ieee754_cosh __P((double));
 extern double __ieee754_fmod __P((double,double));
 extern double __ieee754_pow __P((double,double));
 extern double __ieee754_lgamma_r __P((double,int *));
 extern double __ieee754_gamma_r __P((double,int *));
 extern double __ieee754_lgamma __P((double));
 extern double __ieee754_gamma __P((double));
 extern double __ieee754_log10 __P((double));
 extern double __ieee754_sinh __P((double));
 extern double __ieee754_hypot __P((double,double));
 extern double __ieee754_j0 __P((double));
 extern double __ieee754_j1 __P((double));
 extern double __ieee754_y0 __P((double));
 extern double __ieee754_y1 __P((double));
 extern double __ieee754_jn __P((int,double));
 extern double __ieee754_yn __P((int,double));
 extern double __ieee754_remainder __P((double,double));
 extern int    __ieee754_rem_pio2 __P((double,double*));
 #ifdef _SCALB_INT
 extern double __ieee754_scalb __P((double,int));
 #else
 extern double __ieee754_scalb __P((double,double));
 #endif
 /* fdlibm kernel function */
-extern double __kernel_standard __P((double,double,int));
+extern double __kernel_sin (double,double,int);
-extern double __kernel_sin __P((double,double,int));
+extern double __kernel_cos (double,double);
-extern double __kernel_cos __P((double,double));
+extern double __kernel_tan (double,double,int);
-extern double __kernel_tan __P((double,double,int));
+extern int    __kernel_rem_pio2 (double*,double*,int,int,int,const int*);
 extern int    __kernel_rem_pio2 __P((double*,double*,int,int,int,const int*));
@@ -15,7 +15,7 @@
 * __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.
@@ -26,14 +26,14 @@
 *	   	cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
 *	   where the remez error is
 *	
- * 	|              2     4     6     8     10    12     14 |     -58
+ *	|              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
+ *	|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))
@@ -48,11 +48,7 @@
 #include "fdlibm.h"
-#ifdef __STDC__
+static const double
 static const double 
 #else
 static double 
 #endif
 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 C1  =  4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
 C2  = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
@@ -61,12 +57,7 @@ C4  = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
 C5  =  2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
 C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
-#ifdef __STDC__
+double __kernel_cos(double x, double y)
 	double __kernel_cos(double x, double y)
 #else
 	double __kernel_cos(x, y)
 	double x,y;
 #endif
 {
 	double a,hz,z,r,qx;
 	int ix;
@@ -76,7 +67,7 @@ C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
 	}
 	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 */
@@ -14,12 +14,12 @@
 /*
 * __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
@@ -28,10 +28,10 @@
 * (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 
+ *		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 
+ *		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]:
@@ -61,15 +61,15 @@
 *
 *	nx	dimension of x[]
 *
- *  	prec	an integer indicating the precision:
+ *	prec	an integer indicating the precision:
 *			0	24  bits (single)
 *			1	53  bits (double)
 *			2	64  bits (extended)
 *			3	113 bits (quad)
 *
 *	ipio2[]
- *		integer array, contains the (24*i)-th to (24*i+23)-th 
+ *		integer array, contains the (24*i)-th to (24*i+23)-th
- *		bit of 2/pi after binary point. The corresponding 
+ *		bit of 2/pi after binary point. The corresponding
 *		floating value is
 *
 *			ipio2[i] * 2^(-24(i+1)).
@@ -80,12 +80,12 @@
 *
 * Here is the description of some local variables:
 *
- * 	jk	jk+1 is the initial number of terms of ipio2[] needed
+ *	jk	jk+1 is the initial number of terms of ipio2[] needed
 *		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 
+ *	jz	local integer variable indicating the number of
- *		terms of ipio2[] used. 
+ *		terms of ipio2[] used.
 *
 *	jx	nx - 1
 *
@@ -98,16 +98,16 @@
 *
 *	jp	jp+1 is the number of terms in PIo2[] needed, jp = jk.
 *
- * 	q[]	double array with integral value, representing the
+ *	q[]	double array with integral value, representing the
 *		24-bits chunk of the product of x and 2/pi.
 *
 *	q0	the corresponding exponent of q[0]. Note that the
 *		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.
 *
@@ -121,25 +121,17 @@
 /*
 * Constants:
- * The hexadecimal values are the intended ones for the following 
+ * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the 
+ * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough 
+ * compiler will convert from decimal to binary accurately enough
 * to produce the hexadecimal values shown.
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
 #else
 static int init_jk[] = {2,3,4,6}; 
 #endif
 #ifdef __STDC__
 static const double PIo2[] = {
 #else
 static double PIo2[] = {
 #endif
  1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
  7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
  5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
@@ -150,22 +142,13 @@ static double PIo2[] = {
  2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
 };
-#ifdef __STDC__
+static const double
 static const double			
 #else
 static double			
 #endif
 zero   = 0.0,
 one    = 1.0,
 two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
 twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
-#ifdef __STDC__
+int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int *ipio2)
 	int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int *ipio2) 
 #else
 	int __kernel_rem_pio2(x,y,e0,nx,prec,ipio2) 	
 	double x[], y[]; int e0,nx,prec; int ipio2[];
 #endif
 {
 	int jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
 	double z,fw,f[20],fq[20],q[20];
@@ -258,7 +241,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)((int)(twon24*z));
 		iq[jz] = (int)(z-two24*fw);
 		jz += 1; q0 += 24;
@@ -283,29 +266,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]; 
+		for (i=jz;i>=0;i--) fw += fq[i];
-		y[0] = (ih==0)? fw: -fw; 
+		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 {
@@ -15,10 +15,10 @@
 * 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]
@@ -26,13 +26,13 @@
 *	   	sin(x) ~ x + S1*x + ... + S6*x
 *	   where
 *	
- * 	|sin(x)         2     4     6     8     10     12  |     -58
+ *	|sin(x)         2     4     6     8     10     12  |     -58
- * 	|----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
+ *	|----- - (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
@@ -41,11 +41,7 @@
 #include "fdlibm.h"
-#ifdef __STDC__
+static const double
 static const double 
 #else
 static double 
 #endif
 half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
 S1  = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
 S2  =  8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
@@ -54,12 +50,7 @@ S4  =  2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
 S5  = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
 S6  =  1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
-#ifdef __STDC__
+double __kernel_sin(double x, double y, int iy)
 	double __kernel_sin(double x, double y, int iy)
 #else
 	double __kernel_sin(x, y, iy)
 	double x,y; int iy;		/* iy=0 if y is zero */
 #endif
 {
 	double z,r,v;
 	int ix;
@@ -1,733 +0,0 @@
 /* @(#)k_standard.c 1.3 95/01/18 */
 /*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 *
 */
 #include "fdlibm.h"
 #include <errno.h>
 #ifndef _USE_WRITE
 #include <stdio.h>			/* fputs(), stderr */
 // #define	WRITE2(u,v)	fputs(u, stderr)
 #else	/* !defined(_USE_WRITE) */
 #include <unistd.h>			/* write */
 // #define	WRITE2(u,v)	write(2, u, v)
 #undef fflush
 #endif	/* !defined(_USE_WRITE) */
 static double zero = 0.0;	/* used as const */
 /*
 * Standard conformance (non-IEEE) on exception cases.
 * Mapping:
 *	1 -- acos(|x|>1)
 *	2 -- asin(|x|>1)
 *	3 -- atan2(+-0,+-0)
 *	4 -- hypot overflow
 *	5 -- cosh overflow
 *	6 -- exp overflow
 *	7 -- exp underflow
 *	8 -- y0(0)
 *	9 -- y0(-ve)
 *	10-- y1(0)
 *	11-- y1(-ve)
 *	12-- yn(0)
 *	13-- yn(-ve)
 *	14-- lgamma(finite) overflow
 *	15-- lgamma(-integer)
 *	16-- log(0)
 *	17-- log(x<0)
 *	18-- log10(0)
 *	19-- log10(x<0)
 *	20-- pow(0.0,0.0)
 *	21-- pow(x,y) overflow
 *	22-- pow(x,y) underflow
 *	23-- pow(0,negative)
 *	24-- pow(neg,non-integral)
 *	25-- sinh(finite) overflow
 *	26-- sqrt(negative)
 *      27-- fmod(x,0)
 *      28-- remainder(x,0)
 *	29-- acosh(x<1)
 *	30-- atanh(|x|>1)
 *	31-- atanh(|x|=1)
 *	32-- scalb overflow
 *	33-- scalb underflow
 *	34-- j0(|x|>X_TLOSS)
 *	35-- y0(x>X_TLOSS)
 *	36-- j1(|x|>X_TLOSS)
 *	37-- y1(x>X_TLOSS)
 *	38-- jn(|x|>X_TLOSS, n)
 *	39-- yn(x>X_TLOSS, n)
 *	40-- gamma(finite) overflow
 *	41-- gamma(-integer)
 *	42-- pow(NaN,0.0)
 */
 #ifdef __STDC__
 	double __kernel_standard(double x, double y, int type)
 #else
 	double __kernel_standard(x,y,type)
 	double x,y; int type;
 #endif
 {
 	struct exception exc;
 #ifndef HUGE_VAL	/* this is the only routine that uses HUGE_VAL */
 #define HUGE_VAL inf
 	double inf = 0.0;
 	__HI(inf) = 0x7ff00000;	/* set inf to infinite */
 #endif
 #ifdef _USE_WRITE
 	(void) fflush(stdout);
 #endif
 	exc.arg1 = x;
 	exc.arg2 = y;
 	switch(type) {
 	    case 1:
 		/* acos(|x|>1) */
 		exc.type = DOMAIN;
 		exc.name = "acos";
 		exc.retval = zero;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if(_LIB_VERSION == _SVID_) {
 		//     // (void) WRITE2("acos: DOMAIN error\n", 19);
 		//   }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 2:
 		/* asin(|x|>1) */
 		exc.type = DOMAIN;
 		exc.name = "asin";
 		exc.retval = zero;
 		// if(_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if(_LIB_VERSION == _SVID_) {
 		//     	// (void) WRITE2("asin: DOMAIN error\n", 19);
 		//   }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 3:
 		/* atan2(+-0,+-0) */
 		exc.arg1 = y;
 		exc.arg2 = x;
 		exc.type = DOMAIN;
 		exc.name = "atan2";
 		exc.retval = zero;
 		// if(_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if(_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("atan2: DOMAIN error\n", 20);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 4:
 		/* hypot(finite,finite) overflow */
 		exc.type = OVERFLOW;
 		exc.name = "hypot";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = HUGE;
 		else
 		  exc.retval = HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = ERANGE;
 		// else if (!matherr(&exc)) {
 		// 	errno = ERANGE;
 		// }
 		break;
 	    case 5:
 		/* cosh(finite) overflow */
 		exc.type = OVERFLOW;
 		exc.name = "cosh";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = HUGE;
 		else
 		  exc.retval = HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = ERANGE;
 		// else if (!matherr(&exc)) {
 		// 	errno = ERANGE;
 		// }
 		break;
 	    case 6:
 		/* exp(finite) overflow */
 		exc.type = OVERFLOW;
 		exc.name = "exp";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = HUGE;
 		else
 		  exc.retval = HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = ERANGE;
 		// else if (!matherr(&exc)) {
 		// 	errno = ERANGE;
 		// }
 		break;
 	    case 7:
 		/* exp(finite) underflow */
 		exc.type = UNDERFLOW;
 		exc.name = "exp";
 		exc.retval = zero;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = ERANGE;
 		// else if (!matherr(&exc)) {
 		// 	errno = ERANGE;
 		// }
 		break;
 	    case 8:
 		/* y0(0) = -inf */
 		exc.type = DOMAIN;	/* should be SING for IEEE */
 		exc.name = "y0";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = -HUGE;
 		else
 		  exc.retval = -HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		 //  if (_LIB_VERSION == _SVID_) {
 			// // (void) WRITE2("y0: DOMAIN error\n", 17);
 		 //      }
 		 //  errno = EDOM;
 		// }
 		break;
 	    case 9:
 		/* y0(x<0) = NaN */
 		exc.type = DOMAIN;
 		exc.name = "y0";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = -HUGE;
 		else
 		  exc.retval = -HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("y0: DOMAIN error\n", 17);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 10:
 		/* y1(0) = -inf */
 		exc.type = DOMAIN;	/* should be SING for IEEE */
 		exc.name = "y1";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = -HUGE;
 		else
 		  exc.retval = -HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("y1: DOMAIN error\n", 17);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 11:
 		/* y1(x<0) = NaN */
 		exc.type = DOMAIN;
 		exc.name = "y1";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = -HUGE;
 		else
 		  exc.retval = -HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("y1: DOMAIN error\n", 17);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 12:
 		/* yn(n,0) = -inf */
 		exc.type = DOMAIN;	/* should be SING for IEEE */
 		exc.name = "yn";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = -HUGE;
 		else
 		  exc.retval = -HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("yn: DOMAIN error\n", 17);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 13:
 		/* yn(x<0) = NaN */
 		exc.type = DOMAIN;
 		exc.name = "yn";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = -HUGE;
 		else
 		  exc.retval = -HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("yn: DOMAIN error\n", 17);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 14:
 		/* lgamma(finite) overflow */
 		exc.type = OVERFLOW;
 		exc.name = "lgamma";
                if (_LIB_VERSION == _SVID_)
                  exc.retval = HUGE;
                else
                  exc.retval = HUGE_VAL;
  //               if (_LIB_VERSION == _POSIX_)
 		// 	errno = ERANGE;
  //               else if (!matherr(&exc)) {
  //                       errno = ERANGE;
 		// }
 		break;
 	    case 15:
 		/* lgamma(-integer) or lgamma(0) */
 		exc.type = SING;
 		exc.name = "lgamma";
                if (_LIB_VERSION == _SVID_)
                  exc.retval = HUGE;
                else
                  exc.retval = HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("lgamma: SING error\n", 19);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 16:
 		/* log(0) */
 		exc.type = SING;
 		exc.name = "log";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = -HUGE;
 		else
 		  exc.retval = -HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = ERANGE;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("log: SING error\n", 16);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 17:
 		/* log(x<0) */
 		exc.type = DOMAIN;
 		exc.name = "log";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = -HUGE;
 		else
 		  exc.retval = -HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("log: DOMAIN error\n", 18);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 18:
 		/* log10(0) */
 		exc.type = SING;
 		exc.name = "log10";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = -HUGE;
 		else
 		  exc.retval = -HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = ERANGE;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("log10: SING error\n", 18);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 19:
 		/* log10(x<0) */
 		exc.type = DOMAIN;
 		exc.name = "log10";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = -HUGE;
 		else
 		  exc.retval = -HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("log10: DOMAIN error\n", 20);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 20:
 		/* pow(0.0,0.0) */
 		/* error only if _LIB_VERSION == _SVID_ */
 		exc.type = DOMAIN;
 		exc.name = "pow";
 		exc.retval = zero;
 		// if (_LIB_VERSION != _SVID_) exc.retval = 1.0;
 		// else if (!matherr(&exc)) {
 		// 	// (void) WRITE2("pow(0,0): DOMAIN error\n", 23);
 		// 	errno = EDOM;
 		// }
 		break;
 	    case 21:
 		/* pow(x,y) overflow */
 		exc.type = OVERFLOW;
 		exc.name = "pow";
 		if (_LIB_VERSION == _SVID_) {
 		  exc.retval = HUGE;
 		  y *= 0.5;
 		  if(x<zero&&rint(y)!=y) exc.retval = -HUGE;
 		} else {
 		  exc.retval = HUGE_VAL;
 		  y *= 0.5;
 		  if(x<zero&&rint(y)!=y) exc.retval = -HUGE_VAL;
 		}
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = ERANGE;
 		// else if (!matherr(&exc)) {
 		// 	errno = ERANGE;
 		// }
 		break;
 	    case 22:
 		/* pow(x,y) underflow */
 		exc.type = UNDERFLOW;
 		exc.name = "pow";
 		exc.retval =  zero;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = ERANGE;
 		// else if (!matherr(&exc)) {
 		// 	errno = ERANGE;
 		// }
 		break;
 	    case 23:
 		/* 0**neg */
 		exc.type = DOMAIN;
 		exc.name = "pow";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = zero;
 		else
 		  exc.retval = -HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("pow(0,neg): DOMAIN error\n", 25);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 24:
 		/* neg**non-integral */
 		exc.type = DOMAIN;
 		exc.name = "pow";
 		if (_LIB_VERSION == _SVID_)
 		    exc.retval = zero;
 		else
 		    exc.retval = zero/zero;	/* X/Open allow NaN */
 		// if (_LIB_VERSION == _POSIX_)
 		//    errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("neg**non-integral: DOMAIN error\n", 32);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 25:
 		/* sinh(finite) overflow */
 		exc.type = OVERFLOW;
 		exc.name = "sinh";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = ( (x>zero) ? HUGE : -HUGE);
 		else
 		  exc.retval = ( (x>zero) ? HUGE_VAL : -HUGE_VAL);
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = ERANGE;
 		// else if (!matherr(&exc)) {
 		// 	errno = ERANGE;
 		// }
 		break;
 	    case 26:
 		/* sqrt(x<0) */
 		exc.type = DOMAIN;
 		exc.name = "sqrt";
 		if (_LIB_VERSION == _SVID_)
 		  exc.retval = zero;
 		else
 		  exc.retval = zero/zero;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("sqrt: DOMAIN error\n", 19);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
            case 27:
                /* fmod(x,0) */
                exc.type = DOMAIN;
                exc.name = "fmod";
                if (_LIB_VERSION == _SVID_)
                    exc.retval = x;
 		else
 		    exc.retval = zero/zero;
                // if (_LIB_VERSION == _POSIX_)
                //   errno = EDOM;
                // else if (!matherr(&exc)) {
                //   if (_LIB_VERSION == _SVID_) {
                //     // (void) WRITE2("fmod:  DOMAIN error\n", 20);
                //   }
                //   errno = EDOM;
                // }
                break;
            case 28:
                /* remainder(x,0) */
                exc.type = DOMAIN;
                exc.name = "remainder";
                exc.retval = zero/zero;
                // if (_LIB_VERSION == _POSIX_)
                //   errno = EDOM;
                // else if (!matherr(&exc)) {
                //   if (_LIB_VERSION == _SVID_) {
                //     // (void) WRITE2("remainder: DOMAIN error\n", 24);
                //   }
                //   errno = EDOM;
                // }
                break;
            case 29:
                /* acosh(x<1) */
                exc.type = DOMAIN;
                exc.name = "acosh";
                exc.retval = zero/zero;
                // if (_LIB_VERSION == _POSIX_)
                //   errno = EDOM;
                // else if (!matherr(&exc)) {
                //   if (_LIB_VERSION == _SVID_) {
                //     // (void) WRITE2("acosh: DOMAIN error\n", 20);
                //   }
                //   errno = EDOM;
                // }
                break;
            case 30:
                /* atanh(|x|>1) */
                exc.type = DOMAIN;
                exc.name = "atanh";
                exc.retval = zero/zero;
                // if (_LIB_VERSION == _POSIX_)
                //   errno = EDOM;
                // else if (!matherr(&exc)) {
                //   if (_LIB_VERSION == _SVID_) {
                //     // (void) WRITE2("atanh: DOMAIN error\n", 20);
                //   }
                //   errno = EDOM;
                // }
                break;
            case 31:
                /* atanh(|x|=1) */
                exc.type = SING;
                exc.name = "atanh";
 		exc.retval = x/zero;	/* sign(x)*inf */
                // if (_LIB_VERSION == _POSIX_)
                //   errno = EDOM;
                // else if (!matherr(&exc)) {
                //   if (_LIB_VERSION == _SVID_) {
                //     // (void) WRITE2("atanh: SING error\n", 18);
                //   }
                //   errno = EDOM;
                // }
                break;
 	    case 32:
 		/* scalb overflow; SVID also returns +-HUGE_VAL */
 		exc.type = OVERFLOW;
 		exc.name = "scalb";
 		exc.retval = x > zero ? HUGE_VAL : -HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = ERANGE;
 		// else if (!matherr(&exc)) {
 		// 	errno = ERANGE;
 		// }
 		break;
 	    case 33:
 		/* scalb underflow */
 		exc.type = UNDERFLOW;
 		exc.name = "scalb";
 		exc.retval = copysign(zero,x);
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = ERANGE;
 		// else if (!matherr(&exc)) {
 		// 	errno = ERANGE;
 		// }
 		break;
 	    case 34:
 		/* j0(|x|>X_TLOSS) */
                exc.type = TLOSS;
                exc.name = "j0";
                exc.retval = zero;
                // if (_LIB_VERSION == _POSIX_)
                //         errno = ERANGE;
                // else if (!matherr(&exc)) {
                //         if (_LIB_VERSION == _SVID_) {
                //                 // (void) WRITE2(exc.name, 2);
                //                 // (void) WRITE2(": TLOSS error\n", 14);
                //         }
                //         errno = ERANGE;
                // }
 		break;
 	    case 35:
 		/* y0(x>X_TLOSS) */
                exc.type = TLOSS;
                exc.name = "y0";
                exc.retval = zero;
                // if (_LIB_VERSION == _POSIX_)
                //         errno = ERANGE;
                // else if (!matherr(&exc)) {
                //         if (_LIB_VERSION == _SVID_) {
                //                 // (void) WRITE2(exc.name, 2);
                //                 // (void) WRITE2(": TLOSS error\n", 14);
                //         }
                //         errno = ERANGE;
                // }
 		break;
 	    case 36:
 		/* j1(|x|>X_TLOSS) */
                exc.type = TLOSS;
                exc.name = "j1";
                exc.retval = zero;
                // if (_LIB_VERSION == _POSIX_)
                //         errno = ERANGE;
                // else if (!matherr(&exc)) {
                //         if (_LIB_VERSION == _SVID_) {
                //                 // (void) WRITE2(exc.name, 2);
                //                 // (void) WRITE2(": TLOSS error\n", 14);
                //         }
                //         errno = ERANGE;
                // }
 		break;
 	    case 37:
 		/* y1(x>X_TLOSS) */
                exc.type = TLOSS;
                exc.name = "y1";
                exc.retval = zero;
                // if (_LIB_VERSION == _POSIX_)
                //         errno = ERANGE;
                // else if (!matherr(&exc)) {
                //         if (_LIB_VERSION == _SVID_) {
                //                 // (void) WRITE2(exc.name, 2);
                //                 // (void) WRITE2(": TLOSS error\n", 14);
                //         }
                //         errno = ERANGE;
                // }
 		break;
 	    case 38:
 		/* jn(|x|>X_TLOSS) */
                exc.type = TLOSS;
                exc.name = "jn";
                exc.retval = zero;
                // if (_LIB_VERSION == _POSIX_)
                //         errno = ERANGE;
                // else if (!matherr(&exc)) {
                //         if (_LIB_VERSION == _SVID_) {
                //                 // (void) WRITE2(exc.name, 2);
                //                 // (void) WRITE2(": TLOSS error\n", 14);
                //         }
                //         errno = ERANGE;
                // }
 		break;
 	    case 39:
 		/* yn(x>X_TLOSS) */
                exc.type = TLOSS;
                exc.name = "yn";
                exc.retval = zero;
                // if (_LIB_VERSION == _POSIX_)
                //         errno = ERANGE;
                // else if (!matherr(&exc)) {
                //         if (_LIB_VERSION == _SVID_) {
                //                 // (void) WRITE2(exc.name, 2);
                //                 // (void) WRITE2(": TLOSS error\n", 14);
                //         }
                //         errno = ERANGE;
                // }
 		break;
 	    case 40:
 		/* gamma(finite) overflow */
 		exc.type = OVERFLOW;
 		exc.name = "gamma";
                if (_LIB_VERSION == _SVID_)
                  exc.retval = HUGE;
                else
                  exc.retval = HUGE_VAL;
    //             if (_LIB_VERSION == _POSIX_)
 		  // errno = ERANGE;
    //             else if (!matherr(&exc)) {
    //               errno = ERANGE;
    //             }
 		break;
 	    case 41:
 		/* gamma(-integer) or gamma(0) */
 		exc.type = SING;
 		exc.name = "gamma";
                if (_LIB_VERSION == _SVID_)
                  exc.retval = HUGE;
                else
                  exc.retval = HUGE_VAL;
 		// if (_LIB_VERSION == _POSIX_)
 		//   errno = EDOM;
 		// else if (!matherr(&exc)) {
 		//   if (_LIB_VERSION == _SVID_) {
 		// 	// (void) WRITE2("gamma: SING error\n", 18);
 		//       }
 		//   errno = EDOM;
 		// }
 		break;
 	    case 42:
 		/* pow(NaN,0.0) */
 		/* error only if _LIB_VERSION == _SVID_ & _XOPEN_ */
 		exc.type = DOMAIN;
 		exc.name = "pow";
 		exc.retval = x;
 		// if (_LIB_VERSION == _IEEE_ ||
 		//     _LIB_VERSION == _POSIX_) exc.retval = 1.0;
 		// else if (!matherr(&exc)) {
 		// 	errno = EDOM;
 		// }
 		// break;
 	}
 	return exc.retval;
 }
@@ -10,7 +10,6 @@
 * ====================================================
 */
 /* INDENT OFF */
 /* __kernel_tan( x, y, k )
 * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
 * Input x is assumed to be bounded by ~pi/4 in magnitude.
@@ -26,9 +25,9 @@
 *	   	tan(x) ~ x + T1*x + ... + T13*x
 *	   where
 *
- * 	        |tan(x)         2     4            26   |     -59.2
+ *	        |tan(x)         2     4            26   |     -59.2
- * 	        |----- - (1+T1*x +T2*x +.... +T13*x    )| <= 2
+ *	        |----- - (1+T1*x +T2*x +.... +T13*x    )| <= 2
- * 	        |  x 					|
+ *	        |  x 					|
 *
 *	   Note: tan(x+y) = tan(x) + tan'(x)*y
 *		          ~ tan(x) + (1+x*x)*y
@@ -68,10 +67,9 @@ static const double xxx[] = {
 #define	pio4	xxx[14]
 #define	pio4lo	xxx[15]
 #define	T	xxx
 /* INDENT ON */
-double
+double __kernel_tan(double x, double y, int iy)
-__kernel_tan(double x, double y, int iy) {
+{
 	double z, r, v, w, s;
 	int ix, hx;
@@ -11,15 +11,15 @@
 * ====================================================
 */
-/* __ieee754_acos(x)
+/* 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)))
+ *	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
@@ -37,13 +37,9 @@
 #include "fdlibm.h"
-#ifdef __STDC__
+static const double
-static const double 
+one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
-#else
+pi  =  3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
 static double 
 #endif
 one=  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 pi =  3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
 pio2_hi =  1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
 pio2_lo =  6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
 pS0 =  1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
@@ -57,12 +53,7 @@ qS2 =  2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
 qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
 qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
-#ifdef __STDC__
+double acos(double x)
 	double __ieee754_acos(double x)
 #else
 	double __ieee754_acos(x)
 	double x;
 #endif
 {
 	double z,p,q,r,w,s,c,df;
 	int hx,ix;
@@ -11,13 +11,13 @@
 * ====================================================
 */
-/* __ieee754_asin(x)
+/* 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)
 *
@@ -44,11 +44,7 @@
 #include "fdlibm.h"
 #ifdef __STDC__
 static const double 
 #else
 static double 
 #endif
 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 huge =  1.000e+300,
 pio2_hi =  1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
@@ -66,12 +62,7 @@ qS2 =  2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
 qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
 qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
-#ifdef __STDC__
+double asin(double x)
 	double __ieee754_asin(double x)
 #else
 	double __ieee754_asin(x)
 	double x;
 #endif
 {
 	double t = 0,w,p,q,c,r,s;
 	int hx,ix;
@@ -80,8 +71,8 @@ qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 	if(ix>= 0x3ff00000) {		/* |x|>= 1 */
 	    if(((ix-0x3ff00000)|__LO(x))==0)
 		    /* asin(1)=+-pi/2 with inexact */
-		return x*pio2_hi+x*pio2_lo;	
+		return x*pio2_hi+x*pio2_lo;
-	    return (x-x)/(x-x);		/* asin(|x|>1) is NaN */   
+	    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*/
@@ -109,6 +100,6 @@ qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 	    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;
 }
@@ -26,41 +26,29 @@
 *      [39/16,INF]   atan(x) = atan(INF) + atan( -1/t )
 *
 * Constants:
- * The hexadecimal values are the intended ones for the following 
+ * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the 
+ * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough 
+ * compiler will convert from decimal to binary accurately enough
 * to produce the hexadecimal values shown.
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 static const double atanhi[] = {
 #else
 static double atanhi[] = {
 #endif
  4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
  7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
  9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
  1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
 };
 #ifdef __STDC__
 static const double atanlo[] = {
 #else
 static double atanlo[] = {
 #endif
  2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
  3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
  1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
  6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
 };
 #ifdef __STDC__
 static const double aT[] = {
 #else
 static double aT[] = {
 #endif
  3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
 -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
  1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
@@ -74,20 +62,11 @@ static double aT[] = {
  1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
 };
-#ifdef __STDC__
+static const double
 	static const double 
 #else
 	static double 
 #endif
 one   = 1.0,
 huge   = 1.0e300;
-#ifdef __STDC__
+double atan(double x)
 	double atan(double x)
 #else
 	double atan(x)
 	double x;
 #endif
 {
 	double w,s1,s2,z;
 	int ix,hx,id;
@@ -109,9 +88,9 @@ huge   = 1.0e300;
 	x = fabs(x);
 	if (ix < 0x3ff30000) {		/* |x| < 1.1875 */
 	    if (ix < 0x3fe60000) {	/* 7/16 <=|x|<11/16 */
-		id = 0; x = (2.0*x-one)/(2.0+x); 
+		id = 0; x = (2.0*x-one)/(2.0+x);
 	    } else {			/* 11/16<=|x|< 19/16 */
-		id = 1; x  = (x-one)/(x+one); 
+		id = 1; x  = (x-one)/(x+one);
 	    }
 	} else {
 	    if (ix < 0x40038000) {	/* |x| < 2.4375 */
@@ -12,10 +12,10 @@
 *
 */
-/* __ieee754_atan2(y,x)
+/* 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,
 *
@@ -33,19 +33,15 @@
 *	ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
 *
 * Constants:
- * The hexadecimal values are the intended ones for the following 
+ * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the 
+ * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough 
+ * compiler will convert from decimal to binary accurately enough
 * to produce the hexadecimal values shown.
 */
 #include "fdlibm.h"
-#ifdef __STDC__
+static const double
 static const double 
 #else
 static double 
 #endif
 tiny  = 1.0e-300,
 zero  = 0.0,
 pi_o_4  = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
@@ -53,12 +49,7 @@ pi_o_2  = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
 pi      = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */
 pi_lo   = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
-#ifdef __STDC__
+double atan2(double y, double x)
 	double __ieee754_atan2(double y, double x)
 #else
 	double __ieee754_atan2(y,x)
 	double  y,x;
 #endif
 {  
 	double z;
 	int k,m,hx,hy,ix,iy;
@@ -22,18 +22,9 @@
 #include "fdlibm.h"
 #ifdef __STDC__
 static const double huge = 1.0e300;
 #else
 static double huge = 1.0e300;
 #endif
-#ifdef __STDC__
+double ceil(double x)
 	double ceil(double x)
 #else
 	double ceil(x)
 	double x;
 #endif
 {
 	int i0,i1,j0;
 	unsigned i,j;
@@ -43,7 +34,7 @@ static double huge = 1.0e300;
 	if(j0<20) {
 	    if(j0<0) { 	/* raise inexact if x != 0 */
 		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
-		    if(i0<0) {i0=0x80000000;i1=0;} 
+		    if(i0<0) {i0=0x80000000;i1=0;}
 		    else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;}
 		}
 	    } else {
@@ -62,7 +53,7 @@ static double huge = 1.0e300;
 	    if((i1&i)==0) return x;	/* x is integral */
 	    if(huge+x>0.0) { 		/* raise inexact flag */
 		if(i0>0) {
-		    if(j0==20) i0+=1; 
+		    if(j0==20) i0+=1;
 		    else {
 			j = i1 + (1<<(52-j0));
 			if(j<i1) i0+=1;	/* got a carry */
@@ -19,12 +19,7 @@
 #include "fdlibm.h"
-#ifdef __STDC__
+double copysign(double x, double y)
 	double copysign(double x, double y)
 #else
 	double copysign(x,y)
 	double x,y;
 #endif
 {
 	__HI(x) = (__HI(x)&0x7fffffff)|(__HI(y)&0x80000000);
        return x;
@@ -20,8 +20,8 @@
 *	__ieee754_rem_pio2	... argument reduction routine
 *
 * Method.
- *      Let S,C and T denote the sin, cos and tan respectively on 
+ *      Let S,C and T denote the sin, cos and tan respectively on
- *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 
+ *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
 *	in [-pi/4 , +pi/4], and let n = k mod 4.
 *	We have
 *
@@ -39,17 +39,12 @@
 *      trig(NaN)    is that NaN;
 *
 * Accuracy:
- *	TRIG(x) returns trig(x) nearly rounded 
+ *	TRIG(x) returns trig(x) nearly rounded
 */
 #include "fdlibm.h"
-#ifdef __STDC__
+double cos(double x)
 	double cos(double x)
 #else
 	double cos(x)
 	double x;
 #endif
 {
 	double y[2],z=0.0;
 	int n, ix;
@@ -10,7 +10,7 @@
 * ====================================================
 */
-/* __ieee754_exp(x)
+/* exp(x)
 * Returns the exponential of x.
 *
 * Method
@@ -18,36 +18,36 @@
 *      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 Remes algorithm on [0,0.34658] to generate 
+ *      We use a special Remes algorithm on [0,0.34658] to generate
- * 	a polynomial of degree 5 to approximate R. The maximum error 
+ *	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)
+ *	(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)
@@ -62,24 +62,20 @@
 *	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 
+ * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the 
+ * 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.
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 static const double
 #else
 static double
 #endif
 one	= 1.0,
 halF[2]	= {0.5,-0.5,},
 huge	= 1.0e+300,
@@ -98,12 +94,7 @@ P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
 P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
-#ifdef __STDC__
+double exp(double x)	/* default IEEE double exp */
 	double __ieee754_exp(double x)	/* default IEEE double exp */
 #else
 	double __ieee754_exp(x)	/* default IEEE double exp */
 	double x;
 #endif
 {
 	double y,hi,lo,c,t;
 	int k = 0,xsb;
@@ -125,7 +116,7 @@ P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
 	}
    /* 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 {
@@ -135,7 +126,7 @@ P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
 		lo = t*ln2LO[0];
 	    }
 	    x  = hi - lo;
-	} 
+	}
 	else if(hx < 0x3e300000)  {	/* when |x|<2**-28 */
 	    if(huge+x>one) return one+x;/* trigger inexact */
 	}
@@ -144,7 +135,7 @@ P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
    /* 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) {
 	    __HI(y) += (k<<20);	/* add k to y's exponent */
@@ -17,12 +17,7 @@
 #include "fdlibm.h"
-#ifdef __STDC__
+double fabs(double x)
 	double fabs(double x)
 #else
 	double fabs(x)
 	double x;
 #endif
 {
 	__HI(x) &= 0x7fffffff;
        return x;
@@ -18,14 +18,9 @@
 #include "fdlibm.h"
-#ifdef __STDC__
+int finite(double x)
 	int finite(double x)
 #else
 	int finite(x)
 	double x;
 #endif
 {
-	int hx; 
+	int hx;
 	hx = __HI(x);
 	return  (unsigned)((hx&0x7fffffff)-0x7ff00000)>>31;
 }
@@ -22,18 +22,9 @@
 #include "fdlibm.h"
 #ifdef __STDC__
 static const double huge = 1.0e300;
 #else
 static double huge = 1.0e300;
 #endif
-#ifdef __STDC__
+double floor(double x)
 	double floor(double x)
 #else
 	double floor(x)
 	double x;
 #endif
 {
 	int i0,i1,j0;
 	unsigned i,j;
@@ -43,7 +34,7 @@ static double huge = 1.0e300;
 	if(j0<20) {
 	    if(j0<0) { 	/* raise inexact if x != 0 */
 		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
-		    if(i0>=0) {i0=i1=0;} 
+		    if(i0>=0) {i0=i1=0;}
 		    else if(((i0&0x7fffffff)|i1)!=0)
 			{ i0=0xbff00000;i1=0;}
 		}
@@ -63,7 +54,7 @@ static double huge = 1.0e300;
 	    if((i1&i)==0) return x;	/* x is integral */
 	    if(huge+x>0.0) { 		/* raise inexact flag */
 		if(i0<0) {
-		    if(j0==20) i0+=1; 
+		    if(j0==20) i0+=1;
 		    else {
 			j = i1+(1<<(52-j0));
 			if(j<i1) i0 +=1 ; 	/* got a carry */
@@ -12,25 +12,18 @@
 */
 /*
- * __ieee754_fmod(x,y)
+ * fmod(x,y)
 * Return x mod y in exact arithmetic
 * Method: shift and subtract
 */
 #include "fdlibm.h"
-#ifdef __STDC__
+static const double
-static const double one = 1.0, Zero[] = {0.0, -0.0,};
+one    = 1.0,
-#else
+Zero[] = {0.0, -0.0,};
 static double one = 1.0, Zero[] = {0.0, -0.0,};
 #endif
-#ifdef __STDC__
+double fmod(double x, double y)
 	double __ieee754_fmod(double x, double y)
 #else
 	double __ieee754_fmod(x,y)
 	double x,y ;
 #endif
 {
 	int n,hx,hy,hz,ix,iy,sx,i;
 	unsigned lx,ly,lz;
@@ -18,17 +18,12 @@
 #include "fdlibm.h"
-#ifdef __STDC__
+int isnan(double x)
 	int isnan(double x)
 #else
 	int isnan(x)
 	double x;
 #endif
 {
 	int hx,lx;
 	hx = (__HI(x)&0x7fffffff);
 	lx = __LO(x);
-	hx |= (unsigned)(lx|(-lx))>>31;	
+	hx |= (unsigned)(lx|(-lx))>>31;
 	hx = 0x7ff00000 - hx;
 	return ((unsigned)(hx))>>31;
 }
@@ -1,35 +0,0 @@
 /* @(#)s_lib_version.c 1.3 95/01/18 */
 /*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */
 /*
 * MACRO for standards
 */
 #include "fdlibm.h"
 /*
 * define and initialize _LIB_VERSION
 */
 #ifdef _POSIX_MODE
 _LIB_VERSION_TYPE _LIB_VERSION = _POSIX_;
 #else
 #ifdef _XOPEN_MODE
 _LIB_VERSION_TYPE _LIB_VERSION = _XOPEN_;
 #else
 #ifdef _SVID3_MODE
 _LIB_VERSION_TYPE _LIB_VERSION = _SVID_;
 #else					/* default _IEEE_MODE */
 _LIB_VERSION_TYPE _LIB_VERSION = _IEEE_;
 #endif
 #endif
 #endif
@@ -11,28 +11,28 @@
 * ====================================================
 */
-/* __ieee754_log(x)
+/* log(x)
 * Return the logrithm of x
 *
- * Method :                  
+ * Method :
- *   1. Argument Reduction: find k and f such that 
+ *   1. Argument Reduction: find k and f such that
- *			x = 2^k * (1+f), 
+ *			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 
+ *      We use a special Reme algorithm on [0,0.1716] to generate
- * 	a polynomial of degree 14 to approximate R The maximum error 
+ *	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
 *	    R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
- *  	(the values of Lg1 to Lg7 are listed in the program)
+ *	(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
@@ -40,14 +40,14 @@
 *		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.
 *
@@ -56,19 +56,15 @@
 *	1 ulp (unit in the last place).
 *
 * Constants:
- * The hexadecimal values are the intended ones for the following 
+ * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the 
+ * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough 
+ * compiler will convert from decimal to binary accurately enough
 * to produce the hexadecimal values shown.
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 static const double
 #else
 static double
 #endif
 ln2_hi  =  6.93147180369123816490e-01,	/* 3fe62e42 fee00000 */
 ln2_lo  =  1.90821492927058770002e-10,	/* 3dea39ef 35793c76 */
 two54   =  1.80143985094819840000e+16,  /* 43500000 00000000 */
@@ -80,14 +76,9 @@ Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
 Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
 Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
-static double zero   =  0.0;
+static const double zero   =  0.0;
-#ifdef __STDC__
+double log(double x)
 	double __ieee754_log(double x)
 #else
 	double __ieee754_log(x)
 	double x;
 #endif
 {
 	double hfsq,f,s,z,R,w,t1,t2,dk;
 	int k,hx,i,j;
@@ -118,14 +109,14 @@ static double zero   =  0.0;
 	    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)); 
+	t1= w*(Lg2+w*(Lg4+w*Lg6));
-	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 
+	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
 	i |= j;
 	R = t2+t1;
 	if(i>0) {
@@ -10,14 +10,14 @@
 * ====================================================
 */
-/* __ieee754_pow(x,y) return x**y
+/* pow(x,y) return x**y
 *
 *		      n
 * Method:  Let x =  2   * (1+f)
 *	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)
 *
@@ -45,23 +45,19 @@
 * 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 
+ * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the 
+ * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough 
+ * compiler will convert from decimal to binary accurately enough
 * to produce the hexadecimal values shown.
 */
 #include "fdlibm.h"
-#ifdef __STDC__
+static const double
 static const double 
 #else
 static double 
 #endif
 bp[] = {1.0, 1.5,},
 dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
 dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
@@ -94,12 +90,7 @@ ivln2    =  1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
 ivln2_h  =  1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
 ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
-#ifdef __STDC__
+double pow(double x, double y)
 	double __ieee754_pow(double x, double y)
 #else
 	double __ieee754_pow(x,y)
 	double x, y;
 #endif
 {
 	double z,ax,z_h,z_l,p_h,p_l;
 	double y1,t1,t2,r,s,t,u,v,w;
@@ -113,12 +104,12 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 	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))) 
+	   iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
-		return x+y;	
+		return x+y;
    /* determine if y is an odd int when x < 0
     * yisint = 0	... y is not an integer
@@ -126,7 +117,7 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
     * 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 */
@@ -137,11 +128,11 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 		    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 */
@@ -156,7 +147,7 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 	    if(hy==0x40000000) return x*x; /* y is  2 */
 	    if(hy==0x3fe00000) {	/* y is  0.5 */
 		if(hx>=0)	/* x >= +0 */
-		return sqrt(x);	
+		return sqrt(x);
 	    }
 	}
@@ -169,13 +160,13 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 		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;
 	    }
 	}
-    
+
 	n = (hx>>31)+1;
    /* (x<0)**(non-int) is NaN */
@@ -193,7 +184,7 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 	/* over/underflow if x is not close to one */
 	    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 
+	/* 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-one;		/* t has 20 trailing zeros */
 	    w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
@@ -225,7 +216,7 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 	    __LO(s_h) = 0;
 	/* t_h=ax+bp[k] High */
 	    t_h = zero;
-	    __HI(t_h)=((ix>>1)|0x20000000)+0x00080000+(k<<18); 
+	    __HI(t_h)=((ix>>1)|0x20000000)+0x00080000+(k<<18);
 	    t_l = ax - (t_h-bp[k]);
 	    s_l = v*((u-s_h*t_h)-s_h*t_l);
 	/* compute log(ax) */
@@ -287,7 +278,7 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 	    n = ((n&0x000fffff)|0x00100000)>>(20-k);
 	    if(j<0) n = -n;
 	    p_h -= t;
-	} 
+	}
 	t = p_l+p_h;
 	__LO(t) = 0;
 	u = t*lg2_h;
@@ -1,84 +0,0 @@
 /* @(#)s_rint.c 1.3 95/01/18 */
 /*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */
 /*
 * rint(x)
 * Return x rounded to integral value according to the prevailing
 * rounding mode.
 * Method:
 *	Using floating addition.
 * Exception:
 *	Inexact flag raised if x not equal to rint(x).
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 static const double
 #else
 static double 
 #endif
 TWO52[2]={
  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
 -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
 };
 #ifdef __STDC__
 	double rint(double x)
 #else
 	double rint(x)
 	double x;
 #endif
 {
 	int i0,j0,sx;
 	unsigned i,i1;
 	double w,t;
 	i0 =  __HI(x);
 	sx = (i0>>31)&1;
 	i1 =  __LO(x);
 	j0 = ((i0>>20)&0x7ff)-0x3ff;
 	if(j0<20) {
 	    if(j0<0) { 	
 		if(((i0&0x7fffffff)|i1)==0) return x;
 		i1 |= (i0&0x0fffff);
 		i0 &= 0xfffe0000;
 		i0 |= ((i1|-i1)>>12)&0x80000;
 		__HI(x)=i0;
 	        w = TWO52[sx]+x;
 	        t =  w-TWO52[sx];
 	        i0 = __HI(t);
 	        __HI(t) = (i0&0x7fffffff)|(sx<<31);
 	        return t;
 	    } else {
 		i = (0x000fffff)>>j0;
 		if(((i0&i)|i1)==0) return x; /* x is integral */
 		i>>=1;
 		if(((i0&i)|i1)!=0) {
 		    if(j0==19) i1 = 0x40000000; else
 		    i0 = (i0&(~i))|((0x20000)>>j0);
 		}
 	    }
 	} else if (j0>51) {
 	    if(j0==0x400) return x+x;	/* inf or NaN */
 	    else return x;		/* x is integral */
 	} else {
 	    i = ((unsigned)(0xffffffff))>>(j0-20);
 	    if((i1&i)==0) return x;	/* x is integral */
 	    i>>=1;
 	    if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
 	}
 	__HI(x) = i0;
 	__LO(x) = i1;
 	w = TWO52[sx]+x;
 	return w-TWO52[sx];
 }
@@ -11,31 +11,22 @@
 * ====================================================
 */
-/* 
+/*
 * scalbn (double x, int n)
- * scalbn(x,n) returns x* 2**n  computed by  exponent  
+ * scalbn(x,n) returns x* 2**n  computed by  exponent
- * manipulation rather than by actually performing an 
+ * manipulation rather than by actually performing an
 * exponentiation or a multiplication.
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 static const double
-#else
+two54  = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
-static double
+twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
 #endif
 two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
 twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
 huge   = 1.0e+300,
 tiny   = 1.0e-300;
-#ifdef __STDC__
+double scalbn (double x, int n)
 	double scalbn (double x, int n)
 #else
 	double scalbn (x,n)
 	double x; int n;
 #endif
 {
 	int  k,hx,lx;
 	hx = __HI(x);
@@ -45,7 +36,7 @@ tiny   = 1.0e-300;
            if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
 	    x *= two54; 
 	    hx = __HI(x);
-	    k = ((hx&0x7ff00000)>>20) - 54; 
+	    k = ((hx&0x7ff00000)>>20) - 54;
            if (n< -50000) return tiny*x; 	/*underflow*/
 	    }
        if (k==0x7ff) return x+x;		/* NaN or Inf */
@@ -1,30 +0,0 @@
 /* @(#)s_significand.c 1.3 95/01/18 */
 /*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */
 /*
 * significand(x) computes just
 * 	scalb(x, (double) -ilogb(x)),
 * for exercising the fraction-part(F) IEEE 754-1985 test vector.
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 	double significand(double x)
 #else
 	double significand(x)
 	double x;
 #endif
 {
 	return __ieee754_scalb(x,(double) -ilogb(x));
 }
@@ -20,8 +20,8 @@
 *	__ieee754_rem_pio2	... argument reduction routine
 *
 * Method.
- *      Let S,C and T denote the sin, cos and tan respectively on 
+ *      Let S,C and T denote the sin, cos and tan respectively on
- *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 
+ *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
 *	in [-pi/4 , +pi/4], and let n = k mod 4.
 *	We have
 *
@@ -39,17 +39,12 @@
 *      trig(NaN)    is that NaN;
 *
 * Accuracy:
- *	TRIG(x) returns trig(x) nearly rounded 
+ *	TRIG(x) returns trig(x) nearly rounded
 */
 #include "fdlibm.h"
-#ifdef __STDC__
+double sin(double x)
 	double sin(double x)
 #else
 	double sin(x)
 	double x;
 #endif
 {
 	double y[2],z=0.0;
 	int n, ix;
@@ -11,15 +11,15 @@
 * ====================================================
 */
-/* __ieee754_sqrt(x)
+/* sqrt(x)
 * Return correctly rounded sqrt.
 *           ------------------------------------------
 *	     |  Use the hardware sqrt if you have one |
 *           ------------------------------------------
- * Method: 
+ * Method:
- *   Bit by bit method using integer arithmetic. (Slow, but portable) 
+ *   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
@@ -28,9 +28,9 @@
 *                                     i+1         2
 *	    s  = 2*q , and	y  =  2   * ( y - q  ).		(1)
 *	     i      i            i                 i
- *                                                        
+ *
- *	To compute q    from q , one checks whether 
+ *	To compute q    from q , one checks whether
- *		    i+1       i                       
+ *		    i+1       i
 *
 *			      -(i+1) 2
 *			(q + 2      ) <= y.			(2)
@@ -45,7 +45,7 @@
 *			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
@@ -57,10 +57,10 @@
 *                         -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.
@@ -70,7 +70,7 @@
 *	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
@@ -83,21 +83,14 @@
 #include "fdlibm.h"
-#ifdef __STDC__
+static const double
-static	const double	one	= 1.0, tiny=1.0e-300;
+one  = 1.0,
-#else
+tiny = 1.0e-300;
 static	double	one	= 1.0, tiny=1.0e-300;
 #endif
-#ifdef __STDC__
+double sqrt(double x)
 	double __ieee754_sqrt(double x)
 #else
 	double __ieee754_sqrt(x)
 	double x;
 #endif
 {
 	double z;
-	int 	sign = (int)0x80000000; 
+	int 	sign = (int)0x80000000;
 	unsigned r,t1,s1,ix1,q1;
 	int ix0,s0,q,m,t,i;
@@ -108,7 +101,7 @@ static	double	one	= 1.0, tiny=1.0e-300;
 	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 */
@@ -142,12 +135,12 @@ static	double	one	= 1.0, tiny=1.0e-300;
 	r = 0x00200000;		/* r = moving bit from right to left */
 	while(r!=0) {
-	    t = s0+r; 
+	    t = s0+r;
-	    if(t<=ix0) { 
+	    if(t<=ix0) {
-		s0   = t+r; 
+		s0   = t+r;
-		ix0 -= t; 
+		ix0 -= t;
-		q   += r; 
+		q   += r;
-	    } 
+	    }
 	    ix0 += ix0 + ((ix1&sign)>>31);
 	    ix1 += ix1;
 	    r>>=1;
@@ -155,9 +148,9 @@ static	double	one	= 1.0, tiny=1.0e-300;
 	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;
@@ -178,7 +171,7 @@ static	double	one	= 1.0, tiny=1.0e-300;
 	        if (q1==(unsigned)0xffffffff) { q1=0; q += 1;}
 		else if (z>one) {
 		    if (q1==(unsigned)0xfffffffe) q+=1;
-		    q1+=2; 
+		    q1+=2;
 		} else
 	            q1 += (q1&1);
 	    }
@@ -195,18 +188,18 @@ static	double	one	= 1.0, tiny=1.0e-300;
 /*
 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.
@@ -216,7 +209,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
 	   ------------------------------------------------------
@@ -224,7 +217,7 @@ A.  sqrt(x) by Newton Iteration
 	   ------------------------------------------------------
 	      msb    lsb  msb				      lsb ...order
- 
+
 	     ------------------------  	     ------------------------
 	x0:  |s|   e    |    f1     |	 x1: |          f2           |
 	     ------------------------  	     ------------------------
@@ -249,7 +242,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
@@ -274,12 +267,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
@@ -310,7 +303,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;
@@ -329,7 +322,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)
@@ -350,9 +343,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)
@@ -361,7 +354,7 @@ 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.
@@ -408,27 +401,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
@@ -445,7 +438,6 @@ 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).
- 
+
 */
@@ -19,8 +19,8 @@
 *	__ieee754_rem_pio2	... argument reduction routine
 *
 * Method.
- *      Let S,C and T denote the sin, cos and tan respectively on 
+ *      Let S,C and T denote the sin, cos and tan respectively on
- *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 
+ *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
 *	in [-pi/4 , +pi/4], and let n = k mod 4.
 *	We have
 *
@@ -38,17 +38,12 @@
 *      trig(NaN)    is that NaN;
 *
 * Accuracy:
- *	TRIG(x) returns trig(x) nearly rounded 
+ *	TRIG(x) returns trig(x) nearly rounded
 */
 #include "fdlibm.h"
-#ifdef __STDC__
+double tan(double x)
 	double tan(double x)
 #else
 	double tan(x)
 	double x;
 #endif
 {
 	double y[2],z=0.0;
 	int n, ix;
@@ -1,82 +0,0 @@
 /* @(#)s_tanh.c 1.3 95/01/18 */
 /*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */
 /* Tanh(x)
 * Return the Hyperbolic Tangent of x
 *
 * Method :
 *				       x    -x
 *				      e  - e
 *	0. tanh(x) is defined to be -----------
 *				       x    -x
 *				      e  + e
 *	1. reduce x to non-negative by tanh(-x) = -tanh(x).
 *	2.  0      <= x <= 2**-55 : tanh(x) := x*(one+x)
 *					        -t
 *	    2**-55 <  x <=  1     : tanh(x) := -----; t = expm1(-2x)
 *					       t + 2
 *						     2
 *	    1      <= x <=  22.0  : tanh(x) := 1-  ----- ; t=expm1(2x)
 *						   t + 2
 *	    22.0   <  x <= INF    : tanh(x) := 1.
 *
 * Special cases:
 *	tanh(NaN) is NaN;
 *	only tanh(0)=0 is exact for finite argument.
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 static const double one=1.0, two=2.0, tiny = 1.0e-300;
 #else
 static double one=1.0, two=2.0, tiny = 1.0e-300;
 #endif
 #ifdef __STDC__
 	double tanh(double x)
 #else
 	double tanh(x)
 	double x;
 #endif
 {
 	double t,z;
 	int jx,ix;
    /* High word of |x|. */
 	jx = __HI(x);
 	ix = jx&0x7fffffff;
    /* x is INF or NaN */
 	if(ix>=0x7ff00000) { 
 	    if (jx>=0) return one/x+one;    /* tanh(+-inf)=+-1 */
 	    else       return one/x-one;    /* tanh(NaN) = NaN */
 	}
    /* |x| < 22 */
 	if (ix < 0x40360000) {		/* |x|<22 */
 	    if (ix<0x3c800000) 		/* |x|<2**-55 */
 		return x*(one+x);    	/* tanh(small) = small */
 	    if (ix>=0x3ff00000) {	/* |x|>=1  */
 		t = expm1(two*fabs(x));
 		z = one - two/(t+two);
 	    } else {
 	        t = expm1(-two*fabs(x));
 	        z= -t/(t+two);
 	    }
    /* |x| > 22, return +-1 */
 	} else {
 	    z = one - tiny;		/* raised inexact flag */
 	}
 	return (jx>=0)? z: -z;
 }
@@ -1,39 +0,0 @@
 /* @(#)w_acos.c 1.3 95/01/18 */
 /*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */
 /*
 * wrap_acos(x)
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 	double acos(double x)		/* wrapper acos */
 #else
 	double acos(x)			/* wrapper acos */
 	double x;
 #endif
 {
 #ifdef _IEEE_LIBM
 	return __ieee754_acos(x);
 #else
 	double z;
 	z = __ieee754_acos(x);
 	if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
 	if(fabs(x)>1.0) {
 	        return __kernel_standard(x,x,1); /* acos(|x|>1) */
 	} else
 	    return z;
 #endif
 }
@@ -1,41 +0,0 @@
 /* @(#)w_asin.c 1.3 95/01/18 */
 /*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 *
 */
 /* 
 * wrapper asin(x)
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 	double asin(double x)		/* wrapper asin */
 #else
 	double asin(x)			/* wrapper asin */
 	double x;
 #endif
 {
 #ifdef _IEEE_LIBM
 	return __ieee754_asin(x);
 #else
 	double z;
 	z = __ieee754_asin(x);
 	if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
 	if(fabs(x)>1.0) {
 	        return __kernel_standard(x,x,2); /* asin(|x|>1) */
 	} else
 	    return z;
 #endif
 }
@@ -1,40 +0,0 @@
 /* @(#)w_atan2.c 1.3 95/01/18 */
 /*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 *
 */
 /* 
 * wrapper atan2(y,x)
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 	double atan2(double y, double x)	/* wrapper atan2 */
 #else
 	double atan2(y,x)			/* wrapper atan2 */
 	double y,x;
 #endif
 {
 #ifdef _IEEE_LIBM
 	return __ieee754_atan2(y,x);
 #else
 	double z;
 	z = __ieee754_atan2(y,x);
 	if(_LIB_VERSION == _IEEE_||isnan(x)||isnan(y)) return z;
 	if(x==0.0&&y==0.0) {
 	        return __kernel_standard(y,x,3); /* atan2(+-0,+-0) */
 	} else
 	    return z;
 #endif
 }
@@ -1,48 +0,0 @@
 /* @(#)w_exp.c 1.4 04/04/22 */
 /*
 * ====================================================
 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
 *
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */
 /* 
 * wrapper exp(x)
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 static const double
 #else
 static double
 #endif
 o_threshold=  7.09782712893383973096e+02,  /* 0x40862E42, 0xFEFA39EF */
 u_threshold= -7.45133219101941108420e+02;  /* 0xc0874910, 0xD52D3051 */
 #ifdef __STDC__
 	double exp(double x)		/* wrapper exp */
 #else
 	double exp(x)			/* wrapper exp */
 	double x;
 #endif
 {
 #ifdef _IEEE_LIBM
 	return __ieee754_exp(x);
 #else
 	double z;
 	z = __ieee754_exp(x);
 	if(_LIB_VERSION == _IEEE_) return z;
 	if(finite(x)) {
 	    if(x>o_threshold)
 	        return __kernel_standard(x,x,6); /* exp overflow */
 	    else if(x<u_threshold)
 	        return __kernel_standard(x,x,7); /* exp underflow */
 	} 
 	return z;
 #endif
 }
@@ -1,39 +0,0 @@
 /* @(#)w_fmod.c 1.3 95/01/18 */
 /*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */
 /*
 * wrapper fmod(x,y)
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 	double fmod(double x, double y)	/* wrapper fmod */
 #else
 	double fmod(x,y)		/* wrapper fmod */
 	double x,y;
 #endif
 {
 #ifdef _IEEE_LIBM
 	return __ieee754_fmod(x,y);
 #else
 	double z;
 	z = __ieee754_fmod(x,y);
 	if(_LIB_VERSION == _IEEE_ ||isnan(y)||isnan(x)) return z;
 	if(y==0.0) {
 	        return __kernel_standard(x,y,27); /* fmod(x,0) */
 	} else
 	    return z;
 #endif
 }
@@ -1,39 +0,0 @@
 /* @(#)w_log.c 1.3 95/01/18 */
 /*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */
 /*
 * wrapper log(x)
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 	double log(double x)		/* wrapper log */
 #else
 	double log(x)			/* wrapper log */
 	double x;
 #endif
 {
 #ifdef _IEEE_LIBM
 	return __ieee754_log(x);
 #else
 	double z;
 	z = __ieee754_log(x);
 	if(_LIB_VERSION == _IEEE_ || isnan(x) || x > 0.0) return z;
 	if(x==0.0)
 	    return __kernel_standard(x,x,16); /* log(0) */
 	else 
 	    return __kernel_standard(x,x,17); /* log(x<0) */
 #endif
 }
@@ -1,59 +0,0 @@
 /* @(#)w_pow.c 1.3 95/01/18 */
 /*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */
 /* 
 * wrapper pow(x,y) return x**y
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 	double pow(double x, double y)	/* wrapper pow */
 #else
 	double pow(x,y)			/* wrapper pow */
 	double x,y;
 #endif
 {
 #ifdef _IEEE_LIBM
 	return  __ieee754_pow(x,y);
 #else
 	double z;
 	z=__ieee754_pow(x,y);
 	if(_LIB_VERSION == _IEEE_|| isnan(y)) return z;
 	if(isnan(x)) {
 	    if(y==0.0) 
 	        return __kernel_standard(x,y,42); /* pow(NaN,0.0) */
 	    else 
 		return z;
 	}
 	if(x==0.0){ 
 	    if(y==0.0)
 	        return __kernel_standard(x,y,20); /* pow(0.0,0.0) */
 	    if(finite(y)&&y<0.0)
 	        return __kernel_standard(x,y,23); /* pow(0.0,negative) */
 	    return z;
 	}
 	if(!finite(z)) {
 	    if(finite(x)&&finite(y)) {
 	        if(isnan(z))
 	            return __kernel_standard(x,y,24); /* pow neg**non-int */
 	        else 
 	            return __kernel_standard(x,y,21); /* pow overflow */
 	    }
 	} 
 	if(z==0.0&&finite(x)&&finite(y))
 	    return __kernel_standard(x,y,22); /* pow underflow */
 	return z;
 #endif
 }
@@ -1,38 +0,0 @@
 /* @(#)w_sqrt.c 1.3 95/01/18 */
 /*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */
 /* 
 * wrapper sqrt(x)
 */
 #include "fdlibm.h"
 #ifdef __STDC__
 	double sqrt(double x)		/* wrapper sqrt */
 #else
 	double sqrt(x)			/* wrapper sqrt */
 	double x;
 #endif
 {
 #ifdef _IEEE_LIBM
 	return __ieee754_sqrt(x);
 #else
 	double z;
 	z = __ieee754_sqrt(x);
 	if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
 	if(x<0.0) {
 	    return __kernel_standard(x,x,26); /* sqrt(negative) */
 	} else
 	    return z;
 #endif
 }