From 532fe1553e4803a0e205616ac713e5ec8014f7e4 Mon Sep 17 00:00:00 2001
From: Ed Schouten <ed@nuxi.nl>
Date: Thu, 12 Feb 2015 13:25:30 +0100
Subject: [PATCH 1/9] Remove the non-standard gamma() function. Use lgamma()
 instead.

---
 include/openlibm_math.h | 16 ++++------------
 src/Make.files          |  4 ++--
 src/common.c            |  2 ++
 src/e_gamma.c           | 36 ------------------------------------
 src/e_gamma_r.c         | 33 ---------------------------------
 src/e_gammaf.c          | 37 -------------------------------------
 src/e_gammaf_r.c        | 34 ----------------------------------
 src/math_private.h      |  4 ----
 8 files changed, 8 insertions(+), 158 deletions(-)
 delete mode 100644 src/e_gamma.c
 delete mode 100644 src/e_gamma_r.c
 delete mode 100644 src/e_gammaf.c
 delete mode 100644 src/e_gammaf_r.c

diff --git a/include/openlibm_math.h b/include/openlibm_math.h
index 8605c8b..a3eb614 100644
--- a/include/openlibm_math.h
+++ b/include/openlibm_math.h
@@ -281,10 +281,6 @@ double	y0(double);
 double	y1(double);
 double	yn(int, double);
 
-#if __XSI_VISIBLE <= 500 || __BSD_VISIBLE
-double	gamma(double);
-#endif
-
 #if __XSI_VISIBLE <= 600 || __BSD_VISIBLE
 double	scalb(double, double);
 #endif
@@ -312,10 +308,9 @@ int	finite(double) __pure2;
 int	isnanf(float) __pure2;
 
 /*
- * Reentrant version of gamma & lgamma; passes signgam back by reference
- * as the second argument; user must allocate space for signgam.
+ * Reentrant version of lgamma; passes signgam back by reference as the
+ * second argument; user must allocate space for signgam.
  */
-double	gamma_r(double, int *);
 double	lgamma_r(double, int *);
 
 /*
@@ -402,7 +397,6 @@ float	fminf(float, float) __pure2;
 #if __BSD_VISIBLE
 float	dremf(float, float);
 int	finitef(float) __pure2;
-float	gammaf(float);
 float	j0f(float);
 float	j1f(float);
 float	jnf(int, float);
@@ -412,11 +406,9 @@ float	y1f(float);
 float	ynf(int, float);
 
 /*
- * Float versions of reentrant version of gamma & lgamma; passes
- * signgam back by reference as the second argument; user must
- * allocate space for signgam.
+ * Float versions of reentrant version of lgamma; passes signgam back by
+ * reference as the second argument; user must allocate space for signgam.
  */
-float	gammaf_r(float, int *);
 float	lgammaf_r(float, int *);
 
 /*
diff --git a/src/Make.files b/src/Make.files
index edd64e4..5c75e93 100644
--- a/src/Make.files
+++ b/src/Make.files
@@ -1,8 +1,8 @@
 $(CUR_SRCS) = common.c \
 	e_acos.c e_acosf.c e_acosh.c e_acoshf.c e_asin.c e_asinf.c \
 	e_atan2.c e_atan2f.c e_atanh.c e_atanhf.c e_cosh.c e_coshf.c e_exp.c \
-	e_expf.c e_fmod.c e_fmodf.c e_gamma.c e_gamma_r.c e_gammaf.c \
-	e_gammaf_r.c e_hypot.c e_hypotf.c e_j0.c e_j0f.c e_j1.c e_j1f.c \
+	e_expf.c e_fmod.c e_fmodf.c \
+	e_hypot.c e_hypotf.c e_j0.c e_j0f.c e_j1.c e_j1f.c \
 	e_jn.c e_jnf.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \
 	e_lgammal.c e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c \
 	e_pow.c e_powf.c e_remainder.c e_remainderf.c e_scalb.c e_scalbf.c \
diff --git a/src/common.c b/src/common.c
index 01df0db..5905c20 100644
--- a/src/common.c
+++ b/src/common.c
@@ -1,5 +1,7 @@
 #include "math_private.h"
 
+DLLEXPORT int isopenlibm(void);
+
 DLLEXPORT int isopenlibm(void) {
     return 1;
 }
diff --git a/src/e_gamma.c b/src/e_gamma.c
deleted file mode 100644
index ea6e7e0..0000000
--- a/src/e_gamma.c
+++ /dev/null
@@ -1,36 +0,0 @@
-
-/* @(#)e_gamma.c 1.3 95/01/18 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunSoft, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
- * is preserved.
- * ====================================================
- *
- */
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/e_gamma.c,v 1.8 2008/02/22 02:30:34 das Exp $");
-
-/* __ieee754_gamma(x)
- * Return the logarithm of the Gamma function of x.
- *
- * Method: call __ieee754_gamma_r
- */
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-DLLEXPORT double
-__ieee754_gamma(double x)
-{
-#ifdef OPENLIBM_ONLY_THREAD_SAFE
-	int signgam;
-#endif
-
-	return __ieee754_gamma_r(x,&signgam);
-}
diff --git a/src/e_gamma_r.c b/src/e_gamma_r.c
deleted file mode 100644
index 36a92a1..0000000
--- a/src/e_gamma_r.c
+++ /dev/null
@@ -1,33 +0,0 @@
-
-/* @(#)e_gamma_r.c 1.3 95/01/18 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunSoft, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
- * is preserved.
- * ====================================================
- *
- */
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/e_gamma_r.c,v 1.8 2008/02/22 02:30:34 das Exp $");
-
-/* __ieee754_gamma_r(x, signgamp)
- * Reentrant version of the logarithm of the Gamma function 
- * with user provide pointer for the sign of Gamma(x). 
- *
- * Method: See __ieee754_lgamma_r
- */
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-DLLEXPORT double
-__ieee754_gamma_r(double x, int *signgamp)
-{
-	return __ieee754_lgamma_r(x,signgamp);
-}
diff --git a/src/e_gammaf.c b/src/e_gammaf.c
deleted file mode 100644
index a4c88ad..0000000
--- a/src/e_gammaf.c
+++ /dev/null
@@ -1,37 +0,0 @@
-/* e_gammaf.c -- float version of e_gamma.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/e_gammaf.c,v 1.7 2008/02/22 02:30:35 das Exp $");
-
-/* __ieee754_gammaf(x)
- * Return the logarithm of the Gamma function of x.
- *
- * Method: call __ieee754_gammaf_r
- */
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-DLLEXPORT float
-__ieee754_gammaf(float x)
-{
-#ifdef OPENLIBM_ONLY_THREAD_SAFE
-	int signgam;
-#endif
-
-	return __ieee754_gammaf_r(x,&signgam);
-}
diff --git a/src/e_gammaf_r.c b/src/e_gammaf_r.c
deleted file mode 100644
index 7041118..0000000
--- a/src/e_gammaf_r.c
+++ /dev/null
@@ -1,34 +0,0 @@
-/* e_gammaf_r.c -- float version of e_gamma_r.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/e_gammaf_r.c,v 1.8 2008/02/22 02:30:35 das Exp $");
-
-/* __ieee754_gammaf_r(x, signgamp)
- * Reentrant version of the logarithm of the Gamma function
- * with user provide pointer for the sign of Gamma(x).
- *
- * Method: See __ieee754_lgammaf_r
- */
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-DLLEXPORT float
-__ieee754_gammaf_r(float x, int *signgamp)
-{
-	return __ieee754_lgammaf_r(x,signgamp);
-}
diff --git a/src/math_private.h b/src/math_private.h
index 1ab2309..7f96cf7 100644
--- a/src/math_private.h
+++ b/src/math_private.h
@@ -280,9 +280,7 @@ irint(double x)
 #define	__ieee754_fmod	fmod
 #define	__ieee754_pow	pow
 #define	__ieee754_lgamma lgamma
-#define	__ieee754_gamma	gamma
 #define	__ieee754_lgamma_r lgamma_r
-#define	__ieee754_gamma_r gamma_r
 #define	__ieee754_log10	log10
 #define	__ieee754_sinh	sinh
 #define	__ieee754_hypot	hypot
@@ -306,9 +304,7 @@ irint(double x)
 #define	__ieee754_fmodf	fmodf
 #define	__ieee754_powf	powf
 #define	__ieee754_lgammaf lgammaf
-#define	__ieee754_gammaf gammaf
 #define	__ieee754_lgammaf_r lgammaf_r
-#define	__ieee754_gammaf_r gammaf_r
 #define	__ieee754_log10f log10f
 #define	__ieee754_log2f log2f
 #define	__ieee754_sinhf	sinhf

From eca4ac8bf1cbaf7986260afcc12cab542c5f9a2b Mon Sep 17 00:00:00 2001
From: Ed Schouten <ed@nuxi.nl>
Date: Thu, 12 Feb 2015 13:28:37 +0100
Subject: [PATCH 2/9] Remove the non-standard significand() function.

