From 4088e3a587b6f84f9b51d5e81332b1a56df4b225 Mon Sep 17 00:00:00 2001
From: Elliott Hughes <enh@google.com>
Date: Thu, 3 Aug 2023 13:33:56 -0700
Subject: [PATCH] Sync upstream FreeBSD libm.

Test: treehugger
Change-Id: I583a3e93821d512c975db34fc1610ffd22445d58
---
 libm/upstream-freebsd/lib/msun/ld128/s_expl.c |  28 +--
 libm/upstream-freebsd/lib/msun/ld128/s_logl.c |  24 +-
 libm/upstream-freebsd/lib/msun/src/e_acos.c   |   4 +-
 libm/upstream-freebsd/lib/msun/src/e_acosf.c  |   2 +-
 libm/upstream-freebsd/lib/msun/src/e_acosh.c  |   8 +-
 libm/upstream-freebsd/lib/msun/src/e_acoshf.c |   8 +-
 libm/upstream-freebsd/lib/msun/src/e_asin.c   |   4 +-
 libm/upstream-freebsd/lib/msun/src/e_asinf.c  |   2 +-
 libm/upstream-freebsd/lib/msun/src/e_atan2.c  |   4 +-
 libm/upstream-freebsd/lib/msun/src/e_atan2f.c |   2 +-
 libm/upstream-freebsd/lib/msun/src/e_atanh.c  |   4 +-
 libm/upstream-freebsd/lib/msun/src/e_atanhf.c |   2 +-
 libm/upstream-freebsd/lib/msun/src/e_cosh.c   |   8 +-
 libm/upstream-freebsd/lib/msun/src/e_coshf.c  |   6 +-
 libm/upstream-freebsd/lib/msun/src/e_exp.c    |   4 +-
 libm/upstream-freebsd/lib/msun/src/e_expf.c   |   2 +-
 libm/upstream-freebsd/lib/msun/src/e_fmod.c   |   4 +-
 libm/upstream-freebsd/lib/msun/src/e_fmodf.c  |   4 +-
 libm/upstream-freebsd/lib/msun/src/e_gamma.c  |   8 +-
 .../upstream-freebsd/lib/msun/src/e_gamma_r.c |   8 +-
 libm/upstream-freebsd/lib/msun/src/e_gammaf.c |   8 +-
 .../lib/msun/src/e_gammaf_r.c                 |   8 +-
 libm/upstream-freebsd/lib/msun/src/e_hypot.c  |   4 +-
 libm/upstream-freebsd/lib/msun/src/e_hypotf.c |   6 +-
 libm/upstream-freebsd/lib/msun/src/e_j0.c     |  10 +-
 libm/upstream-freebsd/lib/msun/src/e_j0f.c    |   8 +-
 libm/upstream-freebsd/lib/msun/src/e_j1.c     |   8 +-
 libm/upstream-freebsd/lib/msun/src/e_j1f.c    |   6 +-
 libm/upstream-freebsd/lib/msun/src/e_jn.c     |  28 +--
 libm/upstream-freebsd/lib/msun/src/e_jnf.c    |  26 +-
 libm/upstream-freebsd/lib/msun/src/e_lgamma.c |   8 +-
 .../lib/msun/src/e_lgamma_r.c                 |  16 +-
 .../upstream-freebsd/lib/msun/src/e_lgammaf.c |   8 +-
 .../lib/msun/src/e_lgammaf_r.c                |  14 +-
 libm/upstream-freebsd/lib/msun/src/e_log.c    |   4 +-
 libm/upstream-freebsd/lib/msun/src/e_log10.c  |   2 +-
 libm/upstream-freebsd/lib/msun/src/e_log10f.c |   2 +-
 libm/upstream-freebsd/lib/msun/src/e_log2.c   |   2 +-
 libm/upstream-freebsd/lib/msun/src/e_log2f.c  |   2 +-
 libm/upstream-freebsd/lib/msun/src/e_logf.c   |   2 +-
 libm/upstream-freebsd/lib/msun/src/e_pow.c    |   4 +-
 libm/upstream-freebsd/lib/msun/src/e_powf.c   |   4 +-
 .../lib/msun/src/e_remainder.c                |   6 +-
 .../lib/msun/src/e_remainderf.c               |   4 +-
 libm/upstream-freebsd/lib/msun/src/e_scalb.c  |   6 +-
 libm/upstream-freebsd/lib/msun/src/e_scalbf.c |   4 +-
 libm/upstream-freebsd/lib/msun/src/e_sinh.c   |   6 +-
 libm/upstream-freebsd/lib/msun/src/e_sinhf.c  |   4 +-
 .../lib/msun/src/math_private.h               | 231 +++++-------------
 libm/upstream-freebsd/lib/msun/src/s_asinh.c  |   6 +-
 libm/upstream-freebsd/lib/msun/src/s_asinhf.c |   6 +-
 libm/upstream-freebsd/lib/msun/src/s_cospi.c  |  20 +-
 libm/upstream-freebsd/lib/msun/src/s_erf.c    |   4 +-
 .../lib/msun/src/s_significand.c              |   2 +-
 .../lib/msun/src/s_significandf.c             |   2 +-
 libm/upstream-freebsd/lib/msun/src/s_sinpi.c  |  13 +-
 56 files changed, 245 insertions(+), 385 deletions(-)

diff --git a/libm/upstream-freebsd/lib/msun/ld128/s_expl.c b/libm/upstream-freebsd/lib/msun/ld128/s_expl.c
index 5fc43802b..0274a8f30 100644
--- a/libm/upstream-freebsd/lib/msun/ld128/s_expl.c
+++ b/libm/upstream-freebsd/lib/msun/ld128/s_expl.c
@@ -65,8 +65,6 @@ expl(long double x)
 	int k;
 	uint16_t hx, ix;
 
-	DOPRINT_START(&x);
-
 	/* Filter out exceptional cases. */
 	u.e = x;
 	hx = u.xbits.expsign;
@@ -74,15 +72,15 @@ expl(long double x)
 	if (ix >= BIAS + 13) {		/* |x| >= 8192 or x is NaN */
 		if (ix == BIAS + LDBL_MAX_EXP) {
 			if (hx & 0x8000)  /* x is -Inf or -NaN */
-				RETURNP(-1 / x);
-			RETURNP(x + x);	/* x is +Inf or +NaN */
+				RETURNF(-1 / x);
+			RETURNF(x + x);	/* x is +Inf or +NaN */
 		}
 		if (x > o_threshold)
-			RETURNP(huge * huge);
+			RETURNF(huge * huge);
 		if (x < u_threshold)
-			RETURNP(tiny * tiny);
+			RETURNF(tiny * tiny);
 	} else if (ix < BIAS - 114) {	/* |x| < 0x1p-114 */
-		RETURN2P(1, x);		/* 1 with inexact iff x != 0 */
+		RETURNF(1 + x);		/* 1 with inexact iff x != 0 */
 	}
 
 	ENTERI();
@@ -210,8 +208,6 @@ expm1l(long double x)
 	int k, n, n2;
 	uint16_t hx, ix;
 
-	DOPRINT_START(&x);
-
 	/* Filter out exceptional cases. */
 	u.e = x;
 	hx = u.xbits.expsign;
@@ -219,11 +215,11 @@ expm1l(long double x)
 	if (ix >= BIAS + 7) {		/* |x| >= 128 or x is NaN */
 		if (ix == BIAS + LDBL_MAX_EXP) {
 			if (hx & 0x8000)  /* x is -Inf or -NaN */
-				RETURNP(-1 / x - 1);
-			RETURNP(x + x);	/* x is +Inf or +NaN */
+				RETURNF(-1 / x - 1);
+			RETURNF(x + x);	/* x is +Inf or +NaN */
 		}
 		if (x > o_threshold)
-			RETURNP(huge * huge);
+			RETURNF(huge * huge);
 		/*
 		 * expm1l() never underflows, but it must avoid
 		 * unrepresentable large negative exponents.  We used a
@@ -232,7 +228,7 @@ expm1l(long double x)
 		 * in the same way as large ones here.
 		 */
 		if (hx & 0x8000)	/* x <= -128 */
-			RETURN2P(tiny, -1);	/* good for x < -114ln2 - eps */
+			RETURNF(tiny - 1);	/* good for x < -114ln2 - eps */
 	}
 
 	ENTERI();
