Remove the unused leftovers of the 4.4BSD libm, which was only used

on non-IEEE platforms.  Since the VAX port was discontinued, all
our remaining architectures use IEEE floating point, as will any
future ones.
ok millert@ tb@

34 files changed, 1 insertions, 6185 deletions

diff --git a/lib/libm/Makefile b/lib/libm/Makefile
index 2921d8237ee..2031f8b0df3 100644
--- a/lib/libm/Makefile
+++ b/lib/libm/Makefile
@@ -1,4 +1,4 @@
-#  $OpenBSD: Makefile,v 1.118 2018/03/10 20:52:58 kettenis Exp $
+#  $OpenBSD: Makefile,v 1.119 2018/07/18 20:21:12 naddy Exp $
 #  $NetBSD: Makefile,v 1.28 1995/11/20 22:06:19 jtc Exp $
 #
 #  @(#)Makefile 5.1beta 93/09/24 
@@ -67,7 +67,6 @@ ARCH_SRCS = e_sqrtl.c
 
 .PATH:	${.CURDIR}/man
 .PATH:	${.CURDIR}/src
-.PATH:	${.CURDIR}/noieee_src
 
 LIB=	m
 COMMON_SRCS = b_exp__D.c b_log__D.c b_tgamma.c \
diff --git a/lib/libm/noieee_src/mathimpl.h b/lib/libm/noieee_src/mathimpl.h
deleted file mode 100644
index 163843ab186..00000000000
--- a/lib/libm/noieee_src/mathimpl.h
+++ /dev/null
@@ -1,45 +0,0 @@
-/*	$OpenBSD: mathimpl.h,v 1.10 2009/04/11 20:03:21 martynas Exp $	*/
-/*	$NetBSD: mathimpl.h,v 1.1 1995/10/10 23:36:31 ragge Exp $	*/
-/*
- * Copyright (c) 1988, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- *	@(#)mathimpl.h	8.1 (Berkeley) 6/4/93
- */
-
-/*
- * Functions internal to the math package, yet not static.
- */
-extern double	__exp__E();
-extern double	__log__L();
-
-struct Double {
-	double a, b;
-};
-
-double __exp__D(double, double);
-struct Double __log__D(double);
diff --git a/lib/libm/noieee_src/n_acosh.c b/lib/libm/noieee_src/n_acosh.c
deleted file mode 100644
index 46b2914339f..00000000000
--- a/lib/libm/noieee_src/n_acosh.c
+++ /dev/null
@@ -1,88 +0,0 @@
-/*	$OpenBSD: n_acosh.c,v 1.9 2009/10/27 23:59:29 deraadt Exp $	*/
-/*	$NetBSD: n_acosh.c,v 1.1 1995/10/10 23:36:33 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* ACOSH(X)
- * RETURN THE INVERSE HYPERBOLIC COSINE OF X
- * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 2/16/85;
- * REVISED BY K.C. NG on 3/6/85, 3/24/85, 4/16/85, 8/17/85.
- *
- * Required system supported functions :
- *	sqrt(x)
- *
- * Required kernel function:
- *	log1p(x) 		...return log(1+x)
- *
- * Method :
- *	Based on
- *		acosh(x) = log [ x + sqrt(x*x-1) ]
- *	we have
- *		acosh(x) := log1p(x)+ln2,	if (x > 1.0E20); else
- *		acosh(x) := log1p( sqrt(x-1) * (sqrt(x-1) + sqrt(x+1)) ) .
- *	These formulae avoid the over/underflow complication.
- *
- * Special cases:
- *	acosh(x) is NaN with signal if x<1.
- *	acosh(NaN) is NaN without signal.
- *
- * Accuracy:
- *	acosh(x) returns the exact inverse hyperbolic cosine of x nearly
- *	rounded. In a test run with 512,000 random arguments on a VAX, the
- *	maximum observed error was 3.30 ulps (units of the last place) at
- *	x=1.0070493753568216 .
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-static const double ln2hi = 6.9314718055829871446E-1;
-static const double ln2lo = 1.6465949582897081279E-12;
-
-double
-acosh(double x)
-{
-	double t,big=1.E20; /* big+1==big */
-
-	if (isnan(x))
-		return (x);
-
-    /* return log1p(x) + log(2) if x is large */
-	if(x>big) {t=log1p(x)+ln2lo; return(t+ln2hi);}
-
-	t=sqrt(x-1.0);
-	return(log1p(t*(t+sqrt(x+1.0))));
-}
diff --git a/lib/libm/noieee_src/n_asincos.c b/lib/libm/noieee_src/n_asincos.c
deleted file mode 100644
index a6fe371d39a..00000000000
--- a/lib/libm/noieee_src/n_asincos.c
+++ /dev/null
@@ -1,173 +0,0 @@
-/*	$OpenBSD: n_asincos.c,v 1.15 2016/09/12 19:47:02 guenther Exp $	*/
-/*	$NetBSD: n_asincos.c,v 1.1 1995/10/10 23:36:34 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* ASIN(X)
- * RETURNS ARC SINE OF X
- * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
- * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
- *
- * Required system supported functions:
- *	copysign(x,y)
- *	sqrt(x)
- *
- * Required kernel function:
- *	atan2(y,x)
- *
- * Method:
- *	asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is
- *		  computed as follows
- *			1-x*x                     if x <  0.5,
- *			2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5.
- *
- * Special cases:
- *	if x is NaN, return x itself;
- *	if |x|>1, return NaN.
- *
- * Accuracy:
- * 1)  If atan2() uses machine PI, then
- *
- *	asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded;
- *	and PI is the exact pi rounded to machine precision (see atan2 for
- *      details):
- *
- *	in decimal:
- *		pi = 3.141592653589793 23846264338327 .....
- *    53 bits   PI = 3.141592653589793 115997963 ..... ,
- *    56 bits   PI = 3.141592653589793 227020265 ..... ,
- *
- *	in hexadecimal:
- *		pi = 3.243F6A8885A308D313198A2E....
- *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps
- *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
- *
- *	In a test run with more than 200,000 random arguments on a VAX, the
- *	maximum observed error in ulps (units in the last place) was
- *	2.06 ulps.      (comparing against (PI/pi)*(exact asin(x)));
- *
- * 2)  If atan2() uses true pi, then
- *
- *	asin(x) returns the exact asin(x) with error below about 2 ulps.
- *
- *	In a test run with more than 1,024,000 random arguments on a VAX, the
- *	maximum observed error in ulps (units in the last place) was
- *      1.99 ulps.
- */
-
-#include <math.h>
-
-#include "mathimpl.h"
-
-double
-asin(double x)
-{
-	double s, t, one = 1.0;
-
-	if (isnan(x))
-		return (x);
-
-	s=copysign(x,one);
-	if(s <= 0.5)
-	    return(atan2(x,sqrt(one-x*x)));
-	else
-	    { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); }
-
-}
-DEF_STD(asin);
-LDBL_CLONE(asin);
-
-/* ACOS(X)
- * RETURNS ARC COS OF X
- * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
- * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
- *
- * Required system supported functions:
- *	copysign(x,y)
- *	sqrt(x)
- *
- * Required kernel function:
- *	atan2(y,x)
- *
- * Method:
- *			      ________
- *                           / 1 - x
- *	acos(x) = 2*atan2(  / -------- , 1 ) .
- *                        \/   1 + x
- *
- * Special cases:
- *	if x is NaN, return x itself;
- *	if |x|>1, return NaN.
- *
- * Accuracy:
- * 1)  If atan2() uses machine PI, then
- *
- *	acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded;
- *	and PI is the exact pi rounded to machine precision (see atan2 for
- *      details):
- *
- *	in decimal:
- *		pi = 3.141592653589793 23846264338327 .....
- *    53 bits   PI = 3.141592653589793 115997963 ..... ,
- *    56 bits   PI = 3.141592653589793 227020265 ..... ,
- *
- *	in hexadecimal:
- *		pi = 3.243F6A8885A308D313198A2E....
- *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps
- *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
- *
- *	In a test run with more than 200,000 random arguments on a VAX, the
- *	maximum observed error in ulps (units in the last place) was
- *	2.07 ulps.      (comparing against (PI/pi)*(exact acos(x)));
- *
- * 2)  If atan2() uses true pi, then
- *
- *	acos(x) returns the exact acos(x) with error below about 2 ulps.
- *
- *	In a test run with more than 1,024,000 random arguments on a VAX, the
- *	maximum observed error in ulps (units in the last place) was
- *	2.15 ulps.
- */
-
-double
-acos(double x)
-{
-	double t, one = 1.0;
-
-	if (isnan(x))
-		return (x);
-
-	if( x != -1.0)
-	    t=atan2(sqrt((one-x)/(one+x)),one);
-	else
-	    t=atan2(one,0.0);	/* t = PI/2 */
-	return(t+t);
-}
-DEF_STD(acos);
-LDBL_UNUSED_CLONE(acos);
diff --git a/lib/libm/noieee_src/n_asinh.c b/lib/libm/noieee_src/n_asinh.c
deleted file mode 100644
index 6307c9356b4..00000000000
--- a/lib/libm/noieee_src/n_asinh.c
+++ /dev/null
@@ -1,89 +0,0 @@
-/*	$OpenBSD: n_asinh.c,v 1.10 2009/10/27 23:59:29 deraadt Exp $	*/
-/*	$NetBSD: n_asinh.c,v 1.1 1995/10/10 23:36:35 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* ASINH(X)
- * RETURN THE INVERSE HYPERBOLIC SINE OF X
- * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 2/16/85;
- * REVISED BY K.C. NG on 3/7/85, 3/24/85, 4/16/85.
- *
- * Required system supported functions :
- *	copysign(x,y)
- *	sqrt(x)
- *
- * Required kernel function:
- *	log1p(x) 		...return log(1+x)
- *
- * Method :
- *	Based on
- *		asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
- *	we have
- *	asinh(x) := x  if  1+x*x=1,
- *		 := sign(x)*(log1p(x)+ln2))	 if sqrt(1+x*x)=x, else
- *		 := sign(x)*log1p(|x| + |x|/(1/|x| + sqrt(1+(1/|x|)^2)) )
- *
- * Accuracy:
- *	asinh(x) returns the exact inverse hyperbolic sine of x nearly rounded.
- *	In a test run with 52,000 random arguments on a VAX, the maximum
- *	observed error was 1.58 ulps (units in the last place).
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-static const double ln2hi = 6.9314718055829871446E-1;
-static const double ln2lo = 1.6465949582897081279E-12;
-
-double
-asinh(double x)
-{
-	double t,s;
-	static const double	small=1.0E-10,	/* fl(1+small*small) == 1 */
-				big  =1.0E20,	/* fl(1+big) == big */
-				one  =1.0   ;
-
-	if (isnan(x))
-		return (x);
-
-	if((t=copysign(x,one))>small)
-	    if(t<big) {
-	     	s=one/t; return(copysign(log1p(t+t/(s+sqrt(one+s*s))),x)); }
-	    else	/* if |x| > big */
-		{s=log1p(t)+ln2lo; return(copysign(s+ln2hi,x));}
-	else	/* if |x| < small */
-	    return(x);
-}
diff --git a/lib/libm/noieee_src/n_atan.c b/lib/libm/noieee_src/n_atan.c
deleted file mode 100644
index 79747e2c170..00000000000
--- a/lib/libm/noieee_src/n_atan.c
+++ /dev/null
@@ -1,93 +0,0 @@
-/*	$OpenBSD: n_atan.c,v 1.14 2016/09/12 19:47:02 guenther Exp $	*/
-/*	$NetBSD: n_atan.c,v 1.1 1995/10/10 23:36:36 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* ATAN(X)
- * RETURNS ARC TANGENT OF X
- * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
- * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
- *
- * Required kernel function:
- *	atan2(y,x)
- *
- * Method:
- *	atan(x) = atan2(x,1.0).
- *
- * Special case:
- *	if x is NaN, return x itself.
- *
- * Accuracy:
- * 1)  If atan2() uses machine PI, then
- *
- *	atan(x) returns (PI/pi) * (the exact arc tangent of x) nearly rounded;
- *	and PI is the exact pi rounded to machine precision (see atan2 for
- *      details):
- *
- *	in decimal:
- *		pi = 3.141592653589793 23846264338327 .....
- *    53 bits   PI = 3.141592653589793 115997963 ..... ,
- *    56 bits   PI = 3.141592653589793 227020265 ..... ,
- *
- *	in hexadecimal:
- *		pi = 3.243F6A8885A308D313198A2E....
- *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps
- *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
- *
- *	In a test run with more than 200,000 random arguments on a VAX, the
- *	maximum observed error in ulps (units in the last place) was
- *	0.86 ulps.      (comparing against (PI/pi)*(exact atan(x))).
- *
- * 2)  If atan2() uses true pi, then
- *
- *	atan(x) returns the exact atan(x) with error below about 2 ulps.
- *
- *	In a test run with more than 1,024,000 random arguments on a VAX, the
- *	maximum observed error in ulps (units in the last place) was
- *	0.85 ulps.
- */
-
-#include <math.h>
-
-double
-atan(double x)
-{
-	double one=1.0;
-	return(atan2(x,one));
-}
-DEF_STD(atan);
-LDBL_CLONE(atan);
-
-float
-atanf(float x)
-{
-	float one=1.0f;
-	return(atan2f(x,one));
-}
-DEF_STD(atanf);
diff --git a/lib/libm/noieee_src/n_atan2.c b/lib/libm/noieee_src/n_atan2.c
deleted file mode 100644
index e99441349bc..00000000000
--- a/lib/libm/noieee_src/n_atan2.c
+++ /dev/null
@@ -1,247 +0,0 @@
-/*	$OpenBSD: n_atan2.c,v 1.20 2016/09/12 19:47:02 guenther Exp $	*/
-/*	$NetBSD: n_atan2.c,v 1.1 1995/10/10 23:36:37 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* ATAN2(Y,X)
- * RETURN ARG (X+iY)
- * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85;
- * REVISED BY K.C. NG on 2/7/85, 2/13/85, 3/7/85, 3/30/85, 6/29/85.
- *
- * Required system supported functions :
- *	copysign(x,y)
- *	scalbn(x,y)
- *	logb(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):
- *		ARG (x+iy) = arctan(y/x)   	   ... if x > 0,
- *		ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
- *	3. According to the integer k=4t+0.25 truncated , t=y/x, the argument
- *	   is further reduced to one of the following intervals and the
- *	   arctangent of y/x is evaluated by the corresponding formula:
- *
- *         [0,7/16]	   atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
- *	   [7/16,11/16]    atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) )
- *	   [11/16.19/16]   atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) )
- *	   [19/16,39/16]   atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) )
- *	   [39/16,INF]     atan(y/x) = atan(INF) + atan( -x/y )
- *
- * Special cases:
- * Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y).
- *
- *	ARG( NAN , (anything) ) is NaN;
- *	ARG( (anything), NaN ) is NaN;
- *	ARG(+(anything but NaN), +-0) is +-0  ;
- *	ARG(-(anything but NaN), +-0) is +-PI ;
- *	ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2;
- *	ARG( +INF,+-(anything but INF and NaN) ) is +-0 ;
- *	ARG( -INF,+-(anything but INF and NaN) ) is +-PI;
- *	ARG( +INF,+-INF ) is +-PI/4 ;
- *	ARG( -INF,+-INF ) is +-3PI/4;
- *	ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2;
- *
- * Accuracy:
- *	atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded,
- *	where
- *
- *	in decimal:
- *		pi = 3.141592653589793 23846264338327 .....
- *    53 bits   PI = 3.141592653589793 115997963 ..... ,
- *    56 bits   PI = 3.141592653589793 227020265 ..... ,
- *
- *	in hexadecimal:
- *		pi = 3.243F6A8885A308D313198A2E....
- *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps
- *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
- *
- *	In a test run with 356,000 random argument on [-1,1] * [-1,1] on a
- *	VAX, the maximum observed error was 1.41 ulps (units of the last place)
- *	compared with (PI/pi)*(the exact ARG(x+iy)).
- *
- * Note:
- *	We use machine PI (the true pi rounded) in place of the actual
- *	value of pi for all the trig and inverse trig functions. In general,
- *	if trig is one of sin, cos, tan, then computed trig(y) returns the
- *	exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig
- *	returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the
- *	trig functions have period PI, and trig(arctrig(x)) returns x for
- *	all critical values x.
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include <math.h>
-
-#include "mathimpl.h"
-
-static const double athfhi = 4.6364760900080611433E-1;
-static const double athflo = 1.9338828231967579916E-19;
-static const double PIo4 = 7.8539816339744830676E-1;
-static const double at1fhi = 9.8279372324732906796E-1;
-static const double at1flo = -3.5540295636764633916E-18;
-static const double PIo2 = 1.5707963267948966135E0;
-static const double PI = 3.1415926535897932270E0;
-static const double a1 = 3.3333333333333473730E-1;
-static const double a2 = -2.0000000000017730678E-1;
-static const double a3 = 1.4285714286694640301E-1;
-static const double a4 = -1.1111111135032672795E-1;
-static const double a5 = 9.0909091380563043783E-2;
-static const double a6 = -7.6922954286089459397E-2;
-static const double a7 = 6.6663180891693915586E-2;
-static const double a8 = -5.8772703698290408927E-2;
-static const double a9 = 5.2170707402812969804E-2;
-static const double a10 = -4.4895863157820361210E-2;
-static const double a11 = 3.3006147437343875094E-2;
-static const double a12 = -1.4614844866464185439E-2;
-
-float
-atan2f(float x, float y)
-{
-	return (float)atan2((double)x, (double)y);
-}
-
-double
-atan2(double y, double x)
-{
-	static const double zero=0, one=1, small=1.0E-9, big=1.0E18;
-	double t,z,signy,signx,hi,lo;
-	int k,m;
-
-    /* if x or y is NAN */
-	if (isnan(x))
-		return (x);
-	if (isnan(y))
-		return (y);
-
-    /* copy down the sign of y and x */
-	signy = copysign(one,y) ;
-	signx = copysign(one,x) ;
-
-    /* if x is 1.0, goto begin */
-	if(x==1) { y=copysign(y,one); t=y; if(isfinite(t)) goto begin;}
-
-    /* when y = 0 */
-	if(y==zero) return((signx==one)?y:copysign(PI,signy));
-
-    /* when x = 0 */
-	if(x==zero) return(copysign(PIo2,signy));
-
-    /* when x is INF */
-	if(!isfinite(x))
-	    if(!isfinite(y))
-		return(copysign((signx==one)?PIo4:3*PIo4,signy));
-	    else
-		return(copysign((signx==one)?zero:PI,signy));
-
-    /* when y is INF */
-	if(!isfinite(y)) return(copysign(PIo2,signy));
-
-    /* compute y/x */
-	x=copysign(x,one);
-	y=copysign(y,one);
-	if((m=(k=logb(y))-logb(x)) > 60) t=big+big;
-	    else if(m < -80 ) t=y/x;
-	    else { t = y/x ; y = scalbn(y,-k); x=scalbn(x,-k); }
-
-    /* begin argument reduction */
-begin:
-	if (t < 2.4375) {
-
-	/* truncate 4(t+1/16) to integer for branching */
-	    k = 4 * (t+0.0625);
-	    switch (k) {
-
-	    /* t is in [0,7/16] */
-	    case 0:
-	    case 1:
-		if (t < small) {
-			if (big + small > 0.0)	/* raise inexact flag */
-				return (copysign((signx>zero)?t:PI-t,signy));
-		}
-
-		hi = zero;  lo = zero;  break;
-
-	    /* t is in [7/16,11/16] */
-	    case 2:
-		hi = athfhi; lo = athflo;
-		z = x+x;
-		t = ( (y+y) - x ) / ( z +  y ); break;
-
-	    /* t is in [11/16,19/16] */
-	    case 3:
-	    case 4:
-		hi = PIo4; lo = zero;
-		t = ( y - x ) / ( x + y ); break;
-
-	    /* t is in [19/16,39/16] */
-	    default:
-		hi = at1fhi; lo = at1flo;
-		z = y-x; y=y+y+y; t = x+x;
-		t = ( (z+z)-x ) / ( t + y ); break;
-	    }
-	}
-	/* end of if (t < 2.4375) */
-
-	else
-	{
-	    hi = PIo2; lo = zero;
-
-	    /* t is in [2.4375, big] */
-	    if (t <= big)  t = - x / y;
-
-	    /* t is in [big, INF] */
-	    else {
-		if (big + small > 0.0)	/* raise inexact flag */
-			t = zero;
-	    }
-	}
-    /* end of argument reduction */
-
-    /* compute atan(t) for t in [-.4375, .4375] */
-	z = t*t;
-#if defined(__vax__)
-	z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+
-			z*(a9+z*(a10+z*(a11+z*a12))))))))))));
-#else	/* defined(__vax__) */
-	z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+
-			z*(a9+z*(a10+z*a11)))))))))));
-#endif	/* defined(__vax__) */
-	z = lo - z; z += t; z += hi;
-
-	return(copysign((signx>zero)?z:PI-z,signy));
-}
-DEF_STD(atan2);
-LDBL_CLONE(atan2);
diff --git a/lib/libm/noieee_src/n_atanh.c b/lib/libm/noieee_src/n_atanh.c
deleted file mode 100644
index 513065aca5e..00000000000
--- a/lib/libm/noieee_src/n_atanh.c
+++ /dev/null
@@ -1,79 +0,0 @@
-/*	$OpenBSD: n_atanh.c,v 1.7 2009/10/27 23:59:29 deraadt Exp $	*/
-/*	$NetBSD: n_atanh.c,v 1.1 1995/10/10 23:36:38 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* ATANH(X)
- * RETURN THE HYPERBOLIC ARC TANGENT OF X
- * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85;
- * REVISED BY K.C. NG on 2/7/85, 3/7/85, 8/18/85.
- *
- * Required kernel function:
- *	log1p(x) 	...return log(1+x)
- *
- * Method:
- *	Return
- *                          1              2x                          x
- *		atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
- *                          2             1 - x                      1 - x
- *
- * Special cases:
- *	atanh(x) is NaN if |x| > 1 with signal;
- *	atanh(NaN) is that NaN with no signal;
- *	atanh(+-1) is +-INF with signal.
- *
- * Accuracy:
- *	atanh(x) returns the exact hyperbolic arc tangent of x nearly rounded.
- *	In a test run with 512,000 random arguments on a VAX, the maximum
- *	observed error was 1.87 ulps (units in the last place) at
- *	x= -3.8962076028810414000e-03.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-#if defined(__vax__)
-#include <errno.h>
-#endif	/* defined(__vax__) */
-
-double
-atanh(double x)
-{
-	double z;
-	z = copysign(0.5,x);
-	x = copysign(x,1.0);
-#if defined(__vax__)
-	if (x == 1.0) {
-	    return(copysign(1.0,z)*infnan(ERANGE));	/* sign(x)*INF */
-	}
-#endif	/* defined(__vax__) */
-	x = x/(1.0-x);
-	return( z*log1p(x+x) );
-}
diff --git a/lib/libm/noieee_src/n_cbrt.c b/lib/libm/noieee_src/n_cbrt.c
deleted file mode 100644
index eb2c1e1eaac..00000000000
--- a/lib/libm/noieee_src/n_cbrt.c
+++ /dev/null
@@ -1,115 +0,0 @@
-/*	$OpenBSD: n_cbrt.c,v 1.8 2013/07/15 04:08:26 espie Exp $	*/
-/*	$NetBSD: n_cbrt.c,v 1.1 1995/10/10 23:36:40 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include <sys/types.h>
-
-/* kahan's cube root (53 bits IEEE double precision)
- * for IEEE machines only
- * coded in C by K.C. Ng, 4/30/85
- *
- * Accuracy:
- *	better than 0.667 ulps according to an error analysis. Maximum
- * error observed was 0.666 ulps in an 1,000,000 random arguments test.
- *
- * Warning: this code is semi machine dependent; the ordering of words in
- * a floating point number must be known in advance. I assume that the
- * long interger at the address of a floating point number will be the
- * leading 32 bits of that floating point number (i.e., sign, exponent,
- * and the 20 most significant bits).
- * On a National machine, it has different ordering; therefore, this code
- * must be compiled with flag -DNATIONAL.
- */
-#if !defined(__vax__)
-
-static const unsigned long
-		     B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
-	             B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
-static const double
-	    C= 19./35.,
-	    D= -864./1225.,
-	    E= 99./70.,
-	    F= 45./28.,
-	    G= 5./14.;
-
-float
-cbrtf(float x)
-{
-	return (float)cbrt((double) x);
-}
-
-double
-cbrt(double x)
-{
-	double r,s,t=0.0,w;
-	unsigned long *px = (unsigned long *) &x,
-	              *pt = (unsigned long *) &t,
-		      mexp,sign;
-	const int n0=0,n1=1;
-
-	mexp=px[n0]&0x7ff00000;
-	if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */
-	if(x==0.0) return(x);		/* cbrt(0) is itself */
-
-	sign=px[n0]&0x80000000; /* sign= sign(x) */
-	px[n0] ^= sign;		/* x=|x| */
-
-
-    /* rough cbrt to 5 bits */
-	if(mexp==0) 		/* subnormal number */
-	  {pt[n0]=0x43500000; 	/* set t= 2**54 */
-	   t*=x; pt[n0]=pt[n0]/3+B2;
-	  }
-	else
-	  pt[n0]=px[n0]/3+B1;
-
-
-    /* new cbrt to 23 bits, may be implemented in single precision */
-	r=t*t/x;
-	s=C+r*t;
-	t*=G+F/(s+E+D/s);
-
-    /* chopped to 20 bits and make it larger than cbrt(x) */
-	pt[n1]=0; pt[n0]+=0x00000001;
-
-
-    /* one step newton iteration to 53 bits with error less than 0.667 ulps */
-	s=t*t;		/* t*t is exact */
-	r=x/s;
-	w=t+t;
-	r=(r-t)/(w+r);	/* r-t is exact */
-	t=t+t*r;
-
-
-    /* retore the sign bit */
-	pt[n0] |= sign;
-	return(t);
-}
-#endif
diff --git a/lib/libm/noieee_src/n_cosh.c b/lib/libm/noieee_src/n_cosh.c
deleted file mode 100644
index 72fe1853972..00000000000
--- a/lib/libm/noieee_src/n_cosh.c
+++ /dev/null
@@ -1,118 +0,0 @@
-/*	$OpenBSD: n_cosh.c,v 1.10 2009/10/27 23:59:29 deraadt Exp $	*/
-/*	$NetBSD: n_cosh.c,v 1.1 1995/10/10 23:36:42 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* COSH(X)
- * RETURN THE HYPERBOLIC COSINE OF X
- * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85;
- * REVISED BY K.C. NG on 2/8/85, 2/23/85, 3/7/85, 3/29/85, 4/16/85.
- *
- * Required system supported functions :
- *	copysign(x,y)
- *	scalbn(x,N)
- *
- * Required kernel function:
- *	exp(x)
- *	exp__E(x,c)	...return exp(x+c)-1-x for |x|<0.3465
- *
- * Method :
- *	1. Replace x by |x|.
- *	2.
- *		                                        [ exp(x) - 1 ]^2
- *	    0        <= x <= 0.3465  :  cosh(x) := 1 + -------------------
- *			       			           2*exp(x)
- *
- *		                                   exp(x) +  1/exp(x)
- *	    0.3465   <= x <= 22      :  cosh(x) := -------------------
- *			       			           2
- *	    22       <= x <= lnovfl  :  cosh(x) := exp(x)/2
- *	    lnovfl   <= x <= lnovfl+log(2)
- *				     :  cosh(x) := exp(x)/2 (avoid overflow)
- *	    log(2)+lnovfl <  x <  INF:  overflow to INF
- *
- *	Note: .3465 is a number near one half of ln2.
- *
- * Special cases:
- *	cosh(x) is x if x is +INF, -INF, or NaN.
- *	only cosh(0)=1 is exact for finite x.
- *
- * Accuracy:
- *	cosh(x) returns the exact hyperbolic cosine of x nearly rounded.
- *	In a test run with 768,000 random arguments on a VAX, the maximum
- *	observed error was 1.23 ulps (units in the last place).
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-static const double mln2hi = 8.8029691931113054792E1;
-static const double mln2lo = -4.9650192275318476525E-16;
-static const double lnovfl = 8.8029691931113053016E1;
-
-#if defined(__vax__)
-static max = 126                      ;
-#else	/* defined(__vax__) */
-static max = 1023                     ;
-#endif	/* defined(__vax__) */
-
-double
-cosh(double x)
-{
-	static const double half=1.0/2.0,
-		one=1.0, small=1.0E-18; /* fl(1+small)==1 */
-	double t;
-
-	if (isnan(x))
-		return (x);
-
-	if((x=copysign(x,one)) <= 22)
-	    if(x<0.3465)
-		if(x<small) return(one+x);
-		else {t=x+__exp__E(x,0.0);x=t+t; return(one+t*t/(2.0+x)); }
-
-	    else /* for x lies in [0.3465,22] */
-	        { t=exp(x); return((t+one/t)*half); }
-
-	if( lnovfl <= x && x <= (lnovfl+0.7))
-        /* for x lies in [lnovfl, lnovfl+ln2], decrease x by ln(2^(max+1))
-         * and return 2^max*exp(x) to avoid unnecessary overflow
-         */
-	    return(scalbn(exp((x-mln2hi)-mln2lo), max));
-
-	else
-	    return(exp(x)*half);	/* for large x,  cosh(x)=exp(x)/2 */
-}
diff --git a/lib/libm/noieee_src/n_erf.c b/lib/libm/noieee_src/n_erf.c
deleted file mode 100644
index 00f2df2483e..00000000000
--- a/lib/libm/noieee_src/n_erf.c
+++ /dev/null
@@ -1,388 +0,0 @@
-/*	$OpenBSD: n_erf.c,v 1.9 2016/09/12 19:47:02 guenther Exp $	*/
-/*	$NetBSD: n_erf.c,v 1.1 1995/10/10 23:36:43 ragge Exp $	*/
-/*-
- * Copyright (c) 1992, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-/* Modified Nov 30, 1992 P. McILROY:
- *	Replaced expansions for x >= 1.25 (error 1.7ulp vs ~6ulp)
- * Replaced even+odd with direct calculation for x < .84375,
- * to avoid destructive cancellation.
- *
- * Performance of erfc(x):
- * In 300000 trials in the range [.83, .84375] the
- * maximum observed error was 3.6ulp.
- *
- * In [.84735,1.25] the maximum observed error was <2.5ulp in
- * 100000 runs in the range [1.2, 1.25].
- *
- * In [1.25,26] (Not including subnormal results)
- * the error is < 1.7ulp.
- */
-
-/* double erf(double x)
- * double erfc(double x)
- *			     x
- *		      2      |\
- *     erf(x)  =  ---------  | exp(-t*t)dt
- *		   sqrt(pi) \|
- *			     0
- *
- *     erfc(x) =  1-erf(x)
- *
- * Method:
- *      1. Reduce x to |x| by erf(-x) = -erf(x)
- *	2. For x in [0, 0.84375]
- *	    erf(x)  = x + x*P(x^2)
- *          erfc(x) = 1 - erf(x)           if x<=0.25
- *                  = 0.5 + ((0.5-x)-x*P)  if x in [0.25,0.84375]
- *	   where
- *			2		 2	  4		  20
- *              P =  P(x ) = (p0 + p1 * x + p2 * x + ... + p10 * x  )
- * 	   is an approximation to (erf(x)-x)/x with precision
- *
- *						 -56.45
- *			| P - (erf(x)-x)/x | <= 2
- *
- *
- *	   Remark. The formula is derived by noting
- *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
- *	   and that
- *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
- *	   is close to one. The interval is chosen because the fixed
- *	   point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
- *	   near 0.6174), and by some experiment, 0.84375 is chosen to
- * 	   guarantee the error is less than one ulp for erf.
- *
- *      3. For x in [0.84375,1.25], let s = x - 1, and
- *         c = 0.84506291151 rounded to single (24 bits)
- *         	erf(x)  = c  + P1(s)/Q1(s)
- *         	erfc(x) = (1-c)  - P1(s)/Q1(s)
- *         	|P1/Q1 - (erf(x)-c)| <= 2**-59.06
- *	   Remark: here we use the taylor series expansion at x=1.
- *		erf(1+s) = erf(1) + s*Poly(s)
- *			 = 0.845.. + P1(s)/Q1(s)
- *	   That is, we use rational approximation to approximate
- *			erf(1+s) - (c = (single)0.84506291151)
- *	   Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
- *	   where
- *		P1(s) = degree 6 poly in s
- *		Q1(s) = degree 6 poly in s
- *
- *	4. For x in [1.25, 2]; [2, 4]
- *         	erf(x)  = 1.0 - tiny
- *		erfc(x)	= (1/x)exp(-x*x-(.5*log(pi) -.5z + R(z)/S(z))
- *
- *	Where z = 1/(x*x), R is degree 9, and S is degree 3;
- *
- *      5. For x in [4,28]
- *         	erf(x)  = 1.0 - tiny
- *		erfc(x)	= (1/x)exp(-x*x-(.5*log(pi)+eps + zP(z))
- *
- *	Where P is degree 14 polynomial in 1/(x*x).
- *
- *      Notes:
- *	   Here 4 and 5 make use of the asymptotic series
- *			  exp(-x*x)
- *		erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) );
- *			  x*sqrt(pi)
- *
- *		where for z = 1/(x*x)
- *		P(z) ~ z/2*(-1 + z*3/2*(1 + z*5/2*(-1 + z*7/2*(1 +...))))
- *
- *	   Thus we use rational approximation to approximate
- *              erfc*x*exp(x*x) ~ 1/sqrt(pi);
- *
- *		The error bound for the target function, G(z) for
- *		the interval
- *		[4, 28]:
- * 		|eps + 1/(z)P(z) - G(z)| < 2**(-56.61)
- *		for [2, 4]:
- *      	|R(z)/S(z) - G(z)|	 < 2**(-58.24)
- *		for [1.25, 2]:
- *		|R(z)/S(z) - G(z)|	 < 2**(-58.12)
- *
- *      6. For inf > x >= 28
- *         	erf(x)  = 1 - tiny  (raise inexact)
- *         	erfc(x) = tiny*tiny (raise underflow)
- *
- *      7. Special cases:
- *         	erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
- *         	erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
- *	   	erfc/erf(NaN) is NaN
- */
-
-#if defined(__vax__)
-#define _IEEE	0
-#define TRUNC(x) (double) (float) (x)
-#else
-#define _IEEE	1
-#define TRUNC(x) *(((int *) &x) + 1) &= 0xf8000000
-#define infnan(x) 0.0
-#endif
-
-static double
-tiny	    = 1e-300,
-half	    = 0.5,
-one	    = 1.0,
-two	    = 2.0,
-c 	    = 8.45062911510467529297e-01, /* (float)0.84506291151 */
-/*
- * Coefficients for approximation to erf in [0,0.84375]
- */
-p0t8 = 1.02703333676410051049867154944018394163280,
-p0 =   1.283791670955125638123339436800229927041e-0001,
-p1 =  -3.761263890318340796574473028946097022260e-0001,
-p2 =   1.128379167093567004871858633779992337238e-0001,
-p3 =  -2.686617064084433642889526516177508374437e-0002,
-p4 =   5.223977576966219409445780927846432273191e-0003,
-p5 =  -8.548323822001639515038738961618255438422e-0004,
-p6 =   1.205520092530505090384383082516403772317e-0004,
-p7 =  -1.492214100762529635365672665955239554276e-0005,
-p8 =   1.640186161764254363152286358441771740838e-0006,
-p9 =  -1.571599331700515057841960987689515895479e-0007,
-p10=   1.073087585213621540635426191486561494058e-0008;
-/*
- * Coefficients for approximation to erf in [0.84375,1.25]
- */
-static double
-pa0 =  -2.362118560752659485957248365514511540287e-0003,
-pa1 =   4.148561186837483359654781492060070469522e-0001,
-pa2 =  -3.722078760357013107593507594535478633044e-0001,
-pa3 =   3.183466199011617316853636418691420262160e-0001,
-pa4 =  -1.108946942823966771253985510891237782544e-0001,
-pa5 =   3.547830432561823343969797140537411825179e-0002,
-pa6 =  -2.166375594868790886906539848893221184820e-0003,
-qa1 =   1.064208804008442270765369280952419863524e-0001,
-qa2 =   5.403979177021710663441167681878575087235e-0001,
-qa3 =   7.182865441419627066207655332170665812023e-0002,
-qa4 =   1.261712198087616469108438860983447773726e-0001,
-qa5 =   1.363708391202905087876983523620537833157e-0002,
-qa6 =   1.198449984679910764099772682882189711364e-0002;
-/*
- * log(sqrt(pi)) for large x expansions.
- * The tail (lsqrtPI_lo) is included in the rational
- * approximations.
-*/
-static double
-   lsqrtPI_hi = .5723649429247000819387380943226;
-/*
- * lsqrtPI_lo = .000000000000000005132975581353913;
- *
- * Coefficients for approximation to erfc in [2, 4]
-*/
-static double
-rb0  =	-1.5306508387410807582e-010,	/* includes lsqrtPI_lo */
-rb1  =	 2.15592846101742183841910806188e-008,
-rb2  =	 6.24998557732436510470108714799e-001,
-rb3  =	 8.24849222231141787631258921465e+000,
-rb4  =	 2.63974967372233173534823436057e+001,
-rb5  =	 9.86383092541570505318304640241e+000,
-rb6  =	-7.28024154841991322228977878694e+000,
-rb7  =	 5.96303287280680116566600190708e+000,
-rb8  =	-4.40070358507372993983608466806e+000,
-rb9  =	 2.39923700182518073731330332521e+000,
-rb10 =	-6.89257464785841156285073338950e-001,
-sb1  =	 1.56641558965626774835300238919e+001,
-sb2  =	 7.20522741000949622502957936376e+001,
-sb3  =	 9.60121069770492994166488642804e+001;
-/*
- * Coefficients for approximation to erfc in [1.25, 2]
-*/
-static double
-rc0  =	-2.47925334685189288817e-007,	/* includes lsqrtPI_lo */
-rc1  =	 1.28735722546372485255126993930e-005,
-rc2  =	 6.24664954087883916855616917019e-001,
-rc3  =	 4.69798884785807402408863708843e+000,
-rc4  =	 7.61618295853929705430118701770e+000,
-rc5  =	 9.15640208659364240872946538730e-001,
-rc6  =	-3.59753040425048631334448145935e-001,
-rc7  =	 1.42862267989304403403849619281e-001,
-rc8  =	-4.74392758811439801958087514322e-002,
-rc9  =	 1.09964787987580810135757047874e-002,
-rc10 =	-1.28856240494889325194638463046e-003,
-sc1  =	 9.97395106984001955652274773456e+000,
-sc2  =	 2.80952153365721279953959310660e+001,
-sc3  =	 2.19826478142545234106819407316e+001;
-/*
- * Coefficients for approximation to  erfc in [4,28]
- */
-static double
-rd0  =	-2.1491361969012978677e-016,	/* includes lsqrtPI_lo */
-rd1  =	-4.99999999999640086151350330820e-001,
-rd2  =	 6.24999999772906433825880867516e-001,
-rd3  =	-1.54166659428052432723177389562e+000,
-rd4  =	 5.51561147405411844601985649206e+000,
-rd5  =	-2.55046307982949826964613748714e+001,
-rd6  =	 1.43631424382843846387913799845e+002,
-rd7  =	-9.45789244999420134263345971704e+002,
-rd8  =	 6.94834146607051206956384703517e+003,
-rd9  =	-5.27176414235983393155038356781e+004,
-rd10 =	 3.68530281128672766499221324921e+005,
-rd11 =	-2.06466642800404317677021026611e+006,
-rd12 =	 7.78293889471135381609201431274e+006,
-rd13 =	-1.42821001129434127360582351685e+007;
-
-double
-erf(double x)
-{
-	double R, S, P, Q, ax, s, y, z, r;
-	if(!isfinite(x)) {		/* erf(nan)=nan */
-	    if (isnan(x))
-		return(x);
-	    return (x > 0 ? one : -one); /* erf(+/-inf)= +/-1 */
-	}
-	if ((ax = x) < 0)
-		ax = - ax;
-	if (ax < .84375) {
-	    if (ax < 3.7e-09) {
-		if (ax < 1.0e-308)
-		    return 0.125*(8.0*x+p0t8*x);  /*avoid underflow */
-		return x + p0*x;
-	    }
-	    y = x*x;
-	    r = y*(p1+y*(p2+y*(p3+y*(p4+y*(p5+
-			y*(p6+y*(p7+y*(p8+y*(p9+y*p10)))))))));
-	    return x + x*(p0+r);
-	}
-	if (ax < 1.25) {		/* 0.84375 <= |x| < 1.25 */
-	    s = fabs(x)-one;
-	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
-	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
-	    if (x>=0)
-		return (c + P/Q);
-	    else
-		return (-c - P/Q);
-	}
-	if (ax >= 6.0) {		/* inf>|x|>=6 */
-	    if (x >= 0.0)
-		return (one-tiny);
-	    else
-		return (tiny-one);
-	}
-    /* 1.25 <= |x| < 6 */
-	z = -ax*ax;
-	s = -one/z;
-	if (ax < 2.0) {
-		R = rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*(rc5+
-			s*(rc6+s*(rc7+s*(rc8+s*(rc9+s*rc10)))))))));
-		S = one+s*(sc1+s*(sc2+s*sc3));
-	} else {
-		R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+
-			s*(rb6+s*(rb7+s*(rb8+s*(rb9+s*rb10)))))))));
-		S = one+s*(sb1+s*(sb2+s*sb3));
-	}
-	y = (R/S -.5*s) - lsqrtPI_hi;
-	z += y;
-	z = exp(z)/ax;
-	if (x >= 0)
-		return (one-z);
-	else
-		return (z-one);
-}
-DEF_STD(erf);
-
-double
-erfc(double x)
-{
-	double R, S, P, Q, s, ax, y, z, r;
-	if (!isfinite(x)) {
-		if (isnan(x))		/* erfc(NaN) = NaN */
-			return(x);
-		else if (x > 0)		/* erfc(+-inf)=0,2 */
-			return 0.0;
-		else
-			return 2.0;
-	}
-	if ((ax = x) < 0)
-		ax = -ax;
-	if (ax < .84375) {			/* |x|<0.84375 */
-	    if (ax < 1.38777878078144568e-17)  	/* |x|<2**-56 */
-		return one-x;
-	    y = x*x;
-	    r = y*(p1+y*(p2+y*(p3+y*(p4+y*(p5+
-			y*(p6+y*(p7+y*(p8+y*(p9+y*p10)))))))));
-	    if (ax < .0625) {  	/* |x|<2**-4 */
-		return (one-(x+x*(p0+r)));
-	    } else {
-		r = x*(p0+r);
-		r += (x-half);
-	        return (half - r);
-	    }
-	}
-	if (ax < 1.25) {		/* 0.84375 <= |x| < 1.25 */
-	    s = ax-one;
-	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
-	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
-	    if (x>=0) {
-	        z  = one-c; return z - P/Q;
-	    } else {
-		z = c+P/Q; return one+z;
-	    }
-	}
-	if (ax >= 28)	/* Out of range */
- 		if (x>0)
-			return (tiny*tiny);
-		else
-			return (two-tiny);
-	z = ax;
-	TRUNC(z);
-	y = z - ax; y *= (ax+z);
-	z *= -z;			/* Here z + y = -x^2 */
-		s = one/(-z-y);		/* 1/(x*x) */
-	if (ax >= 4) {			/* 6 <= ax */
-		R = s*(rd1+s*(rd2+s*(rd3+s*(rd4+s*(rd5+
-			s*(rd6+s*(rd7+s*(rd8+s*(rd9+s*(rd10
-			+s*(rd11+s*(rd12+s*rd13))))))))))));
-		y += rd0;
-	} else if (ax >= 2) {
-		R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+
-			s*(rb6+s*(rb7+s*(rb8+s*(rb9+s*rb10)))))))));
-		S = one+s*(sb1+s*(sb2+s*sb3));
-		y += R/S;
-		R = -.5*s;
-	} else {
-		R = rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*(rc5+
-			s*(rc6+s*(rc7+s*(rc8+s*(rc9+s*rc10)))))))));
-		S = one+s*(sc1+s*(sc2+s*sc3));
-		y += R/S;
-		R = -.5*s;
-	}
-	/* return exp(-x^2 - lsqrtPI_hi + R + y)/x;	*/
-	s = ((R + y) - lsqrtPI_hi) + z;
-	y = (((z-s) - lsqrtPI_hi) + R) + y;
-	r = __exp__D(s, y)/x;
-	if (x>0)
-		return r;
-	else
-		return two-r;
-}
-DEF_STD(erfc);
diff --git a/lib/libm/noieee_src/n_exp.c b/lib/libm/noieee_src/n_exp.c
deleted file mode 100644
index 986ff59b8f5..00000000000
--- a/lib/libm/noieee_src/n_exp.c
+++ /dev/null
@@ -1,174 +0,0 @@
-/*	$OpenBSD: n_exp.c,v 1.12 2016/09/12 19:47:02 guenther Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* EXP(X)
- * RETURN THE EXPONENTIAL OF X
- * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
- * CODED IN C BY K.C. NG, 1/19/85;
- * REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86.
- *
- * Required system supported functions:
- *	scalbn(x,n)
- *	copysign(x,y)
- *	isfinite(x)
- *
- * Method:
- *	1. Argument Reduction: given the input x, find r and integer k such
- *	   that
- *	                   x = k*ln2 + r,  |r| <= 0.5*ln2 .
- *	   r will be represented as r := z+c for better accuracy.
- *
- *	2. Compute exp(r) by
- *
- *		exp(r) = 1 + r + r*R1/(2-R1),
- *	   where
- *		R1 = x - x^2*(p1+x^2*(p2+x^2*(p3+x^2*(p4+p5*x^2)))).
- *
- *	3. exp(x) = 2^k * exp(r) .
- *
- * Special cases:
- *	exp(INF) is INF, exp(NaN) is NaN;
- *	exp(-INF)=  0;
- *	for finite argument, only exp(0)=1 is exact.
- *
- * Accuracy:
- *	exp(x) returns the exponential of x nearly rounded. In a test run
- *	with 1,156,000 random arguments on a VAX, the maximum observed
- *	error was 0.869 ulps (units in the last place).
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-static const double ln2hi = 6.9314718055829871446E-1;
-static const double ln2lo = 1.6465949582897081279E-12;
-static const double lnhuge = 9.4961163736712506989E1;
-static const double lntiny = -9.5654310917272452386E1;
-static const double invln2 = 1.4426950408889634148E0;
-static const double p1 = 1.6666666666666602251E-1;
-static const double p2 = -2.7777777777015591216E-3;
-static const double p3 = 6.6137563214379341918E-5;
-static const double p4 = -1.6533902205465250480E-6;
-static const double p5 = 4.1381367970572387085E-8;
-
-double
-exp(double x)
-{
-	double z, hi, lo, c;
-	int k;
-
-	if (isnan(x))
-		return (x);
-
-	if( x <= lnhuge ) {
-		if( x >= lntiny ) {
-
-		    /* argument reduction : x --> x - k*ln2 */
-
-			k=invln2*x+copysign(0.5,x);	/* k=NINT(x/ln2) */
-
-		    /* express x-k*ln2 as hi-lo and let x=hi-lo rounded */
-
-			hi=x-k*ln2hi;
-			x=hi-(lo=k*ln2lo);
-
-		    /* return 2^k*[1+x+x*c/(2+c)]  */
-			z=x*x;
-			c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
-			return  scalbn(1.0+(hi-(lo-(x*c)/(2.0-c))),k);
-
-		}
-		/* end of x > lntiny */
-
-		else
-		     /* exp(-big#) underflows to zero */
-		     if(isfinite(x))  return(scalbn(1.0,-5000));
-
-		     /* exp(-INF) is zero */
-		     else return(0.0);
-	}
-	/* end of x < lnhuge */
-
-	else
-	/* exp(INF) is INF, exp(+big#) overflows to INF */
-	    return( isfinite(x) ?  scalbn(1.0,5000)  : x);
-}
-DEF_STD(exp);
-
-/* returns exp(r = x + c) for |c| < |x| with no overlap.  */
-
-double
-__exp__D(double x, double c)
-{
-	double z, hi, lo;
-	int k;
-
-	if (isnan(x))
-		return (x);
-
-	if ( x <= lnhuge ) {
-		if ( x >= lntiny ) {
-
-		    /* argument reduction : x --> x - k*ln2 */
-			z = invln2*x;
-			k = z + copysign(.5, x);
-
-		    /* express (x+c)-k*ln2 as hi-lo and let x=hi-lo rounded */
-
-			hi=(x-k*ln2hi);			/* Exact. */
-			x= hi - (lo = k*ln2lo-c);
-		    /* return 2^k*[1+x+x*c/(2+c)]  */
-			z=x*x;
-			c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
-			c = (x*c)/(2.0-c);
-
-			return  scalbn(1.+(hi-(lo - c)), k);
-		}
-		/* end of x > lntiny */
-
-		else
-		     /* exp(-big#) underflows to zero */
-		     if(isfinite(x))  return(scalbn(1.0,-5000));
-
-		     /* exp(-INF) is zero */
-		     else return(0.0);
-	}
-	/* end of x < lnhuge */
-
-	else
-	/* exp(INF) is INF, exp(+big#) overflows to INF */
-	    return( isfinite(x) ?  scalbn(1.0,5000)  : x);
-}
diff --git a/lib/libm/noieee_src/n_exp__E.c b/lib/libm/noieee_src/n_exp__E.c
deleted file mode 100644
index ccdda138c45..00000000000
--- a/lib/libm/noieee_src/n_exp__E.c
+++ /dev/null
@@ -1,124 +0,0 @@
-/*	$OpenBSD: n_exp__E.c,v 1.12 2009/10/27 23:59:29 deraadt Exp $	*/
-/*	$NetBSD: n_exp__E.c,v 1.1 1995/10/10 23:36:45 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* exp__E(x,c)
- * ASSUMPTION: c << x  SO THAT  fl(x+c)=x.
- * (c is the correction term for x)
- * exp__E RETURNS
- *
- *			 /  exp(x+c) - 1 - x ,  1E-19 < |x| < .3465736
- *       exp__E(x,c) = 	|
- *			 \  0 ,  |x| < 1E-19.
- *
- * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
- * KERNEL FUNCTION OF EXP, EXPM1, POW FUNCTIONS
- * CODED IN C BY K.C. NG, 1/31/85;
- * REVISED BY K.C. NG on 3/16/85, 4/16/85.
- *
- * Required system supported function:
- *	copysign(x,y)
- *
- * Method:
- *	1. Rational approximation. Let r=x+c.
- *	   Based on
- *                                   2 * sinh(r/2)
- *                exp(r) - 1 =   ----------------------   ,
- *                               cosh(r/2) - sinh(r/2)
- *	   exp__E(r) is computed using
- *                   x*x            (x/2)*W - ( Q - ( 2*P  + x*P ) )
- *                   --- + (c + x*[---------------------------------- + c ])
- *                    2                          1 - W
- * 	   where  P := p1*x^2 + p2*x^4,
- *	          Q := q1*x^2 + q2*x^4 (for 56 bits precision, add q3*x^6)
- *	          W := x/2-(Q-x*P),
- *
- *	   (See the listing below for the values of p1,p2,q1,q2,q3. The poly-
- *	    nomials P and Q may be regarded as the approximations to sinh
- *	    and cosh :
- *		sinh(r/2) =  r/2 + r * P  ,  cosh(r/2) =  1 + Q . )
- *
- *         The coefficients were obtained by a special Remes algorithm.
- *
- * Approximation error:
- *
- *   |	exp(x) - 1			   |        2**(-57),  (IEEE double)
- *   | ------------  -  (exp__E(x,0)+x)/x  |  <=
- *   |	     x			           |	    2**(-69).  (VAX D)
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-static const double p1 = 1.5150724356786683059E-2;
-static const double p2 = 6.3112487873718332688E-5;
-static const double q1 = 1.1363478204690669916E-1;
-static const double q2 = 1.2624568129896839182E-3;
-static const double q3 = 1.5021856115869022674E-6;
-
-double
-__exp__E(double x, double c)
-{
-	static const double zero=0.0, one=1.0, half=1.0/2.0, small=1.0E-19;
-	double z,p,q,xp,xh,w;
-	if(copysign(x,one)>small) {
-           z = x*x  ;
-	   p = z*( p1 +z* p2 );
-#if defined(__vax__)
-           q = z*( q1 +z*( q2 +z* q3 ));
-#else	/* defined(__vax__) */
-           q = z*( q1 +z*  q2 );
-#endif	/* defined(__vax__) */
-           xp= x*p     ;
-	   xh= x*half  ;
-           w = xh-(q-xp)  ;
-	   p = p+p;
-	   c += x*((xh*w-(q-(p+xp)))/(one-w)+c);
-	   return(z*half+c);
-	}
-	/* end of |x| > small */
-
-	else {
-	    if(x != zero) {
-		if (one + small >= 1.0)	/* raise the inexact flag */
-			return(copysign(zero,x));
-	    }
-	    else
-		return(copysign(zero,x));
-	}
-
-	/* NOTREACHED */
-}
diff --git a/lib/libm/noieee_src/n_expm1.c b/lib/libm/noieee_src/n_expm1.c
deleted file mode 100644
index 9494c41f3e4..00000000000
--- a/lib/libm/noieee_src/n_expm1.c
+++ /dev/null
@@ -1,149 +0,0 @@
-/*	$OpenBSD: n_expm1.c,v 1.13 2016/09/12 04:39:47 guenther Exp $	*/
-/*	$NetBSD: n_expm1.c,v 1.1 1995/10/10 23:36:46 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* EXPM1(X)
- * RETURN THE EXPONENTIAL OF X MINUS ONE
- * DOUBLE PRECISION (IEEE 53 BITS, VAX D FORMAT 56 BITS)
- * CODED IN C BY K.C. NG, 1/19/85;
- * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/21/85, 4/16/85.
- *
- * Required system supported functions:
- *	scalbn(x,n)
- *	copysign(x,y)
- *	isfinite(x)
- *
- * Kernel function:
- *	exp__E(x,c)
- *
- * Method:
- *	1. Argument Reduction: given the input x, find r and integer k such
- *	   that
- *	                   x = k*ln2 + r,  |r| <= 0.5*ln2 .
- *	   r will be represented as r := z+c for better accuracy.
- *
- *	2. Compute EXPM1(r)=exp(r)-1 by
- *
- *			EXPM1(r=z+c) := z + exp__E(z,c)
- *
- *	3. EXPM1(x) =  2^k * ( EXPM1(r) + 1-2^-k ).
- *
- * 	Remarks:
- *	   1. When k=1 and z < -0.25, we use the following formula for
- *	      better accuracy:
- *			EXPM1(x) = 2 * ( (z+0.5) + exp__E(z,c) )
- *	   2. To avoid rounding error in 1-2^-k where k is large, we use
- *			EXPM1(x) = 2^k * { [z+(exp__E(z,c)-2^-k )] + 1 }
- *	      when k>56.
- *
- * Special cases:
- *	EXPM1(INF) is INF, EXPM1(NaN) is NaN;
- *	EXPM1(-INF)= -1;
- *	for finite argument, only EXPM1(0)=0 is exact.
- *
- * Accuracy:
- *	EXPM1(x) returns the exact (exp(x)-1) nearly rounded. In a test run with
- *	1,166,000 random arguments on a VAX, the maximum observed error was
- *	.872 ulps (units of the last place).
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-static const double ln2hi = 6.9314718055829871446E-1;
-static const double ln2lo = 1.6465949582897081279E-12;
-static const double lnhuge = 9.4961163736712506989E1;
-static const double invln2 = 1.4426950408889634148E0;
-
-double
-expm1(double x)
-{
-	static const double one=1.0, half=1.0/2.0, tiny=1e-300;
-	double  z,hi,lo,c;
-	int k;
-#if defined(__vax__)
-	static prec=56;
-#else	/* defined(__vax__) */
-	static prec=53;
-#endif	/* defined(__vax__) */
-
-	if (isnan(x))
-		return (x);
-
-	if( x <= lnhuge ) {
-		if( x >= -40.0 ) {
-
-		    /* argument reduction : x - k*ln2 */
-			k= invln2 *x+copysign(0.5,x);	/* k=NINT(x/ln2) */
-			hi=x-k*ln2hi ;
-			z=hi-(lo=k*ln2lo);
-			c=(hi-z)-lo;
-
-			if(k==0) return(z+__exp__E(z,c));
-			if(k==1)
-			    if(z< -0.25)
-				{x=z+half;x +=__exp__E(z,c); return(x+x);}
-			    else
-				{z+=__exp__E(z,c); x=half+z; return(x+x);}
-		    /* end of k=1 */
-
-			else {
-			    if(k<=prec)
-			      { x=one-scalbn(one,-k); z += __exp__E(z,c);}
-			    else if(k<100)
-			      { x = __exp__E(z,c)-scalbn(one,-k); x+=z; z=one;}
-			    else
-			      { x = __exp__E(z,c)+z; z=one;}
-
-			    return (scalbn(x+z,k));
-			}
-		}
-		/* end of x > lnunfl */
-
-		else
-		     /* expm1(-big#) rounded to -1 (inexact) */
-		     if(isfinite(x))
-			return(tiny-one);
-
-		     /* expm1(-INF) is -1 */
-		     else return(-one);
-	}
-	/* end of x < lnhuge */
-
-	else
-	/*  expm1(INF) is INF, expm1(+big#) overflows to INF */
-	    return( isfinite(x) ?  scalbn(one,5000) : x);
-}
diff --git a/lib/libm/noieee_src/n_floor.c b/lib/libm/noieee_src/n_floor.c
deleted file mode 100644
index 6c3888fee02..00000000000
--- a/lib/libm/noieee_src/n_floor.c
+++ /dev/null
@@ -1,174 +0,0 @@
-/*	$OpenBSD: n_floor.c,v 1.21 2016/09/12 19:47:02 guenther Exp $	*/
-/*	$NetBSD: n_floor.c,v 1.1 1995/10/10 23:36:48 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include <math.h>
-
-#include "mathimpl.h"
-
-static const double L = 36028797018963968.0E0;	/* 2**55 */
-static const float F = 8388608E0f;		/* 2**23 */
-
-/*
- * floor(x) := the largest integer no larger than x;
- * ceil(x) := -floor(-x), for all real x.
- *
- * Note: Inexact will be signaled if x is not an integer, as is
- *	customary for IEEE 754.  No other signal can be emitted.
- */
-double
-floor(double x)
-{
-	volatile double y;
-
-	if (isnan(x) ||	x >= L)		/* already an even integer */
-		return x;
-	else if (x < (double)0)
-		return -ceil(-x);
-	else {			/* now 0 <= x < L */
-		y = L+x;		/* destructive store must be forced */
-		y -= L;			/* an integer, and |x-y| < 1 */
-		return x < y ? y-(double)1 : y;
-	}
-}
-DEF_STD(floor);
-LDBL_CLONE(floor);
-
-double
-ceil(double x)
-{
-	volatile double y;
-
-	if (isnan(x) ||	x >= L)		/* already an even integer */
-		return x;
-	else if (x < (double)0)
-		return -floor(-x);
-	else {			/* now 0 <= x < L */
-		y = L+x;		/* destructive store must be forced */
-		y -= L;			/* an integer, and |x-y| < 1 */
-		return x > y ? y+(double)1 : y;
-	}
-}
-DEF_STD(ceil);
-LDBL_UNUSED_CLONE(ceil);
-
-float
-floorf(float x)
-{
-	volatile float y;
-
-	if (isnan(x) || x >= F)		/* already an even integer */
-		return x;
-	else if (x < (float)0)
-		return -ceilf(-x);
-	else {			/* now 0 <= x < F */
-		y = F+x;		/* destructive store must be forced */
-		y -= F;			/* an integer, and |x-y| < 1 */
-		return x < y ? y-(float)1 : y;
-	}
-}
-DEF_STD(floorf);
-
-float
-ceilf(float x)
-{
-	volatile float y;
-
-	if (isnan(x) || x >= F)		/* already an even integer */
-		return x;
-	else if (x < (float)0)
-		return -floorf(-x);
-	else {			/* now 0 <= x < F */
-		y = F+x;		/* destructive store must be forced */
-		y -= F;			/* an integer, and |x-y| < 1 */
-		return x > y ? y+(float)1 : y;
-	}
-}
-DEF_STD(ceilf);
-
-/*
- * algorithm for rint(x) in pseudo-pascal form ...
- *
- * real rint(x): real x;
- *	... delivers integer nearest x in direction of prevailing rounding
- *	... mode
- * const	L = (last consecutive integer)/2
- * 	  = 2**55; for VAX D
- * 	  = 2**52; for IEEE 754 Double
- * real	s,t;
- * begin
- * 	if isnan(x) then return x;		... NaN
- * 	if |x| >= L then return x;		... already an integer
- * 	s := copysign(L,x);
- * 	t := x + s;				... = (x+s) rounded to integer
- * 	return t - s
- * end;
- *
- * Note: Inexact will be signaled if x is not an integer, as is
- *	customary for IEEE 754.  No other signal can be emitted.
- */
-double
-rint(double x)
-{
-	double s;
-	volatile double t;
-	const double one = 1.0;
-
-	if (isnan(x))
-		return (x);
-
-	if (copysign(x, one) >= L)	/* already an integer */
-		return (x);
-
-	s = copysign(L,x);
-	t = x + s;				/* x+s rounded to integer */
-	return (t - s);
-}
-DEF_STD(rint);
-LDBL_CLONE(rint);
-
-float
-rintf(float x)
-{
-	float s;
-	volatile float t;
-	const float one = 1.0f;
-
-	if (isnan(x))
-		return (x);
-
-	if (copysignf(x, one) >= F)	/* already an integer */
-		return (x);
-
-	s = copysignf(F,x);
-	t = x + s;				/* x+s rounded to integer */
-	return (t - s);
-}
-DEF_STD(rintf);
diff --git a/lib/libm/noieee_src/n_fmod.c b/lib/libm/noieee_src/n_fmod.c
deleted file mode 100644
index f43d15dce0c..00000000000
--- a/lib/libm/noieee_src/n_fmod.c
+++ /dev/null
@@ -1,80 +0,0 @@
-/*	$OpenBSD: n_fmod.c,v 1.8 2016/09/12 04:39:47 guenther Exp $	*/
-/*	$NetBSD: n_fmod.c,v 1.1 1995/10/10 23:36:49 ragge Exp $	*/
-/*
- * Copyright (c) 1989, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-/* fmod.c
- *
- * SYNOPSIS
- *
- *    #include <math.h>
- *    double fmod(double x, double y)
- *
- * DESCRIPTION
- *
- *    The fmod function computes the floating-point remainder of x/y.
- *
- * RETURNS
- *
- *    The fmod function returns the value x-i*y, for some integer i
- * such that, if y is nonzero, the result has the same sign as x and
- * magnitude less than the magnitude of y.
- *
- * On a VAX or CCI,
- *
- *    fmod(x,0) traps/faults on floating-point divided-by-zero.
- *
- * On IEEE-754 conforming machines with "isnan()" primitive,
- *
- *    fmod(x,0), fmod(INF,y) are invalid operations and NaN is returned.
- *
- */
-
-double
-fmod(double x, double y)
-{
-	int ir,iy;
-	double r,w;
-
-	if (y == (double)0 || isnan(y) || !isfinite(x))
-	    return (x*y)/(x*y);
-
-	r = fabs(x);
-	y = fabs(y);
-	(void)frexp(y,&iy);
-	while (r >= y) {
-		(void)frexp(r,&ir);
-		w = ldexp(y,ir-iy);
-		r -= w <= r ? w : w*(double)0.5;
-	}
-	return x >= (double)0 ? r : -r;
-}
diff --git a/lib/libm/noieee_src/n_hypot.c b/lib/libm/noieee_src/n_hypot.c
deleted file mode 100644
index 6653e956844..00000000000
--- a/lib/libm/noieee_src/n_hypot.c
+++ /dev/null
@@ -1,195 +0,0 @@
-/*	$OpenBSD: n_hypot.c,v 1.6 2016/09/12 19:47:02 guenther Exp $	*/
-/*	$NetBSD: n_cabs.c,v 1.1 1995/10/10 23:36:39 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* HYPOT(X,Y)
- * RETURN THE SQUARE ROOT OF X^2 + Y^2  WHERE Z=X+iY
- * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 11/28/84;
- * REVISED BY K.C. NG, 7/12/85.
- *
- * Required system supported functions :
- *	copysign(x,y)
- *	isfinite(x)
- *	scalbn(x,N)
- *	sqrt(x)
- *
- * Method :
- *	1. replace x by |x| and y by |y|, and swap x and
- *	   y if y > x (hence x is never smaller than y).
- *	2. Hypot(x,y) is computed by:
- *	   Case I, x/y > 2
- *
- *				       y
- *		hypot = x + -----------------------------
- *			 		    2
- *			    sqrt ( 1 + [x/y]  )  +  x/y
- *
- *	   Case II, x/y <= 2
- *				                   y
- *		hypot = x + --------------------------------------------------
- *				          		     2
- *				     			[x/y]   -  2
- *			   (sqrt(2)+1) + (x-y)/y + -----------------------------
- *			 		    			  2
- *			    			  sqrt ( 1 + [x/y]  )  + sqrt(2)
- *
- *
- *
- * Special cases:
- *	hypot(x,y) is INF if x or y is +INF or -INF; else
- *	hypot(x,y) is NAN if x or y is NAN.
- *
- * Accuracy:
- * 	hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units
- *	in the last place). See Kahan's "Interval Arithmetic Options in the
- *	Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics
- *      1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate
- *	code follows in	comments.) In a test run with 500,000 random arguments
- *	on a VAX, the maximum observed error was .959 ulps.
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-static const double r2p1hi = 2.4142135623730950345E0;
-static const double r2p1lo = 1.4349369327986523769E-17;
-static const double sqrt2 = 1.4142135623730950622E0;
-
-float
-hypotf(float x, float y)
-{
-	return (float)hypot((double) x, (double) y);
-}
-
-double
-hypot(double x, double y)
-{
-	static const double zero=0, one=1,
-		      small=1.0E-18;	/* fl(1+small)==1 */
-	static const ibig=30;	/* fl(1+2**(2*ibig))==1 */
-	double t,r;
-	int exp;
-
-	if(isfinite(x))
-	    if(isfinite(y))
-	    {
-		x=copysign(x,one);
-		y=copysign(y,one);
-		if(y > x)
-		    { t=x; x=y; y=t; }
-		if(x == zero) return(zero);
-		if(y == zero) return(x);
-		exp= logb(x);
-		if (exp - (int)logb(y) > ibig) {
-			if (one + small >= 1.0)	/* raise inexact flag */
-				return(x);	/* return |x| */
-		}
-
-	    /* start computing sqrt(x^2 + y^2) */
-		r=x-y;
-		if(r>y) { 	/* x/y > 2 */
-		    r=x/y;
-		    r=r+sqrt(one+r*r); }
-		else {		/* 1 <= x/y <= 2 */
-		    r/=y; t=r*(r+2.0);
-		    r+=t/(sqrt2+sqrt(2.0+t));
-		    r+=r2p1lo; r+=r2p1hi; }
-
-		r=y/r;
-		return(x+r);
-
-	    }
-
-	    else if(isinf(y))		/* y is +-INF */
-		     return(copysign(y,one));
-	    else
-		     return(y);		/* y is NaN and x is finite */
-
-	else if(isinf(x))		/* x is +-INF */
-	         return (copysign(x,one));
-	else if(isfinite(y))
-	         return(x);		/* x is NaN, y is finite */
-	else if (isnan(y))
-		return (y);
-	else return(copysign(y,one));	/* y is INF */
-}
-DEF_STD(hypot);
-LDBL_CLONE(hypot);
-
-/* A faster but less accurate version of cabs(x,y) */
-#if 0
-double
-hypot(double x, double y)
-{
-	static const double zero=0, one=1;
-		      small=1.0E-18;	/* fl(1+small)==1 */
-	static const ibig=30;	/* fl(1+2**(2*ibig))==1 */
-	double temp;
-	int exp;
-
-	if(isfinite(x))
-	    if(isfinite(y))
-	    {
-		x=copysign(x,one);
-		y=copysign(y,one);
-		if(y > x)
-		    { temp=x; x=y; y=temp; }
-		if(x == zero) return(zero);
-		if(y == zero) return(x);
-		exp= logb(x);
-		x=scalbn(x,-exp);
-		if (exp - (int)logb(y) > ibig) {
-			if (one + small >= 1.0)		/* raise inexact flag */
-				return(scalbn(x,exp));	/* return |x| */
-		}
-		else y=scalbn(y,-exp);
-		return(scalbn(sqrt(x*x+y*y),exp));
-	    }
-
-	    else if(isinf(y))		/* y is +-INF */
-		     return(copysign(y,one));
-	    else
-		     return(y);		/* y is NaN and x is finite */
-
-	else if(isinf(x))		/* x is +-INF */
-	         return (copysign(x,one));
-	else if(isfinite(y))
-	         return(x);		/* x is NaN, y is finite */
-	else if(isnan(y)) return(y);	/* x and y is NaN */
-	else return(copysign(y,one));	/* y is INF */
-}
-#endif
diff --git a/lib/libm/noieee_src/n_j0.c b/lib/libm/noieee_src/n_j0.c
deleted file mode 100644
index 6274b17cde2..00000000000
--- a/lib/libm/noieee_src/n_j0.c
+++ /dev/null
@@ -1,431 +0,0 @@
-/*	$OpenBSD: n_j0.c,v 1.8 2016/09/12 04:39:47 guenther Exp $	*/
-/*	$NetBSD: n_j0.c,v 1.1 1995/10/10 23:36:52 ragge Exp $	*/
-/*-
- * Copyright (c) 1992, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/*
- * 16 December 1992
- * Minor modifications by Peter McIlroy to adapt non-IEEE architecture.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1992 by Sun Microsystems, Inc.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- *
- * ******************* WARNING ********************
- * This is an alpha version of SunPro's FDLIBM (Freely
- * Distributable Math Library) for IEEE double precision
- * arithmetic. FDLIBM is a basic math library written
- * in C that runs on machines that conform to IEEE
- * Standard 754/854. This alpha version is distributed
- * for testing purpose. Those who use this software
- * should report any bugs to
- *
- *		fdlibm-comments@sunpro.eng.sun.com
- *
- * -- K.C. Ng, Oct 12, 1992
- * ************************************************
- */
-
-/* double j0(double x), y0(double 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 - ...
- *	2. Reduce x to |x| since j0(x)=j0(-x),  and
- *	   for x in (0,2)
- *		j0(x) = 1-z/4+ z^2*R0/S0,  where z = x*x;
- *	   (precision:  |j0-1+z/4-z^2R0/S0 |<2**-63.67 )
- *	   for x in (2,inf)
- * 		j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
- * 	   where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
- *	   as follow:
- *		cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
- *			= 1/sqrt(2) * (cos(x) + sin(x))
- *		sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
- *			= 1/sqrt(2) * (sin(x) - cos(x))
- * 	   (To avoid cancellation, use
- *		sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
- * 	    to compute the worse one.)
- *
- *	3 Special cases
- *		j0(nan)= nan
- *		j0(0) = 1
- *		j0(inf) = 0
- *
- * Method -- y0(x):
- *	1. For x<2.
- *	   Since
- *		y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
- *	   therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
- *	   We use the following function to approximate y0,
- *		y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
- *	   where
- *		U(z) = u0 + u1*z + ... + u6*z^6
- *		V(z) = 1  + v1*z + ... + v4*z^4
- *	   with absolute approximation error bounded by 2**-72.
- *	   Note: For tiny x, U/V = u0 and j0(x)~1, hence
- *		y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
- *	2. For x>=2.
- * 		y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
- * 	   where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
- *	   by the method mentioned above.
- *	3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
- */
-
-#include <math.h>
-#include <float.h>
-#include <errno.h>
-
-#if defined(__vax__)
-#define _IEEE	0
-#else
-#define _IEEE	1
-#define infnan(x) (0.0)
-#endif
-
-static double pzero(double), qzero(double);
-
-static double
-huge 	= 1e300,
-zero    = 0.0,
-one	= 1.0,
-invsqrtpi= 5.641895835477562869480794515607725858441e-0001,
-tpi	= 0.636619772367581343075535053490057448,
- 		/* R0/S0 on [0, 2.00] */
-r02 =   1.562499999999999408594634421055018003102e-0002,
-r03 =  -1.899792942388547334476601771991800712355e-0004,
-r04 =   1.829540495327006565964161150603950916854e-0006,
-r05 =  -4.618326885321032060803075217804816988758e-0009,
-s01 =   1.561910294648900170180789369288114642057e-0002,
-s02 =   1.169267846633374484918570613449245536323e-0004,
-s03 =   5.135465502073181376284426245689510134134e-0007,
-s04 =   1.166140033337900097836930825478674320464e-0009;
-
-double
-j0(double x)
-{
-	double z, s,c,ss,cc,r,u,v;
-
-	if (!isfinite(x))
-		if (_IEEE) return one/(x*x);
-		else return (0);
-	x = fabs(x);
-	if (x >= 2.0) {	/* |x| >= 2.0 */
-		s = sin(x);
-		c = cos(x);
-		ss = s-c;
-		cc = s+c;
-		if (x < .5 * DBL_MAX) {  /* make sure x+x not overflow */
-		    z = -cos(x+x);
-		    if ((s*c)<zero) cc = z/ss;
-		    else 	    ss = z/cc;
-		}
-	/*
-	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
-	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
-	 */
-		if (_IEEE && x> 6.80564733841876927e+38) /* 2^129 */
-			z = (invsqrtpi*cc)/sqrt(x);
-		else {
-		    u = pzero(x); v = qzero(x);
-		    z = invsqrtpi*(u*cc-v*ss)/sqrt(x);
-		}
-		return z;
-	}
-	if (x < 1.220703125e-004) {		   /* |x| < 2**-13 */
-	    if (huge+x > one) {			   /* raise inexact if x != 0 */
-	        if (x < 7.450580596923828125e-009) /* |x|<2**-27 */
-			return one;
-	        else return (one - 0.25*x*x);
-	    }
-	}
-	z = x*x;
-	r =  z*(r02+z*(r03+z*(r04+z*r05)));
-	s =  one+z*(s01+z*(s02+z*(s03+z*s04)));
-	if (x < one) {			/* |x| < 1.00 */
-	    return (one + z*(-0.25+(r/s)));
-	} else {
-	    u = 0.5*x;
-	    return ((one+u)*(one-u)+z*(r/s));
-	}
-}
-
-static double
-u00 =  -7.380429510868722527422411862872999615628e-0002,
-u01 =   1.766664525091811069896442906220827182707e-0001,
-u02 =  -1.381856719455968955440002438182885835344e-0002,
-u03 =   3.474534320936836562092566861515617053954e-0004,
-u04 =  -3.814070537243641752631729276103284491172e-0006,
-u05 =   1.955901370350229170025509706510038090009e-0008,
-u06 =  -3.982051941321034108350630097330144576337e-0011,
-v01 =   1.273048348341237002944554656529224780561e-0002,
-v02 =   7.600686273503532807462101309675806839635e-0005,
-v03 =   2.591508518404578033173189144579208685163e-0007,
-v04 =   4.411103113326754838596529339004302243157e-0010;
-
-double
-y0(double x)
-{
-	double z, s, c, ss, cc, u, v;
-    /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
-	if (!isfinite(x))
-		if (_IEEE)
-			return (one/(x+x*x));
-		else
-			return (0);
-        if (x == 0)
-		if (_IEEE)	return (-one/zero);
-		else		return(infnan(-ERANGE));
-        if (x<0)
-		if (_IEEE)	return (zero/zero);
-		else		return (infnan(EDOM));
-        if (x >= 2.00) {	/* |x| >= 2.0 */
-        /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
-         * where x0 = x-pi/4
-         *      Better formula:
-         *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
-         *                      =  1/sqrt(2) * (sin(x) + cos(x))
-         *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
-         *                      =  1/sqrt(2) * (sin(x) - cos(x))
-         * To avoid cancellation, use
-         *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
-         * to compute the worse one.
-         */
-                s = sin(x);
-                c = cos(x);
-                ss = s-c;
-                cc = s+c;
-	/*
-	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
-	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
-	 */
-                if (x < .5 * DBL_MAX) {  /* make sure x+x not overflow */
-                    z = -cos(x+x);
-                    if ((s*c)<zero) cc = z/ss;
-                    else            ss = z/cc;
-                }
-                if (_IEEE && x > 6.80564733841876927e+38) /* > 2^129 */
-			z = (invsqrtpi*ss)/sqrt(x);
-                else {
-                    u = pzero(x); v = qzero(x);
-                    z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
-                }
-                return z;
-	}
-	if (x <= 7.450580596923828125e-009) {		/* x < 2**-27 */
-	    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*(j0(x)*log(x)));
-}
-
-/* The asymptotic expansions of pzero is
- *	1 - 9/128 s^2 + 11025/98304 s^4 - ...,	where s = 1/x.
- * For x >= 2, We approximate pzero by
- * 	pzero(x) = 1 + (R/S)
- * where  R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
- * 	  S = 1 + ps0*s^2 + ... + ps4*s^10
- * and
- *	| pzero(x)-1-R/S | <= 2  ** ( -60.26)
- */
-static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-   0.0,
-  -7.031249999999003994151563066182798210142e-0002,
-  -8.081670412753498508883963849859423939871e+0000,
-  -2.570631056797048755890526455854482662510e+0002,
-  -2.485216410094288379417154382189125598962e+0003,
-  -5.253043804907295692946647153614119665649e+0003,
-};
-static double ps8[5] = {
-   1.165343646196681758075176077627332052048e+0002,
-   3.833744753641218451213253490882686307027e+0003,
-   4.059785726484725470626341023967186966531e+0004,
-   1.167529725643759169416844015694440325519e+0005,
-   4.762772841467309430100106254805711722972e+0004,
-};
-
-static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-  -1.141254646918944974922813501362824060117e-0011,
-  -7.031249408735992804117367183001996028304e-0002,
-  -4.159610644705877925119684455252125760478e+0000,
-  -6.767476522651671942610538094335912346253e+0001,
-  -3.312312996491729755731871867397057689078e+0002,
-  -3.464333883656048910814187305901796723256e+0002,
-};
-static double ps5[5] = {
-   6.075393826923003305967637195319271932944e+0001,
-   1.051252305957045869801410979087427910437e+0003,
-   5.978970943338558182743915287887408780344e+0003,
-   9.625445143577745335793221135208591603029e+0003,
-   2.406058159229391070820491174867406875471e+0003,
-};
-
-static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-  -2.547046017719519317420607587742992297519e-0009,
-  -7.031196163814817199050629727406231152464e-0002,
-  -2.409032215495295917537157371488126555072e+0000,
-  -2.196597747348830936268718293366935843223e+0001,
-  -5.807917047017375458527187341817239891940e+0001,
-  -3.144794705948885090518775074177485744176e+0001,
-};
-static double ps3[5] = {
-   3.585603380552097167919946472266854507059e+0001,
-   3.615139830503038919981567245265266294189e+0002,
-   1.193607837921115243628631691509851364715e+0003,
-   1.127996798569074250675414186814529958010e+0003,
-   1.735809308133357510239737333055228118910e+0002,
-};
-
-static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-  -8.875343330325263874525704514800809730145e-0008,
-  -7.030309954836247756556445443331044338352e-0002,
-  -1.450738467809529910662233622603401167409e+0000,
-  -7.635696138235277739186371273434739292491e+0000,
-  -1.119316688603567398846655082201614524650e+0001,
-  -3.233645793513353260006821113608134669030e+0000,
-};
-static double ps2[5] = {
-   2.222029975320888079364901247548798910952e+0001,
-   1.362067942182152109590340823043813120940e+0002,
-   2.704702786580835044524562897256790293238e+0002,
-   1.538753942083203315263554770476850028583e+0002,
-   1.465761769482561965099880599279699314477e+0001,
-};
-
-static double pzero(double x)
-{
-	double *p,*q,z,r,s;
-	if (x >= 8.00)			   {p = pr8; q= ps8;}
-	else if (x >= 4.54545211791992188) {p = pr5; q= ps5;}
-	else if (x >= 2.85714149475097656) {p = pr3; q= ps3;}
-	else if (x >= 2.00)		   {p = pr2; q= ps2;}
-	z = one/(x*x);
-	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
-	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
-	return one+ r/s;
-}
-
-
-/* For x >= 8, the asymptotic expansions of qzero is
- *	-1/8 s + 75/1024 s^3 - ..., where s = 1/x.
- * We approximate pzero by
- * 	qzero(x) = s*(-1.25 + (R/S))
- * where  R = qr0 + qr1*s^2 + qr2*s^4 + ... + qr5*s^10
- * 	  S = 1 + qs0*s^2 + ... + qs5*s^12
- * and
- *	| qzero(x)/s +1.25-R/S | <= 2  ** ( -61.22)
- */
-static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-   0.0,
-   7.324218749999350414479738504551775297096e-0002,
-   1.176820646822526933903301695932765232456e+0001,
-   5.576733802564018422407734683549251364365e+0002,
-   8.859197207564685717547076568608235802317e+0003,
-   3.701462677768878501173055581933725704809e+0004,
-};
-static double qs8[6] = {
-   1.637760268956898345680262381842235272369e+0002,
-   8.098344946564498460163123708054674227492e+0003,
-   1.425382914191204905277585267143216379136e+0005,
-   8.033092571195144136565231198526081387047e+0005,
-   8.405015798190605130722042369969184811488e+0005,
-  -3.438992935378666373204500729736454421006e+0005,
-};
-
-static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-   1.840859635945155400568380711372759921179e-0011,
-   7.324217666126847411304688081129741939255e-0002,
-   5.835635089620569401157245917610984757296e+0000,
-   1.351115772864498375785526599119895942361e+0002,
-   1.027243765961641042977177679021711341529e+0003,
-   1.989977858646053872589042328678602481924e+0003,
-};
-static double qs5[6] = {
-   8.277661022365377058749454444343415524509e+0001,
-   2.077814164213929827140178285401017305309e+0003,
-   1.884728877857180787101956800212453218179e+0004,
-   5.675111228949473657576693406600265778689e+0004,
-   3.597675384251145011342454247417399490174e+0004,
-  -5.354342756019447546671440667961399442388e+0003,
-};
-
-static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-   4.377410140897386263955149197672576223054e-0009,
-   7.324111800429115152536250525131924283018e-0002,
-   3.344231375161707158666412987337679317358e+0000,
-   4.262184407454126175974453269277100206290e+0001,
-   1.708080913405656078640701512007621675724e+0002,
-   1.667339486966511691019925923456050558293e+0002,
-};
-static double qs3[6] = {
-   4.875887297245871932865584382810260676713e+0001,
-   7.096892210566060535416958362640184894280e+0002,
-   3.704148226201113687434290319905207398682e+0003,
-   6.460425167525689088321109036469797462086e+0003,
-   2.516333689203689683999196167394889715078e+0003,
-  -1.492474518361563818275130131510339371048e+0002,
-};
-
-static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-   1.504444448869832780257436041633206366087e-0007,
-   7.322342659630792930894554535717104926902e-0002,
-   1.998191740938159956838594407540292600331e+0000,
-   1.449560293478857407645853071687125850962e+0001,
-   3.166623175047815297062638132537957315395e+0001,
-   1.625270757109292688799540258329430963726e+0001,
-};
-static double qs2[6] = {
-   3.036558483552191922522729838478169383969e+0001,
-   2.693481186080498724211751445725708524507e+0002,
-   8.447837575953201460013136756723746023736e+0002,
-   8.829358451124885811233995083187666981299e+0002,
-   2.126663885117988324180482985363624996652e+0002,
-  -5.310954938826669402431816125780738924463e+0000,
-};
-
-static double qzero(double x)
-{
-	double *p,*q, s,r,z;
-	if (x >= 8.00)			   {p = qr8; q= qs8;}
-	else if (x >= 4.54545211791992188) {p = qr5; q= qs5;}
-	else if (x >= 2.85714149475097656) {p = qr3; q= qs3;}
-	else if (x >= 2.00)		   {p = qr2; q= qs2;}
-	z = one/(x*x);
-	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
-	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
-	return (-.125 + r/s)/x;
-}
diff --git a/lib/libm/noieee_src/n_j1.c b/lib/libm/noieee_src/n_j1.c
deleted file mode 100644
index 0bb2181b68b..00000000000
--- a/lib/libm/noieee_src/n_j1.c
+++ /dev/null
@@ -1,439 +0,0 @@
-/*	$OpenBSD: n_j1.c,v 1.8 2016/09/12 04:39:47 guenther Exp $	*/
-/*	$NetBSD: n_j1.c,v 1.1 1995/10/10 23:36:53 ragge Exp $	*/
-/*-
- * Copyright (c) 1992, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/*
- * 16 December 1992
- * Minor modifications by Peter McIlroy to adapt non-IEEE architecture.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1992 by Sun Microsystems, Inc.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- *
- * ******************* WARNING ********************
- * This is an alpha version of SunPro's FDLIBM (Freely
- * Distributable Math Library) for IEEE double precision
- * arithmetic. FDLIBM is a basic math library written
- * in C that runs on machines that conform to IEEE
- * Standard 754/854. This alpha version is distributed
- * for testing purpose. Those who use this software
- * should report any bugs to
- *
- *		fdlibm-comments@sunpro.eng.sun.com
- *
- * -- K.C. Ng, Oct 12, 1992
- * ************************************************
- */
-
-/* double j1(double x), y1(double 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 - ...
- *	2. Reduce x to |x| since j1(x)=-j1(-x),  and
- *	   for x in (0,2)
- *		j1(x) = x/2 + x*z*R0/S0,  where z = x*x;
- *	   (precision:  |j1/x - 1/2 - R0/S0 |<2**-61.51 )
- *	   for x in (2,inf)
- * 		j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
- * 		y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
- * 	   where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
- *	   as follows:
- *		cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
- *			=  1/sqrt(2) * (sin(x) - cos(x))
- *		sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
- *			= -1/sqrt(2) * (sin(x) + cos(x))
- * 	   (To avoid cancellation, use
- *		sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
- * 	    to compute the worse one.)
- *
- *	3 Special cases
- *		j1(nan)= nan
- *		j1(0) = 0
- *		j1(inf) = 0
- *
- * Method -- y1(x):
- *	1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
- *	2. For x<2.
- *	   Since
- *		y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
- *	   therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
- *	   We use the following function to approximate y1,
- *		y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
- *	   where for x in [0,2] (abs err less than 2**-65.89)
- *		U(z) = u0 + u1*z + ... + u4*z^4
- *		V(z) = 1  + v1*z + ... + v5*z^5
- *	   Note: For tiny x, 1/x dominate y1 and hence
- *		y1(tiny) = -2/pi/tiny, (choose tiny<2**-54)
- *	3. For x>=2.
- * 		y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
- * 	   where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
- *	   by method mentioned above.
- */
-
-#include <math.h>
-#include <float.h>
-#include <errno.h>
-
-#if defined(__vax__)
-#define _IEEE	0
-#else
-#define _IEEE	1
-#define infnan(x) (0.0)
-#endif
-
-static double pone(), qone();
-
-static double
-huge    = 1e300,
-zero    = 0.0,
-one	= 1.0,
-invsqrtpi= 5.641895835477562869480794515607725858441e-0001,
-tpi	= 0.636619772367581343075535053490057448,
-
-	/* R0/S0 on [0,2] */
-r00 =  -6.250000000000000020842322918309200910191e-0002,
-r01 =   1.407056669551897148204830386691427791200e-0003,
-r02 =  -1.599556310840356073980727783817809847071e-0005,
-r03 =   4.967279996095844750387702652791615403527e-0008,
-s01 =   1.915375995383634614394860200531091839635e-0002,
-s02 =   1.859467855886309024045655476348872850396e-0004,
-s03 =   1.177184640426236767593432585906758230822e-0006,
-s04 =   5.046362570762170559046714468225101016915e-0009,
-s05 =   1.235422744261379203512624973117299248281e-0011;
-
-#define two_129	6.80564733841876926e+038	/* 2^129 */
-#define two_m54	5.55111512312578270e-017	/* 2^-54 */
-
-double
-j1(double x)
-{
-	double z, s,c,ss,cc,r,u,v,y;
-	y = fabs(x);
-	if (!isfinite(x))		/* Inf or NaN */
-		if (isnan(x))
-			return(x);
-		else
-			return (copysign(x, zero));
-	y = fabs(x);
-	if (y >= 2)			/* |x| >= 2.0 */
-	{
-		s = sin(y);
-		c = cos(y);
-		ss = -s-c;
-		cc = s-c;
-		if (y < .5*DBL_MAX) {  	/* make sure y+y not overflow */
-		    z = cos(y+y);
-		    if ((s*c)<zero) cc = z/ss;
-		    else 	    ss = z/cc;
-		}
-	/*
-	 * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
-	 * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
-	 */
-#if !defined(__vax__)
-		if (y > two_129)	 /* x > 2^129 */
-			z = (invsqrtpi*cc)/sqrt(y);
-		else
-#endif /* !defined(__vax__) */
-		{
-		    u = pone(y); v = qone(y);
-		    z = invsqrtpi*(u*cc-v*ss)/sqrt(y);
-		}
-		if (x < 0) return -z;
-		else  	 return  z;
-	}
-	if (y < 7.450580596923828125e-009) {	/* |x|<2**-27 */
-	    if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */
-	}
-	z = x*x;
-	r =  z*(r00+z*(r01+z*(r02+z*r03)));
-	s =  one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
-	r *= x;
-	return (x*0.5+r/s);
-}
-
-static double u0[5] = {
-  -1.960570906462389484206891092512047539632e-0001,
-   5.044387166398112572026169863174882070274e-0002,
-  -1.912568958757635383926261729464141209569e-0003,
-   2.352526005616105109577368905595045204577e-0005,
-   -9.190991580398788465315411784276789663849e-0008,
-};
-static double v0[5] = {
-   1.991673182366499064031901734535479833387e-0002,
-   2.025525810251351806268483867032781294682e-0004,
-   1.356088010975162198085369545564475416398e-0006,
-   6.227414523646214811803898435084697863445e-0009,
-   1.665592462079920695971450872592458916421e-0011,
-};
-
-double
-y1(double x)
-{
-	double z, s, c, ss, cc, u, v;
-    /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
-	if (!isfinite(x))
-		if (x < 0)
-			return(zero/zero);
-		else if (x > 0)
-			return (0);
-		else
-			return(x);
-	if (x <= 0) {
-		if (_IEEE && x == 0) return -one/zero;
-		else if(x == 0) return(infnan(-ERANGE));
-		else if(_IEEE) return (zero/zero);
-		else return(infnan(EDOM));
-	}
-        if (x >= 2)			 /* |x| >= 2.0 */
-	{
-                s = sin(x);
-                c = cos(x);
-                ss = -s-c;
-                cc = s-c;
-		if (x < .5 * DBL_MAX)	/* make sure x+x not overflow */
-		{
-                    z = cos(x+x);
-                    if ((s*c)>zero) cc = z/ss;
-                    else            ss = z/cc;
-                }
-        /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
-         * where x0 = x-3pi/4
-         *      Better formula:
-         *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
-         *                      =  1/sqrt(2) * (sin(x) - cos(x))
-         *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
-         *                      = -1/sqrt(2) * (cos(x) + sin(x))
-         * To avoid cancellation, use
-         *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
-         * to compute the worse one.
-         */
-                if (_IEEE && x>two_129)
-			z = (invsqrtpi*ss)/sqrt(x);
-                else {
-                    u = pone(x); v = qone(x);
-                    z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
-                }
-                return z;
-        }
-        if (x <= two_m54) {    /* x < 2**-54 */
-            return (-tpi/x);
-        }
-        z = x*x;
-        u = u0[0]+z*(u0[1]+z*(u0[2]+z*(u0[3]+z*u0[4])));
-        v = one+z*(v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4]))));
-        return (x*(u/v) + tpi*(j1(x)*log(x)-one/x));
-}
-
-/* For x >= 8, the asymptotic expansions of pone is
- *	1 + 15/128 s^2 - 4725/2^15 s^4 - ...,	where s = 1/x.
- * We approximate pone by
- * 	pone(x) = 1 + (R/S)
- * where  R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
- * 	  S = 1 + ps0*s^2 + ... + ps4*s^10
- * and
- *	| pone(x)-1-R/S | <= 2  ** ( -60.06)
- */
-
-static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-   0.0,
-   1.171874999999886486643746274751925399540e-0001,
-   1.323948065930735690925827997575471527252e+0001,
-   4.120518543073785433325860184116512799375e+0002,
-   3.874745389139605254931106878336700275601e+0003,
-   7.914479540318917214253998253147871806507e+0003,
-};
-static double ps8[5] = {
-   1.142073703756784104235066368252692471887e+0002,
-   3.650930834208534511135396060708677099382e+0003,
-   3.695620602690334708579444954937638371808e+0004,
-   9.760279359349508334916300080109196824151e+0004,
-   3.080427206278887984185421142572315054499e+0004,
-};
-
-static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-   1.319905195562435287967533851581013807103e-0011,
-   1.171874931906140985709584817065144884218e-0001,
-   6.802751278684328781830052995333841452280e+0000,
-   1.083081829901891089952869437126160568246e+0002,
-   5.176361395331997166796512844100442096318e+0002,
-   5.287152013633375676874794230748055786553e+0002,
-};
-static double ps5[5] = {
-   5.928059872211313557747989128353699746120e+0001,
-   9.914014187336144114070148769222018425781e+0002,
-   5.353266952914879348427003712029704477451e+0003,
-   7.844690317495512717451367787640014588422e+0003,
-   1.504046888103610723953792002716816255382e+0003,
-};
-
-static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-   3.025039161373736032825049903408701962756e-0009,
-   1.171868655672535980750284752227495879921e-0001,
-   3.932977500333156527232725812363183251138e+0000,
-   3.511940355916369600741054592597098912682e+0001,
-   9.105501107507812029367749771053045219094e+0001,
-   4.855906851973649494139275085628195457113e+0001,
-};
-static double ps3[5] = {
-   3.479130950012515114598605916318694946754e+0001,
-   3.367624587478257581844639171605788622549e+0002,
-   1.046871399757751279180649307467612538415e+0003,
-   8.908113463982564638443204408234739237639e+0002,
-   1.037879324396392739952487012284401031859e+0002,
-};
-
-static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-   1.077108301068737449490056513753865482831e-0007,
-   1.171762194626833490512746348050035171545e-0001,
-   2.368514966676087902251125130227221462134e+0000,
-   1.224261091482612280835153832574115951447e+0001,
-   1.769397112716877301904532320376586509782e+0001,
-   5.073523125888185399030700509321145995160e+0000,
-};
-static double ps2[5] = {
-   2.143648593638214170243114358933327983793e+0001,
-   1.252902271684027493309211410842525120355e+0002,
-   2.322764690571628159027850677565128301361e+0002,
-   1.176793732871470939654351793502076106651e+0002,
-   8.364638933716182492500902115164881195742e+0000,
-};
-
-static double pone(double x)
-{
-	double *p,*q,z,r,s;
-	if (x >= 8.0) 			   {p = pr8; q= ps8;}
-	else if (x >= 4.54545211791992188) {p = pr5; q= ps5;}
-	else if (x >= 2.85714149475097656) {p = pr3; q= ps3;}
-	else /* if (x >= 2.0) */	   {p = pr2; q= ps2;}
-	z = one/(x*x);
-	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
-	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
-	return (one + r/s);
-}
-
-
-/* For x >= 8, the asymptotic expansions of qone is
- *	3/8 s - 105/1024 s^3 - ..., where s = 1/x.
- * We approximate pone by
- * 	qone(x) = s*(0.375 + (R/S))
- * where  R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
- * 	  S = 1 + qs1*s^2 + ... + qs6*s^12
- * and
- *	| qone(x)/s -0.375-R/S | <= 2  ** ( -61.13)
- */
-
-static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-   0.0,
-  -1.025390624999927207385863635575804210817e-0001,
-  -1.627175345445899724355852152103771510209e+0001,
-  -7.596017225139501519843072766973047217159e+0002,
-  -1.184980667024295901645301570813228628541e+0004,
-  -4.843851242857503225866761992518949647041e+0004,
-};
-static double qs8[6] = {
-   1.613953697007229231029079421446916397904e+0002,
-   7.825385999233484705298782500926834217525e+0003,
-   1.338753362872495800748094112937868089032e+0005,
-   7.196577236832409151461363171617204036929e+0005,
-   6.666012326177764020898162762642290294625e+0005,
-  -2.944902643038346618211973470809456636830e+0005,
-};
-
-static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-  -2.089799311417640889742251585097264715678e-0011,
-  -1.025390502413754195402736294609692303708e-0001,
-  -8.056448281239359746193011295417408828404e+0000,
-  -1.836696074748883785606784430098756513222e+0002,
-  -1.373193760655081612991329358017247355921e+0003,
-  -2.612444404532156676659706427295870995743e+0003,
-};
-static double qs5[6] = {
-   8.127655013843357670881559763225310973118e+0001,
-   1.991798734604859732508048816860471197220e+0003,
-   1.746848519249089131627491835267411777366e+0004,
-   4.985142709103522808438758919150738000353e+0004,
-   2.794807516389181249227113445299675335543e+0004,
-  -4.719183547951285076111596613593553911065e+0003,
-};
-
-static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-  -5.078312264617665927595954813341838734288e-0009,
-  -1.025378298208370901410560259001035577681e-0001,
-  -4.610115811394734131557983832055607679242e+0000,
-  -5.784722165627836421815348508816936196402e+0001,
-  -2.282445407376317023842545937526967035712e+0002,
-  -2.192101284789093123936441805496580237676e+0002,
-};
-static double qs3[6] = {
-   4.766515503237295155392317984171640809318e+0001,
-   6.738651126766996691330687210949984203167e+0002,
-   3.380152866795263466426219644231687474174e+0003,
-   5.547729097207227642358288160210745890345e+0003,
-   1.903119193388108072238947732674639066045e+0003,
-  -1.352011914443073322978097159157678748982e+0002,
-};
-
-static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-  -1.783817275109588656126772316921194887979e-0007,
-  -1.025170426079855506812435356168903694433e-0001,
-  -2.752205682781874520495702498875020485552e+0000,
-  -1.966361626437037351076756351268110418862e+0001,
-  -4.232531333728305108194363846333841480336e+0001,
-  -2.137192117037040574661406572497288723430e+0001,
-};
-static double qs2[6] = {
-   2.953336290605238495019307530224241335502e+0001,
-   2.529815499821905343698811319455305266409e+0002,
-   7.575028348686454070022561120722815892346e+0002,
-   7.393932053204672479746835719678434981599e+0002,
-   1.559490033366661142496448853793707126179e+0002,
-  -4.959498988226281813825263003231704397158e+0000,
-};
-
-static double qone(double x)
-{
-	double *p,*q, s,r,z;
-	if (x >= 8.0)			   {p = qr8; q= qs8;}
-	else if (x >= 4.54545211791992188) {p = qr5; q= qs5;}
-	else if (x >= 2.85714149475097656) {p = qr3; q= qs3;}
-	else /* if (x >= 2.0) */	   {p = qr2; q= qs2;}
-	z = one/(x*x);
-	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
-	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
-	return (.375 + r/s)/x;
-}
diff --git a/lib/libm/noieee_src/n_jn.c b/lib/libm/noieee_src/n_jn.c
deleted file mode 100644
index d8c6111f54f..00000000000
--- a/lib/libm/noieee_src/n_jn.c
+++ /dev/null
@@ -1,306 +0,0 @@
-/*	$OpenBSD: n_jn.c,v 1.8 2016/09/12 04:39:47 guenther Exp $	*/
-/*	$NetBSD: n_jn.c,v 1.1 1995/10/10 23:36:54 ragge Exp $	*/
-/*-
- * Copyright (c) 1992, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/*
- * 16 December 1992
- * Minor modifications by Peter McIlroy to adapt non-IEEE architecture.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1992 by Sun Microsystems, Inc.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- *
- * ******************* WARNING ********************
- * This is an alpha version of SunPro's FDLIBM (Freely
- * Distributable Math Library) for IEEE double precision
- * arithmetic. FDLIBM is a basic math library written
- * in C that runs on machines that conform to IEEE
- * Standard 754/854. This alpha version is distributed
- * for testing purpose. Those who use this software
- * should report any bugs to
- *
- *		fdlibm-comments@sunpro.eng.sun.com
- *
- * -- K.C. Ng, Oct 12, 1992
- * ************************************************
- */
-
-/*
- * jn(int n, double x), yn(int n, double x)
- * floating point Bessel's function of the 1st and 2nd kind
- * of order n
- *
- * Special cases:
- *	y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
- *	y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
- * Note 2. About jn(n,x), yn(n,x)
- *	For n=0, j0(x) is called,
- *	for n=1, j1(x) is called,
- *	for n<x, forward recursion us used starting
- *	from values of j0(x) and j1(x).
- *	for n>x, a continued fraction approximation to
- *	j(n,x)/j(n-1,x) is evaluated and then backward
- *	recursion is used starting from a supposed value
- *	for j(n,x). The resulting value of j(0,x) is
- *	compared with the actual value to correct the
- *	supposed value of j(n,x).
- *
- *	yn(n,x) is similar in all respects, except
- *	that forward recursion is used for all
- *	values of n>1.
- *
- */
-
-#include <math.h>
-#include <float.h>
-#include <errno.h>
-
-#if defined(__vax__)
-#define _IEEE	0
-#else
-#define _IEEE	1
-#define infnan(x) (0.0)
-#endif
-
-static double
-invsqrtpi= 5.641895835477562869480794515607725858441e-0001,
-two  = 2.0,
-zero = 0.0,
-one  = 1.0;
-
-double
-jn(int n, double x)
-{
-	int i, sgn;
-	double a, b, temp;
-	double z, w;
-
-    /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
-     * Thus, J(-n,x) = J(n,-x)
-     */
-    /* if J(n,NaN) is NaN */
-	if (_IEEE && isnan(x)) return x+x;
-	if (n<0){
-		n = -n;
-		x = -x;
-	}
-	if (n==0) return(j0(x));
-	if (n==1) return(j1(x));
-	sgn = (n&1)&(x < zero);		/* even n -- 0, odd n -- sign(x) */
-	x = fabs(x);
-	if (x == 0 || !isfinite (x)) 	/* if x is 0 or inf */
-	    b = zero;
-	else if ((double) n <= x) {
-			/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
-	    if (_IEEE && x >= 8.148143905337944345e+090) {
-					/* x >= 2**302 */
-    /* (x >> n**2)
-     *	    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-     *	    Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-     *	    Let s=sin(x), c=cos(x),
-     *		xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
-     *
-     *		   n	sin(xn)*sqt2	cos(xn)*sqt2
-     *		----------------------------------
-     *		   0	 s-c		 c+s
-     *		   1	-s-c 		-c+s
-     *		   2	-s+c		-c-s
-     *		   3	 s+c		 c-s
-     */
-		switch(n&3) {
-		    case 0: temp =  cos(x)+sin(x); break;
-		    case 1: temp = -cos(x)+sin(x); break;
-		    case 2: temp = -cos(x)-sin(x); break;
-		    case 3: temp =  cos(x)-sin(x); break;
-		}
-		b = invsqrtpi*temp/sqrt(x);
-	    } else {
-	        a = j0(x);
-	        b = j1(x);
-	        for(i=1;i<n;i++){
-		    temp = b;
-		    b = b*((double)(i+i)/x) - a; /* avoid underflow */
-		    a = temp;
-	        }
-	    }
-	} else {
-	    if (x < 1.86264514923095703125e-009) { /* x < 2**-29 */
-    /* x is tiny, return the first Taylor expansion of J(n,x)
-     * J(n,x) = 1/n!*(x/2)^n  - ...
-     */
-		if (n > 33)	/* underflow */
-		    b = zero;
-		else {
-		    temp = x*0.5; b = temp;
-		    for (a=one,i=2;i<=n;i++) {
-			a *= (double)i;		/* a = n! */
-			b *= temp;		/* b = (x/2)^n */
-		    }
-		    b = b/a;
-		}
-	    } else {
-		/* use backward recurrence */
-		/* 			x      x^2      x^2
-		 *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
-		 *			2n  - 2(n+1) - 2(n+2)
-		 *
-		 * 			1      1        1
-		 *  (for large x)   =  ----  ------   ------   .....
-		 *			2n   2(n+1)   2(n+2)
-		 *			-- - ------ - ------ -
-		 *			 x     x         x
-		 *
-		 * Let w = 2n/x and h=2/x, then the above quotient
-		 * is equal to the continued fraction:
-		 *		    1
-		 *	= -----------------------
-		 *		       1
-		 *	   w - -----------------
-		 *			  1
-		 * 	        w+h - ---------
-		 *		       w+2h - ...
-		 *
-		 * To determine how many terms needed, let
-		 * Q(0) = w, Q(1) = w(w+h) - 1,
-		 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
-		 * When Q(k) > 1e4	good for single
-		 * When Q(k) > 1e9	good for double
-		 * When Q(k) > 1e17	good for quadruple
-		 */
-	    /* determine k */
-		double t,v;
-		double q0,q1,h,tmp; int k,m;
-		w  = (n+n)/(double)x; h = 2.0/(double)x;
-		q0 = w;  z = w+h; q1 = w*z - 1.0; k=1;
-		while (q1<1.0e9) {
-			k += 1; z += h;
-			tmp = z*q1 - q0;
-			q0 = q1;
-			q1 = tmp;
-		}
-		m = n+n;
-		for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
-		a = t;
-		b = one;
-		/*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
-		 *  Hence, if n*(log(2n/x)) > ...
-		 *  single 8.8722839355e+01
-		 *  double 7.09782712893383973096e+02
-		 *  long double 1.1356523406294143949491931077970765006170e+04
-		 *  then recurrent value may overflow and the result will
-		 *  likely underflow to zero
-		 */
-		tmp = n;
-		v = two/x;
-		tmp = tmp*log(fabs(v*tmp));
-	    	for (i=n-1;i>0;i--){
-		        temp = b;
-		        b = ((i+i)/x)*b - a;
-		        a = temp;
-		    /* scale b to avoid spurious overflow */
-#			if defined(__vax__)
-#				define BMAX 1e13
-#			else
-#				define BMAX 1e100
-#			endif /* defined(__vax__) */
-			if (b > BMAX) {
-				a /= b;
-				t /= b;
-				b = one;
-			}
-		}
-	    	b = (t*j0(x)/b);
-	    }
-	}
-	return ((sgn == 1) ? -b : b);
-}
-
-double
-yn(int n, double x)
-{
-	int i, sign;
-	double a, b, temp;
-
-    /* Y(n,NaN), Y(n, x < 0) is NaN */
-	if (x <= 0 || isnan(x))
-		if (_IEEE && x < 0) return zero/zero;
-		else if (x < 0)     return (infnan(EDOM));
-		else if (_IEEE)     return -one/zero;
-		else		    return(infnan(-ERANGE));
-	else if (!isfinite(x)) return(0);
-	sign = 1;
-	if (n<0){
-		n = -n;
-		sign = 1 - ((n&1)<<2);
-	}
-	if (n == 0) return(y0(x));
-	if (n == 1) return(sign*y1(x));
-	if(_IEEE && x >= 8.148143905337944345e+090) { /* x > 2**302 */
-    /* (x >> n**2)
-     *	    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-     *	    Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-     *	    Let s=sin(x), c=cos(x),
-     *		xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
-     *
-     *		   n	sin(xn)*sqt2	cos(xn)*sqt2
-     *		----------------------------------
-     *		   0	 s-c		 c+s
-     *		   1	-s-c 		-c+s
-     *		   2	-s+c		-c-s
-     *		   3	 s+c		 c-s
-     */
-		switch (n&3) {
-		    case 0: temp =  sin(x)-cos(x); break;
-		    case 1: temp = -sin(x)-cos(x); break;
-		    case 2: temp = -sin(x)+cos(x); break;
-		    case 3: temp =  sin(x)+cos(x); break;
-		}
-		b = invsqrtpi*temp/sqrt(x);
-	} else {
-	    a = y0(x);
-	    b = y1(x);
-	/* quit if b is -inf */
-	    for (i = 1; i < n && !isfinite(b); i++){
-		temp = b;
-		b = ((double)(i+i)/x)*b - a;
-		a = temp;
-	    }
-	}
-	if (!_IEEE && !isfinite(b))
-		return (infnan(-sign * ERANGE));
-	return ((sign > 0) ? b : -b);
-}
diff --git a/lib/libm/noieee_src/n_lgamma.c b/lib/libm/noieee_src/n_lgamma.c
deleted file mode 100644
index 44a1922e000..00000000000
--- a/lib/libm/noieee_src/n_lgamma.c
+++ /dev/null
@@ -1,327 +0,0 @@
-/*	$OpenBSD: n_lgamma.c,v 1.10 2016/09/12 04:39:47 guenther Exp $	*/
-/*-
- * Copyright (c) 1992, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/*
- * Coded by Peter McIlroy, Nov 1992;
- *
- * The financial support of UUNET Communications Services is greatfully
- * acknowledged.
- */
-
-#include <math.h>
-#include <errno.h>
-
-#include "mathimpl.h"
-
-/* Log gamma function.
- * Error:  x > 0 error < 1.3ulp.
- *	   x > 4, error < 1ulp.
- *	   x > 9, error < .6ulp.
- * 	   x < 0, all bets are off. (When G(x) ~ 1, log(G(x)) ~ 0)
- * Method:
- *	x > 6:
- *		Use the asymptotic expansion (Stirling's Formula)
- *	0 < x < 6:
- *		Use tgamma(x+1) = x*tgamma(x) for argument reduction.
- *		Use rational approximation in
- *		the range 1.2, 2.5
- *		Two approximations are used, one centered at the
- *		minimum to ensure monotonicity; one centered at 2
- *		to maintain small relative error.
- *	x < 0:
- *		Use the reflection formula,
- *		G(1-x)G(x) = PI/sin(PI*x)
- * Special values:
- *	non-positive integer	returns +Inf.
- *	NaN			returns NaN
-*/
-#if defined(__vax__)
-#define _IEEE		0
-/* double and float have same size exponent field */
-#define TRUNC(x)	x = (double) (float) (x)
-#else
-static int endian;
-#define _IEEE		1
-#define TRUNC(x)	*(((int *) &x) + endian) &= 0xf8000000
-#define infnan(x)	0.0
-#endif
-
-static double small_lgam(double);
-static double large_lgam(double);
-static double neg_lgam(double);
-static const double one = 1.0;
-extern int signgam;
-
-#define UNDERFL (1e-1020 * 1e-1020)
-
-#define LEFT	(1.0 - (x0 + .25))
-#define RIGHT	(x0 - .218)
-/*
- * Constants for approximation in [1.244,1.712]
- */
-#define x0	0.461632144968362356785
-#define x0_lo	-.000000000000000015522348162858676890521
-#define a0_hi	-0.12148629128932952880859
-#define a0_lo	.0000000007534799204229502
-#define r0	-2.771227512955130520e-002
-#define r1	-2.980729795228150847e-001
-#define r2	-3.257411333183093394e-001
-#define r3	-1.126814387531706041e-001
-#define r4	-1.129130057170225562e-002
-#define r5	-2.259650588213369095e-005
-#define s0	 1.714457160001714442e+000
-#define s1	 2.786469504618194648e+000
-#define s2	 1.564546365519179805e+000
-#define s3	 3.485846389981109850e-001
-#define s4	 2.467759345363656348e-002
-/*
- * Constants for approximation in [1.71, 2.5]
-*/
-#define a1_hi	4.227843350984671344505727574870e-01
-#define a1_lo	4.670126436531227189e-18
-#define p0	3.224670334241133695662995251041e-01
-#define p1	3.569659696950364669021382724168e-01
-#define p2	1.342918716072560025853732668111e-01
-#define p3	1.950702176409779831089963408886e-02
-#define p4	8.546740251667538090796227834289e-04
-#define q0	1.000000000000000444089209850062e+00
-#define q1	1.315850076960161985084596381057e+00
-#define q2	6.274644311862156431658377186977e-01
-#define q3	1.304706631926259297049597307705e-01
-#define q4	1.102815279606722369265536798366e-02
-#define q5	2.512690594856678929537585620579e-04
-#define q6	-1.003597548112371003358107325598e-06
-/*
- * Stirling's Formula, adjusted for equal-ripple. x in [6,Inf].
-*/
-#define lns2pi	.418938533204672741780329736405
-#define pb0	 8.33333333333333148296162562474e-02
-#define pb1	-2.77777777774548123579378966497e-03
-#define pb2	 7.93650778754435631476282786423e-04
-#define pb3	-5.95235082566672847950717262222e-04
-#define pb4	 8.41428560346653702135821806252e-04
-#define pb5	-1.89773526463879200348872089421e-03
-#define pb6	 5.69394463439411649408050664078e-03
-#define pb7	-1.44705562421428915453880392761e-02
-
-__pure double
-lgamma(double x)
-{
-	double r;
-
-	int signgam = 1;
-#if _IEEE
-	endian = ((*(int *) &one)) ? 1 : 0;
-#endif
-
-	if (!isfinite(x))
-		if (_IEEE)
-			return (x+x);
-		else return (infnan(EDOM));
-
-	if (x > 6 + RIGHT) {
-		r = large_lgam(x);
-		return (r);
-	} else if (x > 1e-16)
-		return (small_lgam(x));
-	else if (x > -1e-16) {
-		if (x < 0) {
-			signgam = -1;
-			x = -x;
-		}
-		return (-log(x));
-	} else
-		return (neg_lgam(x));
-}
-
-float
-lgammaf(float x)
-{
-	return lgamma(x);
-}
-
-/*
- * The gamma() function performs identically to lgamma(), including
- * the use of signgam.
- */
-
-double
-gamma(double x)
-{
-	return lgamma(x);
-}
-
-float
-gammaf(float x)
-{
-	return lgammaf(x);
-}
-
-static double
-large_lgam(double x)
-{
-	double z, p, x1;
-	struct Double t, u, v;
-	u = __log__D(x);
-	u.a -= 1.0;
-	if (x > 1e15) {
-		v.a = x - 0.5;
-		TRUNC(v.a);
-		v.b = (x - v.a) - 0.5;
-		t.a = u.a*v.a;
-		t.b = x*u.b + v.b*u.a;
-		if (_IEEE == 0 && !isfinite(t.a))
-			return(infnan(ERANGE));
-		return(t.a + t.b);
-	}
-	x1 = 1./x;
-	z = x1*x1;
-	p = pb0+z*(pb1+z*(pb2+z*(pb3+z*(pb4+z*(pb5+z*(pb6+z*pb7))))));
-					/* error in approximation = 2.8e-19 */
-
-	p = p*x1;			/* error < 2.3e-18 absolute */
-					/* 0 < p < 1/64 (at x = 5.5) */
-	v.a = x = x - 0.5;
-	TRUNC(v.a);			/* truncate v.a to 26 bits. */
-	v.b = x - v.a;
-	t.a = v.a*u.a;			/* t = (x-.5)*(log(x)-1) */
-	t.b = v.b*u.a + x*u.b;
-	t.b += p; t.b += lns2pi;	/* return t + lns2pi + p */
-	return (t.a + t.b);
-}
-
-static double
-small_lgam(double x)
-{
-	int x_int;
-	double y, z, t, r = 0, p, q, hi, lo;
-	struct Double rr;
-	x_int = (x + .5);
-	y = x - x_int;
-	if (x_int <= 2 && y > RIGHT) {
-		t = y - x0;
-		y--; x_int++;
-		goto CONTINUE;
-	} else if (y < -LEFT) {
-		t = y +(1.0-x0);
-CONTINUE:
-		z = t - x0_lo;
-		p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5))));
-		q = s0+z*(s1+z*(s2+z*(s3+z*s4)));
-		r = t*(z*(p/q) - x0_lo);
-		t = .5*t*t;
-		z = 1.0;
-		switch (x_int) {
-		case 6:	z  = (y + 5);		/* FALLTHROUGH */
-		case 5:	z *= (y + 4);		/* FALLTHROUGH */
-		case 4:	z *= (y + 3);		/* FALLTHROUGH */
-		case 3:	z *= (y + 2);
-			rr = __log__D(z);
-			rr.b += a0_lo; rr.a += a0_hi;
-			return(((r+rr.b)+t+rr.a));
-		case 2: return(((r+a0_lo)+t)+a0_hi);
-		case 0: r -= log1p(x);	/* FALLTHROUGH */
-		default: rr = __log__D(x);
-			rr.a -= a0_hi; rr.b -= a0_lo;
-			return(((r - rr.b) + t) - rr.a);
-		}
-	} else {
-		p = p0+y*(p1+y*(p2+y*(p3+y*p4)));
-		q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6)))));
-		p = p*(y/q);
-		t = (double)(float) y;
-		z = y-t;
-		hi = (double)(float) (p+a1_hi);
-		lo = a1_hi - hi; lo += p; lo += a1_lo;
-		r = lo*y + z*hi;	/* q + r = y*(a0+p/q) */
-		q = hi*t;
-		z = 1.0;
-		switch (x_int) {
-		case 6:	z  = (y + 5);		/* FALLTHROUGH */
-		case 5:	z *= (y + 4);		/* FALLTHROUGH */
-		case 4:	z *= (y + 3);		/* FALLTHROUGH */
-		case 3:	z *= (y + 2);		/* FALLTHROUGH */
-			rr = __log__D(z);
-			r += rr.b; r += q;
-			return(rr.a + r);
-		case 2:	return (q+ r);
-		case 0: rr = __log__D(x);
-			r -= rr.b; r -= log1p(x);
-			r += q; r-= rr.a;
-			return(r);
-		default: rr = __log__D(x);
-			r -= rr.b;
-			q -= rr.a;
-			return (r+q);
-		}
-	}
-}
-
-static double
-neg_lgam(double x)
-{
-	int xi;
-	double y, z, zero = 0.0;
-
-	/* avoid destructive cancellation as much as possible */
-	if (x > -170) {
-		xi = x;
-		if (xi == x)
-			if (_IEEE)
-				return(one/zero);
-			else
-				return(infnan(ERANGE));
-		y = tgamma(x);
-		if (y < 0) {
-			y = -y;
-			signgam = -1;
-		}
-		return (log(y));
-	}
-	z = floor(x + .5);
-	if (z == x) {		/* convention: G(-(integer)) -> +Inf */
-		if (_IEEE)
-			return (one/zero);
-		else
-			return (infnan(ERANGE));
-	}
-	y = .5*ceil(x);
-	if (y == ceil(y))
-		signgam = -1;
-	x = -x;
-	z = fabs(x + z);	/* 0 < z <= .5 */
-	if (z < .25)
-		z = sin(M_PI*z);
-	else
-		z = cos(M_PI*(0.5-z));
-	z = log(M_PI/(z*x));
-	y = large_lgam(x);
-	return (z - y);
-}
diff --git a/lib/libm/noieee_src/n_log.c b/lib/libm/noieee_src/n_log.c
deleted file mode 100644
index ce5a9f40801..00000000000
--- a/lib/libm/noieee_src/n_log.c
+++ /dev/null
@@ -1,482 +0,0 @@
-/*	$OpenBSD: n_log.c,v 1.10 2016/09/12 19:47:02 guenther Exp $	*/
-/*
- * Copyright (c) 1992, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include <math.h>
-#include <errno.h>
-
-#include "mathimpl.h"
-
-/* Table-driven natural logarithm.
- *
- * This code was derived, with minor modifications, from:
- *	Peter Tang, "Table-Driven Implementation of the
- *	Logarithm in IEEE Floating-Point arithmetic." ACM Trans.
- *	Math Software, vol 16. no 4, pp 378-400, Dec 1990).
- *
- * Calculates log(2^m*F*(1+f/F)), |f/j| <= 1/256,
- * where F = j/128 for j an integer in [0, 128].
- *
- * log(2^m) = log2_hi*m + log2_tail*m
- * since m is an integer, the dominant term is exact.
- * m has at most 10 digits (for subnormal numbers),
- * and log2_hi has 11 trailing zero bits.
- *
- * log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table.h
- * logF_hi[] + 512 is exact.
- *
- * log(1+f/F) = 2*f/(2*F + f) + 1/12 * (2*f/(2*F + f))**3 + ...
- * the leading term is calculated to extra precision in two
- * parts, the larger of which adds exactly to the dominant
- * m and F terms.
- * There are two cases:
- *	1. when m, j are non-zero (m | j), use absolute
- *	   precision for the leading term.
- *	2. when m = j = 0, |1-x| < 1/256, and log(x) ~= (x-1).
- *	   In this case, use a relative precision of 24 bits.
- * (This is done differently in the original paper)
- *
- * Special cases:
- *	0	return signalling -Inf
- *	neg	return signalling NaN
- *	+Inf	return +Inf
-*/
-
-#if defined(__vax__)
-#define _IEEE		0
-#define TRUNC(x)	x = (double) (float) (x)
-#else
-#define _IEEE		1
-#define endian		(((*(int *) &one)) ? 1 : 0)
-#define TRUNC(x)	*(((int *) &x) + endian) &= 0xf8000000
-#define infnan(x)	0.0
-#endif
-
-#define N 128
-
-/* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128.
- * Used for generation of extend precision logarithms.
- * The constant 35184372088832 is 2^45, so the divide is exact.
- * It ensures correct reading of logF_head, even for inaccurate
- * decimal-to-binary conversion routines.  (Everybody gets the
- * right answer for integers less than 2^53.)
- * Values for log(F) were generated using error < 10^-57 absolute
- * with the bc -l package.
-*/
-static const double	A1 = 	  .08333333333333178827;
-static const double	A2 = 	  .01250000000377174923;
-static const double	A3 =	 .002232139987919447809;
-static const double	A4 =	.0004348877777076145742;
-
-static const double logF_head[N+1] = {
-	0.,
-	.007782140442060381246,
-	.015504186535963526694,
-	.023167059281547608406,
-	.030771658666765233647,
-	.038318864302141264488,
-	.045809536031242714670,
-	.053244514518837604555,
-	.060624621816486978786,
-	.067950661908525944454,
-	.075223421237524235039,
-	.082443669210988446138,
-	.089612158689760690322,
-	.096729626458454731618,
-	.103796793681567578460,
-	.110814366340264314203,
-	.117783035656430001836,
-	.124703478501032805070,
-	.131576357788617315236,
-	.138402322859292326029,
-	.145182009844575077295,
-	.151916042025732167530,
-	.158605030176659056451,
-	.165249572895390883786,
-	.171850256926518341060,
-	.178407657472689606947,
-	.184922338493834104156,
-	.191394852999565046047,
-	.197825743329758552135,
-	.204215541428766300668,
-	.210564769107350002741,
-	.216873938300523150246,
-	.223143551314024080056,
-	.229374101064877322642,
-	.235566071312860003672,
-	.241719936886966024758,
-	.247836163904594286577,
-	.253915209980732470285,
-	.259957524436686071567,
-	.265963548496984003577,
-	.271933715484010463114,
-	.277868451003087102435,
-	.283768173130738432519,
-	.289633292582948342896,
-	.295464212893421063199,
-	.301261330578199704177,
-	.307025035294827830512,
-	.312755710004239517729,
-	.318453731118097493890,
-	.324119468654316733591,
-	.329753286372579168528,
-	.335355541920762334484,
-	.340926586970454081892,
-	.346466767346100823488,
-	.351976423156884266063,
-	.357455888922231679316,
-	.362905493689140712376,
-	.368325561158599157352,
-	.373716409793814818840,
-	.379078352934811846353,
-	.384411698910298582632,
-	.389716751140440464951,
-	.394993808240542421117,
-	.400243164127459749579,
-	.405465108107819105498,
-	.410659924985338875558,
-	.415827895143593195825,
-	.420969294644237379543,
-	.426084395310681429691,
-	.431173464818130014464,
-	.436236766774527495726,
-	.441274560805140936281,
-	.446287102628048160113,
-	.451274644139630254358,
-	.456237433481874177232,
-	.461175715122408291790,
-	.466089729924533457960,
-	.470979715219073113985,
-	.475845904869856894947,
-	.480688529345570714212,
-	.485507815781602403149,
-	.490303988045525329653,
-	.495077266798034543171,
-	.499827869556611403822,
-	.504556010751912253908,
-	.509261901790523552335,
-	.513945751101346104405,
-	.518607764208354637958,
-	.523248143765158602036,
-	.527867089620485785417,
-	.532464798869114019908,
-	.537041465897345915436,
-	.541597282432121573947,
-	.546132437597407260909,
-	.550647117952394182793,
-	.555141507540611200965,
-	.559615787935399566777,
-	.564070138285387656651,
-	.568504735352689749561,
-	.572919753562018740922,
-	.577315365035246941260,
-	.581691739635061821900,
-	.586049045003164792433,
-	.590387446602107957005,
-	.594707107746216934174,
-	.599008189645246602594,
-	.603290851438941899687,
-	.607555250224322662688,
-	.611801541106615331955,
-	.616029877215623855590,
-	.620240409751204424537,
-	.624433288012369303032,
-	.628608659422752680256,
-	.632766669570628437213,
-	.636907462236194987781,
-	.641031179420679109171,
-	.645137961373620782978,
-	.649227946625615004450,
-	.653301272011958644725,
-	.657358072709030238911,
-	.661398482245203922502,
-	.665422632544505177065,
-	.669430653942981734871,
-	.673422675212350441142,
-	.677398823590920073911,
-	.681359224807238206267,
-	.685304003098281100392,
-	.689233281238557538017,
-	.693147180560117703862
-};
-
-static const double logF_tail[N+1] = {
-	0.,
-	-.00000000000000543229938420049,
-	 .00000000000000172745674997061,
-	-.00000000000001323017818229233,
-	-.00000000000001154527628289872,
-	-.00000000000000466529469958300,
-	 .00000000000005148849572685810,
-	-.00000000000002532168943117445,
-	-.00000000000005213620639136504,
-	-.00000000000001819506003016881,
-	 .00000000000006329065958724544,
-	 .00000000000008614512936087814,
-	-.00000000000007355770219435028,
-	 .00000000000009638067658552277,
-	 .00000000000007598636597194141,
-	 .00000000000002579999128306990,
-	-.00000000000004654729747598444,
-	-.00000000000007556920687451336,
-	 .00000000000010195735223708472,
-	-.00000000000017319034406422306,
-	-.00000000000007718001336828098,
-	 .00000000000010980754099855238,
-	-.00000000000002047235780046195,
-	-.00000000000008372091099235912,
-	 .00000000000014088127937111135,
-	 .00000000000012869017157588257,
-	 .00000000000017788850778198106,
-	 .00000000000006440856150696891,
-	 .00000000000016132822667240822,
-	-.00000000000007540916511956188,
-	-.00000000000000036507188831790,
-	 .00000000000009120937249914984,
-	 .00000000000018567570959796010,
-	-.00000000000003149265065191483,
-	-.00000000000009309459495196889,
-	 .00000000000017914338601329117,
-	-.00000000000001302979717330866,
-	 .00000000000023097385217586939,
-	 .00000000000023999540484211737,
-	 .00000000000015393776174455408,
-	-.00000000000036870428315837678,
-	 .00000000000036920375082080089,
-	-.00000000000009383417223663699,
-	 .00000000000009433398189512690,
-	 .00000000000041481318704258568,
-	-.00000000000003792316480209314,
-	 .00000000000008403156304792424,
-	-.00000000000034262934348285429,
-	 .00000000000043712191957429145,
-	-.00000000000010475750058776541,
-	-.00000000000011118671389559323,
-	 .00000000000037549577257259853,
-	 .00000000000013912841212197565,
-	 .00000000000010775743037572640,
-	 .00000000000029391859187648000,
-	-.00000000000042790509060060774,
-	 .00000000000022774076114039555,
-	 .00000000000010849569622967912,
-	-.00000000000023073801945705758,
-	 .00000000000015761203773969435,
-	 .00000000000003345710269544082,
-	-.00000000000041525158063436123,
-	 .00000000000032655698896907146,
-	-.00000000000044704265010452446,
-	 .00000000000034527647952039772,
-	-.00000000000007048962392109746,
-	 .00000000000011776978751369214,
-	-.00000000000010774341461609578,
-	 .00000000000021863343293215910,
-	 .00000000000024132639491333131,
-	 .00000000000039057462209830700,
-	-.00000000000026570679203560751,
-	 .00000000000037135141919592021,
-	-.00000000000017166921336082431,
-	-.00000000000028658285157914353,
-	-.00000000000023812542263446809,
-	 .00000000000006576659768580062,
-	-.00000000000028210143846181267,
-	 .00000000000010701931762114254,
-	 .00000000000018119346366441110,
-	 .00000000000009840465278232627,
-	-.00000000000033149150282752542,
-	-.00000000000018302857356041668,
-	-.00000000000016207400156744949,
-	 .00000000000048303314949553201,
-	-.00000000000071560553172382115,
-	 .00000000000088821239518571855,
-	-.00000000000030900580513238244,
-	-.00000000000061076551972851496,
-	 .00000000000035659969663347830,
-	 .00000000000035782396591276383,
-	-.00000000000046226087001544578,
-	 .00000000000062279762917225156,
-	 .00000000000072838947272065741,
-	 .00000000000026809646615211673,
-	-.00000000000010960825046059278,
-	 .00000000000002311949383800537,
-	-.00000000000058469058005299247,
-	-.00000000000002103748251144494,
-	-.00000000000023323182945587408,
-	-.00000000000042333694288141916,
-	-.00000000000043933937969737844,
-	 .00000000000041341647073835565,
-	 .00000000000006841763641591466,
-	 .00000000000047585534004430641,
-	 .00000000000083679678674757695,
-	-.00000000000085763734646658640,
-	 .00000000000021913281229340092,
-	-.00000000000062242842536431148,
-	-.00000000000010983594325438430,
-	 .00000000000065310431377633651,
-	-.00000000000047580199021710769,
-	-.00000000000037854251265457040,
-	 .00000000000040939233218678664,
-	 .00000000000087424383914858291,
-	 .00000000000025218188456842882,
-	-.00000000000003608131360422557,
-	-.00000000000050518555924280902,
-	 .00000000000078699403323355317,
-	-.00000000000067020876961949060,
-	 .00000000000016108575753932458,
-	 .00000000000058527188436251509,
-	-.00000000000035246757297904791,
-	-.00000000000018372084495629058,
-	 .00000000000088606689813494916,
-	 .00000000000066486268071468700,
-	 .00000000000063831615170646519,
-	 .00000000000025144230728376072,
-	-.00000000000017239444525614834
-};
-
-double
-log(double x)
-{
-	int m, j;
-	double F, f, g, q, u, u2, v, zero = 0.0, one = 1.0;
-	volatile double u1;
-
-	/* Catch special cases */
-	if (x <= 0)
-		if (_IEEE && x == zero)	/* log(0) = -Inf */
-			return (-one/zero);
-		else if (_IEEE)		/* log(neg) = NaN */
-			return (zero/zero);
-		else if (x == zero)	/* NOT REACHED IF _IEEE */
-			return (infnan(-ERANGE));
-		else
-			return (infnan(EDOM));
-	else if (!isfinite(x))
-		if (_IEEE)		/* x = NaN, Inf */
-			return (x+x);
-		else
-			return (infnan(ERANGE));
-
-	/* Argument reduction: 1 <= g < 2; x/2^m = g;	*/
-	/* y = F*(1 + f/F) for |f| <= 2^-8		*/
-
-	m = logb(x);
-	g = ldexp(x, -m);
-	if (_IEEE && m == -1022) {
-		j = logb(g);
-		m += j;
-		g = ldexp(g, -j);
-	}
-	j = N*(g-1) + .5;
-	F = (1.0/N) * j + 1;	/* F*128 is an integer in [128, 512] */
-	f = g - F;
-
-	/* Approximate expansion for log(1+f/F) ~= u + q */
-	g = 1/(2*F+f);
-	u = 2*f*g;
-	v = u*u;
-	q = u*v*(A1 + v*(A2 + v*(A3 + v*A4)));
-
-    /* case 1: u1 = u rounded to 2^-43 absolute.  Since u < 2^-8,
-     * 	       u1 has at most 35 bits, and F*u1 is exact, as F has < 8 bits.
-     *         It also adds exactly to |m*log2_hi + log_F_head[j] | < 750
-    */
-	if (m | j) {
-		u1 = u + 513;
-		u1 -= 513;
-	}
-
-    /* case 2:	|1-x| < 1/256. The m- and j- dependent terms are zero;
-     * 		u1 = u to 24 bits.
-    */
-	else {
-		u1 = u;
-		TRUNC(u1);
-	}
-	u2 = (2.0*(f - F*u1) - u1*f) * g;
-			/* u1 + u2 = 2f/(2F+f) to extra precision.	*/
-
-	/* log(x) = log(2^m*F*(1+f/F)) =				*/
-	/* (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q);	*/
-	/* (exact) + (tiny)						*/
-
-	u1 += m*logF_head[N] + logF_head[j];		/* exact */
-	u2 = (u2 + logF_tail[j]) + q;			/* tiny */
-	u2 += logF_tail[N]*m;
-	return (u1 + u2);
-}
-DEF_STD(log);
-
-/*
- * Extra precision variant, returning struct {double a, b;};
- * log(x) = a+b to 63 bits, with a rounded to 26 bits.
- */
-struct Double
-__log__D(double x)
-{
-	int m, j;
-	double F, f, g, q, u, v, u2, one = 1.0;
-	volatile double u1;
-	struct Double r;
-
-	/* Argument reduction: 1 <= g < 2; x/2^m = g;	*/
-	/* y = F*(1 + f/F) for |f| <= 2^-8		*/
-
-	m = logb(x);
-	g = ldexp(x, -m);
-	if (_IEEE && m == -1022) {
-		j = logb(g);
-		m += j;
-		g = ldexp(g, -j);
-	}
-	j = N*(g-1) + .5;
-	F = (1.0/N) * j + 1;
-	f = g - F;
-
-	g = 1/(2*F+f);
-	u = 2*f*g;
-	v = u*u;
-	q = u*v*(A1 + v*(A2 + v*(A3 + v*A4)));
-	if (m | j) {
-		u1 = u + 513;
-		u1 -= 513;
-	}
-	else {
-		u1 = u;
-		TRUNC(u1);
-	}
-	u2 = (2.0*(f - F*u1) - u1*f) * g;
-
-	u1 += m*logF_head[N] + logF_head[j];
-
-	u2 +=  logF_tail[j]; u2 += q;
-	u2 += logF_tail[N]*m;
-	r.a = u1 + u2;			/* Only difference is here */
-	TRUNC(r.a);
-	r.b = (u1 - r.a) + u2;
-	return (r);
-}
diff --git a/lib/libm/noieee_src/n_log10.c b/lib/libm/noieee_src/n_log10.c
deleted file mode 100644
index 5c0126a0691..00000000000
--- a/lib/libm/noieee_src/n_log10.c
+++ /dev/null
@@ -1,83 +0,0 @@
-/*	$OpenBSD: n_log10.c,v 1.8 2009/10/27 23:59:29 deraadt Exp $	*/
-/*	$NetBSD: n_log10.c,v 1.1 1995/10/10 23:36:58 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* LOG10(X)
- * RETURN THE BASE 10 LOGARITHM OF x
- * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/20/85;
- * REVISED BY K.C. NG on 1/23/85, 3/7/85, 4/16/85.
- *
- * Required kernel function:
- *	log(x)
- *
- * Method :
- *			     log(x)
- *		log10(x) = ---------  or  [1/log(10)]*log(x)
- *			    log(10)
- *
- *    Note:
- *	  [log(10)]   rounded to 56 bits has error  .0895  ulps,
- *	  [1/log(10)] rounded to 53 bits has error  .198   ulps;
- *	  therefore, for better accuracy, in VAX D format, we divide
- *	  log(x) by log(10), but in IEEE Double format, we multiply
- *	  log(x) by [1/log(10)].
- *
- * Special cases:
- *	log10(x) is NaN with signal if x < 0;
- *	log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
- *	log10(NaN) is that NaN with no signal.
- *
- * Accuracy:
- *	log10(X) returns the exact log10(x) nearly rounded. In a test run
- *	with 1,536,000 random arguments on a VAX, the maximum observed
- *	error was 1.74 ulps (units in the last place).
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-static const double ln10hi = 2.3025850929940456790E0;
-
-double
-log10(double x)
-{
-#if defined(__vax__)
-	return(log(x)/ln10hi);
-#else	/* defined(__vax__) */
-	return(ivln10*log(x));
-#endif	/* defined(__vax__) */
-}
diff --git a/lib/libm/noieee_src/n_log1p.c b/lib/libm/noieee_src/n_log1p.c
deleted file mode 100644
index 154cb1c4901..00000000000
--- a/lib/libm/noieee_src/n_log1p.c
+++ /dev/null
@@ -1,154 +0,0 @@
-/*	$OpenBSD: n_log1p.c,v 1.12 2016/09/12 04:39:47 guenther Exp $	*/
-/*	$NetBSD: n_log1p.c,v 1.1 1995/10/10 23:37:00 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* LOG1P(x)
- * RETURN THE LOGARITHM OF 1+x
- * DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/19/85;
- * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/24/85, 4/16/85.
- *
- * Required system supported functions:
- *	scalbn(x,n)
- *	copysign(x,y)
- *	logb(x)
- *	isfinite(x)
- *
- * Required kernel function:
- *	log__L(z)
- *
- * Method :
- *	1. Argument Reduction: find k and f such that
- *			1+x  = 2^k * (1+f),
- *	   where  sqrt(2)/2 < 1+f < sqrt(2) .
- *
- *	2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
- *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
- *	   log(1+f) is computed by
- *
- *	     		log(1+f) = 2s + s*log__L(s*s)
- *	   where
- *		log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...)))
- *
- *	   See log__L() for the values of the coefficients.
- *
- *	3. Finally,  log(1+x) = k*ln2 + log(1+f).
- *
- *	Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers
- *		   n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last
- *		   20 bits (for VAX D format), or the last 21 bits ( for IEEE
- *		   double) is 0. This ensures n*ln2hi is exactly representable.
- *		2. In step 1, f may not be representable. A correction term c
- *	 	   for f is computed. It follows that the correction term for
- *		   f - t (the leading term of log(1+f) in step 2) is c-c*x. We
- *		   add this correction term to n*ln2lo to attenuate the error.
- *
- *
- * Special cases:
- *	log1p(x) is NaN with signal if x < -1; log1p(NaN) is NaN with no signal;
- *	log1p(INF) is +INF; log1p(-1) is -INF with signal;
- *	only log1p(0)=0 is exact for finite argument.
- *
- * Accuracy:
- *	log1p(x) returns the exact log(1+x) nearly rounded. In a test run
- *	with 1,536,000 random arguments on a VAX, the maximum observed
- *	error was .846 ulps (units in the last place).
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include <errno.h>
-#include "math.h"
-#include "mathimpl.h"
-
-static const double ln2hi = 6.9314718055829871446E-1;
-static const double ln2lo = 1.6465949582897081279E-12;
-static const double sqrt2 = 1.4142135623730950622E0;
-
-double
-log1p(double x)
-{
-	static const double zero=0.0, negone= -1.0, one=1.0,
-		      half=1.0/2.0, small=1.0E-20;   /* 1+small == 1 */
-	double z,s,t,c;
-	int k;
-
-	if (isnan(x))
-		return (x);
-
-	if(isfinite(x)) {
-	   if( x > negone ) {
-
-	   /* argument reduction */
-	      if(copysign(x,one)<small) return(x);
-	      k=logb(one+x); z=scalbn(x,-k); t=scalbn(one,-k);
-	      if(z+t >= sqrt2 )
-		  { k += 1 ; z *= half; t *= half; }
-	      t += negone; x = z + t;
-	      c = (t-x)+z ;		/* correction term for x */
-
- 	   /* compute log(1+x)  */
-              s = x/(2+x); t = x*x*half;
-	      c += (k*ln2lo-c*x);
-	      z = c+s*(t+__log__L(s*s));
-	      x += (z - t) ;
-
-	      return(k*ln2hi+x);
-	   }
-	/* end of if (x > negone) */
-
-	    else {
-#if defined(__vax__)
-		if ( x == negone )
-		    return (infnan(-ERANGE));	/* -INF */
-		else
-		    return (infnan(EDOM));	/* NaN */
-#else	/* defined(__vax__) */
-		/* x = -1, return -INF with signal */
-		if ( x == negone ) return( negone/zero );
-
-		/* negative argument for log, return NaN with signal */
-	        else return ( zero / zero );
-#endif	/* defined(__vax__) */
-	    }
-	}
-    /* end of if (isfinite(x)) */
-
-    /* log(-INF) is NaN */
-	else if(x<0)
-	     return(zero/zero);
-
-    /* log(+INF) is INF */
-	else return(x);
-}
diff --git a/lib/libm/noieee_src/n_log__L.c b/lib/libm/noieee_src/n_log__L.c
deleted file mode 100644
index 84fa2be7c8c..00000000000
--- a/lib/libm/noieee_src/n_log__L.c
+++ /dev/null
@@ -1,86 +0,0 @@
-/*	$OpenBSD: n_log__L.c,v 1.9 2009/10/27 23:59:29 deraadt Exp $	*/
-/*	$NetBSD: n_log__L.c,v 1.1 1995/10/10 23:37:01 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* log__L(Z)
- *		LOG(1+X) - 2S			       X
- * RETURN      ---------------  WHERE Z = S*S,  S = ------- , 0 <= Z <= .0294...
- *		      S				     2 + X
- *
- * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)
- * KERNEL FUNCTION FOR LOG; TO BE USED IN LOG1P, LOG, AND POW FUNCTIONS
- * CODED IN C BY K.C. NG, 1/19/85;
- * REVISED BY K.C. Ng, 2/3/85, 4/16/85.
- *
- * Method :
- *	1. Polynomial approximation: let s = x/(2+x).
- *	   Based on log(1+x) = log(1+s) - log(1-s)
- *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
- *
- *	   (log(1+x) - 2s)/s is computed by
- *
- *	       z*(L1 + z*(L2 + z*(... (L7 + z*L8)...)))
- *
- *	   where z=s*s. (See the listing below for Lk's values.) The
- *	   coefficients are obtained by a special Remes algorithm.
- *
- * Accuracy:
- *	Assuming no rounding error, the maximum magnitude of the approximation
- *	error (absolute) is 2**(-58.49) for IEEE double, and 2**(-63.63)
- *	for VAX D format.
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-static const double L1 = 6.6666666666666703212E-1;
-static const double L2 = 3.9999999999970461961E-1;
-static const double L3 = 2.8571428579395698188E-1;
-static const double L4 = 2.2222221233634724402E-1;
-static const double L5 = 1.8181879517064680057E-1;
-static const double L6 = 1.5382888777946145467E-1;
-static const double L7 = 1.3338356561139403517E-1;
-static const double L8 = 1.2500000000000000000E-1;
-
-double
-__log__L(double z)
-{
-#if defined(__vax__)
-    return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*(L7+z*L8))))))));
-#else	/* defined(__vax__) */
-    return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7)))))));
-#endif	/* defined(__vax__) */
-}
diff --git a/lib/libm/noieee_src/n_lrint.c b/lib/libm/noieee_src/n_lrint.c
deleted file mode 100644
index a0feab5313f..00000000000
--- a/lib/libm/noieee_src/n_lrint.c
+++ /dev/null
@@ -1,57 +0,0 @@
-/*	$OpenBSD: n_lrint.c,v 1.2 2016/09/12 19:47:02 guenther Exp $ */
-/*
- * Copyright (c) 2013 Marc Espie.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE OPENBSD PROJECT AND CONTRIBUTORS
- * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
- * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
- * A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE OPENBSD
- * PROJECT OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
- * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
- * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
- * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
-#include <math.h>
-
-long
-lrint(double x)
-{
-	return (long)rint(x);
-}
-DEF_STD(lrint);
-LDBL_UNUSED_CLONE(lrint);
-
-long
-lrintf(float x)
-{
-	return (long)rintf(x);
-}
-DEF_STD(lrintf);
-
-long long
-llrint(double x)
-{
-	return (long long)rint(x);
-}
-DEF_STD(llrint);
-LDBL_UNUSED_CLONE(llrint);
-
-long long
-llrintf(float x)
-{
-	return (long long)rintf(x);
-}
-DEF_STD(llrintf);
diff --git a/lib/libm/noieee_src/n_nan.c b/lib/libm/noieee_src/n_nan.c
deleted file mode 100644
index 3b0cb57a4ba..00000000000
--- a/lib/libm/noieee_src/n_nan.c
+++ /dev/null
@@ -1,29 +0,0 @@
-/*	$OpenBSD: n_nan.c,v 1.1 2009/03/28 16:16:30 martynas Exp $	*/
-
-/*
- * Written by Martynas Venckus.  Public domain
- */
-
-#include <math.h>
-
-/* ARGSUSED */
-double
-nan(const char *tagp)
-{
-	return (0);
-}
-
-/* ARGSUSED */
-float
-nanf(const char *tagp)
-{
-	return (0);
-}
-
-/* ARGSUSED */
-long double
-nanl(const char *tagp)
-{
-	return (0);
-}
-
diff --git a/lib/libm/noieee_src/n_pow.c b/lib/libm/noieee_src/n_pow.c
deleted file mode 100644
index eeb06e6c845..00000000000
--- a/lib/libm/noieee_src/n_pow.c
+++ /dev/null
@@ -1,205 +0,0 @@
-/*	$OpenBSD: n_pow.c,v 1.13 2016/09/12 04:39:47 guenther Exp $	*/
-/*	$NetBSD: n_pow.c,v 1.1 1995/10/10 23:37:02 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* POW(X,Y)
- * RETURN X**Y
- * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85;
- * REVISED BY K.C. NG on 7/10/85.
- * KERNEL pow_P() REPLACED BY P. McILROY 7/22/92.
- * Required system supported functions:
- *      scalbn(x,n)
- *      logb(x)
- *	copysign(x,y)
- *	isfinite(x)
- *	remainder(x,y)
- *
- * Required kernel functions:
- *	exp__D(a,c)			exp(a + c) for |a| << |c|
- *	struct d_double dlog(x)		r.a + r.b, |r.b| < |r.a|
- *
- * Method
- *	1. Compute and return log(x) in three pieces:
- *		log(x) = n*ln2 + hi + lo,
- *	   where n is an integer.
- *	2. Perform y*log(x) by simulating muti-precision arithmetic and
- *	   return the answer in three pieces:
- *		y*log(x) = m*ln2 + hi + lo,
- *	   where m is an integer.
- *	3. Return x**y = exp(y*log(x))
- *		= 2^m * ( exp(hi+lo) ).
- *
- * Special cases:
- *	(anything) ** 0  is 1 ;
- *	(anything) ** 1  is itself;
- *	(anything) ** NaN is NaN;
- *	NaN ** (anything except 0) is NaN;
- *	+(anything > 1) ** +INF is +INF;
- *	-(anything > 1) ** +INF is NaN;
- *	+-(anything > 1) ** -INF is +0;
- *	+-(anything < 1) ** +INF is +0;
- *	+(anything < 1) ** -INF is +INF;
- *	-(anything < 1) ** -INF is NaN;
- *	+-1 ** +-INF is NaN and signal INVALID;
- *	+0 ** +(anything except 0, NaN)  is +0;
- *	-0 ** +(anything except 0, NaN, odd integer)  is +0;
- *	+0 ** -(anything except 0, NaN)  is +INF and signal DIV-BY-ZERO;
- *	-0 ** -(anything except 0, NaN, odd integer)  is +INF with signal;
- *	-0 ** (odd integer) = -( +0 ** (odd integer) );
- *	+INF ** +(anything except 0,NaN) is +INF;
- *	+INF ** -(anything except 0,NaN) is +0;
- *	-INF ** (odd integer) = -( +INF ** (odd integer) );
- *	-INF ** (even integer) = ( +INF ** (even integer) );
- *	-INF ** -(anything except integer,NaN) is NaN with signal;
- *	-(x=anything) ** (k=integer) is (-1)**k * (x ** k);
- *	-(anything except 0) ** (non-integer) is NaN with signal;
- *
- * Accuracy:
- *	pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX,
- *	and a Zilog Z8000,
- *			pow(integer,integer)
- *	always returns the correct integer provided it is representable.
- *	In a test run with 100,000 random arguments with 0 < x, y < 20.0
- *	on a VAX, the maximum observed error was 1.79 ulps (units in the
- *	last place).
- *
- * Constants :
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include <errno.h>
-#include <math.h>
-
-#include "mathimpl.h"
-
-#if defined(__vax__)
-#define TRUNC(x)	x = (double) (float) x
-#define _IEEE		0
-#else
-#define _IEEE		1
-#define endian		(((*(int *) &one)) ? 1 : 0)
-#define TRUNC(x) 	*(((int *) &x)+endian) &= 0xf8000000
-#define infnan(x)	0.0
-#endif		/* defined(__vax__) */
-
-static const double zero=0.0, one=1.0, two=2.0, negone= -1.0;
-
-static double pow_P(double, double);
-
-double
-pow(double x, double y)
-{
-	double t;
-	if (y==zero)
-		return (one);
-	else if (y==one || isnan(x))
-		return (x);		/* if x is NaN or y=1 */
-	else if (isnan(y))		/* if y is NaN */
-		return (y);
-	else if (!isfinite(y))		/* if y is INF */
-		if ((t=fabs(x))==one)	/* +-1 ** +-INF is NaN */
-			return (y - y);
-		else if (t>one)
-			return ((y<0)? zero : ((x<zero)? y-y : y));
-		else
-			return ((y>0)? zero : ((x<0)? y-y : -y));
-	else if (y==two)
-		return (x*x);
-	else if (y==negone)
-		return (one/x);
-    /* x > 0, x == +0 */
-	else if (copysign(one, x) == one)
-		return (pow_P(x, y));
-
-    /* sign(x)= -1 */
-	/* if y is an even integer */
-	else if ( (t=remainder(y,two)) == zero)
-		return (pow_P(-x, y));
-
-	/* if y is an odd integer */
-	else if (copysign(t,one) == one)
-		return (-pow_P(-x, y));
-
-	/* Henceforth y is not an integer */
-	else if (x==zero)	/* x is -0 */
-		return ((y>zero)? -x : one/(-x));
-	else if (_IEEE)
-		return (zero/zero);
-	else
-		return (infnan(EDOM));
-}
-/* kernel function for x >= 0 */
-static double
-pow_P(double x, double y)
-{
-	struct Double s, t;
-	double  huge = 1e300, tiny = 1e-300;
-
-	if (x == zero)
-		if (y > zero)
-			return (zero);
-		else if (_IEEE)
-			return (huge*huge);
-		else
-			return (infnan(ERANGE));
-	if (x == one)
-		return (one);
-	if (!isfinite(x))
-		if (y < zero)
-			return (zero);
-		else if (_IEEE)
-			return (huge*huge);
-		else
-			return (infnan(ERANGE));
-	if (y >= 7e18)		/* infinity */
-		if (x < 1)
-			return(tiny*tiny);
-		else if (_IEEE)
-			return (huge*huge);
-		else
-			return (infnan(ERANGE));
-
-	/* Return exp(y*log(x)), using simulated extended */
-	/* precision for the log and the multiply.	  */
-
-	s = __log__D(x);
-	t.a = y;
-	TRUNC(t.a);
-	t.b = y - t.a;
-	t.b = s.b*y + t.b*s.a;
-	t.a *= s.a;
-	s.a = t.a + t.b;
-	s.b = (t.a - s.a) + t.b;
-	return (__exp__D(s.a, s.b));
-}
diff --git a/lib/libm/noieee_src/n_sincos.c b/lib/libm/noieee_src/n_sincos.c
deleted file mode 100644
index a3da96c6b9e..00000000000
--- a/lib/libm/noieee_src/n_sincos.c
+++ /dev/null
@@ -1,110 +0,0 @@
-/*	$OpenBSD: n_sincos.c,v 1.16 2016/09/12 19:47:02 guenther Exp $	*/
-/*	$NetBSD: n_sincos.c,v 1.1 1995/10/10 23:37:04 ragge Exp $	*/
-/*
- * Copyright (c) 1987, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include <math.h>
-
-#include "mathimpl.h"
-
-float
-sinf(float x)
-{
-	return (float)sin((double) x);
-}
-DEF_STD(sinf);
-
-double
-sin(double x)
-{
-	double a,c,z;
-
-        if(!isfinite(x))	/* sin(NaN) and sin(INF) must be NaN */
-		return x-x;
-	x=remainder(x,PI2);	/* reduce x into [-PI,PI] */
-	a=copysign(x,one);
-	if (a >= PIo4) {
-		if(a >= PI3o4)		/* ... in [3PI/4,PI] */
-			x = copysign((a = PI-a),x);
-		else {			/* ... in [PI/4,3PI/4]  */
-			a = PIo2-a;		/* rtn. sign(x)*C(PI/2-|x|) */
-			z = a*a;
-			c = cos__C(z);
-			z *= half;
-			a = (z >= thresh ? half-((z-half)-c) : one-(z-c));
-			return copysign(a,x);
-		}
-	}
-
-	if (a < small) {		/* rtn. S(x) */
-		big+a;
-		return x;
-	}
-	return x+x*sin__S(x*x);
-}
-DEF_STD(sin);
-LDBL_CLONE(sin);
-
-float
-cosf(float x)
-{
-	return (float)cos((double) x);
-}
-
-double
-cos(double x)
-{
-	double a,c,z,s = 1.0;
-
-	if(!isfinite(x))	/* cos(NaN) and cos(INF) must be NaN */
-		return x-x;
-	x=remainder(x,PI2);	/* reduce x into [-PI,PI] */
-	a=copysign(x,one);
-	if (a >= PIo4) {
-		if (a >= PI3o4) {	/* ... in [3PI/4,PI] */
-			a = PI-a;
-			s = negone;
-		}
-		else {			/* ... in [PI/4,3PI/4] */
-			a = PIo2-a;
-			return a+a*sin__S(a*a);	/* rtn. S(PI/2-|x|) */
-		}
-	}
-	if (a < small) {
-		big+a;
-		return s;		/* rtn. s*C(a) */
-	}
-	z = a*a;
-	c = cos__C(z);
-	z *= half;
-	a = (z >= thresh ? half-((z-half)-c) : one-(z-c));
-	return copysign(a,s);
-}
-END_STD(cos);
-LDBL_CLONE(cos);
diff --git a/lib/libm/noieee_src/n_sinh.c b/lib/libm/noieee_src/n_sinh.c
deleted file mode 100644
index 1fb55addabe..00000000000
--- a/lib/libm/noieee_src/n_sinh.c
+++ /dev/null
@@ -1,106 +0,0 @@
-/*	$OpenBSD: n_sinh.c,v 1.10 2009/10/27 23:59:29 deraadt Exp $	*/
-/*	$NetBSD: n_sinh.c,v 1.1 1995/10/10 23:37:05 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/* SINH(X)
- * RETURN THE HYPERBOLIC SINE OF X
- * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85;
- * REVISED BY K.C. NG on 2/8/85, 3/7/85, 3/24/85, 4/16/85.
- *
- * Required system supported functions :
- *	copysign(x,y)
- *	scalbn(x,N)
- *
- * Required kernel functions:
- *	expm1(x)	...return exp(x)-1
- *
- * Method :
- *	1. reduce x to non-negative by sinh(-x) = - sinh(x).
- *	2.
- *
- *	                                      expm1(x) + expm1(x)/(expm1(x)+1)
- *	    0 <= x <= lnovfl     : sinh(x) := --------------------------------
- *			       		                      2
- *     lnovfl <= x <= lnovfl+ln2 : sinh(x) := expm1(x)/2 (avoid overflow)
- * lnovfl+ln2 <  x <  INF        :  overflow to INF
- *
- *
- * Special cases:
- *	sinh(x) is x if x is +INF, -INF, or NaN.
- *	only sinh(0)=0 is exact for finite argument.
- *
- * Accuracy:
- *	sinh(x) returns the exact hyperbolic sine of x nearly rounded. In
- *	a test run with 1,024,000 random arguments on a VAX, the maximum
- *	observed error was 1.93 ulps (units in the last place).
- *
- * Constants:
- * The hexadecimal values are the intended ones for the following constants.
- * The decimal values may be used, provided that the compiler will convert
- * from decimal to binary accurately enough to produce the hexadecimal values
- * shown.
- */
-
-#include "math.h"
-#include "mathimpl.h"
-
-static const double mln2hi = 8.8029691931113054792E1;
-static const double mln2lo = -4.9650192275318476525E-16;
-static const double lnovfl = 8.8029691931113053016E1;
-
-#if defined(__vax__)
-static max = 126                      ;
-#else	/* defined(__vax__) */
-static max = 1023                     ;
-#endif	/* defined(__vax__) */
-
-double
-sinh(double x)
-{
-	static const double  one=1.0, half=1.0/2.0 ;
-	double t, sign;
-
-	if (isnan(x))
-		return (x);
-
-	sign=copysign(one,x);
-	x=copysign(x,one);
-	if(x<lnovfl)
-	    {t=expm1(x); return(copysign((t+t/(one+t))*half,sign));}
-
-	else if(x <= lnovfl+0.7)
-		/* subtract x by ln(2^(max+1)) and return 2^max*exp(x)
-	    		to avoid unnecessary overflow */
-	    return(copysign(scalbn(one+expm1((x-mln2hi)-mln2lo),max),sign));
-
-	else  /* sinh(+-INF) = +-INF, sinh(+-big no.) overflow to +-INF */
-	    return( expm1(x)*sign );
-}
diff --git a/lib/libm/noieee_src/n_support.c b/lib/libm/noieee_src/n_support.c
deleted file mode 100644
index 7f310f3339f..00000000000
--- a/lib/libm/noieee_src/n_support.c
+++ /dev/null
@@ -1,505 +0,0 @@
-/*	$OpenBSD: n_support.c,v 1.25 2016/09/12 19:47:02 guenther Exp $	*/
-/*	$NetBSD: n_support.c,v 1.1 1995/10/10 23:37:06 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/*
- * Some IEEE standard 754 recommended functions and remainder and sqrt for
- * supporting the C elementary functions.
- ******************************************************************************
- * WARNING:
- *      These codes are developed (in double) to support the C elementary
- * functions temporarily. They are not universal, and some of them are very
- * slow (in particular, remainder and sqrt is extremely inefficient). Each
- * computer system should have its implementation of these functions using
- * its own assembler.
- ******************************************************************************
- *
- * IEEE 754 required operations:
- *     remainder(x,p)
- *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
- *              nearest x/y; in half way case, choose the even one.
- *     sqrt(x)
- *              returns the square root of x correctly rounded according to
- *		the rounding mod.
- *
- * IEEE 754 recommended functions:
- * (a) copysign(x,y)
- *              returns x with the sign of y.
- * (b) scalbn(x,N)
- *              returns  x * (2**N), for integer values N.
- * (c) logb(x)
- *              returns the unbiased exponent of x, a signed integer in
- *              double precision, except that logb(0) is -INF, logb(INF)
- *              is +INF, and logb(NAN) is that NAN.
- *
- *
- * CODED IN C BY K.C. NG, 11/25/84;
- * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
- */
-
-#include <math.h>
-
-#include "mathimpl.h"
-
-#if defined(__vax__)      /* VAX D format */
-#include <errno.h>
-    static const unsigned short msign=0x7fff, mexp =0x7f80 ;
-    static const short  prep1=57, gap=7, bias=129           ;
-    static const double novf=1.7E38, nunf=3.0E-39;
-#else	/* defined(__vax__) */
-    static const unsigned short msign=0x7fff, mexp =0x7ff0  ;
-    static const short prep1=54, gap=4, bias=1023           ;
-    static const double novf=1.7E308, nunf=3.0E-308;
-#endif	/* defined(__vax__) */
-
-static const double zero = 0.0;
-
-double
-scalbn(double x, int N)
-{
-        int k;
-        unsigned short *px=(unsigned short *) &x;
-
-        if( x == zero )  return(x);
-
-#if defined(__vax__)
-        if( (k= *px & mexp ) != ~msign ) {
-            if (N < -260)
-		return(nunf*nunf);
-	    else if (N > 260) {
-		return(copysign(infnan(ERANGE),x));
-	    }
-#else	/* defined(__vax__) */
-        if( (k= *px & mexp ) != mexp ) {
-            if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
-            if( k == 0 ) {
-                 x *= scalbn(1.0,(int)prep1);  N -= prep1; return(scalbn(x,N));}
-#endif	/* defined(__vax__) */
-
-            if((k = (k>>gap)+ N) > 0 )
-                if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
-                else x=novf+novf;               /* overflow */
-            else
-                if( k > -prep1 )
-                                        /* gradual underflow */
-                    {*px=(*px&~mexp)|(short)(1<<gap); x *= scalbn(1.0,k-1);}
-                else
-                return(nunf*nunf);
-            }
-        return(x);
-}
-DEF_STD(scalbn);
-LDBL_CLONE(scalbn);
-
-double
-copysign(double x, double y)
-{
-        unsigned short  *px=(unsigned short *) &x,
-                        *py=(unsigned short *) &y;
-
-#if defined(__vax__)
-        if ( (*px & mexp) == 0 ) return(x);
-#endif	/* defined(__vax__) */
-
-        *px = ( *px & msign ) | ( *py & ~msign );
-        return(x);
-}
-DEF_STD(copysign);
-LDBL_CLONE(copysign);
-
-double
-logb(double x)
-{
-
-        short *px=(short *) &x, k;
-
-#if defined(__vax__)
-        return (int)(((*px&mexp)>>gap)-bias);
-#else	/* defined(__vax__) */
-        if( (k= *px & mexp ) != mexp )
-            if ( k != 0 )
-                return ( (k>>gap) - bias );
-            else if( x != zero)
-                return ( -1022.0 );
-            else
-                return(-(1.0/zero));
-        else if(isnan(x))
-            return(x);
-        else
-            {*px &= msign; return(x);}
-#endif	/* defined(__vax__) */
-}
-DEF_STD(logb);
-LDBL_UNUSED_CLONE(logb);
-
-double
-remainder(double x, double p)
-{
-        short sign;
-        double hp, dp, tmp;
-        unsigned short k;
-        unsigned short
-              *px=(unsigned short *) &x  ,
-              *pp=(unsigned short *) &p  ,
-              *pd=(unsigned short *) &dp ,
-              *pt=(unsigned short *) &tmp;
-
-        *pp &= msign ;
-
-#if defined(__vax__)
-        if( ( *px & mexp ) == ~msign )	/* is x a reserved operand? */
-#else	/* defined(__vax__) */
-        if( ( *px & mexp ) == mexp )
-#endif	/* defined(__vax__) */
-		return  (x-p)-(x-p);	/* create nan if x is inf */
-	if (p == zero) {
-#if defined(__vax__)
-		return(infnan(EDOM));
-#else	/* defined(__vax__) */
-		return zero/zero;
-#endif	/* defined(__vax__) */
-	}
-
-#if defined(__vax__)
-        if( ( *pp & mexp ) == ~msign )	/* is p a reserved operand? */
-#else	/* defined(__vax__) */
-        if( ( *pp & mexp ) == mexp )
-#endif	/* defined(__vax__) */
-		{ if (p != p) return p; else return x;}
-
-        else  if ( ((*pp & mexp)>>gap) <= 1 )
-                /* subnormal p, or almost subnormal p */
-            { double b; b=scalbn(1.0,(int)prep1);
-              p *= b; x = remainder(x,p); x *= b; return(remainder(x,p)/b);}
-        else  if ( p >= novf/2)
-            { p /= 2 ; x /= 2; return(remainder(x,p)*2);}
-        else
-            {
-                dp=p+p; hp=p/2;
-                sign= *px & ~msign ;
-                *px &= msign       ;
-                while ( x > dp )
-                    {
-                        k=(*px & mexp) - (*pd & mexp) ;
-                        tmp = dp ;
-                        *pt += k ;
-
-#if defined(__vax__)
-                        if( x < tmp ) *pt -= 128 ;
-#else	/* defined(__vax__) */
-                        if( x < tmp ) *pt -= 16 ;
-#endif	/* defined(__vax__) */
-
-                        x -= tmp ;
-                    }
-                if ( x > hp )
-                    { x -= p ;  if ( x >= hp ) x -= p ; }
-
-#if defined(__vax__)
-		if (x)
-#endif	/* defined(__vax__) */
-			*px ^= sign;
-                return( x);
-
-            }
-}
-DEF_STD(remainder);
-
-/* The drem() function is a deprecated alias for remainder(). */
-
-double
-drem(double x, double p)
-{
-	return remainder(x, p);
-}
-
-float
-sqrtf(float x)
-{
-	return (float)sqrt((double) x);
-}
-DEF_STD(sqrtf);
-
-double
-sqrt(double x)
-{
-        double q, s, b, r;
-        double t;
-        int m, n, i;
-#if defined(__vax__)
-        int k=54;
-#else	/* defined(__vax__) */
-        int k=51;
-#endif	/* defined(__vax__) */
-
-    /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
-        if(isnan(x) || x == zero) return(x);
-
-    /* sqrt(negative) is invalid */
-        if(x<zero) {
-#if defined(__vax__)
-		return (infnan(EDOM));	/* NaN */
-#else	/* defined(__vax__) */
-		return(zero/zero);
-#endif	/* defined(__vax__) */
-	}
-
-    /* sqrt(INF) is INF */
-        if(!isfinite(x)) return(x);
-
-    /* scale x to [1,4) */
-        n=logb(x);
-        x=scalbn(x,-n);
-        if((m=logb(x))!=0) x=scalbn(x,-m);       /* subnormal number */
-        m += n;
-        n = m/2;
-        if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
-
-    /* generate sqrt(x) bit by bit (accumulating in q) */
-            q=1.0; s=4.0; x -= 1.0; r=1;
-            for(i=1;i<=k;i++) {
-                t=s+1; x *= 4; r /= 2;
-                if(t<=x) {
-                    s = t+t+2;
-                    x -= t;
-                    q += r;
-                }
-                else
-                    s *= 2;
-                }
-
-    /* generate the last bit and determine the final rounding */
-            r/=2; x *= 4;
-            if(x==zero) goto end;
-	    if (100+r >= 100) {			/* trigger inexact flag */
-            if(s<x) {
-                q+=r; x -=s; s += 2; s *= 2; x *= 4;
-                t = (x-s)-5;
-                b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
-                b=1.0+r/4;   if(b>1.0) t=1;	/* b>1 : Round-to-(+INF) */
-                if(t>=0) q+=r; }	      /* else: Round-to-nearest */
-            else {
-                s *= 2; x *= 4;
-                t = (x-s)-1;
-                b=1.0+3*r/4; if(b==1.0) goto end;
-                b=1.0+r/4;   if(b>1.0) t=1;
-                if(t>=0) q+=r; }
-	    }
-
-end:        return(scalbn(q,n));
-}
-DEF_STD(sqrt);
-LDBL_CLONE(sqrt);
-
-#if 0
-/* REMAINDER(X,Y)
- * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
- * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * INTENDED FOR ASSEMBLY LANGUAGE
- * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
- *
- * Warning: this code should not get compiled in unless ALL of
- * the following machine-dependent routines are supplied.
- *
- * Required machine dependent functions (not on a VAX):
- *     swapINX(i): save inexact flag and reset it to "i"
- *     swapENI(e): save inexact enable and reset it to "e"
- */
-
-double
-remainder(double x, double y)
-{
-	static const n0=0, n1=1, n2=2, n3=3;
-    	static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
-	double hy, y1, t, t1;
-	short k;
-	long n;
-	int i, e;
-	unsigned short xexp,yexp, *px  =(unsigned short *) &x  ,
-	      		nx,nf,	  *py  =(unsigned short *) &y  ,
-	      		sign,	  *pt  =(unsigned short *) &t  ,
-	      			  *pt1 =(unsigned short *) &t1 ;
-
-	xexp = px[n0] & mexp ;	/* exponent of x */
-	yexp = py[n0] & mexp ;	/* exponent of y */
-	sign = px[n0] &0x8000;	/* sign of x     */
-
-/* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */
-	if(isnan(x)) return(x); if(isnan(y)) return(y);	/* x or y is NaN */
-	if( xexp == mexp )   return(zero/zero);		/* x is INF */
-	if(y==zero) return(y/y);
-
-/* save the inexact flag and inexact enable in i and e respectively
- * and reset them to zero
- */
-	i=swapINX(0);	e=swapENI(0);
-
-/* subnormal number */
-	nx=0;
-	if (yexp == 0) {
-		t = 1.0;
-		pt[n0] += m57;
-		y *= t;
-		nx = m57;
-	}
-
-/* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */
-	if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
-
-	nf=nx;
-	py[n0] &= 0x7fff;
-	px[n0] &= 0x7fff;
-
-/* mask off the least significant 27 bits of y */
-	t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t;
-
-/* LOOP: argument reduction on x whenever x > y */
-loop:
-	while ( x > y )
-	{
-	    t=y;
-	    t1=y1;
-	    xexp=px[n0]&mexp;	  /* exponent of x */
-	    k=xexp-yexp-m25;
-	    if(k>0) 	/* if x/y >= 2**26, scale up y so that x/y < 2**26 */
-		{pt[n0]+=k;pt1[n0]+=k;}
-	    n=x/t; x=(x-n*t1)-n*(t-t1);
-	}
-    /* end while (x > y) */
-
-	if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
-
-/* final adjustment */
-
-	hy=y/2.0;
-	if(x>hy||((x==hy)&&n%2==1)) x-=y;
-	px[n0] ^= sign;
-	if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
-
-/* restore inexact flag and inexact enable */
-	swapINX(i); swapENI(e);
-
-	return(x);
-}
-DEF_STD(remainder);
-#endif
-
-#if 0
-/* SQRT
- * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
- * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
- * CODED IN C BY K.C. NG, 3/22/85.
- *
- * Warning: this code should not get compiled in unless ALL of
- * the following machine-dependent routines are supplied.
- *
- * Required machine dependent functions:
- *     swapINX(i)  ...return the status of INEXACT flag and reset it to "i"
- *     swapRM(r)   ...return the current Rounding Mode and reset it to "r"
- *     swapENI(e)  ...return the status of inexact enable and reset it to "e"
- *     addc(t)     ...perform t=t+1 regarding t as a 64 bit unsigned integer
- *     subc(t)     ...perform t=t-1 regarding t as a 64 bit unsigned integer
- */
-
-static const unsigned long table[] = {
-0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
-58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
-21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };
-
-double
-newsqrt(double x)
-{
-        double y, z, t, addc(), subc();
-	double const b54=134217728.*134217728.; /* b54=2**54 */
-        long mx, scalx;
-	long const mexp=0x7ff00000;
-        int i, j, r, e, swapINX(), swapRM(), swapENI();
-        unsigned long *py=(unsigned long *) &y   ,
-                      *pt=(unsigned long *) &t   ,
-                      *px=(unsigned long *) &x   ;
-        const int n0=0, n1=1;
-/* Rounding Mode:  RN ...round-to-nearest
- *                 RZ ...round-towards 0
- *                 RP ...round-towards +INF
- *		   RM ...round-towards -INF
- */
-        const int RN=0, RZ=1, RP=2, RM=3;
-				/* machine dependent: work on a Zilog Z8070
-                                 * and a National 32081 & 16081
-                                 */
-
-/* exceptions */
-	if(isnan(x) || x == 0.0) return(x);	/* sqrt(NaN) is NaN,
-						   sqrt(+-0) = +-0 */
-	if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
-        if((mx=px[n0]&mexp)==mexp) return(x);  /* sqrt(+INF) is +INF */
-
-/* save, reset, initialize */
-        e=swapENI(0);   /* ...save and reset the inexact enable */
-        i=swapINX(0);   /* ...save INEXACT flag */
-        r=swapRM(RN);   /* ...save and reset the Rounding Mode to RN */
-        scalx=0;
-
-/* subnormal number, scale up x to x*2**54 */
-        if(mx==0) {x *= b54 ; scalx-=0x01b00000;}
-
-/* scale x to avoid intermediate over/underflow:
- * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
-        if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
-        if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}
-
-/* magic initial approximation to almost 8 sig. bits */
-        py[n0]=(px[n0]>>1)+0x1ff80000;
-        py[n0]=py[n0]-table[(py[n0]>>15)&31];
-
-/* Heron's rule once with correction to improve y to almost 18 sig. bits */
-        t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
-
-/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
-        t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y;
-        t=z/(t+x) ;  pt[n0]+=0x00100000; y+=t;
-
-/* twiddle last bit to force y correctly rounded */
-        swapRM(RZ);     /* ...set Rounding Mode to round-toward-zero */
-        swapINX(0);     /* ...clear INEXACT flag */
-        swapENI(e);     /* ...restore inexact enable status */
-        t=x/y;          /* ...chopped quotient, possibly inexact */
-        j=swapINX(i);   /* ...read and restore inexact flag */
-        if(j==0) { if(t==y) goto end; else t=subc(t); }  /* ...t=t-ulp */
-        x=b54+0.1;      /* ..trigger inexact flag, sqrt(x) is inexact */
-        if(r==RN) t=addc(t);            /* ...t=t+ulp */
-        else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */
-        y=y+t;                          /* ...chopped sum */
-        py[n0]=py[n0]-0x00100000;       /* ...correctly rounded sqrt(x) */
-end:    py[n0]=py[n0]+scalx;            /* ...scale back y */
-        swapRM(r);                      /* ...restore Rounding Mode */
-        return(y);
-}
-#endif
diff --git a/lib/libm/noieee_src/n_tan.c b/lib/libm/noieee_src/n_tan.c
deleted file mode 100644
index 971ebe7a89f..00000000000
--- a/lib/libm/noieee_src/n_tan.c
+++ /dev/null
@@ -1,74 +0,0 @@
-/*	$OpenBSD: n_tan.c,v 1.16 2016/09/12 19:47:02 guenther Exp $	*/
-/*	$NetBSD: n_tan.c,v 1.1 1995/10/10 23:37:07 ragge Exp $	*/
-/*
- * Copyright (c) 1987, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include <math.h>
-
-#include "mathimpl.h"
-
-float
-tanf(float x)
-{
-	return (float)tan((double) x);
-}
-
-double
-tan(double x)
-{
-	double a,z,ss,cc,c;
-	int k;
-
-	if(!isfinite(x))	/* tan(NaN) and tan(INF) must be NaN */
-		return x-x;
-	x = remainder(x,PI);	/* reduce x into [-PI/2, PI/2] */
-	a = copysign(x,one);	/* ... = abs(x) */
-	if (a >= PIo4) {
-		k = 1;
-		x = copysign(PIo2-a,x);
-	}
-	else {
-		k = 0;
-		if (a < small) {
-			big+a;
-			return x;
-		}
-	}
-	z = x*x;
-	cc = cos__C(z);
-	ss = sin__S(z);
-	z *= half;			/* Next get c = cos(x) accurately */
-	c = (z >= thresh ? half-((z-half)-cc) : one-(z-cc));
-	if (k == 0)
-		return x+(x*(z-(cc-ss)))/c;	/* ... sin/cos */
-	else
-		return c/(x+x*ss);		/* ... cos/sin */
-}
-DEF_STD(tan);
-LDBL_UNUSED_CLONE(tan);
diff --git a/lib/libm/noieee_src/n_tanh.c b/lib/libm/noieee_src/n_tanh.c
deleted file mode 100644
index 3db4a516e63..00000000000
--- a/lib/libm/noieee_src/n_tanh.c
+++ /dev/null
@@ -1,95 +0,0 @@
-/*	$OpenBSD: n_tanh.c,v 1.10 2016/09/12 04:39:47 guenther Exp $	*/
-/*	$NetBSD: n_tanh.c,v 1.1 1995/10/10 23:37:08 ragge Exp $	*/
-/*
- * Copyright (c) 1985, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#include "math.h"
-
-/* TANH(X)
- * RETURN THE HYPERBOLIC TANGENT OF X
- * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85;
- * REVISED BY K.C. NG on 2/8/85, 2/11/85, 3/7/85, 3/24/85.
- *
- * Required system supported functions :
- *	copysign(x,y)
- *	isfinite(x)
- *
- * Required kernel function:
- *	expm1(x)	...exp(x)-1
- *
- * Method :
- *	1. reduce x to non-negative by tanh(-x) = - tanh(x).
- *	2.
- *	    0      <  x <=  1.e-10 :  tanh(x) := x
- *					          -expm1(-2x)
- *	    1.e-10 <  x <=  1      :  tanh(x) := --------------
- *					         expm1(-2x) + 2
- *							  2
- *	    1      <= x <=  22.0   :  tanh(x) := 1 -  ---------------
- *						      expm1(2x) + 2
- *	    22.0   <  x <= INF     :  tanh(x) := 1.
- *
- *	Note: 22 was chosen so that fl(1.0+2/(expm1(2*22)+2)) == 1.
- *
- * Special cases:
- *	tanh(NaN) is NaN;
- *	only tanh(0)=0 is exact for finite argument.
- *
- * Accuracy:
- *	tanh(x) returns the exact hyperbolic tangent of x nealy rounded.
- *	In a test run with 1,024,000 random arguments on a VAX, the maximum
- *	observed error was 2.22 ulps (units in the last place).
- */
-
-double
-tanh(double x)
-{
-	static double one=1.0, two=2.0, small = 1.0e-10, big = 1.0e10;
-	double t, sign;
-
-	if (isnan(x))
-		return (x);
-
-	sign=copysign(one,x);
-	x=copysign(x,one);
-	if(x < 22.0)
-	    if( x > one )
-		return(copysign(one-two/(expm1(x+x)+two),sign));
-	    else if ( x > small )
-		{t= -expm1(-(x+x)); return(copysign(t/(two-t),sign));}
-	    else {		/* raise the INEXACT flag for non-zero x */
-		t = big + x;
-		return(copysign(x,sign) - (t-(big+x)));
-	    }
-	else if(isfinite(x))
-	    return (sign+1.0E-37); /* raise the INEXACT flag */
-	else
-	    return(sign);	/* x is +- INF */
-}
diff --git a/lib/libm/noieee_src/n_tgamma.c b/lib/libm/noieee_src/n_tgamma.c
deleted file mode 100644
index 1bd2f13bbee..00000000000
--- a/lib/libm/noieee_src/n_tgamma.c
+++ /dev/null
@@ -1,363 +0,0 @@
-/*	$OpenBSD: n_tgamma.c,v 1.5 2016/09/12 04:39:47 guenther Exp $	*/
-/*-
- * Copyright (c) 1992, 1993
- *	The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-/*
- * This code by P. McIlroy, Oct 1992;
- *
- * The financial support of UUNET Communications Services is greatfully
- * acknowledged.
- */
-
-#include <math.h>
-#include "mathimpl.h"
-#include <errno.h>
-
-/* METHOD:
- * x < 0: Use reflection formula, G(x) = pi/(sin(pi*x)*x*G(x))
- *	At negative integers, return +Inf and raise overflow.
- *
- * x < 6.5:
- *	Use argument reduction G(x+1) = xG(x) to reach the
- *	range [1.066124,2.066124].  Use a rational
- *	approximation centered at the minimum (x0+1) to
- *	ensure monotonicity.
- *
- * x >= 6.5: Use the asymptotic approximation (Stirling's formula)
- *	adjusted for equal-ripples:
- *
- *	log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + 1/x*P(1/(x*x))
- *
- *	Keep extra precision in multiplying (x-.5)(log(x)-1), to
- *	avoid premature round-off.
- *
- * Special values:
- *	-Inf:			return 0 and raise invalid;
- *	negative integer:	return +Inf and raise overflow;
- *	other x ~< -177.79:	return 0 and raise underflow;
- *	+-0:			return +Inf and raise overflow;
- *	finite x ~> 171.63:	return +Inf and raise div-by-0;
- *	+Inf:			return +Inf and raise div-by-0;
- *	NaN:			return 0 and raise invalid.
- *
- * Accuracy: tgamma(x) is accurate to within
- *	x > 0:  error provably < 0.9ulp.
- *	Maximum observed in 1,000,000 trials was .87ulp.
- *	x < 0:
- *	Maximum observed error < 4ulp in 1,000,000 trials.
- */
-
-static double neg_gam(double);
-static double small_gam(double);
-static double smaller_gam(double);
-static struct Double large_gam(double);
-static struct Double ratfun_gam(double, double);
-
-/*
- * Rational approximation, A0 + x*x*P(x)/Q(x), on the interval
- * [1.066.., 2.066..] accurate to 4.25e-19.
- */
-#define LEFT -.3955078125	/* left boundary for rat. approx */
-#define x0 .461632144968362356785	/* xmin - 1 */
-
-#define a0_hi 0.88560319441088874992
-#define a0_lo -.00000000000000004996427036469019695
-#define P0	 6.21389571821820863029017800727e-01
-#define P1	 2.65757198651533466104979197553e-01
-#define P2	 5.53859446429917461063308081748e-03
-#define P3	 1.38456698304096573887145282811e-03
-#define P4	 2.40659950032711365819348969808e-03
-#define Q0	 1.45019531250000000000000000000e+00
-#define Q1	 1.06258521948016171343454061571e+00
-#define Q2	-2.07474561943859936441469926649e-01
-#define Q3	-1.46734131782005422506287573015e-01
-#define Q4	 3.07878176156175520361557573779e-02
-#define Q5	 5.12449347980666221336054633184e-03
-#define Q6	-1.76012741431666995019222898833e-03
-#define Q7	 9.35021023573788935372153030556e-05
-#define Q8	 6.13275507472443958924745652239e-06
-/*
- * Constants for large x approximation (x in [6, Inf])
- * (Accurate to 2.8*10^-19 absolute)
- */
-#define lns2pi_hi 0.418945312500000
-#define lns2pi_lo -.000006779295327258219670263595
-#define Pa0	 8.33333333333333148296162562474e-02
-#define Pa1	-2.77777777774548123579378966497e-03
-#define Pa2	 7.93650778754435631476282786423e-04
-#define Pa3	-5.95235082566672847950717262222e-04
-#define Pa4	 8.41428560346653702135821806252e-04
-#define Pa5	-1.89773526463879200348872089421e-03
-#define Pa6	 5.69394463439411649408050664078e-03
-#define Pa7	-1.44705562421428915453880392761e-02
-
-static const double zero = 0., one = 1.0, tiny = 1e-300;
-
-/*
- * TRUNC sets trailing bits in a floating-point number to zero.
- * is a temporary variable.
- */
-#if defined(__vax__)
-#define _IEEE		0
-#define TRUNC(x)	x = (double) (float) (x)
-#else
-static int endian;
-#define _IEEE		1
-#define TRUNC(x)	*(((int *) &x) + endian) &= 0xf8000000
-#define infnan(x)	0.0
-#endif
-
-double
-tgamma(double x)
-{
-	struct Double u;
-#if _IEEE
-	endian = (*(int *) &one) ? 1 : 0;
-#endif
-
-	if (x >= 6) {
-		if(x > 171.63)
-			if (_IEEE)
-				return (x/zero);
-			else
-				return (infnan(ERANGE));
-		u = large_gam(x);
-		return(__exp__D(u.a, u.b));
-	} else if (x >= 1.0 + LEFT + x0)
-		return (small_gam(x));
-	else if (x > 1.e-17)
-		return (smaller_gam(x));
-	else if (x > -1.e-17) {
-		if (x == 0.0) {
-			if (!_IEEE)
-				return (infnan(ERANGE));
-		} else {
-			u.a = one - tiny;	/* raise inexact */
-		}
-		return (one/x);
-	} else if (!isfinite(x)) {
-		if (_IEEE)		/* x = NaN, -Inf */
-			return (x - x);
-		else
-			return (infnan(EDOM));
-	 } else
-		return (neg_gam(x));
-}
-
-/*
- * We simply call tgamma() rather than bloating the math library
- * with a float-optimized version of it.  The reason is that tgammaf()
- * is essentially useless, since the function is superexponential
- * and floats have very limited range.  -- das@freebsd.org
- */
-
-float
-tgammaf(float x)
-{
-	return tgamma(x);
-}
-
-/*
- * Accurate to max(ulp(1/128) absolute, 2^-66 relative) error.
- */
-
-static struct Double
-large_gam(double x)
-{
-	double z, p;
-	struct Double t, u, v;
-
-	z = one/(x*x);
-	p = Pa0+z*(Pa1+z*(Pa2+z*(Pa3+z*(Pa4+z*(Pa5+z*(Pa6+z*Pa7))))));
-	p = p/x;
-
-	u = __log__D(x);
-	u.a -= one;
-	v.a = (x -= .5);
-	TRUNC(v.a);
-	v.b = x - v.a;
-	t.a = v.a*u.a;			/* t = (x-.5)*(log(x)-1) */
-	t.b = v.b*u.a + x*u.b;
-	/* return t.a + t.b + lns2pi_hi + lns2pi_lo + p */
-	t.b += lns2pi_lo; t.b += p;
-	u.a = lns2pi_hi + t.b; u.a += t.a;
-	u.b = t.a - u.a;
-	u.b += lns2pi_hi; u.b += t.b;
-	return (u);
-}
-
-/*
- * Good to < 1 ulp.  (provably .90 ulp; .87 ulp on 1,000,000 runs.)
- * It also has correct monotonicity.
- */
-
-static double
-small_gam(double x)
-{
-	double y, ym1, t;
-	struct Double yy, r;
-	y = x - one;
-	ym1 = y - one;
-	if (y <= 1.0 + (LEFT + x0)) {
-		yy = ratfun_gam(y - x0, 0);
-		return (yy.a + yy.b);
-	}
-	r.a = y;
-	TRUNC(r.a);
-	yy.a = r.a - one;
-	y = ym1;
-	yy.b = r.b = y - yy.a;
-	/* Argument reduction: G(x+1) = x*G(x) */
-	for (ym1 = y-one; ym1 > LEFT + x0; y = ym1--, yy.a--) {
-		t = r.a*yy.a;
-		r.b = r.a*yy.b + y*r.b;
-		r.a = t;
-		TRUNC(r.a);
-		r.b += (t - r.a);
-	}
-	/* Return r*tgamma(y). */
-	yy = ratfun_gam(y - x0, 0);
-	y = r.b*(yy.a + yy.b) + r.a*yy.b;
-	y += yy.a*r.a;
-	return (y);
-}
-
-/*
- * Good on (0, 1+x0+LEFT].  Accurate to 1ulp.
- */
-
-static double
-smaller_gam(double x)
-{
-	double t, d;
-	struct Double r, xx;
-	if (x < x0 + LEFT) {
-		t = x;
-		TRUNC(t);
-		d = (t+x)*(x-t);
-		t *= t;
-		xx.a = (t + x);
-		TRUNC(xx.a);
-		xx.b = x - xx.a; xx.b += t; xx.b += d;
-		t = (one-x0); t += x;
-		d = (one-x0); d -= t; d += x;
-		x = xx.a + xx.b;
-	} else {
-		xx.a =  x;
-		TRUNC(xx.a);
-		xx.b = x - xx.a;
-		t = x - x0;
-		d = (-x0 -t); d += x;
-	}
-	r = ratfun_gam(t, d);
-	d = r.a/x;
-	TRUNC(d);
-	r.a -= d*xx.a; r.a -= d*xx.b; r.a += r.b;
-	return (d + r.a/x);
-}
-
-/*
- * returns (z+c)^2 * P(z)/Q(z) + a0
- */
-
-static struct Double
-ratfun_gam(double z, double c)
-{
-	double p, q;
-	struct Double r, t;
-
-	q = Q0 +z*(Q1+z*(Q2+z*(Q3+z*(Q4+z*(Q5+z*(Q6+z*(Q7+z*Q8)))))));
-	p = P0 + z*(P1 + z*(P2 + z*(P3 + z*P4)));
-
-	/* return r.a + r.b = a0 + (z+c)^2*p/q, with r.a truncated to 26 bits. */
-	p = p/q;
-	t.a = z;
-	TRUNC(t.a);			/* t ~= z + c */
-	t.b = (z - t.a) + c;
-	t.b *= (t.a + z);
-	q = (t.a *= t.a);		/* t = (z+c)^2 */
-	TRUNC(t.a);
-	t.b += (q - t.a);
-	r.a = p;
-	TRUNC(r.a);			/* r = P/Q */
-	r.b = p - r.a;
-	t.b = t.b*p + t.a*r.b + a0_lo;
-	t.a *= r.a;			/* t = (z+c)^2*(P/Q) */
-	r.a = t.a + a0_hi;
-	TRUNC(r.a);
-	r.b = ((a0_hi-r.a) + t.a) + t.b;
-	return (r);			/* r = a0 + t */
-}
-
-static double
-neg_gam(double x)
-{
-	int sgn = 1;
-	struct Double lg, lsine;
-	double y, z;
-
-	y = ceil(x);
-	if (y == x)		/* Negative integer. */
-		if (_IEEE)
-			return ((x - x) / zero);
-		else
-			return (infnan(ERANGE));
-	z = y - x;
-	if (z > 0.5)
-		z = one - z;
-	y = 0.5 * y;
-	if (y == ceil(y))
-		sgn = -1;
-	if (z < .25)
-		z = sin(M_PI*z);
-	else
-		z = cos(M_PI*(0.5-z));
-	/* Special case: G(1-x) = Inf; G(x) may be nonzero. */
-	if (x < -170) {
-		if (x < -190)
-			return ((double)sgn*tiny*tiny);
-		y = one - x;		/* exact: 128 < |x| < 255 */
-		lg = large_gam(y);
-		lsine = __log__D(M_PI/z);	/* = TRUNC(log(u)) + small */
-		lg.a -= lsine.a;		/* exact (opposite signs) */
-		lg.b -= lsine.b;
-		y = -(lg.a + lg.b);
-		z = (y + lg.a) + lg.b;
-		y = __exp__D(y, z);
-		if (sgn < 0) y = -y;
-		return (y);
-	}
-	y = one-x;
-	if (one-y == x)
-		y = tgamma(y);
-	else		/* 1-x is inexact */
-		y = -x*tgamma(-x);
-	if (sgn < 0) y = -y;
-	return (M_PI / (y*z));
-}