---
 i387/Make.files         |  7 +++----
 i387/s_significand.S    | 21 ---------------------
 i387/s_significandf.S   | 22 ----------------------
 include/openlibm_math.h | 10 ----------
 src/Make.files          |  2 +-
 src/s_significand.c     | 30 ------------------------------
 src/s_significandf.c    | 27 ---------------------------
 7 files changed, 4 insertions(+), 115 deletions(-)
 delete mode 100644 i387/s_significand.S
 delete mode 100644 i387/s_significandf.S
 delete mode 100644 src/s_significand.c
 delete mode 100644 src/s_significandf.c

diff --git a/i387/Make.files b/i387/Make.files
index 3189957..abfb4fe 100644
--- a/i387/Make.files
+++ b/i387/Make.files
@@ -1,8 +1,7 @@
 $(CUR_SRCS) = e_exp.S e_fmod.S e_log.S e_log10.S \
 	    e_remainder.S e_sqrt.S s_ceil.S s_copysign.S \
 	    s_finite.S s_floor.S s_llrint.S s_logb.S s_lrint.S \
-	    s_remquo.S s_rint.S  s_significand.S s_tan.S \
-	    s_trunc.S
+	    s_remquo.S s_rint.S  s_tan.S s_trunc.S
 
 ifneq ($(OS), WINNT)
 $(CUR_SRCS) += s_scalbn.S s_scalbnf.S s_scalbnl.S
@@ -12,11 +11,11 @@ endif
 $(CUR_SRCS)+= e_log10f.S e_logf.S e_remainderf.S \
 	    e_sqrtf.S s_ceilf.S s_copysignf.S s_floorf.S \
 	    s_llrintf.S s_logbf.S s_lrintf.S \
-	    s_remquof.S s_rintf.S s_significandf.S s_truncf.S
+	    s_remquof.S s_rintf.S s_truncf.S
 
 # long double counterparts
 $(CUR_SRCS)+= e_remainderl.S e_sqrtl.S s_ceill.S s_copysignl.S \
 	    s_floorl.S s_llrintl.S \
 	    s_logbl.S s_lrintl.S s_remquol.S s_rintl.S s_truncl.S
 
-$(CUR_SRCS)+= fenv.c
\ No newline at end of file
+$(CUR_SRCS)+= fenv.c
diff --git a/i387/s_significand.S b/i387/s_significand.S
deleted file mode 100644
index b8552fe..0000000
--- a/i387/s_significand.S
+++ /dev/null
@@ -1,21 +0,0 @@
-/*
- * Written by:
- *	J.T. Conklin (jtc@netbsd.org)
- *	Public domain.
- */
-
-#include <i387/bsd_asm.h>
-//__FBSDID("$FreeBSD: src/lib/msun/i387/s_significand.S,v 1.10 2011/01/07 16:13:12 kib Exp $")
-
-ENTRY(significand)
-	fldl	4(%esp)
-	fxtract
-	fstp	%st(1)
-	ret
-END(significand)
-
-	
-/* Enable stack protection */
-#if defined(__linux__) && defined(__ELF__)
-.section .note.GNU-stack,"",%progbits
-#endif
diff --git a/i387/s_significandf.S b/i387/s_significandf.S
deleted file mode 100644
index eca1679..0000000
--- a/i387/s_significandf.S
+++ /dev/null
@@ -1,22 +0,0 @@
-/*
- * Written by J.T. Conklin <jtc@netbsd.org>.
- * Public domain.
- */
-
-#include <i387/bsd_asm.h>
-
-//__FBSDID("$FreeBSD: src/lib/msun/i387/s_significandf.S,v 1.3 2011/01/07 16:13:12 kib Exp $");
-/* RCSID("$NetBSD: s_significandf.S,v 1.3 1995/05/09 00:24:07 jtc Exp $") */
-
-ENTRY(significandf)
-	flds	4(%esp)
-	fxtract
-	fstp	%st(1)
-	ret
-END(significandf)
-
-	
-/* Enable stack protection */
-#if defined(__linux__) && defined(__ELF__)
-.section .note.GNU-stack,"",%progbits
-#endif
diff --git a/include/openlibm_math.h b/include/openlibm_math.h
index a3eb614..2113390 100644
--- a/include/openlibm_math.h
+++ b/include/openlibm_math.h
@@ -317,11 +317,6 @@ double	lgamma_r(double, int *);
  * Single sine/cosine function.
  */
 void	sincos(double, double *, double *);
-
-/*
- * IEEE Test Vector
- */
-double	significand(double);
 #endif /* __BSD_VISIBLE */
 
 /* float versions of ANSI/POSIX functions */
@@ -415,11 +410,6 @@ float	lgammaf_r(float, int *);
  * Single sine/cosine function.
  */
 void	sincosf(float, float *, float *);
-
-/*
- * float version of IEEE Test Vector
- */
-float	significandf(float);
 #endif	/* __BSD_VISIBLE */
 
 /*
diff --git a/src/Make.files b/src/Make.files
index 5c75e93..dc930e9 100644
--- a/src/Make.files
+++ b/src/Make.files
@@ -28,7 +28,7 @@ $(CUR_SRCS) = common.c \
 	s_nexttowardf.c s_remquo.c s_remquof.c \
 	s_rint.c s_rintf.c s_round.c s_roundf.c s_roundl.c \
 	s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
-	s_signgam.c s_significand.c s_significandf.c s_sin.c s_sincos.c \
+	s_signgam.c s_sin.c s_sincos.c \
 	s_sinf.c s_sincosf.c s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c \
 	s_trunc.c s_truncf.c s_cpow.c  s_cpowf.c \
 	w_cabs.c w_cabsf.c w_drem.c w_dremf.c
diff --git a/src/s_significand.c b/src/s_significand.c
deleted file mode 100644
index 61045bc..0000000
--- a/src/s_significand.c
+++ /dev/null
@@ -1,30 +0,0 @@
-/* @(#)s_signif.c 5.1 93/09/24 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/s_significand.c,v 1.10 2008/02/22 02:30:35 das Exp $");
-
-/*
- * significand(x) computes just
- * 	scalb(x, (double) -ilogb(x)),
- * for exercising the fraction-part(F) IEEE 754-1985 test vector.
- */
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-DLLEXPORT double
-significand(double x)
-{
-	return __ieee754_scalb(x,(double) -ilogb(x));
-}
diff --git a/src/s_significandf.c b/src/s_significandf.c
deleted file mode 100644
index 0ccb9de..0000000
--- a/src/s_significandf.c
+++ /dev/null
@@ -1,27 +0,0 @@
-/* s_significandf.c -- float version of s_significand.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/s_significandf.c,v 1.8 2008/02/22 02:30:36 das Exp $");
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-DLLEXPORT float
-significandf(float x)
-{
-	return __ieee754_scalbf(x,(float) -ilogbf(x));
-}

From ae1e0c309a85ef48822ef184e3a67f5150d436c1 Mon Sep 17 00:00:00 2001
From: Ed Schouten <ed@nuxi.nl>
Date: Thu, 12 Feb 2015 13:29:38 +0100
Subject: [PATCH 3/9] Remove the non-standard drem() function.

---
 include/openlibm_math.h |  1 -
 src/Make.files          |  2 +-
 src/w_drem.c            | 16 ----------------
 src/w_dremf.c           | 17 -----------------
 4 files changed, 1 insertion(+), 35 deletions(-)
 delete mode 100644 src/w_drem.c
 delete mode 100644 src/w_dremf.c

diff --git a/include/openlibm_math.h b/include/openlibm_math.h
index 2113390..faa30de 100644
--- a/include/openlibm_math.h
+++ b/include/openlibm_math.h
@@ -303,7 +303,6 @@ double	trunc(double);
  * BSD math library entry points
  */
 #if __BSD_VISIBLE
-double	drem(double, double);
 int	finite(double) __pure2;
 int	isnanf(float) __pure2;
 
diff --git a/src/Make.files b/src/Make.files
index dc930e9..8c8bf07 100644
--- a/src/Make.files
+++ b/src/Make.files
@@ -31,7 +31,7 @@ $(CUR_SRCS) = common.c \
 	s_signgam.c s_sin.c s_sincos.c \
 	s_sinf.c s_sincosf.c s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c \
 	s_trunc.c s_truncf.c s_cpow.c  s_cpowf.c \
-	w_cabs.c w_cabsf.c w_drem.c w_dremf.c
+	w_cabs.c w_cabsf.c
 
 ifneq ($(OS), WINNT)
 $(CUR_SRCS) += s_nan.c
diff --git a/src/w_drem.c b/src/w_drem.c
deleted file mode 100644
index 8fe000a..0000000
--- a/src/w_drem.c
+++ /dev/null
@@ -1,16 +0,0 @@
-/*
- * drem() wrapper for remainder().
- *
- * Written by J.T. Conklin, <jtc@wimsey.com>
- * Placed into the Public Domain, 1994.
- */
-
-#include <openlibm_math.h>
-#include "math_private.h"
-
-DLLEXPORT double
-drem(x, y)
-	double x, y;
-{
-	return remainder(x, y);
-}
diff --git a/src/w_dremf.c b/src/w_dremf.c
deleted file mode 100644
index faf3d2a..0000000
--- a/src/w_dremf.c
+++ /dev/null
@@ -1,17 +0,0 @@
-/*
- * dremf() wrapper for remainderf().
- *
- * Written by J.T. Conklin, <jtc@wimsey.com>
- * Placed into the Public Domain, 1994.
- */
-/* $FreeBSD: src/lib/msun/src/w_dremf.c,v 1.3 2004/07/28 05:53:18 kan Exp $ */
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-DLLEXPORT float
-dremf(float x, float y)
-{
-	return remainderf(x, y);
-}

