Updated math code from the uClibc 0.9.33 release
authorSam Lantinga <slouken@libsdl.org>
Sat, 04 Nov 2017 15:53:19 -0700
changeset 1168348bcba563d9c
parent 11682 b26412a89fbb
child 11684 eccdf37f8996
Updated math code from the uClibc 0.9.33 release
src/libm/e_atan2.c
src/libm/e_log.c
src/libm/e_pow.c
src/libm/e_rem_pio2.c
src/libm/e_sqrt.c
src/libm/k_cos.c
src/libm/k_rem_pio2.c
src/libm/k_sin.c
src/libm/k_tan.c
src/libm/math_libm.h
src/libm/math_private.h
src/libm/s_atan.c
src/libm/s_copysign.c
src/libm/s_cos.c
src/libm/s_fabs.c
src/libm/s_floor.c
src/libm/s_scalbn.c
src/libm/s_sin.c
     1.1 --- a/src/libm/e_atan2.c	Sat Nov 04 15:34:14 2017 -0700
     1.2 +++ b/src/libm/e_atan2.c	Sat Nov 04 15:53:19 2017 -0700
     1.3 @@ -57,8 +57,8 @@
     1.4  	ix = hx&0x7fffffff;
     1.5  	EXTRACT_WORDS(hy,ly,y);
     1.6  	iy = hy&0x7fffffff;
     1.7 -	if(((ix|((lx|-(int32_t)lx)>>31))>0x7ff00000)||
     1.8 -	   ((iy|((ly|-(int32_t)ly)>>31))>0x7ff00000))	/* x or y is NaN */
     1.9 +	if(((ix|((lx|-lx)>>31))>0x7ff00000)||
    1.10 +	   ((iy|((ly|-ly)>>31))>0x7ff00000))	/* x or y is NaN */
    1.11  	   return x+y;
    1.12  	if(((hx-0x3ff00000)|lx)==0) return atan(y);   /* x=1.0 */
    1.13  	m = ((hy>>31)&1)|((hx>>30)&2);	/* 2*sign(x)+sign(y) */
    1.14 @@ -81,8 +81,8 @@
    1.15  		switch(m) {
    1.16  		    case 0: return  pi_o_4+tiny;/* atan(+INF,+INF) */
    1.17  		    case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
    1.18 -		    case 2: return  3.0*pi_o_4+tiny;/* atan(+INF,-INF) */
    1.19 -		    case 3: return -3.0*pi_o_4-tiny;/* atan(-INF,-INF) */
    1.20 +		    case 2: return  3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
    1.21 +		    case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
    1.22  		}
    1.23  	    } else {
    1.24  		switch(m) {
    1.25 @@ -114,3 +114,21 @@
    1.26  	    	    return  (z-pi_lo)-pi;/* atan(-,-) */
    1.27  	}
    1.28  }
    1.29 +
    1.30 +/*
    1.31 + * wrapper atan2(y,x)
    1.32 + */
    1.33 +#ifndef _IEEE_LIBM
    1.34 +double atan2(double y, double x)
    1.35 +{
    1.36 +	double z = __ieee754_atan2(y, x);
    1.37 +	if (_LIB_VERSION == _IEEE_ || isnan(x) || isnan(y))
    1.38 +		return z;
    1.39 +	if (x == 0.0 && y == 0.0)
    1.40 +		return __kernel_standard(y,x,3); /* atan2(+-0,+-0) */
    1.41 +	return z;
    1.42 +}
    1.43 +#else
    1.44 +strong_alias(__ieee754_atan2, atan2)
    1.45 +#endif
    1.46 +libm_hidden_def(atan2)
     2.1 --- a/src/libm/e_log.c	Sat Nov 04 15:34:14 2017 -0700
     2.2 +++ b/src/libm/e_log.c	Sat Nov 04 15:53:19 2017 -0700
     2.3 @@ -1,4 +1,3 @@
     2.4 -/* @(#)e_log.c 5.1 93/09/24 */
     2.5  /*
     2.6   * ====================================================
     2.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     2.8 @@ -10,11 +9,6 @@
     2.9   * ====================================================
    2.10   */
    2.11  
    2.12 -#if defined(LIBM_SCCS) && !defined(lint)
    2.13 -static const char rcsid[] =
    2.14 -    "$NetBSD: e_log.c,v 1.8 1995/05/10 20:45:49 jtc Exp $";
    2.15 -#endif
    2.16 -
    2.17  /* __ieee754_log(x)
    2.18   * Return the logrithm of x
    2.19   *
    2.20 @@ -69,99 +63,85 @@
    2.21  #include "math_libm.h"
    2.22  #include "math_private.h"
    2.23  
    2.24 -#ifdef __STDC__
    2.25  static const double
    2.26 +ln2_hi  =  6.93147180369123816490e-01,	/* 3fe62e42 fee00000 */
    2.27 +ln2_lo  =  1.90821492927058770002e-10,	/* 3dea39ef 35793c76 */
    2.28 +two54   =  1.80143985094819840000e+16,  /* 43500000 00000000 */
    2.29 +Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
    2.30 +Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
    2.31 +Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
    2.32 +Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
    2.33 +Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
    2.34 +Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
    2.35 +Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
    2.36 +
    2.37 +static const double zero   =  0.0;
    2.38 +
    2.39 +double attribute_hidden __ieee754_log(double x)
    2.40 +{
    2.41 +	double hfsq,f,s,z,R,w,t1,t2,dk;
    2.42 +	int32_t k,hx,i,j;
    2.43 +	u_int32_t lx;
    2.44 +
    2.45 +	EXTRACT_WORDS(hx,lx,x);
    2.46 +
    2.47 +	k=0;
    2.48 +	if (hx < 0x00100000) {			/* x < 2**-1022  */
    2.49 +	    if (((hx&0x7fffffff)|lx)==0)
    2.50 +		return -two54/zero;		/* log(+-0)=-inf */
    2.51 +	    if (hx<0) return (x-x)/zero;	/* log(-#) = NaN */
    2.52 +	    k -= 54; x *= two54; /* subnormal number, scale up x */
    2.53 +	    GET_HIGH_WORD(hx,x);
    2.54 +	}
    2.55 +	if (hx >= 0x7ff00000) return x+x;
    2.56 +	k += (hx>>20)-1023;
    2.57 +	hx &= 0x000fffff;
    2.58 +	i = (hx+0x95f64)&0x100000;
    2.59 +	SET_HIGH_WORD(x,hx|(i^0x3ff00000));	/* normalize x or x/2 */
    2.60 +	k += (i>>20);
    2.61 +	f = x-1.0;
    2.62 +	if((0x000fffff&(2+hx))<3) {	/* |f| < 2**-20 */
    2.63 +	    if(f==zero) {if(k==0) return zero;  else {dk=(double)k;
    2.64 +				 return dk*ln2_hi+dk*ln2_lo;}
    2.65 +	    }
    2.66 +	    R = f*f*(0.5-0.33333333333333333*f);
    2.67 +	    if(k==0) return f-R; else {dk=(double)k;
    2.68 +	    	     return dk*ln2_hi-((R-dk*ln2_lo)-f);}
    2.69 +	}
    2.70 + 	s = f/(2.0+f);
    2.71 +	dk = (double)k;
    2.72 +	z = s*s;
    2.73 +	i = hx-0x6147a;
    2.74 +	w = z*z;
    2.75 +	j = 0x6b851-hx;
    2.76 +	t1= w*(Lg2+w*(Lg4+w*Lg6));
    2.77 +	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
    2.78 +	i |= j;
    2.79 +	R = t2+t1;
    2.80 +	if(i>0) {
    2.81 +	    hfsq=0.5*f*f;
    2.82 +	    if(k==0) return f-(hfsq-s*(hfsq+R)); else
    2.83 +		     return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
    2.84 +	} else {
    2.85 +	    if(k==0) return f-s*(f-R); else
    2.86 +		     return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
    2.87 +	}
    2.88 +}
    2.89 +
    2.90 +/*
    2.91 + * wrapper log(x)
    2.92 + */
    2.93 +#ifndef _IEEE_LIBM
    2.94 +double log(double x)
    2.95 +{
    2.96 +	double z = __ieee754_log(x);
    2.97 +	if (_LIB_VERSION == _IEEE_ || isnan(x) || x > 0.0)
    2.98 +		return z;
    2.99 +	if (x == 0.0)
   2.100 +		return __kernel_standard(x, x, 16); /* log(0) */
   2.101 +	return __kernel_standard(x, x, 17); /* log(x<0) */
   2.102 +}
   2.103  #else
   2.104 -static double
   2.105 +strong_alias(__ieee754_log, log)
   2.106  #endif
   2.107 -  ln2_hi = 6.93147180369123816490e-01,  /* 3fe62e42 fee00000 */
   2.108 -    ln2_lo = 1.90821492927058770002e-10,        /* 3dea39ef 35793c76 */
   2.109 -    two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
   2.110 -    Lg1 = 6.666666666666735130e-01,     /* 3FE55555 55555593 */
   2.111 -    Lg2 = 3.999999999940941908e-01,     /* 3FD99999 9997FA04 */
   2.112 -    Lg3 = 2.857142874366239149e-01,     /* 3FD24924 94229359 */
   2.113 -    Lg4 = 2.222219843214978396e-01,     /* 3FCC71C5 1D8E78AF */
   2.114 -    Lg5 = 1.818357216161805012e-01,     /* 3FC74664 96CB03DE */
   2.115 -    Lg6 = 1.531383769920937332e-01,     /* 3FC39A09 D078C69F */
   2.116 -    Lg7 = 1.479819860511658591e-01;     /* 3FC2F112 DF3E5244 */
   2.117 -
   2.118 -#ifdef __STDC__
   2.119 -static const double zero = 0.0;
   2.120 -#else
   2.121 -static double zero = 0.0;
   2.122 -#endif
   2.123 -
   2.124 -#ifdef __STDC__
   2.125 -double attribute_hidden
   2.126 -__ieee754_log(double x)
   2.127 -#else
   2.128 -double attribute_hidden
   2.129 -__ieee754_log(x)
   2.130 -     double x;
   2.131 -#endif
   2.132 -{
   2.133 -    double hfsq, f, s, z, R, w, t1, t2, dk;
   2.134 -    int32_t k, hx, i, j;
   2.135 -    u_int32_t lx;
   2.136 -
   2.137 -    EXTRACT_WORDS(hx, lx, x);
   2.138 -
   2.139 -    k = 0;
   2.140 -    if (hx < 0x00100000) {      /* x < 2**-1022  */
   2.141 -        if (((hx & 0x7fffffff) | lx) == 0)
   2.142 -            return -two54 / zero;       /* log(+-0)=-inf */
   2.143 -        if (hx < 0)
   2.144 -            return (x - x) / zero;      /* log(-#) = NaN */
   2.145 -        k -= 54;
   2.146 -        x *= two54;             /* subnormal number, scale up x */
   2.147 -        GET_HIGH_WORD(hx, x);
   2.148 -    }
   2.149 -    if (hx >= 0x7ff00000)
   2.150 -        return x + x;
   2.151 -    k += (hx >> 20) - 1023;
   2.152 -    hx &= 0x000fffff;
   2.153 -    i = (hx + 0x95f64) & 0x100000;
   2.154 -    SET_HIGH_WORD(x, hx | (i ^ 0x3ff00000));    /* normalize x or x/2 */
   2.155 -    k += (i >> 20);
   2.156 -    f = x - 1.0;
   2.157 -    if ((0x000fffff & (2 + hx)) < 3) {  /* |f| < 2**-20 */
   2.158 -        if (f == zero) {
   2.159 -            if (k == 0)
   2.160 -                return zero;
   2.161 -            else {
   2.162 -                dk = (double) k;
   2.163 -                return dk * ln2_hi + dk * ln2_lo;
   2.164 -            }
   2.165 -        }
   2.166 -        R = f * f * (0.5 - 0.33333333333333333 * f);
   2.167 -        if (k == 0)
   2.168 -            return f - R;
   2.169 -        else {
   2.170 -            dk = (double) k;
   2.171 -            return dk * ln2_hi - ((R - dk * ln2_lo) - f);
   2.172 -        }
   2.173 -    }
   2.174 -    s = f / (2.0 + f);
   2.175 -    dk = (double) k;
   2.176 -    z = s * s;
   2.177 -    i = hx - 0x6147a;
   2.178 -    w = z * z;
   2.179 -    j = 0x6b851 - hx;
   2.180 -    t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
   2.181 -    t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
   2.182 -    i |= j;
   2.183 -    R = t2 + t1;
   2.184 -    if (i > 0) {
   2.185 -        hfsq = 0.5 * f * f;
   2.186 -        if (k == 0)
   2.187 -            return f - (hfsq - s * (hfsq + R));
   2.188 -        else
   2.189 -            return dk * ln2_hi - ((hfsq - (s * (hfsq + R) + dk * ln2_lo)) -
   2.190 -                                  f);
   2.191 -    } else {
   2.192 -        if (k == 0)
   2.193 -            return f - s * (f - R);
   2.194 -        else
   2.195 -            return dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f);
   2.196 -    }
   2.197 -}
   2.198 +libm_hidden_def(log)
     3.1 --- a/src/libm/e_pow.c	Sat Nov 04 15:34:14 2017 -0700
     3.2 +++ b/src/libm/e_pow.c	Sat Nov 04 15:53:19 2017 -0700
     3.3 @@ -1,4 +1,3 @@
     3.4 -/* @(#)e_pow.c 5.1 93/09/24 */
     3.5  /*
     3.6   * ====================================================
     3.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     3.8 @@ -10,10 +9,6 @@
     3.9   * ====================================================
    3.10   */
    3.11  
    3.12 -#if defined(LIBM_SCCS) && !defined(lint)
    3.13 -static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $";
    3.14 -#endif
    3.15 -
    3.16  /* __ieee754_pow(x,y) return x**y
    3.17   *
    3.18   *		      n
    3.19 @@ -26,25 +21,26 @@
    3.20   *	3. Return x**y = 2**n*exp(y'*log2)
    3.21   *
    3.22   * Special cases:
    3.23 - *	1.  (anything) ** 0  is 1
    3.24 - *	2.  (anything) ** 1  is itself
    3.25 - *	3.  (anything) ** NAN is NAN
    3.26 - *	4.  NAN ** (anything except 0) is NAN
    3.27 - *	5.  +-(|x| > 1) **  +INF is +INF
    3.28 - *	6.  +-(|x| > 1) **  -INF is +0
    3.29 - *	7.  +-(|x| < 1) **  +INF is +0
    3.30 - *	8.  +-(|x| < 1) **  -INF is +INF
    3.31 - *	9.  +-1         ** +-INF is NAN
    3.32 - *	10. +0 ** (+anything except 0, NAN)               is +0
    3.33 - *	11. -0 ** (+anything except 0, NAN, odd integer)  is +0
    3.34 - *	12. +0 ** (-anything except 0, NAN)               is +INF
    3.35 - *	13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
    3.36 - *	14. -0 ** (odd integer) = -( +0 ** (odd integer) )
    3.37 - *	15. +INF ** (+anything except 0,NAN) is +INF
    3.38 - *	16. +INF ** (-anything except 0,NAN) is +0
    3.39 - *	17. -INF ** (anything)  = -0 ** (-anything)
    3.40 - *	18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
    3.41 - *	19. (-anything except 0 and inf) ** (non-integer) is NAN
    3.42 + *	1.  +-1 ** anything  is 1.0
    3.43 + *	2.  +-1 ** +-INF     is 1.0
    3.44 + *	3.  (anything) ** 0  is 1
    3.45 + *	4.  (anything) ** 1  is itself
    3.46 + *	5.  (anything) ** NAN is NAN
    3.47 + *	6.  NAN ** (anything except 0) is NAN
    3.48 + *	7.  +-(|x| > 1) **  +INF is +INF
    3.49 + *	8.  +-(|x| > 1) **  -INF is +0
    3.50 + *	9.  +-(|x| < 1) **  +INF is +0
    3.51 + *	10  +-(|x| < 1) **  -INF is +INF
    3.52 + *	11. +0 ** (+anything except 0, NAN)               is +0
    3.53 + *	12. -0 ** (+anything except 0, NAN, odd integer)  is +0
    3.54 + *	13. +0 ** (-anything except 0, NAN)               is +INF
    3.55 + *	14. -0 ** (-anything except 0, NAN, odd integer)  is +INF
    3.56 + *	15. -0 ** (odd integer) = -( +0 ** (odd integer) )
    3.57 + *	16. +INF ** (+anything except 0,NAN) is +INF
    3.58 + *	17. +INF ** (-anything except 0,NAN) is +0
    3.59 + *	18. -INF ** (anything)  = -0 ** (-anything)
    3.60 + *	19. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
    3.61 + *	20. (-anything except 0 and inf) ** (non-integer) is NAN
    3.62   *
    3.63   * Accuracy:
    3.64   *	pow(x,y) returns x**y nearly rounded. In particular
    3.65 @@ -62,281 +58,281 @@
    3.66  #include "math_libm.h"
    3.67  #include "math_private.h"
    3.68  
    3.69 -libm_hidden_proto(scalbn)
    3.70 -    libm_hidden_proto(fabs)
    3.71 -#ifdef __STDC__
    3.72 -     static const double
    3.73 +static const double
    3.74 +bp[] = {1.0, 1.5,},
    3.75 +dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
    3.76 +dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
    3.77 +zero    =  0.0,
    3.78 +one	=  1.0,
    3.79 +two	=  2.0,
    3.80 +two53	=  9007199254740992.0,	/* 0x43400000, 0x00000000 */
    3.81 +huge	=  1.0e300,
    3.82 +tiny    =  1.0e-300,
    3.83 +	/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
    3.84 +L1  =  5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
    3.85 +L2  =  4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
    3.86 +L3  =  3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
    3.87 +L4  =  2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
    3.88 +L5  =  2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
    3.89 +L6  =  2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
    3.90 +P1   =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
    3.91 +P2   = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
    3.92 +P3   =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
    3.93 +P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
    3.94 +P5   =  4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
    3.95 +lg2  =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
    3.96 +lg2_h  =  6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
    3.97 +lg2_l  = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
    3.98 +ovt =  8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
    3.99 +cp    =  9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
   3.100 +cp_h  =  9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
   3.101 +cp_l  = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
   3.102 +ivln2    =  1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
   3.103 +ivln2_h  =  1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
   3.104 +ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
   3.105 +
   3.106 +double attribute_hidden __ieee754_pow(double x, double y)
   3.107 +{
   3.108 +	double z,ax,z_h,z_l,p_h,p_l;
   3.109 +	double y1,t1,t2,r,s,t,u,v,w;
   3.110 +	int32_t i,j,k,yisint,n;
   3.111 +	int32_t hx,hy,ix,iy;
   3.112 +	u_int32_t lx,ly;
   3.113 +
   3.114 +	EXTRACT_WORDS(hx,lx,x);
   3.115 +    /* x==1: 1**y = 1 (even if y is NaN) */
   3.116 +	if (hx==0x3ff00000 && lx==0) {
   3.117 +		return x;
   3.118 +	}
   3.119 +	ix = hx&0x7fffffff;
   3.120 +
   3.121 +	EXTRACT_WORDS(hy,ly,y);
   3.122 +	iy = hy&0x7fffffff;
   3.123 +
   3.124 +    /* y==zero: x**0 = 1 */
   3.125 +	if((iy|ly)==0) return one;
   3.126 +
   3.127 +    /* +-NaN return x+y */
   3.128 +	if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
   3.129 +	   iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
   3.130 +		return x+y;
   3.131 +
   3.132 +    /* determine if y is an odd int when x < 0
   3.133 +     * yisint = 0	... y is not an integer
   3.134 +     * yisint = 1	... y is an odd int
   3.135 +     * yisint = 2	... y is an even int
   3.136 +     */
   3.137 +	yisint  = 0;
   3.138 +	if(hx<0) {
   3.139 +	    if(iy>=0x43400000) yisint = 2; /* even integer y */
   3.140 +	    else if(iy>=0x3ff00000) {
   3.141 +		k = (iy>>20)-0x3ff;	   /* exponent */
   3.142 +		if(k>20) {
   3.143 +		    j = ly>>(52-k);
   3.144 +		    if((j<<(52-k))==ly) yisint = 2-(j&1);
   3.145 +		} else if(ly==0) {
   3.146 +		    j = iy>>(20-k);
   3.147 +		    if((j<<(20-k))==iy) yisint = 2-(j&1);
   3.148 +		}
   3.149 +	    }
   3.150 +	}
   3.151 +
   3.152 +    /* special value of y */
   3.153 +	if(ly==0) {
   3.154 +	    if (iy==0x7ff00000) {       /* y is +-inf */
   3.155 +	        if (((ix-0x3ff00000)|lx)==0)
   3.156 +		    return one;	        /* +-1**+-inf is 1 (yes, weird rule) */
   3.157 +	        if (ix >= 0x3ff00000)   /* (|x|>1)**+-inf = inf,0 */
   3.158 +		    return (hy>=0) ? y : zero;
   3.159 +	        /* (|x|<1)**-,+inf = inf,0 */
   3.160 +		return (hy<0) ? -y : zero;
   3.161 +	    }
   3.162 +	    if(iy==0x3ff00000) {	/* y is  +-1 */
   3.163 +		if(hy<0) return one/x; else return x;
   3.164 +	    }
   3.165 +	    if(hy==0x40000000) return x*x; /* y is  2 */
   3.166 +	    if(hy==0x3fe00000) {	/* y is  0.5 */
   3.167 +		if(hx>=0)	/* x >= +0 */
   3.168 +		    return __ieee754_sqrt(x);
   3.169 +	    }
   3.170 +	}
   3.171 +
   3.172 +	ax   = fabs(x);
   3.173 +    /* special value of x */
   3.174 +	if(lx==0) {
   3.175 +	    if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
   3.176 +		z = ax;			/*x is +-0,+-inf,+-1*/
   3.177 +		if(hy<0) z = one/z;	/* z = (1/|x|) */
   3.178 +		if(hx<0) {
   3.179 +		    if(((ix-0x3ff00000)|yisint)==0) {
   3.180 +			z = (z-z)/(z-z); /* (-1)**non-int is NaN */
   3.181 +		    } else if(yisint==1)
   3.182 +			z = -z;		/* (x<0)**odd = -(|x|**odd) */
   3.183 +		}
   3.184 +		return z;
   3.185 +	    }
   3.186 +	}
   3.187 +
   3.188 +    /* (x<0)**(non-int) is NaN */
   3.189 +	if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
   3.190 +
   3.191 +    /* |y| is huge */
   3.192 +	if(iy>0x41e00000) { /* if |y| > 2**31 */
   3.193 +	    if(iy>0x43f00000){	/* if |y| > 2**64, must o/uflow */
   3.194 +		if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
   3.195 +		if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
   3.196 +	    }
   3.197 +	/* over/underflow if x is not close to one */
   3.198 +	    if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
   3.199 +	    if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
   3.200 +	/* now |1-x| is tiny <= 2**-20, suffice to compute
   3.201 +	   log(x) by x-x^2/2+x^3/3-x^4/4 */
   3.202 +	    t = x-1;		/* t has 20 trailing zeros */
   3.203 +	    w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
   3.204 +	    u = ivln2_h*t;	/* ivln2_h has 21 sig. bits */
   3.205 +	    v = t*ivln2_l-w*ivln2;
   3.206 +	    t1 = u+v;
   3.207 +	    SET_LOW_WORD(t1,0);
   3.208 +	    t2 = v-(t1-u);
   3.209 +	} else {
   3.210 +	    double s2,s_h,s_l,t_h,t_l;
   3.211 +	    n = 0;
   3.212 +	/* take care subnormal number */
   3.213 +	    if(ix<0x00100000)
   3.214 +		{ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }
   3.215 +	    n  += ((ix)>>20)-0x3ff;
   3.216 +	    j  = ix&0x000fffff;
   3.217 +	/* determine interval */
   3.218 +	    ix = j|0x3ff00000;		/* normalize ix */
   3.219 +	    if(j<=0x3988E) k=0;		/* |x|<sqrt(3/2) */
   3.220 +	    else if(j<0xBB67A) k=1;	/* |x|<sqrt(3)   */
   3.221 +	    else {k=0;n+=1;ix -= 0x00100000;}
   3.222 +	    SET_HIGH_WORD(ax,ix);
   3.223 +
   3.224 +	/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
   3.225 +	    u = ax-bp[k];		/* bp[0]=1.0, bp[1]=1.5 */
   3.226 +	    v = one/(ax+bp[k]);
   3.227 +	    s = u*v;
   3.228 +	    s_h = s;
   3.229 +	    SET_LOW_WORD(s_h,0);
   3.230 +	/* t_h=ax+bp[k] High */
   3.231 +	    t_h = zero;
   3.232 +	    SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18));
   3.233 +	    t_l = ax - (t_h-bp[k]);
   3.234 +	    s_l = v*((u-s_h*t_h)-s_h*t_l);
   3.235 +	/* compute log(ax) */
   3.236 +	    s2 = s*s;
   3.237 +	    r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
   3.238 +	    r += s_l*(s_h+s);
   3.239 +	    s2  = s_h*s_h;
   3.240 +	    t_h = 3.0+s2+r;
   3.241 +	    SET_LOW_WORD(t_h,0);
   3.242 +	    t_l = r-((t_h-3.0)-s2);
   3.243 +	/* u+v = s*(1+...) */
   3.244 +	    u = s_h*t_h;
   3.245 +	    v = s_l*t_h+t_l*s;
   3.246 +	/* 2/(3log2)*(s+...) */
   3.247 +	    p_h = u+v;
   3.248 +	    SET_LOW_WORD(p_h,0);
   3.249 +	    p_l = v-(p_h-u);
   3.250 +	    z_h = cp_h*p_h;		/* cp_h+cp_l = 2/(3*log2) */
   3.251 +	    z_l = cp_l*p_h+p_l*cp+dp_l[k];
   3.252 +	/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
   3.253 +	    t = (double)n;
   3.254 +	    t1 = (((z_h+z_l)+dp_h[k])+t);
   3.255 +	    SET_LOW_WORD(t1,0);
   3.256 +	    t2 = z_l-(((t1-t)-dp_h[k])-z_h);
   3.257 +	}
   3.258 +
   3.259 +	s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
   3.260 +	if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0)
   3.261 +	    s = -one;/* (-ve)**(odd int) */
   3.262 +
   3.263 +    /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
   3.264 +	y1  = y;
   3.265 +	SET_LOW_WORD(y1,0);
   3.266 +	p_l = (y-y1)*t1+y*t2;
   3.267 +	p_h = y1*t1;
   3.268 +	z = p_l+p_h;
   3.269 +	EXTRACT_WORDS(j,i,z);
   3.270 +	if (j>=0x40900000) {				/* z >= 1024 */
   3.271 +	    if(((j-0x40900000)|i)!=0)			/* if z > 1024 */
   3.272 +		return s*huge*huge;			/* overflow */
   3.273 +	    else {
   3.274 +		if(p_l+ovt>z-p_h) return s*huge*huge;	/* overflow */
   3.275 +	    }
   3.276 +	} else if((j&0x7fffffff)>=0x4090cc00 ) {	/* z <= -1075 */
   3.277 +	    if(((j-0xc090cc00)|i)!=0) 		/* z < -1075 */
   3.278 +		return s*tiny*tiny;		/* underflow */
   3.279 +	    else {
   3.280 +		if(p_l<=z-p_h) return s*tiny*tiny;	/* underflow */
   3.281 +	    }
   3.282 +	}
   3.283 +    /*
   3.284 +     * compute 2**(p_h+p_l)
   3.285 +     */
   3.286 +	i = j&0x7fffffff;
   3.287 +	k = (i>>20)-0x3ff;
   3.288 +	n = 0;
   3.289 +	if(i>0x3fe00000) {		/* if |z| > 0.5, set n = [z+0.5] */
   3.290 +	    n = j+(0x00100000>>(k+1));
   3.291 +	    k = ((n&0x7fffffff)>>20)-0x3ff;	/* new k for n */
   3.292 +	    t = zero;
   3.293 +	    SET_HIGH_WORD(t,n&~(0x000fffff>>k));
   3.294 +	    n = ((n&0x000fffff)|0x00100000)>>(20-k);
   3.295 +	    if(j<0) n = -n;
   3.296 +	    p_h -= t;
   3.297 +	}
   3.298 +	t = p_l+p_h;
   3.299 +	SET_LOW_WORD(t,0);
   3.300 +	u = t*lg2_h;
   3.301 +	v = (p_l-(t-p_h))*lg2+t*lg2_l;
   3.302 +	z = u+v;
   3.303 +	w = v-(z-u);
   3.304 +	t  = z*z;
   3.305 +	t1  = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
   3.306 +	r  = (z*t1)/(t1-two)-(w+z*w);
   3.307 +	z  = one-(r-z);
   3.308 +	GET_HIGH_WORD(j,z);
   3.309 +	j += (n<<20);
   3.310 +	if((j>>20)<=0) z = scalbn(z,n);	/* subnormal output */
   3.311 +	else SET_HIGH_WORD(z,j);
   3.312 +	return s*z;
   3.313 +}
   3.314 +
   3.315 +/*
   3.316 + * wrapper pow(x,y) return x**y
   3.317 + */
   3.318 +#ifndef _IEEE_LIBM
   3.319 +double pow(double x, double y)
   3.320 +{
   3.321 +	double z = __ieee754_pow(x, y);
   3.322 +	if (_LIB_VERSION == _IEEE_|| isnan(y))
   3.323 +		return z;
   3.324 +	if (isnan(x)) {
   3.325 +		if (y == 0.0)
   3.326 +			return __kernel_standard(x, y, 42); /* pow(NaN,0.0) */
   3.327 +		return z;
   3.328 +	}
   3.329 +	if (x == 0.0) {
   3.330 +		if (y == 0.0)
   3.331 +	    		return __kernel_standard(x, y, 20); /* pow(0.0,0.0) */
   3.332 +		if (isfinite(y) && y < 0.0)
   3.333 +			return __kernel_standard(x,y,23); /* pow(0.0,negative) */
   3.334 +		return z;
   3.335 +	}
   3.336 +	if (!isfinite(z)) {
   3.337 +		if (isfinite(x) && isfinite(y)) {
   3.338 +			if (isnan(z))
   3.339 +				return __kernel_standard(x, y, 24); /* pow neg**non-int */
   3.340 +			return __kernel_standard(x, y, 21); /* pow overflow */
   3.341 +		}
   3.342 +	}
   3.343 +	if (z == 0.0 && isfinite(x) && isfinite(y))
   3.344 +		return __kernel_standard(x, y, 22); /* pow underflow */
   3.345 +	return z;
   3.346 +}
   3.347  #else
   3.348 -     static double
   3.349 +strong_alias(__ieee754_pow, pow)
   3.350  #endif
   3.351 -       bp[] = { 1.0, 1.5, }, dp_h[] = {
   3.352 -     0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
   3.353 -
   3.354 -         dp_l[] = {
   3.355 -     0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
   3.356 -
   3.357 -         zero = 0.0, one = 1.0, two = 2.0, two53 = 9007199254740992.0,  /* 0x43400000, 0x00000000 */
   3.358 -         huge_val = 1.0e300, tiny = 1.0e-300,
   3.359 -         /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
   3.360 -         L1 = 5.99999999999994648725e-01,       /* 0x3FE33333, 0x33333303 */
   3.361 -         L2 = 4.28571428578550184252e-01,       /* 0x3FDB6DB6, 0xDB6FABFF */
   3.362 -         L3 = 3.33333329818377432918e-01,       /* 0x3FD55555, 0x518F264D */
   3.363 -         L4 = 2.72728123808534006489e-01,       /* 0x3FD17460, 0xA91D4101 */
   3.364 -         L5 = 2.30660745775561754067e-01,       /* 0x3FCD864A, 0x93C9DB65 */
   3.365 -         L6 = 2.06975017800338417784e-01,       /* 0x3FCA7E28, 0x4A454EEF */
   3.366 -         P1 = 1.66666666666666019037e-01,       /* 0x3FC55555, 0x5555553E */
   3.367 -         P2 = -2.77777777770155933842e-03,      /* 0xBF66C16C, 0x16BEBD93 */
   3.368 -         P3 = 6.61375632143793436117e-05,       /* 0x3F11566A, 0xAF25DE2C */
   3.369 -         P4 = -1.65339022054652515390e-06,      /* 0xBEBBBD41, 0xC5D26BF1 */
   3.370 -         P5 = 4.13813679705723846039e-08,       /* 0x3E663769, 0x72BEA4D0 */
   3.371 -         lg2 = 6.93147180559945286227e-01,      /* 0x3FE62E42, 0xFEFA39EF */
   3.372 -         lg2_h = 6.93147182464599609375e-01,    /* 0x3FE62E43, 0x00000000 */
   3.373 -         lg2_l = -1.90465429995776804525e-09,   /* 0xBE205C61, 0x0CA86C39 */
   3.374 -         ovt = 8.0085662595372944372e-0017,     /* -(1024-log2(ovfl+.5ulp)) */
   3.375 -         cp = 9.61796693925975554329e-01,       /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
   3.376 -         cp_h = 9.61796700954437255859e-01,     /* 0x3FEEC709, 0xE0000000 =(float)cp */
   3.377 -         cp_l = -7.02846165095275826516e-09,    /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h */
   3.378 -         ivln2 = 1.44269504088896338700e+00,    /* 0x3FF71547, 0x652B82FE =1/ln2 */
   3.379 -         ivln2_h = 1.44269502162933349609e+00,  /* 0x3FF71547, 0x60000000 =24b 1/ln2 */
   3.380 -         ivln2_l = 1.92596299112661746887e-08;  /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail */
   3.381 -
   3.382 -#ifdef __STDC__
   3.383 -     double attribute_hidden __ieee754_pow(double x, double y)
   3.384 -#else
   3.385 -     double attribute_hidden __ieee754_pow(x, y)
   3.386 -     double x, y;
   3.387 -#endif
   3.388 -     {
   3.389 -         double z, ax, z_h, z_l, p_h, p_l;
   3.390 -         double y1, t1, t2, r, s, t, u, v, w;
   3.391 -         int32_t i, j, k, yisint, n;
   3.392 -         int32_t hx, hy, ix, iy;
   3.393 -         u_int32_t lx, ly;
   3.394 -
   3.395 -         EXTRACT_WORDS(hx, lx, x);
   3.396 -         EXTRACT_WORDS(hy, ly, y);
   3.397 -         ix = hx & 0x7fffffff;
   3.398 -         iy = hy & 0x7fffffff;
   3.399 -
   3.400 -         /* y==zero: x**0 = 1 */
   3.401 -         if ((iy | ly) == 0)
   3.402 -             return one;
   3.403 -
   3.404 -         /* +-NaN return x+y */
   3.405 -         if (ix > 0x7ff00000 || ((ix == 0x7ff00000) && (lx != 0)) ||
   3.406 -             iy > 0x7ff00000 || ((iy == 0x7ff00000) && (ly != 0)))
   3.407 -             return x + y;
   3.408 -
   3.409 -         /* determine if y is an odd int when x < 0
   3.410 -          * yisint = 0       ... y is not an integer
   3.411 -          * yisint = 1       ... y is an odd int
   3.412 -          * yisint = 2       ... y is an even int
   3.413 -          */
   3.414 -         yisint = 0;
   3.415 -         if (hx < 0) {
   3.416 -             if (iy >= 0x43400000)
   3.417 -                 yisint = 2;    /* even integer y */
   3.418 -             else if (iy >= 0x3ff00000) {
   3.419 -                 k = (iy >> 20) - 0x3ff;        /* exponent */
   3.420 -                 if (k > 20) {
   3.421 -                     j = ly >> (52 - k);
   3.422 -                     if ((j << (52 - k)) == ly)
   3.423 -                         yisint = 2 - (j & 1);
   3.424 -                 } else if (ly == 0) {
   3.425 -                     j = iy >> (20 - k);
   3.426 -                     if ((j << (20 - k)) == iy)
   3.427 -                         yisint = 2 - (j & 1);
   3.428 -                 }
   3.429 -             }
   3.430 -         }
   3.431 -
   3.432 -         /* special value of y */
   3.433 -         if (ly == 0) {
   3.434 -             if (iy == 0x7ff00000) {    /* y is +-inf */
   3.435 -                 if (((ix - 0x3ff00000) | lx) == 0)
   3.436 -                     return y - y;      /* inf**+-1 is NaN */
   3.437 -                 else if (ix >= 0x3ff00000)     /* (|x|>1)**+-inf = inf,0 */
   3.438 -                     return (hy >= 0) ? y : zero;
   3.439 -                 else           /* (|x|<1)**-,+inf = inf,0 */
   3.440 -                     return (hy < 0) ? -y : zero;
   3.441 -             }
   3.442 -             if (iy == 0x3ff00000) {    /* y is  +-1 */
   3.443 -                 if (hy < 0)
   3.444 -                     return one / x;
   3.445 -                 else
   3.446 -                     return x;
   3.447 -             }
   3.448 -             if (hy == 0x40000000)
   3.449 -                 return x * x;  /* y is  2 */
   3.450 -             if (hy == 0x3fe00000) {    /* y is  0.5 */
   3.451 -                 if (hx >= 0)   /* x >= +0 */
   3.452 -                     return __ieee754_sqrt(x);
   3.453 -             }
   3.454 -         }
   3.455 -
   3.456 -         ax = fabs(x);
   3.457 -         /* special value of x */
   3.458 -         if (lx == 0) {
   3.459 -             if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) {
   3.460 -                 z = ax;        /* x is +-0,+-inf,+-1 */
   3.461 -                 if (hy < 0)
   3.462 -                     z = one / z;       /* z = (1/|x|) */
   3.463 -                 if (hx < 0) {
   3.464 -                     if (((ix - 0x3ff00000) | yisint) == 0) {
   3.465 -                         z = (z - z) / (z - z); /* (-1)**non-int is NaN */
   3.466 -                     } else if (yisint == 1)
   3.467 -                         z = -z;        /* (x<0)**odd = -(|x|**odd) */
   3.468 -                 }
   3.469 -                 return z;
   3.470 -             }
   3.471 -         }
   3.472 -
   3.473 -         /* (x<0)**(non-int) is NaN */
   3.474 -         if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0)
   3.475 -             return (x - x) / (x - x);
   3.476 -
   3.477 -         /* |y| is huge */
   3.478 -         if (iy > 0x41e00000) { /* if |y| > 2**31 */
   3.479 -             if (iy > 0x43f00000) {     /* if |y| > 2**64, must o/uflow */
   3.480 -                 if (ix <= 0x3fefffff)
   3.481 -                     return (hy < 0) ? huge_val * huge_val : tiny * tiny;
   3.482 -                 if (ix >= 0x3ff00000)
   3.483 -                     return (hy > 0) ? huge_val * huge_val : tiny * tiny;
   3.484 -             }
   3.485 -             /* over/underflow if x is not close to one */
   3.486 -             if (ix < 0x3fefffff)
   3.487 -                 return (hy < 0) ? huge_val * huge_val : tiny * tiny;
   3.488 -             if (ix > 0x3ff00000)
   3.489 -                 return (hy > 0) ? huge_val * huge_val : tiny * tiny;
   3.490 -             /* now |1-x| is tiny <= 2**-20, suffice to compute
   3.491 -                log(x) by x-x^2/2+x^3/3-x^4/4 */
   3.492 -             t = x - 1;         /* t has 20 trailing zeros */
   3.493 -             w = (t * t) * (0.5 - t * (0.3333333333333333333333 - t * 0.25));
   3.494 -             u = ivln2_h * t;   /* ivln2_h has 21 sig. bits */
   3.495 -             v = t * ivln2_l - w * ivln2;
   3.496 -             t1 = u + v;
   3.497 -             SET_LOW_WORD(t1, 0);
   3.498 -             t2 = v - (t1 - u);
   3.499 -         } else {
   3.500 -             double s2, s_h, s_l, t_h, t_l;
   3.501 -             n = 0;
   3.502 -             /* take care subnormal number */
   3.503 -             if (ix < 0x00100000) {
   3.504 -                 ax *= two53;
   3.505 -                 n -= 53;
   3.506 -                 GET_HIGH_WORD(ix, ax);
   3.507 -             }
   3.508 -             n += ((ix) >> 20) - 0x3ff;
   3.509 -             j = ix & 0x000fffff;
   3.510 -             /* determine interval */
   3.511 -             ix = j | 0x3ff00000;       /* normalize ix */
   3.512 -             if (j <= 0x3988E)
   3.513 -                 k = 0;         /* |x|<sqrt(3/2) */
   3.514 -             else if (j < 0xBB67A)
   3.515 -                 k = 1;         /* |x|<sqrt(3)   */
   3.516 -             else {
   3.517 -                 k = 0;
   3.518 -                 n += 1;
   3.519 -                 ix -= 0x00100000;
   3.520 -             }
   3.521 -             SET_HIGH_WORD(ax, ix);
   3.522 -
   3.523 -             /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
   3.524 -             u = ax - bp[k];    /* bp[0]=1.0, bp[1]=1.5 */
   3.525 -             v = one / (ax + bp[k]);
   3.526 -             s = u * v;
   3.527 -             s_h = s;
   3.528 -             SET_LOW_WORD(s_h, 0);
   3.529 -             /* t_h=ax+bp[k] High */
   3.530 -             t_h = zero;
   3.531 -             SET_HIGH_WORD(t_h,
   3.532 -                           ((ix >> 1) | 0x20000000) + 0x00080000 + (k << 18));
   3.533 -             t_l = ax - (t_h - bp[k]);
   3.534 -             s_l = v * ((u - s_h * t_h) - s_h * t_l);
   3.535 -             /* compute log(ax) */
   3.536 -             s2 = s * s;
   3.537 -             r = s2 * s2 * (L1 +
   3.538 -                            s2 * (L2 +
   3.539 -                                  s2 * (L3 +
   3.540 -                                        s2 * (L4 + s2 * (L5 + s2 * L6)))));
   3.541 -             r += s_l * (s_h + s);
   3.542 -             s2 = s_h * s_h;
   3.543 -             t_h = 3.0 + s2 + r;
   3.544 -             SET_LOW_WORD(t_h, 0);
   3.545 -             t_l = r - ((t_h - 3.0) - s2);
   3.546 -             /* u+v = s*(1+...) */
   3.547 -             u = s_h * t_h;
   3.548 -             v = s_l * t_h + t_l * s;
   3.549 -             /* 2/(3log2)*(s+...) */
   3.550 -             p_h = u + v;
   3.551 -             SET_LOW_WORD(p_h, 0);
   3.552 -             p_l = v - (p_h - u);
   3.553 -             z_h = cp_h * p_h;  /* cp_h+cp_l = 2/(3*log2) */
   3.554 -             z_l = cp_l * p_h + p_l * cp + dp_l[k];
   3.555 -             /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
   3.556 -             t = (double) n;
   3.557 -             t1 = (((z_h + z_l) + dp_h[k]) + t);
   3.558 -             SET_LOW_WORD(t1, 0);
   3.559 -             t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
   3.560 -         }
   3.561 -
   3.562 -         s = one;               /* s (sign of result -ve**odd) = -1 else = 1 */
   3.563 -         if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
   3.564 -             s = -one;          /* (-ve)**(odd int) */
   3.565 -
   3.566 -         /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
   3.567 -         y1 = y;
   3.568 -         SET_LOW_WORD(y1, 0);
   3.569 -         p_l = (y - y1) * t1 + y * t2;
   3.570 -         p_h = y1 * t1;
   3.571 -         z = p_l + p_h;
   3.572 -         EXTRACT_WORDS(j, i, z);
   3.573 -         if (j >= 0x40900000) { /* z >= 1024 */
   3.574 -             if (((j - 0x40900000) | i) != 0)   /* if z > 1024 */
   3.575 -                 return s * huge_val * huge_val;        /* overflow */
   3.576 -             else {
   3.577 -                 if (p_l + ovt > z - p_h)
   3.578 -                     return s * huge_val * huge_val;    /* overflow */
   3.579 -             }
   3.580 -         } else if ((j & 0x7fffffff) >= 0x4090cc00) {   /* z <= -1075 */
   3.581 -             if (((j - 0xc090cc00) | i) != 0)   /* z < -1075 */
   3.582 -                 return s * tiny * tiny;        /* underflow */
   3.583 -             else {
   3.584 -                 if (p_l <= z - p_h)
   3.585 -                     return s * tiny * tiny;    /* underflow */
   3.586 -             }
   3.587 -         }
   3.588 -         /*
   3.589 -          * compute 2**(p_h+p_l)
   3.590 -          */
   3.591 -         i = j & 0x7fffffff;
   3.592 -         k = (i >> 20) - 0x3ff;
   3.593 -         n = 0;
   3.594 -         if (i > 0x3fe00000) {  /* if |z| > 0.5, set n = [z+0.5] */
   3.595 -             n = j + (0x00100000 >> (k + 1));
   3.596 -             k = ((n & 0x7fffffff) >> 20) - 0x3ff;      /* new k for n */
   3.597 -             t = zero;
   3.598 -             SET_HIGH_WORD(t, n & ~(0x000fffff >> k));
   3.599 -             n = ((n & 0x000fffff) | 0x00100000) >> (20 - k);
   3.600 -             if (j < 0)
   3.601 -                 n = -n;
   3.602 -             p_h -= t;
   3.603 -         }
   3.604 -         t = p_l + p_h;
   3.605 -         SET_LOW_WORD(t, 0);
   3.606 -         u = t * lg2_h;
   3.607 -         v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
   3.608 -         z = u + v;
   3.609 -         w = v - (z - u);
   3.610 -         t = z * z;
   3.611 -         t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
   3.612 -         r = (z * t1) / (t1 - two) - (w + z * w);
   3.613 -         z = one - (r - z);
   3.614 -         GET_HIGH_WORD(j, z);
   3.615 -         j += (n << 20);
   3.616 -         if ((j >> 20) <= 0)
   3.617 -             z = scalbn(z, n);  /* subnormal output */
   3.618 -         else
   3.619 -             SET_HIGH_WORD(z, j);
   3.620 -         return s * z;
   3.621 -     }
   3.622 +libm_hidden_def(pow)
     4.1 --- a/src/libm/e_rem_pio2.c	Sat Nov 04 15:34:14 2017 -0700
     4.2 +++ b/src/libm/e_rem_pio2.c	Sat Nov 04 15:53:19 2017 -0700
     4.3 @@ -1,4 +1,3 @@
     4.4 -/* @(#)e_rem_pio2.c 5.1 93/09/24 */
     4.5  /*
     4.6   * ====================================================
     4.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     4.8 @@ -10,11 +9,6 @@
     4.9   * ====================================================
    4.10   */
    4.11  
    4.12 -#if defined(LIBM_SCCS) && !defined(lint)
    4.13 -static const char rcsid[] =
    4.14 -    "$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $";
    4.15 -#endif
    4.16 -
    4.17  /* __ieee754_rem_pio2(x,y)
    4.18   *
    4.19   * return the remainder of x rem pi/2 in y[0]+y[1]
    4.20 @@ -24,42 +18,30 @@
    4.21  #include "math_libm.h"
    4.22  #include "math_private.h"
    4.23  
    4.24 -#include "SDL_assert.h"
    4.25 -
    4.26 -libm_hidden_proto(fabs)
    4.27 -
    4.28  /*
    4.29   * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
    4.30   */
    4.31 -#ifdef __STDC__
    4.32 -     static const int32_t two_over_pi[] = {
    4.33 -#else
    4.34 -     static int32_t two_over_pi[] = {
    4.35 -#endif
    4.36 -         0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
    4.37 -         0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
    4.38 -         0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
    4.39 -         0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
    4.40 -         0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
    4.41 -         0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
    4.42 -         0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
    4.43 -         0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
    4.44 -         0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
    4.45 -         0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
    4.46 -         0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
    4.47 -     };
    4.48 +static const int32_t two_over_pi[] = {
    4.49 +0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
    4.50 +0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
    4.51 +0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
    4.52 +0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
    4.53 +0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
    4.54 +0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
    4.55 +0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
    4.56 +0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
    4.57 +0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
    4.58 +0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
    4.59 +0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
    4.60 +};
    4.61  
    4.62 -#ifdef __STDC__
    4.63  static const int32_t npio2_hw[] = {
    4.64 -#else
    4.65 -static int32_t npio2_hw[] = {
    4.66 -#endif
    4.67 -    0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
    4.68 -    0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
    4.69 -    0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
    4.70 -    0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
    4.71 -    0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
    4.72 -    0x404858EB, 0x404921FB,
    4.73 +0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
    4.74 +0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
    4.75 +0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
    4.76 +0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
    4.77 +0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
    4.78 +0x404858EB, 0x404921FB,
    4.79  };
    4.80  
    4.81  /*
    4.82 @@ -72,137 +54,108 @@
    4.83   * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
    4.84   */
    4.85  
    4.86 -#ifdef __STDC__
    4.87  static const double
    4.88 -#else
    4.89 -static double
    4.90 -#endif
    4.91 -  zero = 0.00000000000000000000e+00,    /* 0x00000000, 0x00000000 */
    4.92 -    half = 5.00000000000000000000e-01,  /* 0x3FE00000, 0x00000000 */
    4.93 -    two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
    4.94 -    invpio2 = 6.36619772367581382433e-01,       /* 0x3FE45F30, 0x6DC9C883 */
    4.95 -    pio2_1 = 1.57079632673412561417e+00,        /* 0x3FF921FB, 0x54400000 */
    4.96 -    pio2_1t = 6.07710050650619224932e-11,       /* 0x3DD0B461, 0x1A626331 */
    4.97 -    pio2_2 = 6.07710050630396597660e-11,        /* 0x3DD0B461, 0x1A600000 */
    4.98 -    pio2_2t = 2.02226624879595063154e-21,       /* 0x3BA3198A, 0x2E037073 */
    4.99 -    pio2_3 = 2.02226624871116645580e-21,        /* 0x3BA3198A, 0x2E000000 */
   4.100 -    pio2_3t = 8.47842766036889956997e-32;       /* 0x397B839A, 0x252049C1 */
   4.101 +zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
   4.102 +half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
   4.103 +two24 =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
   4.104 +invpio2 =  6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
   4.105 +pio2_1  =  1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
   4.106 +pio2_1t =  6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
   4.107 +pio2_2  =  6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
   4.108 +pio2_2t =  2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
   4.109 +pio2_3  =  2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
   4.110 +pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
   4.111  
   4.112 -#ifdef __STDC__
   4.113 -int32_t attribute_hidden
   4.114 -__ieee754_rem_pio2(double x, double *y)
   4.115 -#else
   4.116 -int32_t attribute_hidden
   4.117 -__ieee754_rem_pio2(x, y)
   4.118 -     double x, y[];
   4.119 -#endif
   4.120 +int32_t attribute_hidden __ieee754_rem_pio2(double x, double *y)
   4.121  {
   4.122 -    double z = 0.0, w, t, r, fn;
   4.123 -    double tx[3];
   4.124 -    int32_t e0, i, j, nx, n, ix, hx;
   4.125 -    u_int32_t low;
   4.126 +	double z=0.0,w,t,r,fn;
   4.127 +	double tx[3];
   4.128 +	int32_t e0,i,j,nx,n,ix,hx;
   4.129 +	u_int32_t low;
   4.130  
   4.131 -    GET_HIGH_WORD(hx, x);       /* high word of x */
   4.132 -    ix = hx & 0x7fffffff;
   4.133 -    if (ix <= 0x3fe921fb) {     /* |x| ~<= pi/4 , no need for reduction */
   4.134 -        y[0] = x;
   4.135 -        y[1] = 0;
   4.136 -        return 0;
   4.137 -    }
   4.138 -    if (ix < 0x4002d97c) {      /* |x| < 3pi/4, special case with n=+-1 */
   4.139 -        if (hx > 0) {
   4.140 -            z = x - pio2_1;
   4.141 -            if (ix != 0x3ff921fb) {     /* 33+53 bit pi is good enough */
   4.142 -                y[0] = z - pio2_1t;
   4.143 -                y[1] = (z - y[0]) - pio2_1t;
   4.144 -            } else {            /* near pi/2, use 33+33+53 bit pi */
   4.145 -                z -= pio2_2;
   4.146 -                y[0] = z - pio2_2t;
   4.147 -                y[1] = (z - y[0]) - pio2_2t;
   4.148 -            }
   4.149 -            return 1;
   4.150 -        } else {                /* negative x */
   4.151 -            z = x + pio2_1;
   4.152 -            if (ix != 0x3ff921fb) {     /* 33+53 bit pi is good enough */
   4.153 -                y[0] = z + pio2_1t;
   4.154 -                y[1] = (z - y[0]) + pio2_1t;
   4.155 -            } else {            /* near pi/2, use 33+33+53 bit pi */
   4.156 -                z += pio2_2;
   4.157 -                y[0] = z + pio2_2t;
   4.158 -                y[1] = (z - y[0]) + pio2_2t;
   4.159 -            }
   4.160 -            return -1;
   4.161 -        }
   4.162 -    }
   4.163 -    if (ix <= 0x413921fb) {     /* |x| ~<= 2^19*(pi/2), medium size */
   4.164 -        t = fabs(x);
   4.165 -        n = (int32_t) (t * invpio2 + half);
   4.166 -        fn = (double) n;
   4.167 -        r = t - fn * pio2_1;
   4.168 -        w = fn * pio2_1t;       /* 1st round good to 85 bit */
   4.169 -        if (n < 32 && ix != npio2_hw[n - 1]) {
   4.170 -            y[0] = r - w;       /* quick check no cancellation */
   4.171 -        } else {
   4.172 -            u_int32_t high;
   4.173 -            j = ix >> 20;
   4.174 -            y[0] = r - w;
   4.175 -            GET_HIGH_WORD(high, y[0]);
   4.176 -            i = j - ((high >> 20) & 0x7ff);
   4.177 -            if (i > 16) {       /* 2nd iteration needed, good to 118 */
   4.178 -                t = r;
   4.179 -                w = fn * pio2_2;
   4.180 -                r = t - w;
   4.181 -                w = fn * pio2_2t - ((t - r) - w);
   4.182 -                y[0] = r - w;
   4.183 -                GET_HIGH_WORD(high, y[0]);
   4.184 -                i = j - ((high >> 20) & 0x7ff);
   4.185 -                if (i > 49) {   /* 3rd iteration need, 151 bits acc */
   4.186 -                    t = r;      /* will cover all possible cases */
   4.187 -                    w = fn * pio2_3;
   4.188 -                    r = t - w;
   4.189 -                    w = fn * pio2_3t - ((t - r) - w);
   4.190 -                    y[0] = r - w;
   4.191 -                }
   4.192 -            }
   4.193 -        }
   4.194 -        y[1] = (r - y[0]) - w;
   4.195 -        if (hx < 0) {
   4.196 -            y[0] = -y[0];
   4.197 -            y[1] = -y[1];
   4.198 -            return -n;
   4.199 -        } else
   4.200 -            return n;
   4.201 -    }
   4.202 +	GET_HIGH_WORD(hx,x);		/* high word of x */
   4.203 +	ix = hx&0x7fffffff;
   4.204 +	if(ix<=0x3fe921fb)   /* |x| ~<= pi/4 , no need for reduction */
   4.205 +	    {y[0] = x; y[1] = 0; return 0;}
   4.206 +	if(ix<0x4002d97c) {  /* |x| < 3pi/4, special case with n=+-1 */
   4.207 +	    if(hx>0) {
   4.208 +		z = x - pio2_1;
   4.209 +		if(ix!=0x3ff921fb) { 	/* 33+53 bit pi is good enough */
   4.210 +		    y[0] = z - pio2_1t;
   4.211 +		    y[1] = (z-y[0])-pio2_1t;
   4.212 +		} else {		/* near pi/2, use 33+33+53 bit pi */
   4.213 +		    z -= pio2_2;
   4.214 +		    y[0] = z - pio2_2t;
   4.215 +		    y[1] = (z-y[0])-pio2_2t;
   4.216 +		}
   4.217 +		return 1;
   4.218 +	    } else {	/* negative x */
   4.219 +		z = x + pio2_1;
   4.220 +		if(ix!=0x3ff921fb) { 	/* 33+53 bit pi is good enough */
   4.221 +		    y[0] = z + pio2_1t;
   4.222 +		    y[1] = (z-y[0])+pio2_1t;
   4.223 +		} else {		/* near pi/2, use 33+33+53 bit pi */
   4.224 +		    z += pio2_2;
   4.225 +		    y[0] = z + pio2_2t;
   4.226 +		    y[1] = (z-y[0])+pio2_2t;
   4.227 +		}
   4.228 +		return -1;
   4.229 +	    }
   4.230 +	}
   4.231 +	if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
   4.232 +	    t  = fabs(x);
   4.233 +	    n  = (int32_t) (t*invpio2+half);
   4.234 +	    fn = (double)n;
   4.235 +	    r  = t-fn*pio2_1;
   4.236 +	    w  = fn*pio2_1t;	/* 1st round good to 85 bit */
   4.237 +	    if(n<32&&ix!=npio2_hw[n-1]) {
   4.238 +		y[0] = r-w;	/* quick check no cancellation */
   4.239 +	    } else {
   4.240 +	        u_int32_t high;
   4.241 +	        j  = ix>>20;
   4.242 +	        y[0] = r-w;
   4.243 +		GET_HIGH_WORD(high,y[0]);
   4.244 +	        i = j-((high>>20)&0x7ff);
   4.245 +	        if(i>16) {  /* 2nd iteration needed, good to 118 */
   4.246 +		    t  = r;
   4.247 +		    w  = fn*pio2_2;
   4.248 +		    r  = t-w;
   4.249 +		    w  = fn*pio2_2t-((t-r)-w);
   4.250 +		    y[0] = r-w;
   4.251 +		    GET_HIGH_WORD(high,y[0]);
   4.252 +		    i = j-((high>>20)&0x7ff);
   4.253 +		    if(i>49)  {	/* 3rd iteration need, 151 bits acc */
   4.254 +		    	t  = r;	/* will cover all possible cases */
   4.255 +		    	w  = fn*pio2_3;
   4.256 +		    	r  = t-w;
   4.257 +		    	w  = fn*pio2_3t-((t-r)-w);
   4.258 +		    	y[0] = r-w;
   4.259 +		    }
   4.260 +		}
   4.261 +	    }
   4.262 +	    y[1] = (r-y[0])-w;
   4.263 +	    if(hx<0) 	{y[0] = -y[0]; y[1] = -y[1]; return -n;}
   4.264 +	    else	 return n;
   4.265 +	}
   4.266      /*
   4.267       * all other (large) arguments
   4.268       */
   4.269 -    if (ix >= 0x7ff00000) {     /* x is inf or NaN */
   4.270 -        y[0] = y[1] = x - x;
   4.271 -        return 0;
   4.272 -    }
   4.273 +	if(ix>=0x7ff00000) {		/* x is inf or NaN */
   4.274 +	    y[0]=y[1]=x-x; return 0;
   4.275 +	}
   4.276      /* set z = scalbn(|x|,ilogb(x)-23) */
   4.277 -    GET_LOW_WORD(low, x);
   4.278 -    SET_LOW_WORD(z, low);
   4.279 -    e0 = (ix >> 20) - 1046;     /* e0 = ilogb(z)-23; */
   4.280 -    SET_HIGH_WORD(z, ix - ((int32_t) (e0 << 20)));
   4.281 -    for (i = 0; i < 2; i++) {
   4.282 -        tx[i] = (double) ((int32_t) (z));
   4.283 -        z = (z - tx[i]) * two24;
   4.284 -    }
   4.285 -    tx[2] = z;
   4.286 -    nx = 3;
   4.287 -
   4.288 -    /* If this assertion ever fires, here's the static analysis that warned about it:
   4.289 -        http://buildbot.libsdl.org/sdl-static-analysis/sdl-macosx-static-analysis/sdl-macosx-static-analysis-1101/report-8c6ccb.html#EndPath */
   4.290 -    SDL_assert((tx[0] != zero) || (tx[1] != zero) || (tx[2] != zero));
   4.291 -
   4.292 -    while (nx && tx[nx - 1] == zero)
   4.293 -        nx--;                   /* skip zero term */
   4.294 -    n = __kernel_rem_pio2(tx, y, e0, nx, 2, two_over_pi);
   4.295 -    if (hx < 0) {
   4.296 -        y[0] = -y[0];
   4.297 -        y[1] = -y[1];
   4.298 -        return -n;
   4.299 -    }
   4.300 -    return n;
   4.301 +	GET_LOW_WORD(low,x);
   4.302 +	SET_LOW_WORD(z,low);
   4.303 +	e0 	= (ix>>20)-1046;	/* e0 = ilogb(z)-23; */
   4.304 +	SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20)));
   4.305 +	for(i=0;i<2;i++) {
   4.306 +		tx[i] = (double)((int32_t)(z));
   4.307 +		z     = (z-tx[i])*two24;
   4.308 +	}
   4.309 +	tx[2] = z;
   4.310 +	nx = 3;
   4.311 +	while(tx[nx-1]==zero) nx--;	/* skip zero term */
   4.312 +	n  =  __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
   4.313 +	if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
   4.314 +	return n;
   4.315  }
     5.1 --- a/src/libm/e_sqrt.c	Sat Nov 04 15:34:14 2017 -0700
     5.2 +++ b/src/libm/e_sqrt.c	Sat Nov 04 15:53:19 2017 -0700
     5.3 @@ -1,4 +1,3 @@
     5.4 -/* @(#)e_sqrt.c 5.1 93/09/24 */
     5.5  /*
     5.6   * ====================================================
     5.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     5.8 @@ -10,11 +9,6 @@
     5.9   * ====================================================
    5.10   */
    5.11  
    5.12 -#if defined(LIBM_SCCS) && !defined(lint)
    5.13 -static const char rcsid[] =
    5.14 -    "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $";
    5.15 -#endif
    5.16 -
    5.17  /* __ieee754_sqrt(x)
    5.18   * Return correctly rounded sqrt.
    5.19   *           ------------------------------------------
    5.20 @@ -88,125 +82,124 @@
    5.21  #include "math_libm.h"
    5.22  #include "math_private.h"
    5.23  
    5.24 -#ifdef __STDC__
    5.25  static const double one = 1.0, tiny = 1.0e-300;
    5.26 -#else
    5.27 -static double one = 1.0, tiny = 1.0e-300;
    5.28 -#endif
    5.29  
    5.30 -#ifdef __STDC__
    5.31 -double attribute_hidden
    5.32 -__ieee754_sqrt(double x)
    5.33 -#else
    5.34 -double attribute_hidden
    5.35 -__ieee754_sqrt(x)
    5.36 -     double x;
    5.37 -#endif
    5.38 +double attribute_hidden __ieee754_sqrt(double x)
    5.39  {
    5.40 -    double z;
    5.41 -    int32_t sign = (int) 0x80000000;
    5.42 -    int32_t ix0, s0, q, m, t, i;
    5.43 -    u_int32_t r, t1, s1, ix1, q1;
    5.44 +	double z;
    5.45 +	int32_t sign = (int)0x80000000;
    5.46 +	int32_t ix0,s0,q,m,t,i;
    5.47 +	u_int32_t r,t1,s1,ix1,q1;
    5.48  
    5.49 -    EXTRACT_WORDS(ix0, ix1, x);
    5.50 +	EXTRACT_WORDS(ix0,ix1,x);
    5.51  
    5.52      /* take care of Inf and NaN */
    5.53 -    if ((ix0 & 0x7ff00000) == 0x7ff00000) {
    5.54 -        return x * x + x;       /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
    5.55 -                                   sqrt(-inf)=sNaN */
    5.56 -    }
    5.57 +	if((ix0&0x7ff00000)==0x7ff00000) {
    5.58 +	    return x*x+x;		/* sqrt(NaN)=NaN, sqrt(+inf)=+inf
    5.59 +					   sqrt(-inf)=sNaN */
    5.60 +	}
    5.61      /* take care of zero */
    5.62 -    if (ix0 <= 0) {
    5.63 -        if (((ix0 & (~sign)) | ix1) == 0)
    5.64 -            return x;           /* sqrt(+-0) = +-0 */
    5.65 -        else if (ix0 < 0)
    5.66 -            return (x - x) / (x - x);   /* sqrt(-ve) = sNaN */
    5.67 -    }
    5.68 +	if(ix0<=0) {
    5.69 +	    if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
    5.70 +	    else if(ix0<0)
    5.71 +		return (x-x)/(x-x);		/* sqrt(-ve) = sNaN */
    5.72 +	}
    5.73      /* normalize x */
    5.74 -    m = (ix0 >> 20);
    5.75 -    if (m == 0) {               /* subnormal x */
    5.76 -        while (ix0 == 0) {
    5.77 -            m -= 21;
    5.78 -            ix0 |= (ix1 >> 11);
    5.79 -            ix1 <<= 21;
    5.80 -        }
    5.81 -        for (i = 0; (ix0 & 0x00100000) == 0; i++)
    5.82 -            ix0 <<= 1;
    5.83 -        m -= i - 1;
    5.84 -        ix0 |= (ix1 >> (32 - i));
    5.85 -        ix1 <<= i;
    5.86 -    }
    5.87 -    m -= 1023;                  /* unbias exponent */
    5.88 -    ix0 = (ix0 & 0x000fffff) | 0x00100000;
    5.89 -    if (m & 1) {                /* odd m, double x to make it even */
    5.90 -        ix0 += ix0 + ((ix1 & sign) >> 31);
    5.91 -        ix1 += ix1;
    5.92 -    }
    5.93 -    m >>= 1;                    /* m = [m/2] */
    5.94 +	m = (ix0>>20);
    5.95 +	if(m==0) {				/* subnormal x */
    5.96 +	    while(ix0==0) {
    5.97 +		m -= 21;
    5.98 +		ix0 |= (ix1>>11); ix1 <<= 21;
    5.99 +	    }
   5.100 +	    for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
   5.101 +	    m -= i-1;
   5.102 +	    ix0 |= (ix1>>(32-i));
   5.103 +	    ix1 <<= i;
   5.104 +	}
   5.105 +	m -= 1023;	/* unbias exponent */
   5.106 +	ix0 = (ix0&0x000fffff)|0x00100000;
   5.107 +	if(m&1){	/* odd m, double x to make it even */
   5.108 +	    ix0 += ix0 + ((ix1&sign)>>31);
   5.109 +	    ix1 += ix1;
   5.110 +	}
   5.111 +	m >>= 1;	/* m = [m/2] */
   5.112  
   5.113      /* generate sqrt(x) bit by bit */
   5.114 -    ix0 += ix0 + ((ix1 & sign) >> 31);
   5.115 -    ix1 += ix1;
   5.116 -    q = q1 = s0 = s1 = 0;       /* [q,q1] = sqrt(x) */
   5.117 -    r = 0x00200000;             /* r = moving bit from right to left */
   5.118 +	ix0 += ix0 + ((ix1&sign)>>31);
   5.119 +	ix1 += ix1;
   5.120 +	q = q1 = s0 = s1 = 0;	/* [q,q1] = sqrt(x) */
   5.121 +	r = 0x00200000;		/* r = moving bit from right to left */
   5.122  
   5.123 -    while (r != 0) {
   5.124 -        t = s0 + r;
   5.125 -        if (t <= ix0) {
   5.126 -            s0 = t + r;
   5.127 -            ix0 -= t;
   5.128 -            q += r;
   5.129 -        }
   5.130 -        ix0 += ix0 + ((ix1 & sign) >> 31);
   5.131 -        ix1 += ix1;
   5.132 -        r >>= 1;
   5.133 -    }
   5.134 +	while(r!=0) {
   5.135 +	    t = s0+r;
   5.136 +	    if(t<=ix0) {
   5.137 +		s0   = t+r;
   5.138 +		ix0 -= t;
   5.139 +		q   += r;
   5.140 +	    }
   5.141 +	    ix0 += ix0 + ((ix1&sign)>>31);
   5.142 +	    ix1 += ix1;
   5.143 +	    r>>=1;
   5.144 +	}
   5.145  
   5.146 -    r = sign;
   5.147 -    while (r != 0) {
   5.148 -        t1 = s1 + r;
   5.149 -        t = s0;
   5.150 -        if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) {
   5.151 -            s1 = t1 + r;
   5.152 -            if (((t1 & sign) == sign) && (s1 & sign) == 0)
   5.153 -                s0 += 1;
   5.154 -            ix0 -= t;
   5.155 -            if (ix1 < t1)
   5.156 -                ix0 -= 1;
   5.157 -            ix1 -= t1;
   5.158 -            q1 += r;
   5.159 -        }
   5.160 -        ix0 += ix0 + ((ix1 & sign) >> 31);
   5.161 -        ix1 += ix1;
   5.162 -        r >>= 1;
   5.163 -    }
   5.164 +	r = sign;
   5.165 +	while(r!=0) {
   5.166 +	    t1 = s1+r;
   5.167 +	    t  = s0;
   5.168 +	    if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
   5.169 +		s1  = t1+r;
   5.170 +		if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
   5.171 +		ix0 -= t;
   5.172 +		if (ix1 < t1) ix0 -= 1;
   5.173 +		ix1 -= t1;
   5.174 +		q1  += r;
   5.175 +	    }
   5.176 +	    ix0 += ix0 + ((ix1&sign)>>31);
   5.177 +	    ix1 += ix1;
   5.178 +	    r>>=1;
   5.179 +	}
   5.180  
   5.181      /* use floating add to find out rounding direction */
   5.182 -    if ((ix0 | ix1) != 0) {
   5.183 -        z = one - tiny;         /* trigger inexact flag */
   5.184 -        if (z >= one) {
   5.185 -            z = one + tiny;
   5.186 -            if (q1 == (u_int32_t) 0xffffffff) {
   5.187 -                q1 = 0;
   5.188 -                q += 1;
   5.189 -            } else if (z > one) {
   5.190 -                if (q1 == (u_int32_t) 0xfffffffe)
   5.191 -                    q += 1;
   5.192 -                q1 += 2;
   5.193 -            } else
   5.194 -                q1 += (q1 & 1);
   5.195 -        }
   5.196 -    }
   5.197 -    ix0 = (q >> 1) + 0x3fe00000;
   5.198 -    ix1 = q1 >> 1;
   5.199 -    if ((q & 1) == 1)
   5.200 -        ix1 |= sign;
   5.201 -    ix0 += (m << 20);
   5.202 -    INSERT_WORDS(z, ix0, ix1);
   5.203 -    return z;
   5.204 +	if((ix0|ix1)!=0) {
   5.205 +	    z = one-tiny; /* trigger inexact flag */
   5.206 +	    if (z>=one) {
   5.207 +	        z = one+tiny;
   5.208 +	        if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
   5.209 +		else if (z>one) {
   5.210 +		    if (q1==(u_int32_t)0xfffffffe) q+=1;
   5.211 +		    q1+=2;
   5.212 +		} else
   5.213 +	            q1 += (q1&1);
   5.214 +	    }
   5.215 +	}
   5.216 +	ix0 = (q>>1)+0x3fe00000;
   5.217 +	ix1 =  q1>>1;
   5.218 +	if ((q&1)==1) ix1 |= sign;
   5.219 +	ix0 += (m <<20);
   5.220 +	INSERT_WORDS(z,ix0,ix1);
   5.221 +	return z;
   5.222  }
   5.223  
   5.224  /*
   5.225 + * wrapper sqrt(x)
   5.226 + */
   5.227 +#ifndef _IEEE_LIBM
   5.228 +double sqrt(double x)
   5.229 +{
   5.230 +	double z = __ieee754_sqrt(x);
   5.231 +	if (_LIB_VERSION == _IEEE_ || isnan(x))
   5.232 +		return z;
   5.233 +	if (x < 0.0)
   5.234 +		return __kernel_standard(x, x, 26); /* sqrt(negative) */
   5.235 +	return z;
   5.236 +}
   5.237 +#else
   5.238 +strong_alias(__ieee754_sqrt, sqrt)
   5.239 +#endif
   5.240 +libm_hidden_def(sqrt)
   5.241 +
   5.242 +
   5.243 +/*
   5.244  Other methods  (use floating-point arithmetic)
   5.245  -------------
   5.246  (This is a copy of a drafted paper by Prof W. Kahan
     6.1 --- a/src/libm/k_cos.c	Sat Nov 04 15:34:14 2017 -0700
     6.2 +++ b/src/libm/k_cos.c	Sat Nov 04 15:53:19 2017 -0700
     6.3 @@ -1,4 +1,3 @@
     6.4 -/* @(#)k_cos.c 5.1 93/09/24 */
     6.5  /*
     6.6   * ====================================================
     6.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     6.8 @@ -10,11 +9,6 @@
     6.9   * ====================================================
    6.10   */
    6.11  
    6.12 -#if defined(LIBM_SCCS) && !defined(lint)
    6.13 -static const char rcsid[] =
    6.14 -    "$NetBSD: k_cos.c,v 1.8 1995/05/10 20:46:22 jtc Exp $";
    6.15 -#endif
    6.16 -
    6.17  /*
    6.18   * __kernel_cos( x,  y )
    6.19   * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
    6.20 @@ -53,48 +47,36 @@
    6.21  #include "math_libm.h"
    6.22  #include "math_private.h"
    6.23  
    6.24 -#ifdef __STDC__
    6.25  static const double
    6.26 -#else
    6.27 -static double
    6.28 -#endif
    6.29 -  one = 1.00000000000000000000e+00,     /* 0x3FF00000, 0x00000000 */
    6.30 -    C1 = 4.16666666666666019037e-02,    /* 0x3FA55555, 0x5555554C */
    6.31 -    C2 = -1.38888888888741095749e-03,   /* 0xBF56C16C, 0x16C15177 */
    6.32 -    C3 = 2.48015872894767294178e-05,    /* 0x3EFA01A0, 0x19CB1590 */
    6.33 -    C4 = -2.75573143513906633035e-07,   /* 0xBE927E4F, 0x809C52AD */
    6.34 -    C5 = 2.08757232129817482790e-09,    /* 0x3E21EE9E, 0xBDB4B1C4 */
    6.35 -    C6 = -1.13596475577881948265e-11;   /* 0xBDA8FAE9, 0xBE8838D4 */
    6.36 +one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
    6.37 +C1  =  4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
    6.38 +C2  = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
    6.39 +C3  =  2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
    6.40 +C4  = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
    6.41 +C5  =  2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
    6.42 +C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
    6.43  
    6.44 -#ifdef __STDC__
    6.45 -double attribute_hidden
    6.46 -__kernel_cos(double x, double y)
    6.47 -#else
    6.48 -double attribute_hidden
    6.49 -__kernel_cos(x, y)
    6.50 -     double x, y;
    6.51 -#endif
    6.52 +double attribute_hidden __kernel_cos(double x, double y)
    6.53  {
    6.54 -    double a, hz, z, r, qx;
    6.55 -    int32_t ix;
    6.56 -    GET_HIGH_WORD(ix, x);
    6.57 -    ix &= 0x7fffffff;           /* ix = |x|'s high word */
    6.58 -    if (ix < 0x3e400000) {      /* if x < 2**27 */
    6.59 -        if (((int) x) == 0)
    6.60 -            return one;         /* generate inexact */
    6.61 -    }
    6.62 -    z = x * x;
    6.63 -    r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6)))));
    6.64 -    if (ix < 0x3FD33333)        /* if |x| < 0.3 */
    6.65 -        return one - (0.5 * z - (z * r - x * y));
    6.66 -    else {
    6.67 -        if (ix > 0x3fe90000) {  /* x > 0.78125 */
    6.68 -            qx = 0.28125;
    6.69 -        } else {
    6.70 -            INSERT_WORDS(qx, ix - 0x00200000, 0);       /* x/4 */
    6.71 -        }
    6.72 -        hz = 0.5 * z - qx;
    6.73 -        a = one - qx;
    6.74 -        return a - (hz - (z * r - x * y));
    6.75 -    }
    6.76 +	double a,hz,z,r,qx;
    6.77 +	int32_t ix;
    6.78 +	GET_HIGH_WORD(ix,x);
    6.79 +	ix &= 0x7fffffff;			/* ix = |x|'s high word*/
    6.80 +	if(ix<0x3e400000) {			/* if x < 2**27 */
    6.81 +	    if(((int)x)==0) return one;		/* generate inexact */
    6.82 +	}
    6.83 +	z  = x*x;
    6.84 +	r  = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
    6.85 +	if(ix < 0x3FD33333) 			/* if |x| < 0.3 */
    6.86 +	    return one - (0.5*z - (z*r - x*y));
    6.87 +	else {
    6.88 +	    if(ix > 0x3fe90000) {		/* x > 0.78125 */
    6.89 +		qx = 0.28125;
    6.90 +	    } else {
    6.91 +	        INSERT_WORDS(qx,ix-0x00200000,0);	/* x/4 */
    6.92 +	    }
    6.93 +	    hz = 0.5*z-qx;
    6.94 +	    a  = one-qx;
    6.95 +	    return a - (hz - (z*r-x*y));
    6.96 +	}
    6.97  }
     7.1 --- a/src/libm/k_rem_pio2.c	Sat Nov 04 15:34:14 2017 -0700
     7.2 +++ b/src/libm/k_rem_pio2.c	Sat Nov 04 15:53:19 2017 -0700
     7.3 @@ -1,4 +1,3 @@
     7.4 -/* @(#)k_rem_pio2.c 5.1 93/09/24 */
     7.5  /*
     7.6   * ====================================================
     7.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     7.8 @@ -10,11 +9,6 @@
     7.9   * ====================================================
    7.10   */
    7.11  
    7.12 -#if defined(LIBM_SCCS) && !defined(lint)
    7.13 -static const char rcsid[] =
    7.14 -    "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $";
    7.15 -#endif
    7.16 -
    7.17  /*
    7.18   * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
    7.19   * double x[],y[]; int e0,nx,prec; int ipio2[];
    7.20 @@ -134,235 +128,171 @@
    7.21  #include "math_libm.h"
    7.22  #include "math_private.h"
    7.23  
    7.24 -#include "SDL_assert.h"
    7.25 +static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
    7.26  
    7.27 -libm_hidden_proto(scalbn)
    7.28 -    libm_hidden_proto(floor)
    7.29 -#ifdef __STDC__
    7.30 -     static const int init_jk[] = { 2, 3, 4, 6 };       /* initial value for jk */
    7.31 -#else
    7.32 -     static int init_jk[] = { 2, 3, 4, 6 };
    7.33 -#endif
    7.34 -
    7.35 -#ifdef __STDC__
    7.36  static const double PIo2[] = {
    7.37 -#else
    7.38 -static double PIo2[] = {
    7.39 -#endif
    7.40 -    1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
    7.41 -    7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
    7.42 -    5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
    7.43 -    3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
    7.44 -    1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
    7.45 -    1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
    7.46 -    2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
    7.47 -    2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
    7.48 +  1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
    7.49 +  7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
    7.50 +  5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
    7.51 +  3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
    7.52 +  1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
    7.53 +  1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
    7.54 +  2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
    7.55 +  2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
    7.56  };
    7.57  
    7.58 -#ifdef __STDC__
    7.59  static const double
    7.60 -#else
    7.61 -static double
    7.62 -#endif
    7.63 -  zero = 0.0, one = 1.0, two24 = 1.67772160000000000000e+07,    /* 0x41700000, 0x00000000 */
    7.64 -    twon24 = 5.96046447753906250000e-08;        /* 0x3E700000, 0x00000000 */
    7.65 +zero   = 0.0,
    7.66 +one    = 1.0,
    7.67 +two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
    7.68 +twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
    7.69  
    7.70 -#ifdef __STDC__
    7.71 -int attribute_hidden
    7.72 -__kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec,
    7.73 -                  const int32_t * ipio2)
    7.74 -#else
    7.75 -int attribute_hidden
    7.76 -__kernel_rem_pio2(x, y, e0, nx, prec, ipio2)
    7.77 -     double x[], y[];
    7.78 -     int e0, nx, prec;
    7.79 -     int32_t ipio2[];
    7.80 -#endif
    7.81 +int attribute_hidden __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int32_t *ipio2)
    7.82  {
    7.83 -    int32_t jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
    7.84 -    double z, fw, f[20], fq[20], q[20];
    7.85 +	int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
    7.86 +	double z,fw,f[20],fq[20],q[20];
    7.87  
    7.88 -    /* initialize jk */
    7.89 -    SDL_assert((prec >= 0) && (prec < SDL_arraysize(init_jk)));
    7.90 -    jk = init_jk[prec];
    7.91 -    SDL_assert((jk >= 2) && (jk <= 6));
    7.92 -    jp = jk;
    7.93 +    /* initialize jk*/
    7.94 +	jk = init_jk[prec];
    7.95 +	jp = jk;
    7.96  
    7.97      /* determine jx,jv,q0, note that 3>q0 */
    7.98 -    SDL_assert(nx > 0);
    7.99 -    jx = nx - 1;
   7.100 -    jv = (e0 - 3) / 24;
   7.101 -    if (jv < 0)
   7.102 -        jv = 0;
   7.103 -    q0 = e0 - 24 * (jv + 1);
   7.104 +	jx =  nx-1;
   7.105 +	jv = (e0-3)/24; if(jv<0) jv=0;
   7.106 +	q0 =  e0-24*(jv+1);
   7.107  
   7.108      /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
   7.109 -    j = jv - jx;
   7.110 -    m = jx + jk;
   7.111 -    for (i = 0; i <= m; i++, j++)
   7.112 -        f[i] = (j < 0) ? zero : (double) ipio2[j];
   7.113 +	j = jv-jx; m = jx+jk;
   7.114 +	for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
   7.115  
   7.116      /* compute q[0],q[1],...q[jk] */
   7.117 -    for (i = 0; i <= jk; i++) {
   7.118 -        for (j = 0, fw = 0.0; j <= jx; j++) {
   7.119 -            const int32_t idx = jx + i - j;
   7.120 -            SDL_assert(idx >= 0);
   7.121 -            SDL_assert(idx < 20);
   7.122 -            SDL_assert(idx <= m);
   7.123 -            fw += x[j] * f[idx];
   7.124 -        }
   7.125 -        q[i] = fw;
   7.126 -    }
   7.127 +	for (i=0;i<=jk;i++) {
   7.128 +	    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
   7.129 +	}
   7.130  
   7.131 -    jz = jk;
   7.132 -  recompute:
   7.133 +	jz = jk;
   7.134 +recompute:
   7.135      /* distill q[] into iq[] reversingly */
   7.136 -    for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
   7.137 -        fw = (double) ((int32_t) (twon24 * z));
   7.138 -        iq[i] = (int32_t) (z - two24 * fw);
   7.139 -        z = q[j - 1] + fw;
   7.140 -    }
   7.141 +	for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
   7.142 +	    fw    =  (double)((int32_t)(twon24* z));
   7.143 +	    iq[i] =  (int32_t)(z-two24*fw);
   7.144 +	    z     =  q[j-1]+fw;
   7.145 +	}
   7.146  
   7.147      /* compute n */
   7.148 -    z = scalbn(z, q0);          /* actual value of z */
   7.149 -    z -= 8.0 * floor(z * 0.125);        /* trim off integer >= 8 */
   7.150 -    n = (int32_t) z;
   7.151 -    z -= (double) n;
   7.152 -    ih = 0;
   7.153 -    if (q0 > 0) {               /* need iq[jz-1] to determine n */
   7.154 -        i = (iq[jz - 1] >> (24 - q0));
   7.155 -        n += i;
   7.156 -        iq[jz - 1] -= i << (24 - q0);
   7.157 -        ih = iq[jz - 1] >> (23 - q0);
   7.158 -    } else if (q0 == 0)
   7.159 -        ih = iq[jz - 1] >> 23;
   7.160 -    else if (z >= 0.5)
   7.161 -        ih = 2;
   7.162 +	z  = scalbn(z,q0);		/* actual value of z */
   7.163 +	z -= 8.0*floor(z*0.125);		/* trim off integer >= 8 */
   7.164 +	n  = (int32_t) z;
   7.165 +	z -= (double)n;
   7.166 +	ih = 0;
   7.167 +	if(q0>0) {	/* need iq[jz-1] to determine n */
   7.168 +	    i  = (iq[jz-1]>>(24-q0)); n += i;
   7.169 +	    iq[jz-1] -= i<<(24-q0);
   7.170 +	    ih = iq[jz-1]>>(23-q0);
   7.171 +	}
   7.172 +	else if(q0==0) ih = iq[jz-1]>>23;
   7.173 +	else if(z>=0.5) ih=2;
   7.174  
   7.175 -    if (ih > 0) {               /* q > 0.5 */
   7.176 -        n += 1;
   7.177 -        carry = 0;
   7.178 -        for (i = 0; i < jz; i++) {      /* compute 1-q */
   7.179 -            j = iq[i];
   7.180 -            if (carry == 0) {
   7.181 -                if (j != 0) {
   7.182 -                    carry = 1;
   7.183 -                    iq[i] = 0x1000000 - j;
   7.184 -                }
   7.185 -            } else
   7.186 -                iq[i] = 0xffffff - j;
   7.187 -        }
   7.188 -        if (q0 > 0) {           /* rare case: chance is 1 in 12 */
   7.189 -            switch (q0) {
   7.190 -            case 1:
   7.191 -                iq[jz - 1] &= 0x7fffff;
   7.192 -                break;
   7.193 -            case 2:
   7.194 -                iq[jz - 1] &= 0x3fffff;
   7.195 -                break;
   7.196 -            }
   7.197 -        }
   7.198 -        if (ih == 2) {
   7.199 -            z = one - z;
   7.200 -            if (carry != 0)
   7.201 -                z -= scalbn(one, q0);
   7.202 -        }
   7.203 -    }
   7.204 +	if(ih>0) {	/* q > 0.5 */
   7.205 +	    n += 1; carry = 0;
   7.206 +	    for(i=0;i<jz ;i++) {	/* compute 1-q */
   7.207 +		j = iq[i];
   7.208 +		if(carry==0) {
   7.209 +		    if(j!=0) {
   7.210 +			carry = 1; iq[i] = 0x1000000- j;
   7.211 +		    }
   7.212 +		} else  iq[i] = 0xffffff - j;
   7.213 +	    }
   7.214 +	    if(q0>0) {		/* rare case: chance is 1 in 12 */
   7.215 +	        switch(q0) {
   7.216 +	        case 1:
   7.217 +	    	   iq[jz-1] &= 0x7fffff; break;
   7.218 +	    	case 2:
   7.219 +	    	   iq[jz-1] &= 0x3fffff; break;
   7.220 +	        }
   7.221 +	    }
   7.222 +	    if(ih==2) {
   7.223 +		z = one - z;
   7.224 +		if(carry!=0) z -= scalbn(one,q0);
   7.225 +	    }
   7.226 +	}
   7.227  
   7.228      /* check if recomputation is needed */
   7.229 -    if (z == zero) {
   7.230 -        j = 0;
   7.231 -        for (i = jz - 1; i >= jk; i--)
   7.232 -            j |= iq[i];
   7.233 -        if (j == 0) {           /* need recomputation */
   7.234 -            for (k = 1; iq[jk - k] == 0; k++);  /* k = no. of terms needed */
   7.235 +	if(z==zero) {
   7.236 +	    j = 0;
   7.237 +	    for (i=jz-1;i>=jk;i--) j |= iq[i];
   7.238 +	    if(j==0) { /* need recomputation */
   7.239 +		for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */
   7.240  
   7.241 -            for (i = jz + 1; i <= jz + k; i++) {        /* add q[jz+1] to q[jz+k] */
   7.242 -                f[jx + i] = (double) ipio2[jv + i];
   7.243 -                for (j = 0, fw = 0.0; j <= jx; j++)
   7.244 -                    fw += x[j] * f[jx + i - j];
   7.245 -                q[i] = fw;
   7.246 -            }
   7.247 -            jz += k;
   7.248 -            goto recompute;
   7.249 -        }
   7.250 -    }
   7.251 +		for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */
   7.252 +		    f[jx+i] = (double) ipio2[jv+i];
   7.253 +		    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
   7.254 +		    q[i] = fw;
   7.255 +		}
   7.256 +		jz += k;
   7.257 +		goto recompute;
   7.258 +	    }
   7.259 +	}
   7.260  
   7.261      /* chop off zero terms */
   7.262 -    if (z == 0.0) {
   7.263 -        jz -= 1;
   7.264 -        q0 -= 24;
   7.265 -        while (iq[jz] == 0) {
   7.266 -            jz--;
   7.267 -            q0 -= 24;
   7.268 -        }
   7.269 -    } else {                    /* break z into 24-bit if necessary */
   7.270 -        z = scalbn(z, -q0);
   7.271 -        if (z >= two24) {
   7.272 -            fw = (double) ((int32_t) (twon24 * z));
   7.273 -            iq[jz] = (int32_t) (z - two24 * fw);
   7.274 -            jz += 1;
   7.275 -            q0 += 24;
   7.276 -            iq[jz] = (int32_t) fw;
   7.277 -        } else
   7.278 -            iq[jz] = (int32_t) z;
   7.279 -    }
   7.280 +	if(z==0.0) {
   7.281 +	    jz -= 1; q0 -= 24;
   7.282 +	    while(iq[jz]==0) { jz--; q0-=24;}
   7.283 +	} else { /* break z into 24-bit if necessary */
   7.284 +	    z = scalbn(z,-q0);
   7.285 +	    if(z>=two24) {
   7.286 +		fw = (double)((int32_t)(twon24*z));
   7.287 +		iq[jz] = (int32_t)(z-two24*fw);
   7.288 +		jz += 1; q0 += 24;
   7.289 +		iq[jz] = (int32_t) fw;
   7.290 +	    } else iq[jz] = (int32_t) z ;
   7.291 +	}
   7.292  
   7.293      /* convert integer "bit" chunk to floating-point value */
   7.294 -    fw = scalbn(one, q0);
   7.295 -    for (i = jz; i >= 0; i--) {
   7.296 -        q[i] = fw * (double) iq[i];
   7.297 -        fw *= twon24;
   7.298 -    }
   7.299 +	fw = scalbn(one,q0);
   7.300 +	for(i=jz;i>=0;i--) {
   7.301 +	    q[i] = fw*(double)iq[i]; fw*=twon24;
   7.302 +	}
   7.303  
   7.304      /* compute PIo2[0,...,jp]*q[jz,...,0] */
   7.305 -    for (i = jz; i >= 0; i--) {
   7.306 -        for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
   7.307 -            fw += PIo2[k] * q[i + k];
   7.308 -        fq[jz - i] = fw;
   7.309 -    }
   7.310 +	for(i=jz;i>=0;i--) {
   7.311 +	    for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
   7.312 +	    fq[jz-i] = fw;
   7.313 +	}
   7.314  
   7.315      /* compress fq[] into y[] */
   7.316 -    switch (prec) {
   7.317 -    case 0:
   7.318 -        fw = 0.0;
   7.319 -        for (i = jz; i >= 0; i--)
   7.320 -            fw += fq[i];
   7.321 -        y[0] = (ih == 0) ? fw : -fw;
   7.322 -        break;
   7.323 -    case 1:
   7.324 -    case 2:
   7.325 -        fw = 0.0;
   7.326 -        for (i = jz; i >= 0; i--)
   7.327 -            fw += fq[i];
   7.328 -        y[0] = (ih == 0) ? fw : -fw;
   7.329 -        fw = fq[0] - fw;
   7.330 -        for (i = 1; i <= jz; i++)
   7.331 -            fw += fq[i];
   7.332 -        y[1] = (ih == 0) ? fw : -fw;
   7.333 -        break;
   7.334 -    case 3:                    /* painful */
   7.335 -        for (i = jz; i > 0; i--) {
   7.336 -            fw = fq[i - 1] + fq[i];
   7.337 -            fq[i] += fq[i - 1] - fw;
   7.338 -            fq[i - 1] = fw;
   7.339 -        }
   7.340 -        for (i = jz; i > 1; i--) {
   7.341 -            fw = fq[i - 1] + fq[i];
   7.342 -            fq[i] += fq[i - 1] - fw;
   7.343 -            fq[i - 1] = fw;
   7.344 -        }
   7.345 -        for (fw = 0.0, i = jz; i >= 2; i--)
   7.346 -            fw += fq[i];
   7.347 -        if (ih == 0) {
   7.348 -            y[0] = fq[0];
   7.349 -            y[1] = fq[1];
   7.350 -            y[2] = fw;
   7.351 -        } else {
   7.352 -            y[0] = -fq[0];
   7.353 -            y[1] = -fq[1];
   7.354 -            y[2] = -fw;
   7.355 -        }
   7.356 -    }
   7.357 -    return n & 7;
   7.358 +	switch(prec) {
   7.359 +	    case 0:
   7.360 +		fw = 0.0;
   7.361 +		for (i=jz;i>=0;i--) fw += fq[i];
   7.362 +		y[0] = (ih==0)? fw: -fw;
   7.363 +		break;
   7.364 +	    case 1:
   7.365 +	    case 2:
   7.366 +		fw = 0.0;
   7.367 +		for (i=jz;i>=0;i--) fw += fq[i];
   7.368 +		y[0] = (ih==0)? fw: -fw;
   7.369 +		fw = fq[0]-fw;
   7.370 +		for (i=1;i<=jz;i++) fw += fq[i];
   7.371 +		y[1] = (ih==0)? fw: -fw;
   7.372 +		break;
   7.373 +	    case 3:	/* painful */
   7.374 +		for (i=jz;i>0;i--) {
   7.375 +		    fw      = fq[i-1]+fq[i];
   7.376 +		    fq[i]  += fq[i-1]-fw;
   7.377 +		    fq[i-1] = fw;
   7.378 +		}
   7.379 +		for (i=jz;i>1;i--) {
   7.380 +		    fw      = fq[i-1]+fq[i];
   7.381 +		    fq[i]  += fq[i-1]-fw;
   7.382 +		    fq[i-1] = fw;
   7.383 +		}
   7.384 +		for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
   7.385 +		if(ih==0) {
   7.386 +		    y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
   7.387 +		} else {
   7.388 +		    y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
   7.389 +		}
   7.390 +	}
   7.391 +	return n&7;
   7.392  }
     8.1 --- a/src/libm/k_sin.c	Sat Nov 04 15:34:14 2017 -0700
     8.2 +++ b/src/libm/k_sin.c	Sat Nov 04 15:53:19 2017 -0700
     8.3 @@ -1,4 +1,3 @@
     8.4 -/* @(#)k_sin.c 5.1 93/09/24 */
     8.5  /*
     8.6   * ====================================================
     8.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     8.8 @@ -10,11 +9,6 @@
     8.9   * ====================================================
    8.10   */
    8.11  
    8.12 -#if defined(LIBM_SCCS) && !defined(lint)
    8.13 -static const char rcsid[] =
    8.14 -    "$NetBSD: k_sin.c,v 1.8 1995/05/10 20:46:31 jtc Exp $";
    8.15 -#endif
    8.16 -
    8.17  /* __kernel_sin( x, y, iy)
    8.18   * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
    8.19   * Input x is assumed to be bounded by ~pi/4 in magnitude.
    8.20 @@ -46,42 +40,26 @@
    8.21  #include "math_libm.h"
    8.22  #include "math_private.h"
    8.23  
    8.24 -#ifdef __STDC__
    8.25  static const double
    8.26 -#else
    8.27 -static double
    8.28 -#endif
    8.29 -  half = 5.00000000000000000000e-01,    /* 0x3FE00000, 0x00000000 */
    8.30 -    S1 = -1.66666666666666324348e-01,   /* 0xBFC55555, 0x55555549 */
    8.31 -    S2 = 8.33333333332248946124e-03,    /* 0x3F811111, 0x1110F8A6 */
    8.32 -    S3 = -1.98412698298579493134e-04,   /* 0xBF2A01A0, 0x19C161D5 */
    8.33 -    S4 = 2.75573137070700676789e-06,    /* 0x3EC71DE3, 0x57B1FE7D */
    8.34 -    S5 = -2.50507602534068634195e-08,   /* 0xBE5AE5E6, 0x8A2B9CEB */
    8.35 -    S6 = 1.58969099521155010221e-10;    /* 0x3DE5D93A, 0x5ACFD57C */
    8.36 +half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
    8.37 +S1  = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
    8.38 +S2  =  8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
    8.39 +S3  = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
    8.40 +S4  =  2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
    8.41 +S5  = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
    8.42 +S6  =  1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
    8.43  
    8.44 -#ifdef __STDC__
    8.45 -double attribute_hidden
    8.46 -__kernel_sin(double x, double y, int iy)
    8.47 -#else
    8.48 -double attribute_hidden
    8.49 -__kernel_sin(x, y, iy)
    8.50 -     double x, y;
    8.51 -     int iy;                    /* iy=0 if y is zero */
    8.52 -#endif
    8.53 +double attribute_hidden __kernel_sin(double x, double y, int iy)
    8.54  {
    8.55 -    double z, r, v;
    8.56 -    int32_t ix;
    8.57 -    GET_HIGH_WORD(ix, x);
    8.58 -    ix &= 0x7fffffff;           /* high word of x */
    8.59 -    if (ix < 0x3e400000) {      /* |x| < 2**-27 */
    8.60 -        if ((int) x == 0)
    8.61 -            return x;
    8.62 -    }                           /* generate inexact */
    8.63 -    z = x * x;
    8.64 -    v = z * x;
    8.65 -    r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6)));
    8.66 -    if (iy == 0)
    8.67 -        return x + v * (S1 + z * r);
    8.68 -    else
    8.69 -        return x - ((z * (half * y - v * r) - y) - v * S1);
    8.70 +	double z,r,v;
    8.71 +	int32_t ix;
    8.72 +	GET_HIGH_WORD(ix,x);
    8.73 +	ix &= 0x7fffffff;			/* high word of x */
    8.74 +	if(ix<0x3e400000)			/* |x| < 2**-27 */
    8.75 +	   {if((int)x==0) return x;}		/* generate inexact */
    8.76 +	z	=  x*x;
    8.77 +	v	=  z*x;
    8.78 +	r	=  S2+z*(S3+z*(S4+z*(S5+z*S6)));
    8.79 +	if(iy==0) return x+v*(S1+z*r);
    8.80 +	else      return x-((z*(half*y-v*r)-y)-v*S1);
    8.81  }
     9.1 --- a/src/libm/k_tan.c	Sat Nov 04 15:34:14 2017 -0700
     9.2 +++ b/src/libm/k_tan.c	Sat Nov 04 15:53:19 2017 -0700
     9.3 @@ -66,7 +66,7 @@
     9.4    2.59073051863633712884e-05, /* 0x3EFB2A70, 0x74BF7AD4 */
     9.5  };
     9.6  
     9.7 -double __kernel_tan(double x, double y, int iy)
     9.8 +double attribute_hidden __kernel_tan(double x, double y, int iy)
     9.9  {
    9.10  	double z,r,v,w,s;
    9.11  	int32_t ix,hx;
    10.1 --- a/src/libm/math_libm.h	Sat Nov 04 15:34:14 2017 -0700
    10.2 +++ b/src/libm/math_libm.h	Sat Nov 04 15:53:19 2017 -0700
    10.3 @@ -28,6 +28,7 @@
    10.4  double SDL_uclibc_cos(double x);         
    10.5  double SDL_uclibc_fabs(double x);        
    10.6  double SDL_uclibc_floor(double x);
    10.7 +double SDL_uclibc_fmod(double x, double y);
    10.8  double SDL_uclibc_log(double x);
    10.9  double SDL_uclibc_pow(double x, double y);    
   10.10  double SDL_uclibc_scalbn(double x, int n);
    11.1 --- a/src/libm/math_private.h	Sat Nov 04 15:34:14 2017 -0700
    11.2 +++ b/src/libm/math_private.h	Sat Nov 04 15:53:19 2017 -0700
    11.3 @@ -21,9 +21,11 @@
    11.4  #include "SDL_endian.h"
    11.5  /* #include <sys/types.h> */
    11.6  
    11.7 +#define _IEEE_LIBM
    11.8  #define attribute_hidden
    11.9  #define libm_hidden_proto(x)
   11.10  #define libm_hidden_def(x)
   11.11 +#define strong_alias(x, y)
   11.12  
   11.13  #ifndef __HAIKU__ /* already defined in a system header. */
   11.14  typedef unsigned int u_int32_t;
   11.15 @@ -35,6 +37,7 @@
   11.16  #define cos             SDL_uclibc_cos
   11.17  #define fabs            SDL_uclibc_fabs
   11.18  #define floor           SDL_uclibc_floor
   11.19 +#define __ieee754_fmod  SDL_uclibc_fmod
   11.20  #define __ieee754_log   SDL_uclibc_log
   11.21  #define __ieee754_pow   SDL_uclibc_pow
   11.22  #define scalbn          SDL_uclibc_scalbn
    12.1 --- a/src/libm/s_atan.c	Sat Nov 04 15:34:14 2017 -0700
    12.2 +++ b/src/libm/s_atan.c	Sat Nov 04 15:53:19 2017 -0700
    12.3 @@ -112,4 +112,3 @@
    12.4  	}
    12.5  }
    12.6  libm_hidden_def(atan)
    12.7 -
    13.1 --- a/src/libm/s_copysign.c	Sat Nov 04 15:34:14 2017 -0700
    13.2 +++ b/src/libm/s_copysign.c	Sat Nov 04 15:53:19 2017 -0700
    13.3 @@ -1,4 +1,3 @@
    13.4 -/* @(#)s_copysign.c 5.1 93/09/24 */
    13.5  /*
    13.6   * ====================================================
    13.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    13.8 @@ -10,11 +9,6 @@
    13.9   * ====================================================
   13.10   */
   13.11  
   13.12 -#if defined(LIBM_SCCS) && !defined(lint)
   13.13 -static const char rcsid[] =
   13.14 -    "$NetBSD: s_copysign.c,v 1.8 1995/05/10 20:46:57 jtc Exp $";
   13.15 -#endif
   13.16 -
   13.17  /*
   13.18   * copysign(double x, double y)
   13.19   * copysign(x,y) returns a value with the magnitude of x and
   13.20 @@ -24,19 +18,12 @@
   13.21  #include "math_libm.h"
   13.22  #include "math_private.h"
   13.23  
   13.24 -libm_hidden_proto(copysign)
   13.25 -#ifdef __STDC__
   13.26 -     double copysign(double x, double y)
   13.27 -#else
   13.28 -     double copysign(x, y)
   13.29 -     double x, y;
   13.30 -#endif
   13.31 +double copysign(double x, double y)
   13.32  {
   13.33 -    u_int32_t hx, hy;
   13.34 -    GET_HIGH_WORD(hx, x);
   13.35 -    GET_HIGH_WORD(hy, y);
   13.36 -    SET_HIGH_WORD(x, (hx & 0x7fffffff) | (hy & 0x80000000));
   13.37 -    return x;
   13.38 +	u_int32_t hx,hy;
   13.39 +	GET_HIGH_WORD(hx,x);
   13.40 +	GET_HIGH_WORD(hy,y);
   13.41 +	SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000));
   13.42 +        return x;
   13.43  }
   13.44 -
   13.45  libm_hidden_def(copysign)
    14.1 --- a/src/libm/s_cos.c	Sat Nov 04 15:34:14 2017 -0700
    14.2 +++ b/src/libm/s_cos.c	Sat Nov 04 15:53:19 2017 -0700
    14.3 @@ -1,4 +1,3 @@
    14.4 -/* @(#)s_cos.c 5.1 93/09/24 */
    14.5  /*
    14.6   * ====================================================
    14.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    14.8 @@ -10,11 +9,6 @@
    14.9   * ====================================================
   14.10   */
   14.11  
   14.12 -#if defined(LIBM_SCCS) && !defined(lint)
   14.13 -static const char rcsid[] =
   14.14 -    "$NetBSD: s_cos.c,v 1.7 1995/05/10 20:47:02 jtc Exp $";
   14.15 -#endif
   14.16 -
   14.17  /* cos(x)
   14.18   * Return cosine function of x.
   14.19   *
   14.20 @@ -49,43 +43,31 @@
   14.21  #include "math_libm.h"
   14.22  #include "math_private.h"
   14.23  
   14.24 -libm_hidden_proto(cos)
   14.25 -#ifdef __STDC__
   14.26 -     double cos(double x)
   14.27 -#else
   14.28 -     double cos(x)
   14.29 -     double x;
   14.30 -#endif
   14.31 +double cos(double x)
   14.32  {
   14.33 -    double y[2], z = 0.0;
   14.34 -    int32_t n, ix;
   14.35 +	double y[2],z=0.0;
   14.36 +	int32_t n, ix;
   14.37  
   14.38      /* High word of x. */
   14.39 -    GET_HIGH_WORD(ix, x);
   14.40 +	GET_HIGH_WORD(ix,x);
   14.41  
   14.42      /* |x| ~< pi/4 */
   14.43 -    ix &= 0x7fffffff;
   14.44 -    if (ix <= 0x3fe921fb)
   14.45 -        return __kernel_cos(x, z);
   14.46 +	ix &= 0x7fffffff;
   14.47 +	if(ix <= 0x3fe921fb) return __kernel_cos(x,z);
   14.48  
   14.49      /* cos(Inf or NaN) is NaN */
   14.50 -    else if (ix >= 0x7ff00000)
   14.51 -        return x - x;
   14.52 +	else if (ix>=0x7ff00000) return x-x;
   14.53  
   14.54      /* argument reduction needed */
   14.55 -    else {
   14.56 -        n = __ieee754_rem_pio2(x, y);
   14.57 -        switch (n & 3) {
   14.58 -        case 0:
   14.59 -            return __kernel_cos(y[0], y[1]);
   14.60 -        case 1:
   14.61 -            return -__kernel_sin(y[0], y[1], 1);
   14.62 -        case 2:
   14.63 -            return -__kernel_cos(y[0], y[1]);
   14.64 -        default:
   14.65 -            return __kernel_sin(y[0], y[1], 1);
   14.66 -        }
   14.67 -    }
   14.68 +	else {
   14.69 +	    n = __ieee754_rem_pio2(x,y);
   14.70 +	    switch(n&3) {
   14.71 +		case 0: return  __kernel_cos(y[0],y[1]);
   14.72 +		case 1: return -__kernel_sin(y[0],y[1],1);
   14.73 +		case 2: return -__kernel_cos(y[0],y[1]);
   14.74 +		default:
   14.75 +		        return  __kernel_sin(y[0],y[1],1);
   14.76 +	    }
   14.77 +	}
   14.78  }
   14.79 -
   14.80  libm_hidden_def(cos)
    15.1 --- a/src/libm/s_fabs.c	Sat Nov 04 15:34:14 2017 -0700
    15.2 +++ b/src/libm/s_fabs.c	Sat Nov 04 15:53:19 2017 -0700
    15.3 @@ -1,4 +1,3 @@
    15.4 -/* @(#)s_fabs.c 5.1 93/09/24 */
    15.5  /*
    15.6   * ====================================================
    15.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    15.8 @@ -10,30 +9,21 @@
    15.9   * ====================================================
   15.10   */
   15.11  
   15.12 -#if defined(LIBM_SCCS) && !defined(lint)
   15.13 -static const char rcsid[] =
   15.14 -    "$NetBSD: s_fabs.c,v 1.7 1995/05/10 20:47:13 jtc Exp $";
   15.15 -#endif
   15.16 -
   15.17  /*
   15.18   * fabs(x) returns the absolute value of x.
   15.19   */
   15.20  
   15.21 +/*#include <features.h>*/
   15.22 +/* Prevent math.h from defining a colliding inline */
   15.23 +#undef __USE_EXTERN_INLINES
   15.24  #include "math_libm.h"
   15.25  #include "math_private.h"
   15.26  
   15.27 -libm_hidden_proto(fabs)
   15.28 -#ifdef __STDC__
   15.29 -     double fabs(double x)
   15.30 -#else
   15.31 -     double fabs(x)
   15.32 -     double x;
   15.33 -#endif
   15.34 +double fabs(double x)
   15.35  {
   15.36 -    u_int32_t high;
   15.37 -    GET_HIGH_WORD(high, x);
   15.38 -    SET_HIGH_WORD(x, high & 0x7fffffff);
   15.39 -    return x;
   15.40 +	u_int32_t high;
   15.41 +	GET_HIGH_WORD(high,x);
   15.42 +	SET_HIGH_WORD(x,high&0x7fffffff);
   15.43 +        return x;
   15.44  }
   15.45 -
   15.46  libm_hidden_def(fabs)
    16.1 --- a/src/libm/s_floor.c	Sat Nov 04 15:34:14 2017 -0700
    16.2 +++ b/src/libm/s_floor.c	Sat Nov 04 15:53:19 2017 -0700
    16.3 @@ -1,4 +1,3 @@
    16.4 -/* @(#)s_floor.c 5.1 93/09/24 */
    16.5  /*
    16.6   * ====================================================
    16.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    16.8 @@ -10,11 +9,6 @@
    16.9   * ====================================================
   16.10   */
   16.11  
   16.12 -#if defined(LIBM_SCCS) && !defined(lint)
   16.13 -static const char rcsid[] =
   16.14 -    "$NetBSD: s_floor.c,v 1.8 1995/05/10 20:47:20 jtc Exp $";
   16.15 -#endif
   16.16 -
   16.17  /*
   16.18   * floor(x)
   16.19   * Return x rounded toward -inf to integral value
   16.20 @@ -24,73 +18,54 @@
   16.21   *	Inexact flag raised if x not equal to floor(x).
   16.22   */
   16.23  
   16.24 +/*#include <features.h>*/
   16.25 +/* Prevent math.h from defining a colliding inline */
   16.26 +#undef __USE_EXTERN_INLINES
   16.27  #include "math_libm.h"
   16.28  #include "math_private.h"
   16.29  
   16.30 -#ifdef __STDC__
   16.31 -static const double huge_val = 1.0e300;
   16.32 -#else
   16.33 -static double huge_val = 1.0e300;
   16.34 -#endif
   16.35 +static const double huge = 1.0e300;
   16.36  
   16.37 -libm_hidden_proto(floor)
   16.38 -#ifdef __STDC__
   16.39 -     double floor(double x)
   16.40 -#else
   16.41 -     double floor(x)
   16.42 -     double x;
   16.43 -#endif
   16.44 +double floor(double x)
   16.45  {
   16.46 -    int32_t i0, i1, j0;
   16.47 -    u_int32_t i, j;
   16.48 -    EXTRACT_WORDS(i0, i1, x);
   16.49 -    j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
   16.50 -    if (j0 < 20) {
   16.51 -        if (j0 < 0) {           /* raise inexact if x != 0 */
   16.52 -            if (huge_val + x > 0.0) {       /* return 0*sign(x) if |x|<1 */
   16.53 -                if (i0 >= 0) {
   16.54 -                    i0 = i1 = 0;
   16.55 -                } else if (((i0 & 0x7fffffff) | i1) != 0) {
   16.56 -                    i0 = 0xbff00000;
   16.57 -                    i1 = 0;
   16.58 -                }
   16.59 -            }
   16.60 -        } else {
   16.61 -            i = (0x000fffff) >> j0;
   16.62 -            if (((i0 & i) | i1) == 0)
   16.63 -                return x;       /* x is integral */
   16.64 -            if (huge_val + x > 0.0) {       /* raise inexact flag */
   16.65 -                if (i0 < 0)
   16.66 -                    i0 += (0x00100000) >> j0;
   16.67 -                i0 &= (~i);
   16.68 -                i1 = 0;
   16.69 -            }
   16.70 -        }
   16.71 -    } else if (j0 > 51) {
   16.72 -        if (j0 == 0x400)
   16.73 -            return x + x;       /* inf or NaN */
   16.74 -        else
   16.75 -            return x;           /* x is integral */
   16.76 -    } else {
   16.77 -        i = ((u_int32_t) (0xffffffff)) >> (j0 - 20);
   16.78 -        if ((i1 & i) == 0)
   16.79 -            return x;           /* x is integral */
   16.80 -        if (huge_val + x > 0.0) {   /* raise inexact flag */
   16.81 -            if (i0 < 0) {
   16.82 -                if (j0 == 20)
   16.83 -                    i0 += 1;
   16.84 -                else {
   16.85 -                    j = i1 + (1 << (52 - j0));
   16.86 -                    if (j < (u_int32_t) i1)
   16.87 -                        i0 += 1;        /* got a carry */
   16.88 -                    i1 = j;
   16.89 -                }
   16.90 -            }
   16.91 -            i1 &= (~i);
   16.92 -        }
   16.93 -    }
   16.94 -    INSERT_WORDS(x, i0, i1);
   16.95 -    return x;
   16.96 +	int32_t i0,i1,j0;
   16.97 +	u_int32_t i,j;
   16.98 +	EXTRACT_WORDS(i0,i1,x);
   16.99 +	j0 = ((i0>>20)&0x7ff)-0x3ff;
  16.100 +	if(j0<20) {
  16.101 +	    if(j0<0) { 	/* raise inexact if x != 0 */
  16.102 +		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
  16.103 +		    if(i0>=0) {i0=i1=0;}
  16.104 +		    else if(((i0&0x7fffffff)|i1)!=0)
  16.105 +			{ i0=0xbff00000;i1=0;}
  16.106 +		}
  16.107 +	    } else {
  16.108 +		i = (0x000fffff)>>j0;
  16.109 +		if(((i0&i)|i1)==0) return x; /* x is integral */
  16.110 +		if(huge+x>0.0) {	/* raise inexact flag */
  16.111 +		    if(i0<0) i0 += (0x00100000)>>j0;
  16.112 +		    i0 &= (~i); i1=0;
  16.113 +		}
  16.114 +	    }
  16.115 +	} else if (j0>51) {
  16.116 +	    if(j0==0x400) return x+x;	/* inf or NaN */
  16.117 +	    else return x;		/* x is integral */
  16.118 +	} else {
  16.119 +	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
  16.120 +	    if((i1&i)==0) return x;	/* x is integral */
  16.121 +	    if(huge+x>0.0) { 		/* raise inexact flag */
  16.122 +		if(i0<0) {
  16.123 +		    if(j0==20) i0+=1;
  16.124 +		    else {
  16.125 +			j = i1+(1<<(52-j0));
  16.126 +			if(j<i1) i0 +=1 ; 	/* got a carry */
  16.127 +			i1=j;
  16.128 +		    }
  16.129 +		}
  16.130 +		i1 &= (~i);
  16.131 +	    }
  16.132 +	}
  16.133 +	INSERT_WORDS(x,i0,i1);
  16.134 +	return x;
  16.135  }
  16.136 -
  16.137  libm_hidden_def(floor)
    17.1 --- a/src/libm/s_scalbn.c	Sat Nov 04 15:34:14 2017 -0700
    17.2 +++ b/src/libm/s_scalbn.c	Sat Nov 04 15:53:19 2017 -0700
    17.3 @@ -1,4 +1,3 @@
    17.4 -/* @(#)s_scalbn.c 5.1 93/09/24 */
    17.5  /*
    17.6   * ====================================================
    17.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    17.8 @@ -10,70 +9,69 @@
    17.9   * ====================================================
   17.10   */
   17.11  
   17.12 -#if defined(LIBM_SCCS) && !defined(lint)
   17.13 -static const char rcsid[] =
   17.14 -    "$NetBSD: s_scalbn.c,v 1.8 1995/05/10 20:48:08 jtc Exp $";
   17.15 -#endif
   17.16 -
   17.17  /*
   17.18 - * scalbn (double x, int n)
   17.19 - * scalbn(x,n) returns x* 2**n  computed by  exponent
   17.20 + * scalbln(double x, long n)
   17.21 + * scalbln(x,n) returns x * 2**n computed by exponent
   17.22   * manipulation rather than by actually performing an
   17.23   * exponentiation or a multiplication.
   17.24   */
   17.25  
   17.26  #include "math_libm.h"
   17.27  #include "math_private.h"
   17.28 +#include <limits.h>
   17.29  
   17.30 -libm_hidden_proto(copysign)
   17.31 -#ifdef __STDC__
   17.32 -     static const double
   17.33 +static const double
   17.34 +two54  = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
   17.35 +twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
   17.36 +huge   = 1.0e+300,
   17.37 +tiny   = 1.0e-300;
   17.38 +
   17.39 +double scalbln(double x, long n)
   17.40 +{
   17.41 +	int32_t k, hx, lx;
   17.42 +
   17.43 +	EXTRACT_WORDS(hx, lx, x);
   17.44 +	k = (hx & 0x7ff00000) >> 20; /* extract exponent */
   17.45 +	if (k == 0) { /* 0 or subnormal x */
   17.46 +		if ((lx | (hx & 0x7fffffff)) == 0)
   17.47 +			return x; /* +-0 */
   17.48 +		x *= two54;
   17.49 +		GET_HIGH_WORD(hx, x);
   17.50 +		k = ((hx & 0x7ff00000) >> 20) - 54;
   17.51 +	}
   17.52 +	if (k == 0x7ff)
   17.53 +		return x + x; /* NaN or Inf */
   17.54 +	k = k + n;
   17.55 +	if (k > 0x7fe)
   17.56 +		return huge * copysign(huge, x); /* overflow */
   17.57 +	if (n < -50000)
   17.58 +		return tiny * copysign(tiny, x); /* underflow */
   17.59 +	if (k > 0) { /* normal result */
   17.60 +		SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
   17.61 +		return x;
   17.62 +	}
   17.63 +	if (k <= -54) {
   17.64 +		if (n > 50000) /* in case integer overflow in n+k */
   17.65 +			return huge * copysign(huge, x); /* overflow */
   17.66 +		return tiny * copysign(tiny, x); /* underflow */
   17.67 +	}
   17.68 +	k += 54; /* subnormal result */
   17.69 +	SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
   17.70 +	return x * twom54;
   17.71 +}
   17.72 +libm_hidden_def(scalbln)
   17.73 +
   17.74 +#if LONG_MAX == INT_MAX
   17.75 +/* strong_alias(scalbln, scalbn) - "error: conflicting types for 'scalbn'"
   17.76 + * because it tries to declare "typeof(scalbln) scalbn;"
   17.77 + * which tries to give "long" parameter to scalbn.
   17.78 + * Doing it by hand:
   17.79 + */
   17.80 +__typeof(scalbn) scalbn __attribute__((alias("scalbln")));
   17.81  #else
   17.82 -     static double
   17.83 +double scalbn(double x, int n)
   17.84 +{
   17.85 +	return scalbln(x, n);
   17.86 +}
   17.87  #endif
   17.88 -       two54 = 1.80143985094819840000e+16,      /* 0x43500000, 0x00000000 */
   17.89 -         twom54 = 5.55111512312578270212e-17,   /* 0x3C900000, 0x00000000 */
   17.90 -         huge_val = 1.0e+300, tiny = 1.0e-300;
   17.91 -
   17.92 -libm_hidden_proto(scalbn)
   17.93 -#ifdef __STDC__
   17.94 -     double scalbn(double x, int n)
   17.95 -#else
   17.96 -     double scalbn(x, n)
   17.97 -     double x;
   17.98 -     int n;
   17.99 -#endif
  17.100 -{
  17.101 -    int32_t k, hx, lx;
  17.102 -    EXTRACT_WORDS(hx, lx, x);
  17.103 -    k = (hx & 0x7ff00000) >> 20;        /* extract exponent */
  17.104 -    if (k == 0) {               /* 0 or subnormal x */
  17.105 -        if ((lx | (hx & 0x7fffffff)) == 0)
  17.106 -            return x;           /* +-0 */
  17.107 -        x *= two54;
  17.108 -        GET_HIGH_WORD(hx, x);
  17.109 -        k = ((hx & 0x7ff00000) >> 20) - 54;
  17.110 -        if (n < -50000)
  17.111 -            return tiny * x;    /* underflow */
  17.112 -    }
  17.113 -    if (k == 0x7ff)
  17.114 -        return x + x;           /* NaN or Inf */
  17.115 -    k = k + n;
  17.116 -    if (k > 0x7fe)
  17.117 -        return huge_val * copysign(huge_val, x);        /* overflow  */
  17.118 -    if (k > 0) {                /* normal result */
  17.119 -        SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
  17.120 -        return x;
  17.121 -    }
  17.122 -    if (k <= -54) {
  17.123 -        if (n > 50000)          /* in case integer overflow in n+k */
  17.124 -            return huge_val * copysign(huge_val, x);    /* overflow */
  17.125 -        else
  17.126 -            return tiny * copysign(tiny, x);    /* underflow */
  17.127 -    }
  17.128 -    k += 54;                    /* subnormal result */
  17.129 -    SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
  17.130 -    return x * twom54;
  17.131 -}
  17.132 -
  17.133  libm_hidden_def(scalbn)
    18.1 --- a/src/libm/s_sin.c	Sat Nov 04 15:34:14 2017 -0700
    18.2 +++ b/src/libm/s_sin.c	Sat Nov 04 15:53:19 2017 -0700
    18.3 @@ -1,4 +1,3 @@
    18.4 -/* @(#)s_sin.c 5.1 93/09/24 */
    18.5  /*
    18.6   * ====================================================
    18.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    18.8 @@ -10,11 +9,6 @@
    18.9   * ====================================================
   18.10   */
   18.11  
   18.12 -#if defined(LIBM_SCCS) && !defined(lint)
   18.13 -static const char rcsid[] =
   18.14 -    "$NetBSD: s_sin.c,v 1.7 1995/05/10 20:48:15 jtc Exp $";
   18.15 -#endif
   18.16 -
   18.17  /* sin(x)
   18.18   * Return sine function of x.
   18.19   *
   18.20 @@ -49,43 +43,31 @@
   18.21  #include "math_libm.h"
   18.22  #include "math_private.h"
   18.23  
   18.24 -libm_hidden_proto(sin)
   18.25 -#ifdef __STDC__
   18.26 -     double sin(double x)
   18.27 -#else
   18.28 -     double sin(x)
   18.29 -     double x;
   18.30 -#endif
   18.31 +double sin(double x)
   18.32  {
   18.33 -    double y[2], z = 0.0;
   18.34 -    int32_t n, ix;
   18.35 +	double y[2],z=0.0;
   18.36 +	int32_t n, ix;
   18.37  
   18.38      /* High word of x. */
   18.39 -    GET_HIGH_WORD(ix, x);
   18.40 +	GET_HIGH_WORD(ix,x);
   18.41  
   18.42      /* |x| ~< pi/4 */
   18.43 -    ix &= 0x7fffffff;
   18.44 -    if (ix <= 0x3fe921fb)
   18.45 -        return __kernel_sin(x, z, 0);
   18.46 +	ix &= 0x7fffffff;
   18.47 +	if(ix <= 0x3fe921fb) return __kernel_sin(x,z,0);
   18.48  
   18.49      /* sin(Inf or NaN) is NaN */
   18.50 -    else if (ix >= 0x7ff00000)
   18.51 -        return x - x;
   18.52 +	else if (ix>=0x7ff00000) return x-x;
   18.53  
   18.54      /* argument reduction needed */
   18.55 -    else {
   18.56 -        n = __ieee754_rem_pio2(x, y);
   18.57 -        switch (n & 3) {
   18.58 -        case 0:
   18.59 -            return __kernel_sin(y[0], y[1], 1);
   18.60 -        case 1:
   18.61 -            return __kernel_cos(y[0], y[1]);
   18.62 -        case 2:
   18.63 -            return -__kernel_sin(y[0], y[1], 1);
   18.64 -        default:
   18.65 -            return -__kernel_cos(y[0], y[1]);
   18.66 -        }
   18.67 -    }
   18.68 +	else {
   18.69 +	    n = __ieee754_rem_pio2(x,y);
   18.70 +	    switch(n&3) {
   18.71 +		case 0: return  __kernel_sin(y[0],y[1],1);
   18.72 +		case 1: return  __kernel_cos(y[0],y[1]);
   18.73 +		case 2: return -__kernel_sin(y[0],y[1],1);
   18.74 +		default:
   18.75 +			return -__kernel_cos(y[0],y[1]);
   18.76 +	    }
   18.77 +	}
   18.78  }
   18.79 -
   18.80  libm_hidden_def(sin)