@@ -244,7 +240,7 @@ expm1l(long double x)
 		if (x < T3) {
 			if (ix < BIAS - 113) {	/* |x| < 0x1p-113 */
 				/* x (rounded) with inexact if x != 0: */
-				RETURNPI(x == 0 ? x :
+				RETURNI(x == 0 ? x :
 				    (0x1p200 * x + fabsl(x)) * 0x1p-200);
 			}
 			q = x * x2 * C3 + x2 * x2 * (C4 + x * (C5 + x * (C6 +
@@ -265,9 +261,9 @@ expm1l(long double x)
 		hx2_hi = x_hi * x_hi / 2;
 		hx2_lo = x_lo * (x + x_hi) / 2;
 		if (ix >= BIAS - 7)
-			RETURN2PI(hx2_hi + x_hi, hx2_lo + x_lo + q);
+			RETURNI((hx2_hi + x_hi) + (hx2_lo + x_lo + q));
 		else
-			RETURN2PI(x, hx2_lo + q + hx2_hi);
+			RETURNI(x + (hx2_lo + q + hx2_hi));
 	}
 
 	/* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
diff --git a/libm/upstream-freebsd/lib/msun/ld128/s_logl.c b/libm/upstream-freebsd/lib/msun/ld128/s_logl.c
index 40a22c0f1..bc538840a 100644
--- a/libm/upstream-freebsd/lib/msun/ld128/s_logl.c
+++ b/libm/upstream-freebsd/lib/msun/ld128/s_logl.c
@@ -573,24 +573,23 @@ log1pl(long double x)
 	int i, k;
 	int16_t ax, hx;
 
-	DOPRINT_START(&x);
 	EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
 	if (hx < 0x3fff) {		/* x < 1, or x neg NaN */
 		ax = hx & 0x7fff;
 		if (ax >= 0x3fff) {	/* x <= -1, or x neg NaN */
 			if (ax == 0x3fff && (lx | llx) == 0)
-				RETURNP(-1 / zero);	/* log1p(-1) = -Inf */
+				RETURNF(-1 / zero);	/* log1p(-1) = -Inf */
 			/* log1p(x < 1, or x NaN) = qNaN: */
-			RETURNP((x - x) / (x - x));
+			RETURNF((x - x) / (x - x));
 		}
 		if (ax <= 0x3f8d) {	/* |x| < 2**-113 */
 			if ((int)x == 0)
-				RETURNP(x);	/* x with inexact if x != 0 */
+				RETURNF(x);	/* x with inexact if x != 0 */
 		}
 		f_hi = 1;
 		f_lo = x;
 	} else if (hx >= 0x7fff) {	/* x +Inf or non-neg NaN */
-		RETURNP(x + x);		/* log1p(Inf or NaN) = Inf or qNaN */
+		RETURNF(x + x);		/* log1p(Inf or NaN) = Inf or qNaN */
 	} else if (hx < 0x40e1) {	/* 1 <= x < 2**226 */
 		f_hi = x;
 		f_lo = 1;
@@ -669,7 +668,7 @@ log1pl(long double x)
 #endif
 
 	_3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
-	RETURN2PI(val_hi, val_lo);
+	RETURNI(val_hi + val_lo);
 }
 
 #ifdef STRUCT_RETURN
@@ -680,7 +679,6 @@ logl(long double x)
 	struct ld r;
 
 	ENTERI();
-	DOPRINT_START(&x);
 	k_logl(x, &r);
 	RETURNSPI(&r);
 }
@@ -708,15 +706,13 @@ log10l(long double x)
 	long double hi, lo;
 
 	ENTERI();
-	DOPRINT_START(&x);
 	k_logl(x, &r);
 	if (!r.lo_set)
-		RETURNPI(r.hi);
+		RETURNI(r.hi);
 	_2sumF(r.hi, r.lo);
 	hi = (float)r.hi;
 	lo = r.lo + (r.hi - hi);
-	RETURN2PI(invln10_hi * hi,
-	    invln10_lo_plus_hi * lo + invln10_lo * hi);
+	RETURNI(invln10_hi * hi + (invln10_lo_plus_hi * lo + invln10_lo * hi));
 }
 
 long double
@@ -726,15 +722,13 @@ log2l(long double x)
 	long double hi, lo;
 
 	ENTERI();
-	DOPRINT_START(&x);
 	k_logl(x, &r);
 	if (!r.lo_set)
-		RETURNPI(r.hi);
+		RETURNI(r.hi);
 	_2sumF(r.hi, r.lo);
 	hi = (float)r.hi;
 	lo = r.lo + (r.hi - hi);
-	RETURN2PI(invln2_hi * hi,
-	    invln2_lo_plus_hi * lo + invln2_lo * hi);
+	RETURNI(invln2_hi * hi + (invln2_lo_plus_hi * lo + invln2_lo * hi));
 }
 
 #endif /* STRUCT_RETURN */
diff --git a/libm/upstream-freebsd/lib/msun/src/e_acos.c b/libm/upstream-freebsd/lib/msun/src/e_acos.c
index 1f6dca5bb..6623355ba 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_acos.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_acos.c
@@ -14,7 +14,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_acos(x)
+/* acos(x)
  * Method :                  
  *	acos(x)  = pi/2 - asin(x)
  *	acos(-x) = pi/2 + asin(x)
@@ -62,7 +62,7 @@ qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
 qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 
 double
-__ieee754_acos(double x)
+acos(double x)
 {
 	double z,p,q,r,w,s,c,df;
 	int32_t hx,ix;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_acosf.c b/libm/upstream-freebsd/lib/msun/src/e_acosf.c
index c9f62cc40..64f1c5afb 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_acosf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_acosf.c
@@ -32,7 +32,7 @@ pS2 = -8.6563630030e-03,
 qS1 = -7.0662963390e-01;
 
 float
-__ieee754_acosf(float x)
+acosf(float x)
 {
 	float z,p,q,r,w,s,c,df;
 	int32_t hx,ix;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_acosh.c b/libm/upstream-freebsd/lib/msun/src/e_acosh.c
index 358c8bd6a..794799582 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_acosh.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_acosh.c
@@ -15,7 +15,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_acosh(x)
+/* acosh(x)
  * Method :
  *	Based on 
  *		acosh(x) = log [ x + sqrt(x*x-1) ]
@@ -39,7 +39,7 @@ one	= 1.0,
 ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
 
 double
-__ieee754_acosh(double x)
+acosh(double x)
 {
 	double t;
 	int32_t hx;
@@ -51,12 +51,12 @@ __ieee754_acosh(double x)
 	    if(hx >=0x7ff00000) {	/* x is inf of NaN */
 	        return x+x;
 	    } else 