From c08bd962f2f00e09e46989c18fcdd8cbbc2caf6d Mon Sep 17 00:00:00 2001
From: Ed Schouten <ed@nuxi.nl>
Date: Thu, 12 Feb 2015 13:35:19 +0100
Subject: [PATCH 4/9] Remove the non-standard finite() function.

---
 i387/Make.files         |  2 +-
 i387/s_finite.S         | 23 -----------------------
 include/openlibm_math.h |  2 --
 src/Make.files          |  1 -
 src/s_finite.c          | 31 -------------------------------
 src/s_finitef.c         | 34 ----------------------------------
 6 files changed, 1 insertion(+), 92 deletions(-)
 delete mode 100644 i387/s_finite.S
 delete mode 100644 src/s_finite.c
 delete mode 100644 src/s_finitef.c

diff --git a/i387/Make.files b/i387/Make.files
index abfb4fe..369c541 100644
--- a/i387/Make.files
+++ b/i387/Make.files
@@ -1,6 +1,6 @@
 $(CUR_SRCS) = e_exp.S e_fmod.S e_log.S e_log10.S \
 	    e_remainder.S e_sqrt.S s_ceil.S s_copysign.S \
-	    s_finite.S s_floor.S s_llrint.S s_logb.S s_lrint.S \
+	    s_floor.S s_llrint.S s_logb.S s_lrint.S \
 	    s_remquo.S s_rint.S  s_tan.S s_trunc.S
 
 ifneq ($(OS), WINNT)
diff --git a/i387/s_finite.S b/i387/s_finite.S
deleted file mode 100644
index 41015b4..0000000
--- a/i387/s_finite.S
+++ /dev/null
@@ -1,23 +0,0 @@
-/*
- * Written by:
- *	J.T. Conklin (jtc@netbsd.org)
- *	Public domain.
- */
-
-#include <i387/bsd_asm.h>
-//__FBSDID("$FreeBSD: src/lib/msun/i387/s_finite.S,v 1.10 2011/01/07 16:13:12 kib Exp $")
-
-ENTRY(finite)
-	movl	8(%esp),%eax
-	andl	$0x7ff00000, %eax
-	cmpl	$0x7ff00000, %eax
-	setneb	%al
-	andl	$0x000000ff, %eax
-	ret
-END(finite)
-
-	
-/* Enable stack protection */
-#if defined(__linux__) && defined(__ELF__)
-.section .note.GNU-stack,"",%progbits
-#endif
diff --git a/include/openlibm_math.h b/include/openlibm_math.h
index faa30de..f786eb5 100644
--- a/include/openlibm_math.h
+++ b/include/openlibm_math.h
@@ -303,7 +303,6 @@ double	trunc(double);
  * BSD math library entry points
  */
 #if __BSD_VISIBLE
-int	finite(double) __pure2;
 int	isnanf(float) __pure2;
 
 /*
@@ -390,7 +389,6 @@ float	fminf(float, float) __pure2;
  */
 #if __BSD_VISIBLE
 float	dremf(float, float);
-int	finitef(float) __pure2;
 float	j0f(float);
 float	j1f(float);
 float	jnf(int, float);
diff --git a/src/Make.files b/src/Make.files
index 8c8bf07..ef79a4a 100644
--- a/src/Make.files
+++ b/src/Make.files
@@ -15,7 +15,6 @@ $(CUR_SRCS) = common.c \
 	s_copysign.c s_copysignf.c s_cos.c s_cosf.c \
 	s_csqrt.c s_csqrtf.c s_erf.c s_erff.c \
 	s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabs.c s_fabsf.c s_fdim.c \
-	s_finite.c s_finitef.c \
 	s_floor.c s_floorf.c s_fma.c s_fmaf.c \
 	s_fmax.c s_fmaxf.c s_fmin.c \
 	s_fminf.c s_fpclassify.c \
diff --git a/src/s_finite.c b/src/s_finite.c
deleted file mode 100644
index 8118140..0000000
--- a/src/s_finite.c
+++ /dev/null
@@ -1,31 +0,0 @@
-/* @(#)s_finite.c 5.1 93/09/24 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/s_finite.c,v 1.9 2008/02/22 02:30:35 das Exp $");
-
-/*
- * finite(x) returns 1 is x is finite, else 0;
- * no branching!
- */
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-DLLEXPORT int
-finite(double x)
-{
-	int32_t hx;
-	GET_HIGH_WORD(hx,x);
-	return (int)((u_int32_t)((hx&0x7fffffff)-0x7ff00000)>>31);
-}
diff --git a/src/s_finitef.c b/src/s_finitef.c
deleted file mode 100644
index b26b0c7..0000000
--- a/src/s_finitef.c
+++ /dev/null
@@ -1,34 +0,0 @@
-/* s_finitef.c -- float version of s_finite.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/s_finitef.c,v 1.7 2008/02/22 02:30:35 das Exp $");
-
-/*
- * finitef(x) returns 1 is x is finite, else 0;
- * no branching!
- */
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-DLLEXPORT int
-finitef(float x)
-{
-	int32_t ix;
-	GET_FLOAT_WORD(ix,x);
-	return (int)((u_int32_t)((ix&0x7fffffff)-0x7f800000)>>31);
-}

From ce4982acf850517075a6ddb2350bb4c016fa4002 Mon Sep 17 00:00:00 2001
From: Ed Schouten <ed@nuxi.nl>
Date: Thu, 12 Feb 2015 13:38:13 +0100
Subject: [PATCH 5/9] Remove scalb(). The scalbn() or scalbln() function should
 be used.

---
 include/openlibm_math.h |  5 ----
 src/Make.files          |  2 +-
 src/e_scalb.c           | 48 --------------------------------
 src/e_scalbf.c          | 44 -----------------------------
 src/math_private.h      |  2 --
 test/libm-test.c        | 61 -----------------------------------------
 6 files changed, 1 insertion(+), 161 deletions(-)
 delete mode 100644 src/e_scalb.c
 delete mode 100644 src/e_scalbf.c

diff --git a/include/openlibm_math.h b/include/openlibm_math.h
index f786eb5..d97475c 100644
--- a/include/openlibm_math.h
+++ b/include/openlibm_math.h
@@ -280,10 +280,6 @@ double	jn(int, double);
 double	y0(double);
 double	y1(double);
 double	yn(int, double);
-
-#if __XSI_VISIBLE <= 600 || __BSD_VISIBLE
-double	scalb(double, double);
-#endif
 #endif /* __BSD_VISIBLE || __XSI_VISIBLE */
 
 #if __BSD_VISIBLE || __ISO_C_VISIBLE >= 1999
@@ -392,7 +388,6 @@ float	dremf(float, float);
 float	j0f(float);
 float	j1f(float);
 float	jnf(int, float);
-float	scalbf(float, float);
 float	y0f(float);
 float	y1f(float);
 float	ynf(int, float);
diff --git a/src/Make.files b/src/Make.files
index ef79a4a..1489aa8 100644
--- a/src/Make.files
+++ b/src/Make.files
@@ -5,7 +5,7 @@ $(CUR_SRCS) = common.c \
 	e_hypot.c e_hypotf.c e_j0.c e_j0f.c e_j1.c e_j1f.c \
 	e_jn.c e_jnf.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \
 	e_lgammal.c e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c \
-	e_pow.c e_powf.c e_remainder.c e_remainderf.c e_scalb.c e_scalbf.c \
+	e_pow.c e_powf.c e_remainder.c e_remainderf.c \
 	e_rem_pio2.c e_rem_pio2f.c \
 	e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c \
 	k_cos.c k_exp.c k_expf.c k_rem_pio2.c k_sin.c k_tan.c \
diff --git a/src/e_scalb.c b/src/e_scalb.c
deleted file mode 100644
index c9a4208..0000000
--- a/src/e_scalb.c
+++ /dev/null
@@ -1,48 +0,0 @@
-
-/* @(#)e_scalb.c 1.3 95/01/18 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunSoft, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
- * is preserved.
- * ====================================================
- */
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/e_scalb.c,v 1.13 2008/02/22 02:30:35 das Exp $");
-
-/*
- * __ieee754_scalb(x, fn) is provide for
- * passing various standard test suite. One 
- * should use scalbn() instead.
- */
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-#ifdef _SCALB_INT
-DLLEXPORT double
-__ieee754_scalb(double x, int fn)
-#else
-DLLEXPORT double
-__ieee754_scalb(double x, double fn)
-#endif
-{
-#ifdef _SCALB_INT
-	return scalbn(x,fn);
-#else
-	if (isnan(x)||isnan(fn)) return x*fn;
-	if (!isfinite(fn)) {
-	    if(fn>0.0) return x*fn;
-	    else       return x/(-fn);
-	}
-	if (rint(fn)!=fn) return (fn-fn)/(fn-fn);
-	if ( fn > 65000.0) return scalbn(x, 65000);
-	if (-fn > 65000.0) return scalbn(x,-65000);
-	return scalbn(x,(int)fn);
-#endif
-}
diff --git a/src/e_scalbf.c b/src/e_scalbf.c
deleted file mode 100644
index 32258f7..0000000
--- a/src/e_scalbf.c
+++ /dev/null
@@ -1,44 +0,0 @@
-/* e_scalbf.c -- float version of e_scalb.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/e_scalbf.c,v 1.13 2008/02/22 02:30:35 das Exp $");
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-#ifdef _SCALB_INT
-DLLEXPORT float
-__ieee754_scalbf(float x, int fn)
-#else
-DLLEXPORT float
-__ieee754_scalbf(float x, float fn)
-#endif
-{
-#ifdef _SCALB_INT
-	return scalbnf(x,fn);
-#else
-	if (isnan(x)||isnan(fn)) return x*fn;
-	if (!isfinite(fn)) {
-	    if(fn>(float)0.0) return x*fn;
-	    else       return x/(-fn);
-	}
-	if (rintf(fn)!=fn) return (fn-fn)/(fn-fn);
-	if ( fn > (float)65000.0) return scalbnf(x, 65000);
-	if (-fn > (float)65000.0) return scalbnf(x,-65000);
-	return scalbnf(x,(int)fn);
-#endif
-}
diff --git a/src/math_private.h b/src/math_private.h
index 7f96cf7..7bfd76b 100644
--- a/src/math_private.h
+++ b/src/math_private.h
@@ -291,7 +291,6 @@ irint(double x)
 #define	__ieee754_jn	jn
 #define	__ieee754_yn	yn
 #define	__ieee754_remainder remainder
-#define	__ieee754_scalb	scalb
 #define	__ieee754_sqrtf	sqrtf
 #define	__ieee754_acosf	acosf
 #define	__ieee754_acoshf acoshf
@@ -316,7 +315,6 @@ irint(double x)
 #define	__ieee754_jnf	jnf
 #define	__ieee754_ynf	ynf
 #define	__ieee754_remainderf remainderf
-#define	__ieee754_scalbf scalbf
 
 /* fdlibm kernel function */
 int	__kernel_rem_pio2(double*,double*,int,int,int);
diff --git a/test/libm-test.c b/test/libm-test.c
index e6955e9..b5903ce 100644
--- a/test/libm-test.c
+++ b/test/libm-test.c
@@ -3807,66 +3807,6 @@ round_test (void)
 }
 
 