-		return __ieee754_log(x)+ln2;	/* acosh(huge)=log(2x) */
+		return log(x)+ln2;	/* acosh(huge)=log(2x) */
 	} else if(((hx-0x3ff00000)|lx)==0) {
 	    return 0.0;			/* acosh(1) = 0 */
 	} else if (hx > 0x40000000) {	/* 2**28 > x > 2 */
 	    t=x*x;
-	    return __ieee754_log(2.0*x-one/(x+sqrt(t-one)));
+	    return log(2.0*x-one/(x+sqrt(t-one)));
 	} else {			/* 1<x<2 */
 	    t = x-one;
 	    return log1p(t+sqrt(2.0*t+t*t));
diff --git a/libm/upstream-freebsd/lib/msun/src/e_acoshf.c b/libm/upstream-freebsd/lib/msun/src/e_acoshf.c
index f529b20d1..781ccf2a1 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_acoshf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_acoshf.c
@@ -24,7 +24,7 @@ one	= 1.0,
 ln2	= 6.9314718246e-01;  /* 0x3f317218 */
 
 float
-__ieee754_acoshf(float x)
+acoshf(float x)
 {
 	float t;
 	int32_t hx;
@@ -35,14 +35,14 @@ __ieee754_acoshf(float x)
 	    if(hx >=0x7f800000) {	/* x is inf of NaN */
 	        return x+x;
 	    } else
-		return __ieee754_logf(x)+ln2;	/* acosh(huge)=log(2x) */
+		return logf(x)+ln2;	/* acosh(huge)=log(2x) */
 	} else if (hx==0x3f800000) {
 	    return 0.0;			/* acosh(1) = 0 */
 	} else if (hx > 0x40000000) {	/* 2**28 > x > 2 */
 	    t=x*x;
-	    return __ieee754_logf((float)2.0*x-one/(x+__ieee754_sqrtf(t-one)));
+	    return logf((float)2.0*x-one/(x+sqrtf(t-one)));
 	} else {			/* 1<x<2 */
 	    t = x-one;
-	    return log1pf(t+__ieee754_sqrtf((float)2.0*t+t*t));
+	    return log1pf(t+sqrtf((float)2.0*t+t*t));
 	}
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/e_asin.c b/libm/upstream-freebsd/lib/msun/src/e_asin.c
index 931b27049..fa180ab57 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_asin.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_asin.c
@@ -14,7 +14,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_asin(x)
+/* asin(x)
  * Method :                  
  *	Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
  *	we approximate asin(x) on [0,0.5] by
@@ -68,7 +68,7 @@ qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
 qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 
 double
-__ieee754_asin(double x)
+asin(double x)
 {
 	double t=0.0,w,p,q,c,r,s;
 	int32_t hx,ix;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_asinf.c b/libm/upstream-freebsd/lib/msun/src/e_asinf.c
index deaabb6a3..db4b9b603 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_asinf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_asinf.c
@@ -32,7 +32,7 @@ static const double
 pio2 =  1.570796326794896558e+00;
 
 float
-__ieee754_asinf(float x)
+asinf(float x)
 {
 	double s;
 	float t,w,p,q;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_atan2.c b/libm/upstream-freebsd/lib/msun/src/e_atan2.c
index 231a1611e..0b2e72102 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_atan2.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_atan2.c
@@ -15,7 +15,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_atan2(y,x)
+/* atan2(y,x)
  * Method :
  *	1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
  *	2. Reduce x to positive by (if x and y are unexceptional): 
@@ -58,7 +58,7 @@ static volatile double
 pi_lo   = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
 
 double
-__ieee754_atan2(double y, double x)
+atan2(double y, double x)
 {
 	double z;
 	int32_t k,m,hx,hy,ix,iy;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_atan2f.c b/libm/upstream-freebsd/lib/msun/src/e_atan2f.c
index 346d76746..4ea001df9 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_atan2f.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_atan2f.c
@@ -30,7 +30,7 @@ static volatile float
 pi_lo   = -8.7422776573e-08; /* 0xb3bbbd2e */
 
 float
-__ieee754_atan2f(float y, float x)
+atan2f(float y, float x)
 {
 	float z;
 	int32_t k,m,hx,hy,ix,iy;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_atanh.c b/libm/upstream-freebsd/lib/msun/src/e_atanh.c
index 422ff2698..41f3bcaca 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_atanh.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_atanh.c
@@ -15,7 +15,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_atanh(x)
+/* atanh(x)
  * Method :
  *    1.Reduced x to positive by atanh(-x) = -atanh(x)
  *    2.For x>=0.5
@@ -42,7 +42,7 @@ static const double one = 1.0, huge = 1e300;
 static const double zero = 0.0;
 
 double
-__ieee754_atanh(double x)
+atanh(double x)
 {
 	double t;
 	int32_t hx,ix;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_atanhf.c b/libm/upstream-freebsd/lib/msun/src/e_atanhf.c
index 4bd6a8f9b..46643beb5 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_atanhf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_atanhf.c
@@ -24,7 +24,7 @@ static const float one = 1.0, huge = 1e30;
 static const float zero = 0.0;
 
 float
-__ieee754_atanhf(float x)
+atanhf(float x)
 {
 	float t;
 	int32_t hx,ix;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_cosh.c b/libm/upstream-freebsd/lib/msun/src/e_cosh.c
index 246b5fbec..071663eb4 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_cosh.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_cosh.c
@@ -14,7 +14,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_cosh(x)
+/* cosh(x)
  * Method : 
  * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
  *	1. Replace x by |x| (cosh(x) = cosh(-x)). 
@@ -43,7 +43,7 @@ __FBSDID("$FreeBSD$");
 static const double one = 1.0, half=0.5, huge = 1.0e300;
 
 double
-__ieee754_cosh(double x)
+cosh(double x)
 {
 	double t,w;
 	int32_t ix;
@@ -65,12 +65,12 @@ __ieee754_cosh(double x)
 
     /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
 	if (ix < 0x40360000) {
-		t = __ieee754_exp(fabs(x));
+		t = exp(fabs(x));
 		return half*t+half/t;
 	}
 
     /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
-	if (ix < 0x40862E42)  return half*__ieee754_exp(fabs(x));
+	if (ix < 0x40862E42)  return half*exp(fabs(x));
 
     /* |x| in [log(maxdouble), overflowthresold] */
 	if (ix<=0x408633CE)
diff --git a/libm/upstream-freebsd/lib/msun/src/e_coshf.c b/libm/upstream-freebsd/lib/msun/src/e_coshf.c
index 95a0d6ee6..1673315d7 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_coshf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_coshf.c
@@ -22,7 +22,7 @@ __FBSDID("$FreeBSD$");
 static const float one = 1.0, half=0.5, huge = 1.0e30;
 
 float
-__ieee754_coshf(float x)
+coshf(float x)
 {
 	float t,w;
 	int32_t ix;
@@ -43,12 +43,12 @@ __ieee754_coshf(float x)
 
     /* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */
 	if (ix < 0x41100000) {
-		t = __ieee754_expf(fabsf(x));
+		t = expf(fabsf(x));
 		return half*t+half/t;
 	}
 
     /* |x| in [9, log(maxfloat)] return half*exp(|x|) */
-	if (ix < 0x42b17217)  return half*__ieee754_expf(fabsf(x));
+	if (ix < 0x42b17217)  return half*expf(fabsf(x));
 
     /* |x| in [log(maxfloat), overflowthresold] */
 	if (ix<=0x42b2d4fc)
diff --git a/libm/upstream-freebsd/lib/msun/src/e_exp.c b/libm/upstream-freebsd/lib/msun/src/e_exp.c
index dd04d8e83..59da39217 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_exp.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_exp.c
@@ -13,7 +13,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_exp(x)
+/* exp(x)
  * Returns the exponential of x.
  *
  * Method
@@ -102,7 +102,7 @@ huge	= 1.0e+300,
 twom1000= 9.33263618503218878990e-302;     /* 2**-1000=0x01700000,0*/
 
 double
-__ieee754_exp(double x)	/* default IEEE double exp */
+exp(double x)	/* default IEEE double exp */
 {
 	double y,hi=0.0,lo=0.0,c,t,twopk;
 	int32_t k=0,xsb;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_expf.c b/libm/upstream-freebsd/lib/msun/src/e_expf.c
index 4903d55c5..620d341be 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_expf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_expf.c
@@ -43,7 +43,7 @@ huge	= 1.0e+30,
 twom100 = 7.8886090522e-31;      /* 2**-100=0x0d800000 */
 
 float
-__ieee754_expf(float x)
+expf(float x)
 {
 	float y,hi=0.0,lo=0.0,c,t,twopk;
 	int32_t k=0,xsb;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_fmod.c b/libm/upstream-freebsd/lib/msun/src/e_fmod.c
index 3a28dc4ff..6d5f5332c 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_fmod.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_fmod.c
@@ -15,7 +15,7 @@
 __FBSDID("$FreeBSD$");
 
 /* 
- * __ieee754_fmod(x,y)
+ * fmod(x,y)
  * Return x mod y in exact arithmetic
  * Method: shift and subtract
  */
@@ -28,7 +28,7 @@ __FBSDID("$FreeBSD$");
 static const double one = 1.0, Zero[] = {0.0, -0.0,};
 
 double
-__ieee754_fmod(double x, double y)
+fmod(double x, double y)
 {
 	int32_t n,hx,hy,hz,ix,iy,sx,i;
 	u_int32_t lx,ly,lz;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_fmodf.c b/libm/upstream-freebsd/lib/msun/src/e_fmodf.c
index 1b6bf36f9..3cef9213c 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_fmodf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_fmodf.c
@@ -17,7 +17,7 @@
 __FBSDID("$FreeBSD$");
 
 /*
- * __ieee754_fmodf(x,y)
+ * fmodf(x,y)
  * Return x mod y in exact arithmetic
  * Method: shift and subtract
  */
@@ -28,7 +28,7 @@ __FBSDID("$FreeBSD$");
 static const float one = 1.0, Zero[] = {0.0, -0.0,};
 
 float
-__ieee754_fmodf(float x, float y)
+fmodf(float x, float y)
 {
 	int32_t n,hx,hy,hz,ix,iy,sx,i;
 
diff --git a/libm/upstream-freebsd/lib/msun/src/e_gamma.c b/libm/upstream-freebsd/lib/msun/src/e_gamma.c
index 28fb5ccba..a13f3e22c 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_gamma.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_gamma.c
@@ -15,10 +15,10 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_gamma(x)
+/* gamma(x)
  * Return the logarithm of the Gamma function of x.
  *
- * Method: call __ieee754_gamma_r
+ * Method: call gamma_r
  */
 
 #include "math.h"
@@ -27,7 +27,7 @@ __FBSDID("$FreeBSD$");
 extern int signgam;
 
 double
-__ieee754_gamma(double x)
+gamma(double x)
 {
-	return __ieee754_gamma_r(x,&signgam);
+	return gamma_r(x,&signgam);
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/e_gamma_r.c b/libm/upstream-freebsd/lib/msun/src/e_gamma_r.c
index 2c423dce8..2d996cabc 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_gamma_r.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_gamma_r.c
@@ -15,18 +15,18 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_gamma_r(x, signgamp)
+/* gamma_r(x, signgamp)
  * Reentrant version of the logarithm of the Gamma function 
  * with user provide pointer for the sign of Gamma(x). 
  *
- * Method: See __ieee754_lgamma_r
+ * Method: See lgamma_r
  */
 
 #include "math.h"
 #include "math_private.h"
 
 double
-__ieee754_gamma_r(double x, int *signgamp)
+gamma_r(double x, int *signgamp)
 {
-	return __ieee754_lgamma_r(x,signgamp);
+	return lgamma_r(x,signgamp);
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/e_gammaf.c b/libm/upstream-freebsd/lib/msun/src/e_gammaf.c
index c1b1668df..563c14822 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_gammaf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_gammaf.c
@@ -16,10 +16,10 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_gammaf(x)
+/* gammaf(x)
  * Return the logarithm of the Gamma function of x.
  *
- * Method: call __ieee754_gammaf_r
+ * Method: call gammaf_r
  */
 
 #include "math.h"
@@ -28,7 +28,7 @@ __FBSDID("$FreeBSD$");
 extern int signgam;
 
 float
-__ieee754_gammaf(float x)
+gammaf(float x)
 {
-	return __ieee754_gammaf_r(x,&signgam);
+	return gammaf_r(x,&signgam);
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/e_gammaf_r.c b/libm/upstream-freebsd/lib/msun/src/e_gammaf_r.c
index 9d7831b55..d7fc2db1e 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_gammaf_r.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_gammaf_r.c
@@ -16,18 +16,18 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_gammaf_r(x, signgamp)
+/* gammaf_r(x, signgamp)
  * Reentrant version of the logarithm of the Gamma function
  * with user provide pointer for the sign of Gamma(x).
  *
- * Method: See __ieee754_lgammaf_r
+ * Method: See lgammaf_r
  */
 
 #include "math.h"
 #include "math_private.h"
 
 float
-__ieee754_gammaf_r(float x, int *signgamp)
+gammaf_r(float x, int *signgamp)
 {
-	return __ieee754_lgammaf_r(x,signgamp);
+	return lgammaf_r(x,signgamp);
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/e_hypot.c b/libm/upstream-freebsd/lib/msun/src/e_hypot.c
index 7c455bb66..8e3f9317b 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_hypot.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_hypot.c
@@ -14,7 +14,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_hypot(x,y)
+/* hypot(x,y)
  *
  * Method :                  
  *	If (assume round-to-nearest) z=x*x+y*y 
@@ -52,7 +52,7 @@ __FBSDID("$FreeBSD$");
 #include "math_private.h"
 
 double
-__ieee754_hypot(double x, double y)
+hypot(double x, double y)
 {
 	double a,b,t1,t2,y1,y2,w;
 	int32_t j,k,ha,hb;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_hypotf.c b/libm/upstream-freebsd/lib/msun/src/e_hypotf.c
index 00610268f..a3b8c8666 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_hypotf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_hypotf.c
@@ -20,7 +20,7 @@ __FBSDID("$FreeBSD$");
 #include "math_private.h"
 
 float
-__ieee754_hypotf(float x, float y)
+hypotf(float x, float y)
 {
 	float a,b,t1,t2,y1,y2,w;
 	int32_t j,k,ha,hb;
@@ -67,14 +67,14 @@ __ieee754_hypotf(float x, float y)
 	if (w>b) {
 	    SET_FLOAT_WORD(t1,ha&0xfffff000);
 	    t2 = a-t1;
-	    w  = __ieee754_sqrtf(t1*t1-(b*(-b)-t2*(a+t1)));
+	    w  = sqrtf(t1*t1-(b*(-b)-t2*(a+t1)));
 	} else {
 	    a  = a+a;
 	    SET_FLOAT_WORD(y1,hb&0xfffff000);
 	    y2 = b - y1;
 	    SET_FLOAT_WORD(t1,(ha+0x00800000)&0xfffff000);
 	    t2 = a - t1;
-	    w  = __ieee754_sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
+	    w  = sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
 	}
 	if(k!=0) {
 	    SET_FLOAT_WORD(t1,(127+k)<<23);
diff --git a/libm/upstream-freebsd/lib/msun/src/e_j0.c b/libm/upstream-freebsd/lib/msun/src/e_j0.c
index 5d862b6f8..c43ab6996 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_j0.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_j0.c
@@ -13,7 +13,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_j0(x), __ieee754_y0(x)
+/* j0(x), y0(x)
  * Bessel function of the first and second kinds of order zero.
  * Method -- j0(x):
  *	1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
@@ -83,7 +83,7 @@ S04  =  1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
 static const double zero = 0, qrtr = 0.25;
 
 double
-__ieee754_j0(double x)
+j0(double x)
 {
 	double z, s,c,ss,cc,r,u,v;
 	int32_t hx,ix;
@@ -143,7 +143,7 @@ v03  =  2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */
 v04  =  4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */
 
 double
-__ieee754_y0(double x)
+y0(double x)
 {
 	double z, s,c,ss,cc,u,v;
 	int32_t hx,ix,lx;
@@ -192,12 +192,12 @@ __ieee754_y0(double x)
                 return z;
 	}
 	if(ix<=0x3e400000) {	/* x < 2**-27 */
-	    return(u00 + tpi*__ieee754_log(x));
+	    return(u00 + tpi*log(x));
 	}
 	z = x*x;
 	u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
 	v = one+z*(v01+z*(v02+z*(v03+z*v04)));
-	return(u/v + tpi*(__ieee754_j0(x)*__ieee754_log(x)));
+	return(u/v + tpi*(j0(x)*log(x)));
 }
 
 /* The asymptotic expansions of pzero is
diff --git a/libm/upstream-freebsd/lib/msun/src/e_j0f.c b/libm/upstream-freebsd/lib/msun/src/e_j0f.c
index 1c5ef4da1..290be04fc 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_j0f.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_j0f.c
@@ -45,7 +45,7 @@ S04  =  1.1661400734e-09; /* 0x30a045e8 */
 static const float zero = 0, qrtr = 0.25;
 
 float
-__ieee754_j0f(float x)
+j0f(float x)
 {
 	float z, s,c,ss,cc,r,u,v;
 	int32_t hx,ix;
@@ -105,7 +105,7 @@ v03  =  2.5915085189e-07, /* 0x348b216c */
 v04  =  4.4111031494e-10; /* 0x2ff280c2 */
 
 float
-__ieee754_y0f(float x)
+y0f(float x)
 {
 	float z, s,c,ss,cc,u,v;
 	int32_t hx,ix;
@@ -147,12 +147,12 @@ __ieee754_y0f(float x)
                 return z;
 	}
 	if(ix<=0x39000000) {	/* x < 2**-13 */
-	    return(u00 + tpi*__ieee754_logf(x));
+	    return(u00 + tpi*logf(x));
 	}
 	z = x*x;
 	u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
 	v = one+z*(v01+z*(v02+z*(v03+z*v04)));
-	return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x)));
+	return(u/v + tpi*(j0f(x)*logf(x)));
 }
 
 /* The asymptotic expansions of pzero is
diff --git a/libm/upstream-freebsd/lib/msun/src/e_j1.c b/libm/upstream-freebsd/lib/msun/src/e_j1.c
index fb4462738..ee3f6fcc6 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_j1.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_j1.c
@@ -13,7 +13,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_j1(x), __ieee754_y1(x)
+/* j1(x), y1(x)
  * Bessel function of the first and second kinds of order zero.
  * Method -- j1(x):
  *	1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
@@ -84,7 +84,7 @@ s05  =  1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */
 static const double zero    = 0.0;
 
 double
-__ieee754_j1(double x)
+j1(double x)
 {
 	double z, s,c,ss,cc,r,u,v,y;
 	int32_t hx,ix;
@@ -140,7 +140,7 @@ static const double V0[5] = {
 };
 
 double
-__ieee754_y1(double x)
+y1(double x)
 {
 	double z, s,c,ss,cc,u,v;
 	int32_t hx,ix,lx;
@@ -190,7 +190,7 @@ __ieee754_y1(double x)
         z = x*x;
         u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
         v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
-        return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x));
+        return(x*(u/v) + tpi*(j1(x)*log(x)-one/x));
 }
 
 /* For x >= 8, the asymptotic expansions of pone is
diff --git a/libm/upstream-freebsd/lib/msun/src/e_j1f.c b/libm/upstream-freebsd/lib/msun/src/e_j1f.c
index c6c45c107..e1f4498b7 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_j1f.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_j1f.c
@@ -46,7 +46,7 @@ s05  =  1.2354227016e-11; /* 0x2d59567e */
 static const float zero    = 0.0;
 
 float
-__ieee754_j1f(float x)
+j1f(float x)
 {
 	float z, s,c,ss,cc,r,u,v,y;
 	int32_t hx,ix;
@@ -102,7 +102,7 @@ static const float V0[5] = {
 };
 
 float
-__ieee754_y1f(float x)
+y1f(float x)
 {
 	float z, s,c,ss,cc,u,v;
 	int32_t hx,ix;
@@ -145,7 +145,7 @@ __ieee754_y1f(float x)
         z = x*x;
         u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
         v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
-        return(x*(u/v) + tpi*(__ieee754_j1f(x)*__ieee754_logf(x)-one/x));
+        return(x*(u/v) + tpi*(j1f(x)*logf(x)-one/x));
 }
 
 /* For x >= 8, the asymptotic expansions of pone is
diff --git a/libm/upstream-freebsd/lib/msun/src/e_jn.c b/libm/upstream-freebsd/lib/msun/src/e_jn.c
index 5aaebd400..6b876ce69 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_jn.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_jn.c
@@ -14,7 +14,7 @@
 __FBSDID("$FreeBSD$");
 
 /*
- * __ieee754_jn(n, x), __ieee754_yn(n, x)
+ * jn(n, x), yn(n, x)
  * floating point Bessel's function of the 1st and 2nd kind
  * of order n
  *
@@ -51,7 +51,7 @@ one   =  1.00000000000000000000e+00; /* 0x3FF00000, 0x00000000 */
 static const double zero  =  0.00000000000000000000e+00;
 
 double
-__ieee754_jn(int n, double x)
+jn(int n, double x)
 {
 	int32_t i,hx,ix,lx, sgn;
 	double a, b, c, s, temp, di;
@@ -69,8 +69,8 @@ __ieee754_jn(int n, double x)
 		x = -x;
 		hx ^= 0x80000000;
 	}
-	if(n==0) return(__ieee754_j0(x));
-	if(n==1) return(__ieee754_j1(x));
+	if(n==0) return(j0(x));
+	if(n==1) return(j1(x));
 	sgn = (n&1)&(hx>>31);	/* even n -- 0, odd n -- sign(x) */
 	x = fabs(x);
 	if((ix|lx)==0||ix>=0x7ff00000) 	/* if x is 0 or inf */
@@ -100,8 +100,8 @@ __ieee754_jn(int n, double x)
 		}
 		b = invsqrtpi*temp/sqrt(x);
 	    } else {
-	        a = __ieee754_j0(x);
-	        b = __ieee754_j1(x);
+	        a = j0(x);
+	        b = j1(x);
 	        for(i=1;i<n;i++){
 		    temp = b;
 		    b = b*((double)(i+i)/x) - a; /* avoid underflow */
@@ -177,7 +177,7 @@ __ieee754_jn(int n, double x)
 		 */
 		tmp = n;
 		v = two/x;
-		tmp = tmp*__ieee754_log(fabs(v*tmp));
+		tmp = tmp*log(fabs(v*tmp));
 		if(tmp<7.09782712893383973096e+02) {
 	    	    for(i=n-1,di=(double)(i+i);i>0;i--){
 		        temp = b;
@@ -201,8 +201,8 @@ __ieee754_jn(int n, double x)
 			}
 	     	    }
 		}
-		z = __ieee754_j0(x);
-		w = __ieee754_j1(x);
+		z = j0(x);
+		w = j1(x);
 		if (fabs(z) >= fabs(w))
 		    b = (t*z/b);
 		else
@@ -213,7 +213,7 @@ __ieee754_jn(int n, double x)
 }
 
 double
-__ieee754_yn(int n, double x)
+yn(int n, double x)
 {
 	int32_t i,hx,ix,lx;
 	int32_t sign;
@@ -232,8 +232,8 @@ __ieee754_yn(int n, double x)
 		n = -n;
 		sign = 1 - ((n&1)<<1);
 	}
-	if(n==0) return(__ieee754_y0(x));
-	if(n==1) return(sign*__ieee754_y1(x));
+	if(n==0) return(y0(x));
+	if(n==1) return(sign*y1(x));
 	if(ix==0x7ff00000) return zero;
 	if(ix>=0x52D00000) { /* x > 2**302 */
     /* (x >> n**2)
@@ -259,8 +259,8 @@ __ieee754_yn(int n, double x)
 		b = invsqrtpi*temp/sqrt(x);
 	} else {
 	    u_int32_t high;
-	    a = __ieee754_y0(x);
-	    b = __ieee754_y1(x);
+	    a = y0(x);
+	    b = y1(x);
 	/* quit if b is -inf */
 	    GET_HIGH_WORD(high,b);
 	    for(i=1;i<n&&high!=0xfff00000;i++){
diff --git a/libm/upstream-freebsd/lib/msun/src/e_jnf.c b/libm/upstream-freebsd/lib/msun/src/e_jnf.c
index 965feeb66..ba586226d 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_jnf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_jnf.c
@@ -32,7 +32,7 @@ one   =  1.0000000000e+00; /* 0x3F800000 */
 static const float zero  =  0.0000000000e+00;
 
 float
-__ieee754_jnf(int n, float x)
+jnf(int n, float x)
 {
 	int32_t i,hx,ix, sgn;
 	float a, b, temp, di;
@@ -50,16 +50,16 @@ __ieee754_jnf(int n, float x)
 		x = -x;
 		hx ^= 0x80000000;
 	}
-	if(n==0) return(__ieee754_j0f(x));
-	if(n==1) return(__ieee754_j1f(x));
+	if(n==0) return(j0f(x));
+	if(n==1) return(j1f(x));
 	sgn = (n&1)&(hx>>31);	/* even n -- 0, odd n -- sign(x) */
 	x = fabsf(x);
 	if(ix==0||ix>=0x7f800000) 	/* if x is 0 or inf */
 	    b = zero;
 	else if((float)n<=x) {
 		/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
-	    a = __ieee754_j0f(x);
-	    b = __ieee754_j1f(x);
+	    a = j0f(x);
+	    b = j1f(x);
 	    for(i=1;i<n;i++){
 		temp = b;
 		b = b*((float)(i+i)/x) - a; /* avoid underflow */
@@ -134,7 +134,7 @@ __ieee754_jnf(int n, float x)
 		 */
 		tmp = n;
 		v = two/x;
-		tmp = tmp*__ieee754_logf(fabsf(v*tmp));
+		tmp = tmp*logf(fabsf(v*tmp));
 		if(tmp<(float)8.8721679688e+01) {
 	    	    for(i=n-1,di=(float)(i+i);i>0;i--){
 		        temp = b;
@@ -158,8 +158,8 @@ __ieee754_jnf(int n, float x)
 			}
 	     	    }
 		}
-		z = __ieee754_j0f(x);
-		w = __ieee754_j1f(x);
+		z = j0f(x);
+		w = j1f(x);
 		if (fabsf(z) >= fabsf(w))
 		    b = (t*z/b);
 		else
@@ -170,7 +170,7 @@ __ieee754_jnf(int n, float x)
 }
 
 float
-__ieee754_ynf(int n, float x)
+ynf(int n, float x)
 {
 	int32_t i,hx,ix,ib;
 	int32_t sign;
@@ -186,12 +186,12 @@ __ieee754_ynf(int n, float x)
 		n = -n;
 		sign = 1 - ((n&1)<<1);
 	}
-	if(n==0) return(__ieee754_y0f(x));
-	if(n==1) return(sign*__ieee754_y1f(x));
+	if(n==0) return(y0f(x));
+	if(n==1) return(sign*y1f(x));
 	if(ix==0x7f800000) return zero;
 
-	a = __ieee754_y0f(x);
-	b = __ieee754_y1f(x);
+	a = y0f(x);
+	b = y1f(x);
 	/* quit if b is -inf */
 	GET_FLOAT_WORD(ib,b);
 	for(i=1;i<n&&ib!=0xff800000;i++){
diff --git a/libm/upstream-freebsd/lib/msun/src/e_lgamma.c b/libm/upstream-freebsd/lib/msun/src/e_lgamma.c
index 43f517560..9c4a30e6e 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_lgamma.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_lgamma.c
@@ -15,10 +15,10 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_lgamma(x)
+/* lgamma(x)
  * Return the logarithm of the Gamma function of x.
  *
- * Method: call __ieee754_lgamma_r
+ * Method: call lgamma_r
  */
 
 #include <float.h>
@@ -29,9 +29,9 @@ __FBSDID("$FreeBSD$");
 extern int signgam;
 
 double
-__ieee754_lgamma(double x)
+lgamma(double x)
 {
-	return __ieee754_lgamma_r(x,&signgam);
+	return lgamma_r(x,&signgam);
 }
 
 #if (LDBL_MANT_DIG == 53)
diff --git a/libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c b/libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c
index 48da493fe..c020b638e 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c
@@ -13,7 +13,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_lgamma_r(x, signgamp)
+/* lgamma_r(x, signgamp)
  * Reentrant version of the logarithm of the Gamma function
  * with user provide pointer for the sign of Gamma(x).
  *
@@ -199,7 +199,7 @@ sin_pi(double x)
 
 
 double
-__ieee754_lgamma_r(double x, int *signgamp)
+lgamma_r(double x, int *signgamp)
 {
 	double nadj,p,p1,p2,p3,q,r,t,w,y,z;
 	int32_t hx;
@@ -217,7 +217,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
 	if(ix<0x3c700000) {	/* |x|<2**-56, return -log(|x|) */
 	    if((ix|lx)==0)
 	        return one/vzero;
-	    return -__ieee754_log(fabs(x));
+	    return -log(fabs(x));
 	}
 
     /* purge negative integers and start evaluation for other x < 0 */
@@ -227,7 +227,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
 		return one/vzero;
 	    t = sin_pi(x);
 	    if(t==zero) return one/vzero; /* -integer */
-	    nadj = __ieee754_log(pi/fabs(t*x));
+	    nadj = log(pi/fabs(t*x));
 	    if(t<zero) *signgamp = -1;
 	    x = -x;
 	}
@@ -237,7 +237,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
     /* for x < 2.0 */
 	else if(ix<0x40000000) {
 	    if(ix<=0x3feccccc) { 	/* lgamma(x) = lgamma(x+1)-log(x) */
-		r = -__ieee754_log(x);
+		r = -log(x);
 		if(ix>=0x3FE76944) {y = one-x; i= 0;}
 		else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;}
 	  	else {y = x; i=2;}
@@ -282,18 +282,18 @@ __ieee754_lgamma_r(double x, int *signgamp)
 	    case 5: z *= (y+4);		/* FALLTHRU */
 	    case 4: z *= (y+3);		/* FALLTHRU */
 	    case 3: z *= (y+2);		/* FALLTHRU */
-		    r += __ieee754_log(z); break;
+		    r += log(z); break;
 	    }
     /* 8.0 <= x < 2**56 */
 	} else if (ix < 0x43700000) {
-	    t = __ieee754_log(x);
+	    t = log(x);
 	    z = one/x;
 	    y = z*z;
 	    w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
 	    r = (x-half)*(t-one)+w;
 	} else
     /* 2**56 <= x <= inf */
-	    r =  x*(__ieee754_log(x)-one);
+	    r =  x*(log(x)-one);
 	if(hx<0) r = nadj - r;
 	return r;
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/e_lgammaf.c b/libm/upstream-freebsd/lib/msun/src/e_lgammaf.c
index 1e2c55273..00a816c35 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_lgammaf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_lgammaf.c
@@ -16,10 +16,10 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_lgammaf(x)
+/* lgammaf(x)
  * Return the logarithm of the Gamma function of x.
  *
- * Method: call __ieee754_lgammaf_r
+ * Method: call lgammaf_r
  */
 
 #include "math.h"
@@ -28,7 +28,7 @@ __FBSDID("$FreeBSD$");
 extern int signgam;
 
 float
-__ieee754_lgammaf(float x)
+lgammaf(float x)
 {
-	return __ieee754_lgammaf_r(x,&signgam);
+	return lgammaf_r(x,&signgam);
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c b/libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c
index 48346c336..fdd23218e 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c
@@ -120,7 +120,7 @@ sin_pif(float x)
 
 
 float
-__ieee754_lgammaf_r(float x, int *signgamp)
+lgammaf_r(float x, int *signgamp)
 {
 	float nadj,p,p1,p2,q,r,t,w,y,z;
 	int32_t hx;
@@ -138,7 +138,7 @@ __ieee754_lgammaf_r(float x, int *signgamp)
 	if(ix<0x32000000) {		/* |x|<2**-27, return -log(|x|) */
 	    if(ix==0)
 	        return one/vzero;
-	    return -__ieee754_logf(fabsf(x));
+	    return -logf(fabsf(x));
 	}
 
     /* purge negative integers and start evaluation for other x < 0 */
@@ -148,7 +148,7 @@ __ieee754_lgammaf_r(float x, int *signgamp)
 		return one/vzero;
 	    t = sin_pif(x);
 	    if(t==zero) return one/vzero; /* -integer */
-	    nadj = __ieee754_logf(pi/fabsf(t*x));
+	    nadj = logf(pi/fabsf(t*x));
 	    if(t<zero) *signgamp = -1;
 	    x = -x;
 	}
@@ -158,7 +158,7 @@ __ieee754_lgammaf_r(float x, int *signgamp)
     /* for x < 2.0 */
 	else if(ix<0x40000000) {
 	    if(ix<=0x3f666666) { 	/* lgamma(x) = lgamma(x+1)-log(x) */
-		r = -__ieee754_logf(x);
+		r = -logf(x);
 		if(ix>=0x3f3b4a20) {y = one-x; i= 0;}
 		else if(ix>=0x3e6d3308) {y= x-(tc-one); i=1;}
 	  	else {y = x; i=2;}
@@ -198,18 +198,18 @@ __ieee754_lgammaf_r(float x, int *signgamp)
 	    case 5: z *= (y+4);		/* FALLTHRU */
 	    case 4: z *= (y+3);		/* FALLTHRU */
 	    case 3: z *= (y+2);		/* FALLTHRU */
-		    r += __ieee754_logf(z); break;
+		    r += logf(z); break;
 	    }
     /* 8.0 <= x < 2**27 */
 	} else if (ix < 0x4d000000) {
-	    t = __ieee754_logf(x);
+	    t = logf(x);
 	    z = one/x;
 	    y = z*z;
 	    w = w0+z*(w1+y*w2);
 	    r = (x-half)*(t-one)+w;
 	} else
     /* 2**27 <= x <= inf */
-	    r =  x*(__ieee754_logf(x)-one);
+	    r =  x*(logf(x)-one);
 	if(hx<0) r = nadj - r;
 	return r;
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/e_log.c b/libm/upstream-freebsd/lib/msun/src/e_log.c
index 68bc1070b..03ce82061 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_log.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_log.c
@@ -14,7 +14,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_log(x)
+/* log(x)
  * Return the logrithm of x
  *
  * Method :                  
@@ -86,7 +86,7 @@ static const double zero   =  0.0;
 static volatile double vzero = 0.0;
 
 double
-__ieee754_log(double x)
+log(double x)
 {
 	double hfsq,f,s,z,R,w,t1,t2,dk;
 	int32_t k,hx,i,j;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_log10.c b/libm/upstream-freebsd/lib/msun/src/e_log10.c
index 3c89ed2d8..595c23808 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_log10.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_log10.c
@@ -39,7 +39,7 @@ static const double zero   =  0.0;
 static volatile double vzero = 0.0;
 
 double
-__ieee754_log10(double x)
+log10(double x)
 {
 	double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2;
 	int32_t i,k,hx;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_log10f.c b/libm/upstream-freebsd/lib/msun/src/e_log10f.c
index 9856df2e7..d0c3a53c8 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_log10f.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_log10f.c
@@ -31,7 +31,7 @@ static const float zero   =  0.0;
 static volatile float vzero = 0.0;
 
 float
-__ieee754_log10f(float x)
+log10f(float x)
 {
 	float f,hfsq,hi,lo,r,y;
 	int32_t i,k,hx;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_log2.c b/libm/upstream-freebsd/lib/msun/src/e_log2.c
index 4766cdb49..10b1c00d5 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_log2.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_log2.c
@@ -39,7 +39,7 @@ static const double zero   =  0.0;
 static volatile double vzero = 0.0;
 
 double
-__ieee754_log2(double x)
+log2(double x)
 {
 	double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
 	int32_t i,k,hx;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_log2f.c b/libm/upstream-freebsd/lib/msun/src/e_log2f.c
index 1794484e6..956f33a34 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_log2f.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_log2f.c
@@ -29,7 +29,7 @@ static const float zero   =  0.0;
 static volatile float vzero = 0.0;
 
 float
-__ieee754_log2f(float x)
+log2f(float x)
 {
 	float f,hfsq,hi,lo,r,y;
 	int32_t i,k,hx;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_logf.c b/libm/upstream-freebsd/lib/msun/src/e_logf.c
index ec3985fcb..68a4d5d88 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_logf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_logf.c
@@ -33,7 +33,7 @@ static const float zero   =  0.0;
 static volatile float vzero = 0.0;
 
 float
-__ieee754_logf(float x)
+logf(float x)
 {
 	float hfsq,f,s,z,R,w,t1,t2,dk;
 	int32_t k,ix,i,j;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_pow.c b/libm/upstream-freebsd/lib/msun/src/e_pow.c
index 52b2f5ceb..adc64c998 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_pow.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_pow.c
@@ -12,7 +12,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_pow(x,y) return x**y
+/* pow(x,y) return x**y
  *
  *		      n
  * Method:  Let x =  2   * (1+f)
@@ -98,7 +98,7 @@ ivln2_h  =  1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
 ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 
 double
-__ieee754_pow(double x, double y)
+pow(double x, double y)
 {
 	double z,ax,z_h,z_l,p_h,p_l;
 	double y1,t1,t2,r,s,t,u,v,w;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_powf.c b/libm/upstream-freebsd/lib/msun/src/e_powf.c
index 122da455f..f5a2c70c7 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_powf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_powf.c
@@ -56,7 +56,7 @@ ivln2_h  =  1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
 ivln2_l  =  7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
 
 float
-__ieee754_powf(float x, float y)
+powf(float x, float y)
 {
 	float z,ax,z_h,z_l,p_h,p_l;
 	float y1,t1,t2,r,s,sn,t,u,v,w;
@@ -108,7 +108,7 @@ __ieee754_powf(float x, float y)
 	if(hy==0x40000000) return x*x; /* y is  2 */
 	if(hy==0x3f000000) {	/* y is  0.5 */
 	    if(hx>=0)	/* x >= +0 */
-	    return __ieee754_sqrtf(x);
+	    return sqrtf(x);
 	}
 
 	ax   = fabsf(x);
diff --git a/libm/upstream-freebsd/lib/msun/src/e_remainder.c b/libm/upstream-freebsd/lib/msun/src/e_remainder.c
index a4ae0b780..13156d8cb 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_remainder.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_remainder.c
@@ -14,7 +14,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_remainder(x,p)
+/* remainder(x,p)
  * Return :                  
  * 	returns  x REM p  =  x - [x/p]*p as if in infinite 
  * 	precise arithmetic, where [x/p] is the (infinite bit) 
@@ -32,7 +32,7 @@ static const double zero = 0.0;
 
 
 double
-__ieee754_remainder(double x, double p)
+remainder(double x, double p)
 {
 	int32_t hx,hp;
 	u_int32_t sx,lx,lp;
@@ -52,7 +52,7 @@ __ieee754_remainder(double x, double p)
 	    return nan_mix_op(x, p, *)/nan_mix_op(x, p, *);
 
 
-	if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p);	/* now x < 2p */
+	if (hp<=0x7fdfffff) x = fmod(x,p+p);	/* now x < 2p */
 	if (((hx-hp)|(lx-lp))==0) return zero*x;
 	x  = fabs(x);
 	p  = fabs(p);
diff --git a/libm/upstream-freebsd/lib/msun/src/e_remainderf.c b/libm/upstream-freebsd/lib/msun/src/e_remainderf.c
index 8004493de..e0dcfd181 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_remainderf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_remainderf.c
@@ -23,7 +23,7 @@ static const float zero = 0.0;
 
 
 float
-__ieee754_remainderf(float x, float p)
+remainderf(float x, float p)
 {
 	int32_t hx,hp;
 	u_int32_t sx;
@@ -42,7 +42,7 @@ __ieee754_remainderf(float x, float p)
 	    return nan_mix_op(x, p, *)/nan_mix_op(x, p, *);
 
 
-	if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p);	/* now x < 2p */
+	if (hp<=0x7effffff) x = fmodf(x,p+p);	/* now x < 2p */
 	if ((hx-hp)==0) return zero*x;
 	x  = fabsf(x);
 	p  = fabsf(p);
diff --git a/libm/upstream-freebsd/lib/msun/src/e_scalb.c b/libm/upstream-freebsd/lib/msun/src/e_scalb.c
index c0a7b5b75..84a68939b 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_scalb.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_scalb.c
@@ -15,7 +15,7 @@
 __FBSDID("$FreeBSD$");
 
 /*
- * __ieee754_scalb(x, fn) is provide for
+ * scalb(x, fn) is provide for
  * passing various standard test suite. One 
  * should use scalbn() instead.
  */
@@ -25,10 +25,10 @@ __FBSDID("$FreeBSD$");
 
 #ifdef _SCALB_INT
 double
-__ieee754_scalb(double x, int fn)
+scalb(double x, int fn)
 #else
 double
-__ieee754_scalb(double x, double fn)
+scalb(double x, double fn)
 #endif
 {
 #ifdef _SCALB_INT
diff --git a/libm/upstream-freebsd/lib/msun/src/e_scalbf.c b/libm/upstream-freebsd/lib/msun/src/e_scalbf.c
index d49e9041f..28483a5e2 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_scalbf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_scalbf.c
@@ -21,10 +21,10 @@ __FBSDID("$FreeBSD$");
 
 #ifdef _SCALB_INT
 float
-__ieee754_scalbf(float x, int fn)
+scalbf(float x, int fn)
 #else
 float
-__ieee754_scalbf(float x, float fn)
+scalbf(float x, float fn)
 #endif
 {
 #ifdef _SCALB_INT
diff --git a/libm/upstream-freebsd/lib/msun/src/e_sinh.c b/libm/upstream-freebsd/lib/msun/src/e_sinh.c
index 6c01f4a3c..9fe89996a 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_sinh.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_sinh.c
@@ -14,7 +14,7 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/* __ieee754_sinh(x)
+/* sinh(x)
  * Method : 
  * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
  *	1. Replace x by |x| (sinh(-x) = -sinh(x)). 
@@ -40,7 +40,7 @@ __FBSDID("$FreeBSD$");
 static const double one = 1.0, shuge = 1.0e307;
 
 double
-__ieee754_sinh(double x)
+sinh(double x)
 {
 	double t,h;
 	int32_t ix,jx;
@@ -64,7 +64,7 @@ __ieee754_sinh(double x)
 	}
 
     /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
-	if (ix < 0x40862E42)  return h*__ieee754_exp(fabs(x));
+	if (ix < 0x40862E42)  return h*exp(fabs(x));
 
     /* |x| in [log(maxdouble), overflowthresold] */
 	if (ix<=0x408633CE)
diff --git a/libm/upstream-freebsd/lib/msun/src/e_sinhf.c b/libm/upstream-freebsd/lib/msun/src/e_sinhf.c
index 1be2dc397..082beb14f 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_sinhf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_sinhf.c
@@ -22,7 +22,7 @@ __FBSDID("$FreeBSD$");
 static const float one = 1.0, shuge = 1.0e37;
 
 float
-__ieee754_sinhf(float x)
+sinhf(float x)
 {
 	float t,h;
 	int32_t ix,jx;
@@ -45,7 +45,7 @@ __ieee754_sinhf(float x)
 	}
 
     /* |x| in [9, logf(maxfloat)] return 0.5*exp(|x|) */
-	if (ix < 0x42b17217)  return h*__ieee754_expf(fabsf(x));
+	if (ix < 0x42b17217)  return h*expf(fabsf(x));
 
     /* |x| in [logf(maxfloat), overflowthresold] */
 	if (ix<=0x42b2d4fc)
diff --git a/libm/upstream-freebsd/lib/msun/src/math_private.h b/libm/upstream-freebsd/lib/msun/src/math_private.h
index df526e71e..a55d97a78 100644
--- a/libm/upstream-freebsd/lib/msun/src/math_private.h
+++ b/libm/upstream-freebsd/lib/msun/src/math_private.h
@@ -624,7 +624,7 @@ rnintf(__float_t x)
  * The complications for extra precision are smaller for rnintl() since it
  * can safely assume that the rounding precision has been increased from
  * its default to FP_PE on x86.  We don't exploit that here to get small
- * optimizations from limiting the rangle to double.  We just need it for
+ * optimizations from limiting the range to double.  We just need it for
  * the magic number to work with long doubles.  ld128 callers should use
  * rnint() instead of this if possible.  ld80 callers should prefer
  * rnintl() since for amd64 this avoids swapping the register set, while
@@ -688,6 +688,59 @@ irintl(long double x)
 }
 #endif
 
+/*
+ * The following are fast floor macros for 0 <= |x| < 0x1p(N-1), where
+ * N is the precision of the type of x. These macros are used in the
+ * half-cycle trignometric functions (e.g., sinpi(x)).
+ */
+#define	FFLOORF(x, j0, ix) do {			\
+	(j0) = (((ix) >> 23) & 0xff) - 0x7f;	\
+	(ix) &= ~(0x007fffff >> (j0));		\
+	SET_FLOAT_WORD((x), (ix));		\
+} while (0)
+
+#define	FFLOOR(x, j0, ix, lx) do {				\
+	(j0) = (((ix) >> 20) & 0x7ff) - 0x3ff;			\
+	if ((j0) < 20) {					\
+		(ix) &= ~(0x000fffff >> (j0));			\
+		(lx) = 0;					\
+	} else {						\
+		(lx) &= ~((uint32_t)0xffffffff >> ((j0) - 20));	\
+	}							\
+	INSERT_WORDS((x), (ix), (lx));				\
+} while (0)
+
+#define	FFLOORL80(x, j0, ix, lx) do {			\
+	j0 = ix - 0x3fff + 1;				\
+	if ((j0) < 32) {				\
+		(lx) = ((lx) >> 32) << 32;		\
+		(lx) &= ~((((lx) << 32)-1) >> (j0));	\
+	} else {					\
+		uint64_t _m;				\
+		_m = (uint64_t)-1 >> (j0);		\
+		if ((lx) & _m) (lx) &= ~_m;		\
+	}						\
+	INSERT_LDBL80_WORDS((x), (ix), (lx));		\
+} while (0)
+
+#define FFLOORL128(x, ai, ar) do {			\
+	union IEEEl2bits u;				\
+	uint64_t m;					\
+	int e;						\
+	u.e = (x);					\
+	e = u.bits.exp - 16383;				\
+	if (e < 48) {					\
+		m = ((1llu << 49) - 1) >> (e + 1);	\
+		u.bits.manh &= ~m;			\
+		u.bits.manl = 0;			\
+	} else {					\
+		m = (uint64_t)-1 >> (e - 48);		\
+		u.bits.manl &= ~m;			\
+	}						\
+	(ai) = u.e;					\
+	(ar) = (x) - (ai);				\
+} while (0)
+
 #ifdef DEBUG
 #if defined(__amd64__) || defined(__i386__)
 #define	breakpoint()	asm("int $3")
@@ -698,191 +751,25 @@ irintl(long double x)
 #endif
 #endif
 
-/* Write a pari script to test things externally. */
-#ifdef DOPRINT
-#include <stdio.h>
-
-#ifndef DOPRINT_SWIZZLE
-#define	DOPRINT_SWIZZLE		0
-#endif
-
-#ifdef DOPRINT_LD80
-
-#define	DOPRINT_START(xp) do {						\
-	uint64_t __lx;							\
-	uint16_t __hx;							\
-									\
-	/* Hack to give more-problematic args. */			\
-	EXTRACT_LDBL80_WORDS(__hx, __lx, *xp);				\
-	__lx ^= DOPRINT_SWIZZLE;					\
-	INSERT_LDBL80_WORDS(*xp, __hx, __lx);				\
-	printf("x = %.21Lg; ", (long double)*xp);			\
-} while (0)
-#define	DOPRINT_END1(v)							\
-	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
-#define	DOPRINT_END2(hi, lo)						\
-	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
-	    (long double)(hi), (long double)(lo))
-
-#elif defined(DOPRINT_D64)
-
-#define	DOPRINT_START(xp) do {						\
-	uint32_t __hx, __lx;						\
-									\
-	EXTRACT_WORDS(__hx, __lx, *xp);					\
-	__lx ^= DOPRINT_SWIZZLE;					\
-	INSERT_WORDS(*xp, __hx, __lx);					\
-	printf("x = %.21Lg; ", (long double)*xp);			\
-} while (0)
-#define	DOPRINT_END1(v)							\
-	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
-#define	DOPRINT_END2(hi, lo)						\
-	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
-	    (long double)(hi), (long double)(lo))
-
-#elif defined(DOPRINT_F32)
-
-#define	DOPRINT_START(xp) do {						\
-	uint32_t __hx;							\
-									\
-	GET_FLOAT_WORD(__hx, *xp);					\
-	__hx ^= DOPRINT_SWIZZLE;					\
-	SET_FLOAT_WORD(*xp, __hx);					\
-	printf("x = %.21Lg; ", (long double)*xp);			\
-} while (0)
-#define	DOPRINT_END1(v)							\
-	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
-#define	DOPRINT_END2(hi, lo)						\
-	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
-	    (long double)(hi), (long double)(lo))
-
-#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
-
-#ifndef DOPRINT_SWIZZLE_HIGH
-#define	DOPRINT_SWIZZLE_HIGH	0
-#endif
-
-#define	DOPRINT_START(xp) do {						\
-	uint64_t __lx, __llx;						\
-	uint16_t __hx;							\
-									\
-	EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp);			\
-	__llx ^= DOPRINT_SWIZZLE;					\
-	__lx ^= DOPRINT_SWIZZLE_HIGH;					\
-	INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx);			\
-	printf("x = %.36Lg; ", (long double)*xp);					\
-} while (0)
-#define	DOPRINT_END1(v)							\
-	printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
-#define	DOPRINT_END2(hi, lo)						\
-	printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n",		\
-	    (long double)(hi), (long double)(lo))
-
-#endif /* DOPRINT_LD80 */
-
-#else /* !DOPRINT */
-#define	DOPRINT_START(xp)
-#define	DOPRINT_END1(v)
-#define	DOPRINT_END2(hi, lo)
-#endif /* DOPRINT */
-
-#define	RETURNP(x) do {			\
-	DOPRINT_END1(x);		\
-	RETURNF(x);			\
-} while (0)
-#define	RETURNPI(x) do {		\
-	DOPRINT_END1(x);		\
-	RETURNI(x);			\
-} while (0)
-#define	RETURN2P(x, y) do {		\
-	DOPRINT_END2((x), (y));		\
-	RETURNF((x) + (y));		\
-} while (0)
-#define	RETURN2PI(x, y) do {		\
-	DOPRINT_END2((x), (y));		\
-	RETURNI((x) + (y));		\
-} while (0)
 #ifdef STRUCT_RETURN
 #define	RETURNSP(rp) do {		\
 	if (!(rp)->lo_set)		\
-		RETURNP((rp)->hi);	\
-	RETURN2P((rp)->hi, (rp)->lo);	\
+		RETURNF((rp)->hi);	\
+	RETURNF((rp)->hi + (rp)->lo);	\
 } while (0)
 #define	RETURNSPI(rp) do {		\
 	if (!(rp)->lo_set)		\
-		RETURNPI((rp)->hi);	\
-	RETURN2PI((rp)->hi, (rp)->lo);	\
+		RETURNI((rp)->hi);	\
+	RETURNI((rp)->hi + (rp)->lo);	\
 } while (0)
 #endif
+
 #define	SUM2P(x, y) ({			\
 	const __typeof (x) __x = (x);	\
 	const __typeof (y) __y = (y);	\
-					\
-	DOPRINT_END2(__x, __y);		\
 	__x + __y;			\
 })
 
-/*
- * ieee style elementary functions
- *
- * We rename functions here to improve other sources' diffability
- * against fdlibm.
- */
-#define	__ieee754_sqrt	sqrt
-#define	__ieee754_acos	acos
-#define	__ieee754_acosh	acosh
-#define	__ieee754_log	log
-#define	__ieee754_log2	log2
-#define	__ieee754_atanh	atanh
-#define	__ieee754_asin	asin
-#define	__ieee754_atan2	atan2
-#define	__ieee754_exp	exp
-#define	__ieee754_cosh	cosh
-#define	__ieee754_fmod	fmod
-#define	__ieee754_pow	pow
-#define	__ieee754_lgamma lgamma
-#define	__ieee754_gamma	gamma
-#define	__ieee754_lgamma_r lgamma_r
-#define	__ieee754_gamma_r gamma_r
-#define	__ieee754_log10	log10
-#define	__ieee754_sinh	sinh
-#define	__ieee754_hypot	hypot
-#define	__ieee754_j0	j0
-#define	__ieee754_j1	j1
-#define	__ieee754_y0	y0
-#define	__ieee754_y1	y1
-#define	__ieee754_jn	jn
-#define	__ieee754_yn	yn
-#define	__ieee754_remainder remainder
-#define	__ieee754_scalb	scalb
-#define	__ieee754_sqrtf	sqrtf
-#define	__ieee754_acosf	acosf
-#define	__ieee754_acoshf acoshf
-#define	__ieee754_logf	logf
-#define	__ieee754_atanhf atanhf
-#define	__ieee754_asinf	asinf
-#define	__ieee754_atan2f atan2f
-#define	__ieee754_expf	expf
-#define	__ieee754_coshf	coshf
-#define	__ieee754_fmodf	fmodf
-#define	__ieee754_powf	powf
-#define	__ieee754_lgammaf lgammaf
-#define	__ieee754_gammaf gammaf
-#define	__ieee754_lgammaf_r lgammaf_r
-#define	__ieee754_gammaf_r gammaf_r
-#define	__ieee754_log10f log10f
-#define	__ieee754_log2f log2f
-#define	__ieee754_sinhf	sinhf
-#define	__ieee754_hypotf hypotf
-#define	__ieee754_j0f	j0f
-#define	__ieee754_j1f	j1f
-#define	__ieee754_y0f	y0f
-#define	__ieee754_y1f	y1f
-#define	__ieee754_jnf	jnf
-#define	__ieee754_ynf	ynf
-#define	__ieee754_remainderf remainderf
-#define	__ieee754_scalbf scalbf
-
 /* fdlibm kernel function */
 int	__kernel_rem_pio2(double*,double*,int,int,int);
 