-static void
-scalb_test (void)
-{
-
-  init_max_error ();
-
-  check_float ("scalb (2.0, 0.5) == NaN plus invalid exception",  FUNC(scalb) (2.0, 0.5), nan_value, 0, 0, INVALID_EXCEPTION);
-  check_float ("scalb (3.0, -2.5) == NaN plus invalid exception",  FUNC(scalb) (3.0, -2.5), nan_value, 0, 0, INVALID_EXCEPTION);
-
-  check_float ("scalb (0, NaN) == NaN",  FUNC(scalb) (0, nan_value), nan_value, 0, 0, 0);
-  check_float ("scalb (1, NaN) == NaN",  FUNC(scalb) (1, nan_value), nan_value, 0, 0, 0);
-
-  check_float ("scalb (1, 0) == 1",  FUNC(scalb) (1, 0), 1, 0, 0, 0);
-  check_float ("scalb (-1, 0) == -1",  FUNC(scalb) (-1, 0), -1, 0, 0, 0);
-
-  check_float ("scalb (0, inf) == NaN plus invalid exception",  FUNC(scalb) (0, plus_infty), nan_value, 0, 0, INVALID_EXCEPTION);
-  check_float ("scalb (-0, inf) == NaN plus invalid exception",  FUNC(scalb) (minus_zero, plus_infty), nan_value, 0, 0, INVALID_EXCEPTION);
-
-  check_float ("scalb (0, 2) == 0",  FUNC(scalb) (0, 2), 0, 0, 0, 0);
-  check_float ("scalb (-0, -4) == -0",  FUNC(scalb) (minus_zero, -4), minus_zero, 0, 0, 0);
-  check_float ("scalb (0, 0) == 0",  FUNC(scalb) (0, 0), 0, 0, 0, 0);
-  check_float ("scalb (-0, 0) == -0",  FUNC(scalb) (minus_zero, 0), minus_zero, 0, 0, 0);
-  check_float ("scalb (0, -1) == 0",  FUNC(scalb) (0, -1), 0, 0, 0, 0);
-  check_float ("scalb (-0, -10) == -0",  FUNC(scalb) (minus_zero, -10), minus_zero, 0, 0, 0);
-  check_float ("scalb (0, -inf) == 0",  FUNC(scalb) (0, minus_infty), 0, 0, 0, 0);
-  check_float ("scalb (-0, -inf) == -0",  FUNC(scalb) (minus_zero, minus_infty), minus_zero, 0, 0, 0);
-
-  check_float ("scalb (inf, -1) == inf",  FUNC(scalb) (plus_infty, -1), plus_infty, 0, 0, 0);
-  check_float ("scalb (-inf, -10) == -inf",  FUNC(scalb) (minus_infty, -10), minus_infty, 0, 0, 0);
-  check_float ("scalb (inf, 0) == inf",  FUNC(scalb) (plus_infty, 0), plus_infty, 0, 0, 0);
-  check_float ("scalb (-inf, 0) == -inf",  FUNC(scalb) (minus_infty, 0), minus_infty, 0, 0, 0);
-  check_float ("scalb (inf, 2) == inf",  FUNC(scalb) (plus_infty, 2), plus_infty, 0, 0, 0);
-  check_float ("scalb (-inf, 100) == -inf",  FUNC(scalb) (minus_infty, 100), minus_infty, 0, 0, 0);
-
-  check_float ("scalb (0.1, -inf) == 0.0",  FUNC(scalb) (0.1L, minus_infty), 0.0, 0, 0, 0);
-  check_float ("scalb (-0.1, -inf) == -0",  FUNC(scalb) (-0.1L, minus_infty), minus_zero, 0, 0, 0);
-
-  check_float ("scalb (1, inf) == inf",  FUNC(scalb) (1, plus_infty), plus_infty, 0, 0, 0);
-  check_float ("scalb (-1, inf) == -inf",  FUNC(scalb) (-1, plus_infty), minus_infty, 0, 0, 0);
-  check_float ("scalb (inf, inf) == inf",  FUNC(scalb) (plus_infty, plus_infty), plus_infty, 0, 0, 0);
-  check_float ("scalb (-inf, inf) == -inf",  FUNC(scalb) (minus_infty, plus_infty), minus_infty, 0, 0, 0);
-
-  check_float ("scalb (inf, -inf) == NaN plus invalid exception",  FUNC(scalb) (plus_infty, minus_infty), nan_value, 0, 0, INVALID_EXCEPTION);
-  check_float ("scalb (-inf, -inf) == NaN plus invalid exception",  FUNC(scalb) (minus_infty, minus_infty), nan_value, 0, 0, INVALID_EXCEPTION);
-
-  check_float ("scalb (NaN, 1) == NaN",  FUNC(scalb) (nan_value, 1), nan_value, 0, 0, 0);
-  check_float ("scalb (1, NaN) == NaN",  FUNC(scalb) (1, nan_value), nan_value, 0, 0, 0);
-  check_float ("scalb (NaN, 0) == NaN",  FUNC(scalb) (nan_value, 0), nan_value, 0, 0, 0);
-  check_float ("scalb (0, NaN) == NaN",  FUNC(scalb) (0, nan_value), nan_value, 0, 0, 0);
-  check_float ("scalb (NaN, inf) == NaN",  FUNC(scalb) (nan_value, plus_infty), nan_value, 0, 0, 0);
-  check_float ("scalb (inf, NaN) == NaN",  FUNC(scalb) (plus_infty, nan_value), nan_value, 0, 0, 0);
-  check_float ("scalb (NaN, NaN) == NaN",  FUNC(scalb) (nan_value, nan_value), nan_value, 0, 0, 0);
-
-  check_float ("scalb (0.8, 4) == 12.8",  FUNC(scalb) (0.8L, 4), 12.8L, 0, 0, 0);
-  check_float ("scalb (-0.854375, 5) == -27.34",  FUNC(scalb) (-0.854375L, 5), -27.34L, 0, 0, 0);
-
-  print_max_error ("scalb", 0, 0);
-}
-
-
 static void
 scalbn_test (void)
 {
@@ -4518,7 +4458,6 @@ main (int argc, char **argv)
   logb_test ();
   modf_test ();
   ilogb_test ();
-  scalb_test ();
   scalbn_test ();
   scalbln_test ();
 

From 8c8693cf7962be7b0e53c491287aeb2fb13ae6f2 Mon Sep 17 00:00:00 2001
From: Ed Schouten <ed@nuxi.nl>
Date: Thu, 12 Feb 2015 13:48:13 +0100
Subject: [PATCH 6/9] Remove [jy][01n]f(). X/Open only standardizes the double
 versions.

---
 include/openlibm_math.h |   8 -
 src/Make.files          |   5 +-
 src/e_j0f.c             | 344 ----------------------------------------
 src/e_j1f.c             | 340 ---------------------------------------
 src/e_jnf.c             | 200 -----------------------
 src/math_private.h      |   6 -
 test/libm-test.c        |   6 +
 7 files changed, 8 insertions(+), 901 deletions(-)
 delete mode 100644 src/e_j0f.c
 delete mode 100644 src/e_j1f.c
 delete mode 100644 src/e_jnf.c

diff --git a/include/openlibm_math.h b/include/openlibm_math.h
index d97475c..4d53605 100644
--- a/include/openlibm_math.h
+++ b/include/openlibm_math.h
@@ -384,14 +384,6 @@ float	fminf(float, float) __pure2;
  * float versions of BSD math library entry points
  */
 #if __BSD_VISIBLE
-float	dremf(float, float);
-float	j0f(float);
-float	j1f(float);
-float	jnf(int, float);
-float	y0f(float);
-float	y1f(float);
-float	ynf(int, float);
-
 /*
  * Float versions of reentrant version of lgamma; passes signgam back by
  * reference as the second argument; user must allocate space for signgam.
diff --git a/src/Make.files b/src/Make.files
index 1489aa8..855c4d3 100644
--- a/src/Make.files
+++ b/src/Make.files
@@ -1,9 +1,8 @@
 $(CUR_SRCS) = common.c \
 	e_acos.c e_acosf.c e_acosh.c e_acoshf.c e_asin.c e_asinf.c \
 	e_atan2.c e_atan2f.c e_atanh.c e_atanhf.c e_cosh.c e_coshf.c e_exp.c \
-	e_expf.c e_fmod.c e_fmodf.c \
-	e_hypot.c e_hypotf.c e_j0.c e_j0f.c e_j1.c e_j1f.c \
-	e_jn.c e_jnf.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \
+	e_expf.c e_fmod.c e_fmodf.c e_hypot.c e_hypotf.c e_j0.c e_j1.c \
+	e_jn.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \
 	e_lgammal.c e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c \
 	e_pow.c e_powf.c e_remainder.c e_remainderf.c \
 	e_rem_pio2.c e_rem_pio2f.c \
diff --git a/src/e_j0f.c b/src/e_j0f.c
deleted file mode 100644
index 74d753b..0000000
--- a/src/e_j0f.c
+++ /dev/null
@@ -1,344 +0,0 @@
-/* e_j0f.c -- float version of e_j0.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include <assert.h>
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/e_j0f.c,v 1.8 2008/02/22 02:30:35 das Exp $");
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-static float pzerof(float), qzerof(float);
-
-static const float
-huge 	= 1e30,
-one	= 1.0,
-invsqrtpi=  5.6418961287e-01, /* 0x3f106ebb */
-tpi      =  6.3661974669e-01, /* 0x3f22f983 */
- 		/* R0/S0 on [0, 2.00] */
-R02  =  1.5625000000e-02, /* 0x3c800000 */
-R03  = -1.8997929874e-04, /* 0xb947352e */
-R04  =  1.8295404516e-06, /* 0x35f58e88 */
-R05  = -4.6183270541e-09, /* 0xb19eaf3c */
-S01  =  1.5619102865e-02, /* 0x3c7fe744 */
-S02  =  1.1692678527e-04, /* 0x38f53697 */
-S03  =  5.1354652442e-07, /* 0x3509daa6 */
-S04  =  1.1661400734e-09; /* 0x30a045e8 */
-
-static const float zero = 0.0;
-
-DLLEXPORT float
-__ieee754_j0f(float x)
-{
-	float z, s,c,ss,cc,r,u,v;
-	int32_t hx,ix;
-
-	GET_FLOAT_WORD(hx,x);
-	ix = hx&0x7fffffff;
-	if(ix>=0x7f800000) return one/(x*x);
-	x = fabsf(x);
-	if(ix >= 0x40000000) {	/* |x| >= 2.0 */
-		s = sinf(x);
-		c = cosf(x);
-		ss = s-c;
-		cc = s+c;
-		if(ix<0x7f000000) {  /* make sure x+x not overflow */
-		    z = -cosf(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(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(x);
-		else {
-		    u = pzerof(x); v = qzerof(x);
-		    z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);
-		}
-		return z;
-	}
-	if(ix<0x39000000) {	/* |x| < 2**-13 */
-	    if(huge+x>one) {	/* raise inexact if x != 0 */
-	        if(ix<0x32000000) return one;	/* |x|<2**-27 */
-	        else 	      return one - (float)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(ix < 0x3F800000) {	/* |x| < 1.00 */
-	    return one + z*((float)-0.25+(r/s));
-	} else {
-	    u = (float)0.5*x;
-	    return((one+u)*(one-u)+z*(r/s));
-	}
-}
-
-static const float
-u00  = -7.3804296553e-02, /* 0xbd9726b5 */
-u01  =  1.7666645348e-01, /* 0x3e34e80d */
-u02  = -1.3818567619e-02, /* 0xbc626746 */
-u03  =  3.4745343146e-04, /* 0x39b62a69 */
-u04  = -3.8140706238e-06, /* 0xb67ff53c */
-u05  =  1.9559013964e-08, /* 0x32a802ba */
-u06  = -3.9820518410e-11, /* 0xae2f21eb */
-v01  =  1.2730483897e-02, /* 0x3c509385 */
-v02  =  7.6006865129e-05, /* 0x389f65e0 */
-v03  =  2.5915085189e-07, /* 0x348b216c */
-v04  =  4.4111031494e-10; /* 0x2ff280c2 */
-
-DLLEXPORT float
-__ieee754_y0f(float x)
-{
-	float z, s,c,ss,cc,u,v;
-	int32_t hx,ix;
-
-	GET_FLOAT_WORD(hx,x);
-        ix = 0x7fffffff&hx;
-    /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
-	if(ix>=0x7f800000) return  one/(x+x*x);
-        if(ix==0) return -one/zero;
-        if(hx<0) return zero/zero;
-        if(ix >= 0x40000000) {  /* |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 = sinf(x);
-                c = cosf(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(ix<0x7f000000) {  /* make sure x+x not overflow */
-                    z = -cosf(x+x);
-                    if ((s*c)<zero) cc = z/ss;
-                    else            ss = z/cc;
-                }
-                if(ix>0x80000000) z = (invsqrtpi*ss)/sqrtf(x);
-                else {
-                    u = pzerof(x); v = qzerof(x);
-                    z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
-                }
-                return z;
-	}
-	if(ix<=0x32000000) {	/* x < 2**-27 */
-	    return(u00 + tpi*__ieee754_logf(x));
-	}
-	z = x*x;
-	u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
-	v = one+z*(v01+z*(v02+z*(v03+z*v04)));
-	return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x)));
-}
-
-/* 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 const float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-  0.0000000000e+00, /* 0x00000000 */
- -7.0312500000e-02, /* 0xbd900000 */
- -8.0816707611e+00, /* 0xc1014e86 */
- -2.5706311035e+02, /* 0xc3808814 */
- -2.4852163086e+03, /* 0xc51b5376 */
- -5.2530439453e+03, /* 0xc5a4285a */
-};
-static const float pS8[5] = {
-  1.1653436279e+02, /* 0x42e91198 */
-  3.8337448730e+03, /* 0x456f9beb */
-  4.0597855469e+04, /* 0x471e95db */
-  1.1675296875e+05, /* 0x47e4087c */
-  4.7627726562e+04, /* 0x473a0bba */
-};
-static const float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
- -1.1412546255e-11, /* 0xad48c58a */
- -7.0312492549e-02, /* 0xbd8fffff */
- -4.1596107483e+00, /* 0xc0851b88 */
- -6.7674766541e+01, /* 0xc287597b */
- -3.3123129272e+02, /* 0xc3a59d9b */
- -3.4643338013e+02, /* 0xc3ad3779 */
-};
-static const float pS5[5] = {
-  6.0753936768e+01, /* 0x42730408 */
-  1.0512523193e+03, /* 0x44836813 */
-  5.9789707031e+03, /* 0x45bad7c4 */
-  9.6254453125e+03, /* 0x461665c8 */
-  2.4060581055e+03, /* 0x451660ee */
-};
-
-static const float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
- -2.5470459075e-09, /* 0xb12f081b */
- -7.0311963558e-02, /* 0xbd8fffb8 */
- -2.4090321064e+00, /* 0xc01a2d95 */
- -2.1965976715e+01, /* 0xc1afba52 */
- -5.8079170227e+01, /* 0xc2685112 */
- -3.1447946548e+01, /* 0xc1fb9565 */
-};
-static const float pS3[5] = {
-  3.5856033325e+01, /* 0x420f6c94 */
-  3.6151397705e+02, /* 0x43b4c1ca */
-  1.1936077881e+03, /* 0x44953373 */
-  1.1279968262e+03, /* 0x448cffe6 */
-  1.7358093262e+02, /* 0x432d94b8 */
-};
-
-static const float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
- -8.8753431271e-08, /* 0xb3be98b7 */
- -7.0303097367e-02, /* 0xbd8ffb12 */
- -1.4507384300e+00, /* 0xbfb9b1cc */
- -7.6356959343e+00, /* 0xc0f4579f */
- -1.1193166733e+01, /* 0xc1331736 */
- -3.2336456776e+00, /* 0xc04ef40d */
-};
-static const float pS2[5] = {
-  2.2220300674e+01, /* 0x41b1c32d */
-  1.3620678711e+02, /* 0x430834f0 */
-  2.7047027588e+02, /* 0x43873c32 */
-  1.5387539673e+02, /* 0x4319e01a */
-  1.4657617569e+01, /* 0x416a859a */
-};
-
-	/* Note: This function is only called for ix>=0x40000000 (see above) */
-	static float pzerof(float x)
-{
-	const float *p,*q;
-	float z,r,s;
-	int32_t ix;
-	GET_FLOAT_WORD(ix,x);
-	ix &= 0x7fffffff;
-        assert(ix>=0x40000000 && ix<=0x48000000);
-	if(ix>=0x41000000)     {p = pR8; q= pS8;}
-	else if(ix>=0x40f71c58){p = pR5; q= pS5;}
-	else if(ix>=0x4036db68){p = pR3; q= pS3;}
-	else                   {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 const float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-  0.0000000000e+00, /* 0x00000000 */
-  7.3242187500e-02, /* 0x3d960000 */
-  1.1768206596e+01, /* 0x413c4a93 */
-  5.5767340088e+02, /* 0x440b6b19 */
-  8.8591972656e+03, /* 0x460a6cca */
-  3.7014625000e+04, /* 0x471096a0 */
-};
-static const float qS8[6] = {
-  1.6377603149e+02, /* 0x4323c6aa */
-  8.0983447266e+03, /* 0x45fd12c2 */
-  1.4253829688e+05, /* 0x480b3293 */
-  8.0330925000e+05, /* 0x49441ed4 */
-  8.4050156250e+05, /* 0x494d3359 */
- -3.4389928125e+05, /* 0xc8a7eb69 */
-};
-
-static const float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-  1.8408595828e-11, /* 0x2da1ec79 */
-  7.3242180049e-02, /* 0x3d95ffff */
-  5.8356351852e+00, /* 0x40babd86 */
-  1.3511157227e+02, /* 0x43071c90 */
-  1.0272437744e+03, /* 0x448067cd */
-  1.9899779053e+03, /* 0x44f8bf4b */
-};
-static const float qS5[6] = {
-  8.2776611328e+01, /* 0x42a58da0 */
-  2.0778142090e+03, /* 0x4501dd07 */
-  1.8847289062e+04, /* 0x46933e94 */
-  5.6751113281e+04, /* 0x475daf1d */
-  3.5976753906e+04, /* 0x470c88c1 */
- -5.3543427734e+03, /* 0xc5a752be */
-};
-
-static const float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
-  4.3774099900e-09, /* 0x3196681b */
-  7.3241114616e-02, /* 0x3d95ff70 */
-  3.3442313671e+00, /* 0x405607e3 */
-  4.2621845245e+01, /* 0x422a7cc5 */
-  1.7080809021e+02, /* 0x432acedf */
-  1.6673394775e+02, /* 0x4326bbe4 */
-};
-static const float qS3[6] = {
-  4.8758872986e+01, /* 0x42430916 */
-  7.0968920898e+02, /* 0x44316c1c */
-  3.7041481934e+03, /* 0x4567825f */
-  6.4604252930e+03, /* 0x45c9e367 */
-  2.5163337402e+03, /* 0x451d4557 */
- -1.4924745178e+02, /* 0xc3153f59 */
-};
-
-static const float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-  1.5044444979e-07, /* 0x342189db */
-  7.3223426938e-02, /* 0x3d95f62a */
-  1.9981917143e+00, /* 0x3fffc4bf */
-  1.4495602608e+01, /* 0x4167edfd */
-  3.1666231155e+01, /* 0x41fd5471 */
-  1.6252708435e+01, /* 0x4182058c */
-};
-static const float qS2[6] = {
-  3.0365585327e+01, /* 0x41f2ecb8 */
-  2.6934811401e+02, /* 0x4386ac8f */
-  8.4478375244e+02, /* 0x44533229 */
-  8.8293585205e+02, /* 0x445cbbe5 */
-  2.1266638184e+02, /* 0x4354aa98 */
- -5.3109550476e+00, /* 0xc0a9f358 */
-};
-
-	/* Note: This function is only called for ix>=0x40000000 (see above) */
-	static float qzerof(float x)
-{
-	const float *p,*q;
-	float s,r,z;
-	int32_t ix;
-	GET_FLOAT_WORD(ix,x);
-	ix &= 0x7fffffff;
-        assert(ix>=0x40000000 && ix<=0x48000000);
-	if(ix>=0x41000000)     {p = qR8; q= qS8;}
-	else if(ix>=0x40f71c58){p = qR5; q= qS5;}
-	else if(ix>=0x4036db68){p = qR3; q= qS3;}
-	else                   {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 (-(float).125 + r/s)/x;
-}
diff --git a/src/e_j1f.c b/src/e_j1f.c
deleted file mode 100644
index afed84f..0000000
--- a/src/e_j1f.c
+++ /dev/null
@@ -1,340 +0,0 @@
-/* e_j1f.c -- float version of e_j1.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include <assert.h>
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/e_j1f.c,v 1.8 2008/02/22 02:30:35 das Exp $");
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-static float ponef(float), qonef(float);
-
-static const float
-huge    = 1e30,
-one	= 1.0,
-invsqrtpi=  5.6418961287e-01, /* 0x3f106ebb */
-tpi      =  6.3661974669e-01, /* 0x3f22f983 */
-	/* R0/S0 on [0,2] */
-r00  = -6.2500000000e-02, /* 0xbd800000 */
-r01  =  1.4070566976e-03, /* 0x3ab86cfd */
-r02  = -1.5995563444e-05, /* 0xb7862e36 */
-r03  =  4.9672799207e-08, /* 0x335557d2 */
-s01  =  1.9153760746e-02, /* 0x3c9ce859 */
-s02  =  1.8594678841e-04, /* 0x3942fab6 */
-s03  =  1.1771846857e-06, /* 0x359dffc2 */
-s04  =  5.0463624390e-09, /* 0x31ad6446 */
-s05  =  1.2354227016e-11; /* 0x2d59567e */
-
-static const float zero    = 0.0;
-
-DLLEXPORT float
-__ieee754_j1f(float x)
-{
-	float z, s,c,ss,cc,r,u,v,y;
-	int32_t hx,ix;
-
-	GET_FLOAT_WORD(hx,x);
-	ix = hx&0x7fffffff;
-	if(ix>=0x7f800000) return one/x;
-	y = fabsf(x);
-	if(ix >= 0x40000000) {	/* |x| >= 2.0 */
-		s = sinf(y);
-		c = cosf(y);
-		ss = -s-c;
-		cc = s-c;
-		if(ix<0x7f000000) {  /* make sure y+y not overflow */
-		    z = cosf(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(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(y);
-		else {
-		    u = ponef(y); v = qonef(y);
-		    z = invsqrtpi*(u*cc-v*ss)/sqrtf(y);
-		}
-		if(hx<0) return -z;
-		else  	 return  z;
-	}
-	if(ix<0x32000000) {	/* |x|<2**-27 */
-	    if(huge+x>one) return (float)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*(float)0.5+r/s);
-}
-
-static const float U0[5] = {
- -1.9605709612e-01, /* 0xbe48c331 */
-  5.0443872809e-02, /* 0x3d4e9e3c */
- -1.9125689287e-03, /* 0xbafaaf2a */
-  2.3525259166e-05, /* 0x37c5581c */
- -9.1909917899e-08, /* 0xb3c56003 */
-};
-static const float V0[5] = {
-  1.9916731864e-02, /* 0x3ca3286a */
-  2.0255257550e-04, /* 0x3954644b */
-  1.3560879779e-06, /* 0x35b602d4 */
-  6.2274145840e-09, /* 0x31d5f8eb */
-  1.6655924903e-11, /* 0x2d9281cf */
-};
-
-DLLEXPORT float
-__ieee754_y1f(float x)
-{
-	float z, s,c,ss,cc,u,v;
-	int32_t hx,ix;
-
-	GET_FLOAT_WORD(hx,x);
-        ix = 0x7fffffff&hx;
-    /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
-	if(ix>=0x7f800000) return  one/(x+x*x);
-        if(ix==0) return -one/zero;
-        if(hx<0) return zero/zero;
-        if(ix >= 0x40000000) {  /* |x| >= 2.0 */
-                s = sinf(x);
-                c = cosf(x);
-                ss = -s-c;
-                cc = s-c;
-                if(ix<0x7f000000) {  /* make sure x+x not overflow */
-                    z = cosf(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(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x);
-                else {
-                    u = ponef(x); v = qonef(x);
-                    z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
-                }
-                return z;
-        }
-        if(ix<=0x24800000) {    /* 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*(__ieee754_j1f(x)*__ieee754_logf(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 const float pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-  0.0000000000e+00, /* 0x00000000 */
-  1.1718750000e-01, /* 0x3df00000 */
-  1.3239480972e+01, /* 0x4153d4ea */
-  4.1205184937e+02, /* 0x43ce06a3 */
-  3.8747453613e+03, /* 0x45722bed */
-  7.9144794922e+03, /* 0x45f753d6 */
-};
-static const float ps8[5] = {
-  1.1420736694e+02, /* 0x42e46a2c */
-  3.6509309082e+03, /* 0x45642ee5 */
-  3.6956207031e+04, /* 0x47105c35 */
-  9.7602796875e+04, /* 0x47bea166 */
-  3.0804271484e+04, /* 0x46f0a88b */
-};
-
-static const float pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
-  1.3199052094e-11, /* 0x2d68333f */
-  1.1718749255e-01, /* 0x3defffff */
-  6.8027510643e+00, /* 0x40d9b023 */
-  1.0830818176e+02, /* 0x42d89dca */
-  5.1763616943e+02, /* 0x440168b7 */
-  5.2871520996e+02, /* 0x44042dc6 */
-};
-static const float ps5[5] = {
-  5.9280597687e+01, /* 0x426d1f55 */
-  9.9140142822e+02, /* 0x4477d9b1 */
-  5.3532670898e+03, /* 0x45a74a23 */
-  7.8446904297e+03, /* 0x45f52586 */
-  1.5040468750e+03, /* 0x44bc0180 */
-};
-
-static const float pr3[6] = {
-  3.0250391081e-09, /* 0x314fe10d */
-  1.1718686670e-01, /* 0x3defffab */
-  3.9329774380e+00, /* 0x407bb5e7 */
-  3.5119403839e+01, /* 0x420c7a45 */
-  9.1055007935e+01, /* 0x42b61c2a */
-  4.8559066772e+01, /* 0x42423c7c */
-};
-static const float ps3[5] = {
-  3.4791309357e+01, /* 0x420b2a4d */
-  3.3676245117e+02, /* 0x43a86198 */
-  1.0468714600e+03, /* 0x4482dbe3 */
-  8.9081134033e+02, /* 0x445eb3ed */
-  1.0378793335e+02, /* 0x42cf936c */
-};
-
-static const float pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
-  1.0771083225e-07, /* 0x33e74ea8 */
-  1.1717621982e-01, /* 0x3deffa16 */
-  2.3685150146e+00, /* 0x401795c0 */
-  1.2242610931e+01, /* 0x4143e1bc */
-  1.7693971634e+01, /* 0x418d8d41 */
-  5.0735230446e+00, /* 0x40a25a4d */
-};
-static const float ps2[5] = {
-  2.1436485291e+01, /* 0x41ab7dec */
-  1.2529022980e+02, /* 0x42fa9499 */
-  2.3227647400e+02, /* 0x436846c7 */
-  1.1767937469e+02, /* 0x42eb5bd7 */
-  8.3646392822e+00, /* 0x4105d590 */
-};
-
-	/* Note: This function is only called for ix>=0x40000000 (see above) */
-	static float ponef(float x)
-{
-	const float *p,*q;
-	float z,r,s;
-        int32_t ix;
-	GET_FLOAT_WORD(ix,x);
-	ix &= 0x7fffffff;
-        assert(ix>=0x40000000 && ix<=0x48000000);
-        if(ix>=0x41000000)     {p = pr8; q= ps8;}
-        else if(ix>=0x40f71c58){p = pr5; q= ps5;}
-        else if(ix>=0x4036db68){p = pr3; q= ps3;}
-        else                   {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 const float qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
-  0.0000000000e+00, /* 0x00000000 */
- -1.0253906250e-01, /* 0xbdd20000 */
- -1.6271753311e+01, /* 0xc1822c8d */
- -7.5960174561e+02, /* 0xc43de683 */
- -1.1849806641e+04, /* 0xc639273a */
- -4.8438511719e+04, /* 0xc73d3683 */
-};
-static const float qs8[6] = {
-  1.6139537048e+02, /* 0x43216537 */
-  7.8253862305e+03, /* 0x45f48b17 */
-  1.3387534375e+05, /* 0x4802bcd6 */
-  7.1965775000e+05, /* 0x492fb29c */
-  6.6660125000e+05, /* 0x4922be94 */
- -2.9449025000e+05, /* 0xc88fcb48 */
-};
-
-static const float qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
- -2.0897993405e-11, /* 0xadb7d219 */
- -1.0253904760e-01, /* 0xbdd1fffe */
- -8.0564479828e+00, /* 0xc100e736 */
- -1.8366960144e+02, /* 0xc337ab6b */
- -1.3731937256e+03, /* 0xc4aba633 */
- -2.6124443359e+03, /* 0xc523471c */
-};
-static const float qs5[6] = {
-  8.1276550293e+01, /* 0x42a28d98 */
-  1.9917987061e+03, /* 0x44f8f98f */
-  1.7468484375e+04, /* 0x468878f8 */
-  4.9851425781e+04, /* 0x4742bb6d */
-  2.7948074219e+04, /* 0x46da5826 */
- -4.7191835938e+03, /* 0xc5937978 */
-};
-
-static const float qr3[6] = {
- -5.0783124372e-09, /* 0xb1ae7d4f */
- -1.0253783315e-01, /* 0xbdd1ff5b */
- -4.6101160049e+00, /* 0xc0938612 */
- -5.7847221375e+01, /* 0xc267638e */
- -2.2824453735e+02, /* 0xc3643e9a */
- -2.1921012878e+02, /* 0xc35b35cb */
-};
-static const float qs3[6] = {
-  4.7665153503e+01, /* 0x423ea91e */
-  6.7386511230e+02, /* 0x4428775e */
-  3.3801528320e+03, /* 0x45534272 */
-  5.5477290039e+03, /* 0x45ad5dd5 */
-  1.9031191406e+03, /* 0x44ede3d0 */
- -1.3520118713e+02, /* 0xc3073381 */
-};
-
-static const float qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
- -1.7838172539e-07, /* 0xb43f8932 */
- -1.0251704603e-01, /* 0xbdd1f475 */
- -2.7522056103e+00, /* 0xc0302423 */
- -1.9663616180e+01, /* 0xc19d4f16 */
- -4.2325313568e+01, /* 0xc2294d1f */
- -2.1371921539e+01, /* 0xc1aaf9b2 */
-};
-static const float qs2[6] = {
-  2.9533363342e+01, /* 0x41ec4454 */
-  2.5298155212e+02, /* 0x437cfb47 */
-  7.5750280762e+02, /* 0x443d602e */
-  7.3939318848e+02, /* 0x4438d92a */
-  1.5594900513e+02, /* 0x431bf2f2 */
- -4.9594988823e+00, /* 0xc09eb437 */
-};
-
-	/* Note: This function is only called for ix>=0x40000000 (see above) */
-	static float qonef(float x)
-{
-	const float *p,*q;
-	float  s,r,z;
-	int32_t ix;
-	GET_FLOAT_WORD(ix,x);
-	ix &= 0x7fffffff;
-        assert(ix>=0x40000000 && ix<=0x48000000);
-	if(ix>=0x40200000)     {p = qr8; q= qs8;}
-	else if(ix>=0x40f71c58){p = qr5; q= qs5;}
-	else if(ix>=0x4036db68){p = qr3; q= qs3;}
-	else                   {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 ((float).375 + r/s)/x;
-}
diff --git a/src/e_jnf.c b/src/e_jnf.c
deleted file mode 100644
index a2b83ee..0000000
--- a/src/e_jnf.c
+++ /dev/null
@@ -1,200 +0,0 @@
-/* e_jnf.c -- float version of e_jn.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include "cdefs-compat.h"
-//__FBSDID("$FreeBSD: src/lib/msun/src/e_jnf.c,v 1.11 2010/11/13 10:54:10 uqs Exp $");
-
-#include <openlibm_math.h>
-
-#include "math_private.h"
-
-static const float
-two   =  2.0000000000e+00, /* 0x40000000 */
-one   =  1.0000000000e+00; /* 0x3F800000 */
-
-static const float zero  =  0.0000000000e+00;
-
-DLLEXPORT float
-__ieee754_jnf(int n, float x)
-{
-	int32_t i,hx,ix, sgn;
-	float a, b, temp, di;
-	float 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)
-     */
-	GET_FLOAT_WORD(hx,x);
-	ix = 0x7fffffff&hx;
-    /* if J(n,NaN) is NaN */
-	if(ix>0x7f800000) return x+x;
-	if(n<0){
-		n = -n;
-		x = -x;
-		hx ^= 0x80000000;
-	}
-	if(n==0) return(__ieee754_j0f(x));
-	if(n==1) return(__ieee754_j1f(x));
-	sgn = (n&1)&(hx>>31);	/* even n -- 0, odd n -- sign(x) */
-	x = fabsf(x);
-	if(ix==0||ix>=0x7f800000) 	/* if x is 0 or inf */
-	    b = zero;
-	else if((float)n<=x) {
-		/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
-	    a = __ieee754_j0f(x);
-	    b = __ieee754_j1f(x);
-	    for(i=1;i<n;i++){
-		temp = b;
-		b = b*((float)(i+i)/x) - a; /* avoid underflow */
-		a = temp;
-	    }
-	} else {
-	    if(ix<0x30800000) {	/* 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*(float)0.5; b = temp;
-		    for (a=one,i=2;i<=n;i++) {
-			a *= (float)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 */
-		float t,v;
-		float q0,q1,h,tmp; int32_t k,m;
-		w  = (n+n)/(float)x; h = (float)2.0/(float)x;
-		q0 = w;  z = w+h; q1 = w*z - (float)1.0; k=1;
-		while(q1<(float)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 is
-		 *  likely underflow to zero
-		 */
-		tmp = n;
-		v = two/x;
-		tmp = tmp*__ieee754_logf(fabsf(v*tmp));
-		if(tmp<(float)8.8721679688e+01) {
-	    	    for(i=n-1,di=(float)(i+i);i>0;i--){
-		        temp = b;
-			b *= di;
-			b  = b/x - a;
-		        a = temp;
-			di -= two;
-	     	    }
-		} else {
-	    	    for(i=n-1,di=(float)(i+i);i>0;i--){
-		        temp = b;
-			b *= di;
-			b  = b/x - a;
-		        a = temp;
-			di -= two;
-		    /* scale b to avoid spurious overflow */
-			if(b>(float)1e10) {
-			    a /= b;
-			    t /= b;
-			    b  = one;
-			}
-	     	    }
-		}
-		z = __ieee754_j0f(x);
-		w = __ieee754_j1f(x);
-		if (fabsf(z) >= fabsf(w))
-		    b = (t*z/b);
-		else
-		    b = (t*w/a);
-	    }
-	}
-	if(sgn==1) return -b; else return b;
-}
-
-DLLEXPORT float
-__ieee754_ynf(int n, float x)
-{
-	int32_t i,hx,ix,ib;
-	int32_t sign;
-	float a, b, temp;
-
-	GET_FLOAT_WORD(hx,x);
-	ix = 0x7fffffff&hx;
-    /* if Y(n,NaN) is NaN */
-	if(ix>0x7f800000) return x+x;
-	if(ix==0) return -one/zero;
-	if(hx<0) return zero/zero;
-	sign = 1;
-	if(n<0){
-		n = -n;
-		sign = 1 - ((n&1)<<1);
-	}
-	if(n==0) return(__ieee754_y0f(x));
-	if(n==1) return(sign*__ieee754_y1f(x));
-	if(ix==0x7f800000) return zero;
-
-	a = __ieee754_y0f(x);
-	b = __ieee754_y1f(x);
-	/* quit if b is -inf */
-	GET_FLOAT_WORD(ib,b);
-	for(i=1;i<n&&ib!=0xff800000;i++){
-	    temp = b;
-	    b = ((float)(i+i)/x)*b - a;
-	    GET_FLOAT_WORD(ib,b);
-	    a = temp;
-	}
-	if(sign>0) return b; else return -b;
-}
diff --git a/src/math_private.h b/src/math_private.h
index 7bfd76b..50745f8 100644
--- a/src/math_private.h
+++ b/src/math_private.h
@@ -308,12 +308,6 @@ irint(double x)
 #define	__ieee754_log2f log2f
 #define	__ieee754_sinhf	sinhf
 #define	__ieee754_hypotf hypotf
-#define	__ieee754_j0f	j0f
-#define	__ieee754_j1f	j1f
-#define	__ieee754_y0f	y0f
-#define	__ieee754_y1f	y1f
-#define	__ieee754_jnf	jnf
-#define	__ieee754_ynf	ynf
 #define	__ieee754_remainderf remainderf
 
 /* fdlibm kernel function */
diff --git a/test/libm-test.c b/test/libm-test.c
index b5903ce..99f4718 100644
--- a/test/libm-test.c
+++ b/test/libm-test.c
@@ -2918,6 +2918,7 @@ isnormal_test (void)
   print_max_error ("isnormal", 0, 0);
 }
 
+#ifdef TEST_DOUBLE
 static void
 j0_test (void)
 {
@@ -3055,6 +3056,7 @@ jn_test (void)
 
   print_max_error ("jn", DELTAjn, 0);
 }
+#endif
 
 
 static void
@@ -4121,6 +4123,7 @@ trunc_test (void)
   print_max_error ("trunc", 0, 0);
 }
 
+#ifdef TEST_DOUBLE
 static void
 y0_test (void)
 {
@@ -4256,6 +4259,7 @@ yn_test (void)
   print_max_error ("yn", DELTAyn, 0);
 
 }
+#endif
 
 
 
@@ -4534,6 +4538,7 @@ main (int argc, char **argv)
   ctanh_test ();
 #endif
 
+#ifdef TEST_DOUBLE
   /* Bessel functions:  */
   j0_test ();
   j1_test ();
@@ -4541,6 +4546,7 @@ main (int argc, char **argv)
   y0_test ();
   y1_test ();
   yn_test ();
+#endif
 
   if (output_ulps)
     fclose (ulps_file);

From 16902a6ce1bb60b37894ab01a93dbb3e89f349b6 Mon Sep 17 00:00:00 2001
From: Ed Schouten <ed@nuxi.nl>
Date: Thu, 12 Feb 2015 13:49:28 +0100
Subject: [PATCH 7/9] Reconnect cbrtl() to the build. It seems to build
 properly.

---
 src/Make.files | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/Make.files b/src/Make.files
index 855c4d3..9109359 100644
--- a/src/Make.files
+++ b/src/Make.files
@@ -52,8 +52,7 @@ $(CUR_SRCS) +=	e_acosl.c e_asinl.c e_atan2l.c e_fmodl.c \
 	s_casinl.c s_ctanl.c \
 	s_cimagl.c s_conjl.c s_creall.c s_cacoshl.c s_catanhl.c s_casinhl.c \
 	s_catanl.c s_csinl.c s_cacosl.c s_cexpl.c s_csinhl.c s_ccoshl.c \
-	s_clogl.c s_ctanhl.c s_ccosl.c
-#	s_cbrtl.c
+	s_clogl.c s_ctanhl.c s_ccosl.c s_cbrtl.c
 endif
 
 # C99 complex functions

From b6ff8bbe910ea4cc0d26650d65d60ecc2184255f Mon Sep 17 00:00:00 2001
From: Ed Schouten <ed@nuxi.nl>
Date: Thu, 12 Feb 2015 13:58:35 +0100
Subject: [PATCH 8/9] Clean up a change that I sneaked in accidentically.

I often build the code with -Wmissing-prototypes to ensure that we don't
accidentically pollute the symbol namespace. If we want to provide a
symbol such as isopenlibm(), make sure we also declare it in
<openlibm_math.h>.
---
 include/openlibm_math.h | 3 +++
 src/common.c            | 4 ++--
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/include/openlibm_math.h b/include/openlibm_math.h
index 4d53605..81d5843 100644
--- a/include/openlibm_math.h
+++ b/include/openlibm_math.h
@@ -194,6 +194,9 @@ extern int signgam;
 #if defined(__cplusplus)
 extern "C" {
 #endif
+/* Symbol present when OpenLibm is used. */
+int isopenlibm(void);
+
 /*
  * ANSI/POSIX
  */
diff --git a/src/common.c b/src/common.c
index 5905c20..373f579 100644
--- a/src/common.c
+++ b/src/common.c
@@ -1,6 +1,6 @@
-#include "math_private.h"
+#include <openlibm_math.h>
 
-DLLEXPORT int isopenlibm(void);
+#include "math_private.h"
 
 DLLEXPORT int isopenlibm(void) {
     return 1;

From d64ea4e5f97fac82e8d05dd4a47e711a8cf8ca2e Mon Sep 17 00:00:00 2001
From: Ed Schouten <ed@nuxi.nl>
Date: Thu, 12 Feb 2015 14:15:15 +0100
Subject: [PATCH 9/9] Make small changes to libm-test.c to improve portability.

- Don't define llabs(). This breaks if llabs() is already a macro, which
  is allowed by the C standard/POSIX. llabs() was introduced in C99, so
  I think we can safely assume it is present on all interesting systems.

- Cast the parameters to fabs() to the floating point type. Clang has
  introduced some interesting warnings that trigger if the arguments to
  fabs*() are not the right type.
---
 test/libm-test.c | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/test/libm-test.c b/test/libm-test.c
index 99f4718..604f13b 100644
--- a/test/libm-test.c
+++ b/test/libm-test.c
@@ -712,7 +712,6 @@ check_longlong (const char *test_name, long long int computed,
 
   test_exceptions (test_name, exceptions);
   noTests++;
-#define llabs(x) (x < 0 ? -x : x)
   if (llabs (diff) <= max_ulp)
     ok = 1;
 
@@ -2531,15 +2530,15 @@ fabs_test (void)
 {
   init_max_error ();
 
-  check_float ("fabs (0) == 0",  FUNC(fabs) (0), 0, 0, 0, 0);
+  check_float ("fabs (0) == 0",  FUNC(fabs) ((FLOAT)0.0), 0, 0, 0, 0);
   check_float ("fabs (-0) == 0",  FUNC(fabs) (minus_zero), 0, 0, 0, 0);
 
   check_float ("fabs (inf) == inf",  FUNC(fabs) (plus_infty), plus_infty, 0, 0, 0);
   check_float ("fabs (-inf) == inf",  FUNC(fabs) (minus_infty), plus_infty, 0, 0, 0);
   check_float ("fabs (NaN) == NaN",  FUNC(fabs) (nan_value), nan_value, 0, 0, 0);
 
-  check_float ("fabs (38.0) == 38.0",  FUNC(fabs) (38.0), 38.0, 0, 0, 0);
-  check_float ("fabs (-e) == e",  FUNC(fabs) (-M_El), M_El, 0, 0, 0);
+  check_float ("fabs (38.0) == 38.0",  FUNC(fabs) ((FLOAT)38.0), 38.0, 0, 0, 0);
+  check_float ("fabs (-e) == e",  FUNC(fabs) ((FLOAT)-M_El), M_El, 0, 0, 0);
 
   print_max_error ("fabs", 0, 0);
 }