diff --git a/libm/upstream-freebsd/lib/msun/src/s_asinh.c b/libm/upstream-freebsd/lib/msun/src/s_asinh.c
index cbb3d463a..a1b9169be 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_asinh.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_asinh.c
@@ -46,13 +46,13 @@ asinh(double x)
 	    if(huge+x>one) return x;	/* return x inexact except 0 */
 	}
 	if(ix>0x41b00000) {	/* |x| > 2**28 */
-	    w = __ieee754_log(fabs(x))+ln2;
+	    w = log(fabs(x))+ln2;
 	} else if (ix>0x40000000) {	/* 2**28 > |x| > 2.0 */
 	    t = fabs(x);
-	    w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t));
+	    w = log(2.0*t+one/(sqrt(x*x+one)+t));
 	} else {		/* 2.0 > |x| > 2**-28 */
 	    t = x*x;
-	    w =log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t)));
+	    w =log1p(fabs(x)+t/(one+sqrt(one+t)));
 	}
 	if(hx>0) return w; else return -w;
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/s_asinhf.c b/libm/upstream-freebsd/lib/msun/src/s_asinhf.c
index c1620dd58..72bcefed9 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_asinhf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_asinhf.c
@@ -36,13 +36,13 @@ asinhf(float x)
 	    if(huge+x>one) return x;	/* return x inexact except 0 */
 	}
 	if(ix>0x4d800000) {	/* |x| > 2**28 */
-	    w = __ieee754_logf(fabsf(x))+ln2;
+	    w = logf(fabsf(x))+ln2;
 	} else if (ix>0x40000000) {	/* 2**28 > |x| > 2.0 */
 	    t = fabsf(x);
-	    w = __ieee754_logf((float)2.0*t+one/(__ieee754_sqrtf(x*x+one)+t));
+	    w = logf((float)2.0*t+one/(sqrtf(x*x+one)+t));
 	} else {		/* 2.0 > |x| > 2**-28 */
 	    t = x*x;
-	    w =log1pf(fabsf(x)+t/(one+__ieee754_sqrtf(one+t)));
+	    w =log1pf(fabsf(x)+t/(one+sqrtf(one+t)));
 	}
 	if(hx>0) return w; else return -w;
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/s_cospi.c b/libm/upstream-freebsd/lib/msun/src/s_cospi.c
index 2e2f92733..f97570dc8 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_cospi.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_cospi.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017 Steven G. Kargl
+ * Copyright (c) 2017, 2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -104,20 +104,10 @@ cospi(double x)
 	}
 
 	if (ix < 0x43300000) {		/* 1 <= |x| < 0x1p52 */
-		/* Determine integer part of ax. */
-		j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
-		if (j0 < 20) {
-			ix &= ~(0x000fffff >> j0);
-			lx = 0;
-		} else {
-			lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
-		}
-		INSERT_WORDS(x, ix, lx);
-
+		FFLOOR(x, j0, ix, lx);	/* Integer part of ax. */
 		ax -= x;
 		EXTRACT_WORDS(ix, lx, ax);
 
-
 		if (ix < 0x3fe00000) {		/* |x| < 0.5 */
 			if (ix < 0x3fd00000)	/* |x| < 0.25 */
 				c = ix == 0 ? 1 : __kernel_cospi(ax);
@@ -143,9 +133,11 @@ cospi(double x)
 		return (vzero / vzero);
 
 	/*
-	 * |x| >= 0x1p52 is always an even integer, so return 1.
+	 * For 0x1p52 <= |x| < 0x1p53 need to determine if x is an even
+	 * or odd integer to return +1 or -1.
+	 * For |x| >= 0x1p53, it is always an even integer, so return 1.
 	 */
-	return (1);
+	return (ix < 0x43400000 ? ((lx & 1) ? -1 : 1) : 1);
 }
 
 #if LDBL_MANT_DIG == 53
diff --git a/libm/upstream-freebsd/lib/msun/src/s_erf.c b/libm/upstream-freebsd/lib/msun/src/s_erf.c
index 5f228474f..ab2dc1905 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_erf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_erf.c
@@ -238,7 +238,7 @@ erf(double x)
 	}
 	z  = x;
 	SET_LOW_WORD(z,0);
-	r  =  __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
+	r  =  exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
 	if(hx>=0) return one-r/x; else return  r/x-one;
 }
 
@@ -297,7 +297,7 @@ erfc(double x)
 	    }
 	    z  = x;
 	    SET_LOW_WORD(z,0);
-	    r  =  __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
+	    r  =  exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S);
 	    if(hx>0) return r/x; else return two-r/x;
 	} else {
 	    if(hx>0) return tiny*tiny; else return two-tiny;
diff --git a/libm/upstream-freebsd/lib/msun/src/s_significand.c b/libm/upstream-freebsd/lib/msun/src/s_significand.c
index 356e3001f..eed80ece1 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_significand.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_significand.c
@@ -25,5 +25,5 @@ __FBSDID("$FreeBSD$");
 double
 significand(double x)
 {
-	return __ieee754_scalb(x,(double) -ilogb(x));
+	return scalb(x,(double) -ilogb(x));
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/s_significandf.c b/libm/upstream-freebsd/lib/msun/src/s_significandf.c
index ad030e239..b33ab7b6f 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_significandf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_significandf.c
@@ -22,5 +22,5 @@ __FBSDID("$FreeBSD$");
 float
 significandf(float x)
 {
-	return __ieee754_scalbf(x,(float) -ilogbf(x));
+	return scalbf(x,(float) -ilogbf(x));
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/s_sinpi.c b/libm/upstream-freebsd/lib/msun/src/s_sinpi.c
index bc3759e56..8b388de86 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_sinpi.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_sinpi.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2017 Steven G. Kargl
+ * Copyright (c) 2017, 2023 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -118,16 +118,7 @@ sinpi(double x)
 	}
 
 	if (ix < 0x43300000) {			/* 1 <= |x| < 0x1p52 */
-		/* Determine integer part of ax. */
-		j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
-		if (j0 < 20) {
-			ix &= ~(0x000fffff >> j0);
-			lx = 0;
-		} else {
-			lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
-		}
-		INSERT_WORDS(x, ix, lx);
-
+		FFLOOR(x, j0, ix, lx);	/* Integer part of ax. */
 		ax -= x;
 		EXTRACT_WORDS(ix, lx, ax);