Replace the ldexp() implementation in libc with a renamed copy of the

scalbn() implementation from libm.  (The two functions are defined to
be identical, but ldexp() lives in libc for backwards compatibility.)
The old ldexp() implementation...
  - was more complicated than this one
  - set errno instead of raising FP exceptions
  - got some corner cases wrong
    (e.g. ldexp(1.0, 2000) in round-to-zero mode)

The new implementation lives in libc/gen instead of
libc/$MACHINE_ARCH/gen, since we don't need N copies of a
machine-independent file.  The amd64 and i386 platforms
retain their fast and correct MD implementations and
override this one.

6 files changed, 123 insertions, 739 deletions

diff --git a/lib/libc/alpha/gen/ldexp.c b/lib/libc/alpha/gen/ldexp.c
deleted file mode 100644
index 78decbb..0000000
--- a/lib/libc/alpha/gen/ldexp.c
+++ /dev/null
@@ -1,137 +0,0 @@
-/*
- * Copyright (c) 1994, 1995 Carnegie-Mellon University.
- * All rights reserved.
- *
- * Author: Chris G. Demetriou
- * 
- * Permission to use, copy, modify and distribute this software and
- * its documentation is hereby granted, provided that both the copyright
- * notice and this permission notice appear in all copies of the
- * software, derivative works or modified versions, and any portions
- * thereof, and that both notices appear in supporting documentation.
- * 
- * CARNEGIE MELLON ALLOWS FREE USE OF THIS SOFTWARE IN ITS "AS IS" 
- * CONDITION.  CARNEGIE MELLON DISCLAIMS ANY LIABILITY OF ANY KIND 
- * FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE USE OF THIS SOFTWARE.
- * 
- * Carnegie Mellon requests users of this software to return to
- *
- *  Software Distribution Coordinator  or  Software.Distribution@CS.CMU.EDU
- *  School of Computer Science
- *  Carnegie Mellon University
- *  Pittsburgh PA 15213-3890
- *
- * any improvements or extensions that they make and grant Carnegie the
- * rights to redistribute these changes.
- *
- *	$NetBSD: ldexp.c,v 1.1 1995/02/10 17:50:24 cgd Exp $
- */
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
-
-#include <sys/types.h>
-#include <machine/ieee.h>
-#include <errno.h>
-#include <math.h>
-
-/*
- * double ldexp(double val, int exp)
- * returns: val * (2**exp)
- */
-double
-ldexp(val, exp)
-	double val;
-	int exp;
-{
-	int oldexp, newexp, mulexp;
-	union doub {
-		double v;
-		struct ieee_double s;
-	} u, mul;
-
-	/*
-	 * If input is zero, or no change, just return input.
-	 * Likewise, if input is Inf or NaN, just return it.
-	 */
-	u.v = val;
-	oldexp = u.s.dbl_exp;
-	if (val == 0 || exp == 0 || oldexp == DBL_EXP_INFNAN)
-		return (val);
-
-	/*
-	 * Compute new exponent and check for over/under flow.
-	 * Underflow, unfortunately, could mean switching to denormal.
-	 * If result out of range, set ERANGE and return 0 if too small
-	 * or Inf if too big, with the same sign as the input value.
-	 */
-	newexp = oldexp + exp;
-	if (newexp >= DBL_EXP_INFNAN) {
-		/* u.s.dbl_sign = val < 0; -- already set */
-		u.s.dbl_exp = DBL_EXP_INFNAN;
-		u.s.dbl_frach = u.s.dbl_fracl = 0;
-		errno = ERANGE;
-		return (u.v);		/* Inf */
-	}
-	if (newexp <= 0) {
-		/*
-		 * The output number is either a denormal or underflows
-		 * (see comments in machine/ieee.h).
-		 */
-		if (newexp <= -DBL_FRACBITS) {
-			/* u.s.dbl_sign = val < 0; -- already set */
-			u.s.dbl_exp = 0;
-			u.s.dbl_frach = u.s.dbl_fracl = 0;
-			errno = ERANGE;
-			return (u.v);		/* zero */
-		}
-		/*
-		 * We are going to produce a denorm.  Our `exp' argument
-		 * might be as small as -2097, and we cannot compute
-		 * 2^-2097, so we may have to do this as many as three
-		 * steps (not just two, as for positive `exp's below).
-		 */
-		mul.v = 0;
-		while (exp <= -DBL_EXP_BIAS) {
-			mul.s.dbl_exp = 1;
-			val *= mul.v;
-			exp += DBL_EXP_BIAS - 1;
-		}
-		mul.s.dbl_exp = exp + DBL_EXP_BIAS;
-		val *= mul.v;
-		return (val);
-	}
-
-	/*
-	 * Newexp is positive.
-	 *
-	 * If oldexp is zero, we are starting with a denorm, and simply
-	 * adjusting the exponent will produce bogus answers.  We need
-	 * to fix that first.
-	 */
-	if (oldexp == 0) {
-		/*
-		 * Multiply by 2^mulexp to make the number normalizable.
-		 * We cannot multiply by more than 2^1023, but `exp'
-		 * argument might be as large as 2046.  A single
-		 * adjustment, however, will normalize the number even
-		 * for huge `exp's, and then we can use exponent
-		 * arithmetic just as for normal `double's.
-		 */
-		mulexp = exp <= DBL_EXP_BIAS ? exp : DBL_EXP_BIAS;
-		mul.v = 0;
-		mul.s.dbl_exp = mulexp + DBL_EXP_BIAS;
-		val *= mul.v;
-		if (mulexp == exp)
-			return (val);
-		u.v = val;
-		newexp -= mulexp;
-	}
-
-	/*
-	 * Both oldexp and newexp are positive; just replace the
-	 * old exponent with the new one.
-	 */
-	u.s.dbl_exp = newexp;
-	return (u.v);
-}
diff --git a/lib/libc/arm/gen/ldexp.c b/lib/libc/arm/gen/ldexp.c
deleted file mode 100644
index 9aed74d..0000000
--- a/lib/libc/arm/gen/ldexp.c
+++ /dev/null
@@ -1,155 +0,0 @@
-/*-
- * Copyright (c) 1999 The NetBSD Foundation, Inc.
- * All rights reserved.
- *
- * This code is derived from software contributed to The NetBSD Foundation
- * by Charles M. Hannum.
- *
- * 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. All advertising materials mentioning features or use of this software
- *    must display the following acknowledgement:
- *        This product includes software developed by the NetBSD
- *        Foundation, Inc. and its contributors.
- * 4. Neither the name of The NetBSD Foundation 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 NETBSD FOUNDATION, INC. 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 FOUNDATION 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/cdefs.h>
-#if 0
-#if defined(LIBC_SCCS) && !defined(lint)
-__RCSID("$NetBSD: ldexp.c,v 1.8 1999/08/30 18:28:26 mycroft Exp $");
-#endif /* LIBC_SCCS and not lint */
-#endif
-__FBSDID("$FreeBSD$");
-
-#include <sys/types.h>
-#include <errno.h>
-#include <math.h>
-#include <machine/ieee.h>
-
-/*
- * Multiply the given value by 2^exp.
- */
-double
-ldexp(val, exp)
-	double val;
-	int exp;
-{
-	int oldexp, newexp;
-	union {
-		double v;
-		struct ieee_double s;
-	} u, mul;
-
-	u.v = val;
-	oldexp = u.s.dbl_exp;
-
-	/*
-	 * If input is zero, Inf or NaN, just return it.
-	 */
-	if (u.v == 0.0 || oldexp == DBL_EXP_INFNAN)
-		return (val);
-
-	if (oldexp == 0) {
-		/*
-		 * u.v is denormal.  We must adjust it so that the exponent
-		 * arithmetic below will work.
-		 */
-		if (exp <= DBL_EXP_BIAS) {
-			/*
-			 * Optimization: if the scaling can be done in a single
-			 * multiply, or underflows, just do it now.
-			 */
-			if (exp <= -DBL_FRACBITS) {
-				errno = ERANGE;
-				return (0.0);
-			}
-			mul.v = 0.0;
-			mul.s.dbl_exp = exp + DBL_EXP_BIAS;
-			u.v *= mul.v;
-			if (u.v == 0.0) {
-				errno = ERANGE;
-				return (0.0);
-			}
-			return (u.v);
-		} else {
-			/*
-			 * We know that exp is very large, and therefore the
-			 * result cannot be denormal (though it may be Inf).
-			 * Shift u.v by just enough to make it normal.
-			 */
-			mul.v = 0.0;
-			mul.s.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS;
-			u.v *= mul.v;
-			exp -= DBL_FRACBITS;
-			oldexp = u.s.dbl_exp;
-		}
-	}
-
-	/*
-	 * u.v is now normalized and oldexp has been adjusted if necessary.
-	 * Calculate the new exponent and check for underflow and overflow.
-	 */
-	newexp = oldexp + exp;
-
-	if (newexp <= 0) {
-		/*
-		 * The output number is either denormal or underflows (see
-		 * comments in machine/ieee.h).
-		 */
-		if (newexp <= -DBL_FRACBITS) {
-			errno = ERANGE;
-			return (0.0);
-		}
-		/*
-		 * Denormalize the result.  We do this with a multiply.  If exp
-		 * is very large, it won't fit in a double, so we have to
-		 * adjust the exponent first.  This is safe because we know
-		 * that u.v is normal at this point.
-		 */
-		if (exp <= -DBL_EXP_BIAS) {
-			u.s.dbl_exp = 1;
-			exp += oldexp - 1;
-		}
-		mul.v = 0.0;
-		mul.s.dbl_exp = exp + DBL_EXP_BIAS;
-		u.v *= mul.v;
-		return (u.v);
-	} else if (newexp >= DBL_EXP_INFNAN) {
-		/*
-		 * The result overflowed; return +/-Inf.
-		 */
-		u.s.dbl_exp = DBL_EXP_INFNAN;
-		u.s.dbl_frach = 0;
-		u.s.dbl_fracl = 0;
-		errno = ERANGE;
-		return (u.v);
-	} else {
-		/*
-		 * The result is normal; just replace the old exponent with the
-		 * new one.
-		 */
-		u.s.dbl_exp = newexp;
-		return (u.v);
-	}
-}
diff --git a/lib/libc/gen/ldexp.c b/lib/libc/gen/ldexp.c
new file mode 100644
index 0000000..887f673
--- /dev/null
+++ b/lib/libc/gen/ldexp.c
@@ -0,0 +1,123 @@
+/* @(#)s_scalbn.c 5.1 93/09/24 */
+/* @(#)fdlibm.h 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 <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <sys/types.h>
+#include <machine/endian.h>
+#include <math.h>
+
+/* Bit fiddling routines copied from msun/src/math_private.h,v 1.15 */
+
+#if BYTE_ORDER == BIG_ENDIAN
+
+typedef union
+{
+  double value;
+  struct
+  {
+    u_int32_t msw;
+    u_int32_t lsw;
+  } parts;
+} ieee_double_shape_type;
+
+#endif
+
+#if BYTE_ORDER == LITTLE_ENDIAN
+
+typedef union
+{
+  double value;
+  struct
+  {
+    u_int32_t lsw;
+    u_int32_t msw;
+  } parts;
+} ieee_double_shape_type;
+
+#endif
+
+/* Get two 32 bit ints from a double.  */
+
+#define EXTRACT_WORDS(ix0,ix1,d)				\
+do {								\
+  ieee_double_shape_type ew_u;					\
+  ew_u.value = (d);						\
+  (ix0) = ew_u.parts.msw;					\
+  (ix1) = ew_u.parts.lsw;					\
+} while (0)
+
+/* Get the more significant 32 bit int from a double.  */
+
+#define GET_HIGH_WORD(i,d)					\
+do {								\
+  ieee_double_shape_type gh_u;					\
+  gh_u.value = (d);						\
+  (i) = gh_u.parts.msw;						\
+} while (0)
+
+/* Set the more significant 32 bits of a double from an int.  */
+
+#define SET_HIGH_WORD(d,v)					\
+do {								\
+  ieee_double_shape_type sh_u;					\
+  sh_u.value = (d);						\
+  sh_u.parts.msw = (v);						\
+  (d) = sh_u.value;						\
+} while (0)
+
+
+static const double
+two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
+huge   = 1.0e+300,
+tiny   = 1.0e-300;
+
+static double
+_copysign(double x, double y)
+{
+	u_int32_t hx,hy;
+	GET_HIGH_WORD(hx,x);
+	GET_HIGH_WORD(hy,y);
+	SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000));
+	return x;
+}
+
+double
+ldexp(double x, int n)
+{
+	int32_t k,hx,lx;
+	EXTRACT_WORDS(hx,lx,x);
+        k = (hx&0x7ff00000)>>20;		/* extract exponent */
+        if (k==0) {				/* 0 or subnormal x */
+            if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
+	    x *= two54;
+	    GET_HIGH_WORD(hx,x);
+	    k = ((hx&0x7ff00000)>>20) - 54;
+            if (n< -50000) return tiny*x; 	/*underflow*/
+	    }
+        if (k==0x7ff) return x+x;		/* NaN or Inf */
+        k = k+n;
+        if (k >  0x7fe) return huge*_copysign(huge,x); /* overflow  */
+        if (k > 0) 				/* normal result */
+	    {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
+        if (k <= -54) {
+            if (n > 50000) 	/* in case integer overflow in n+k */
+		return huge*_copysign(huge,x);	/*overflow*/
+	    else return tiny*_copysign(tiny,x); 	/*underflow*/
+	}
+        k += 54;				/* subnormal result */
+	SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
+        return x*twom54;
+}
diff --git a/lib/libc/ia64/gen/ldexp.c b/lib/libc/ia64/gen/ldexp.c
deleted file mode 100644
index c431417..0000000
--- a/lib/libc/ia64/gen/ldexp.c
+++ /dev/null
@@ -1,137 +0,0 @@
-/*	$NetBSD: ldexp.c,v 1.1 1995/02/10 17:50:24 cgd Exp $	*/
-
-/*
- * Copyright (c) 1994, 1995 Carnegie-Mellon University.
- * All rights reserved.
- *
- * Author: Chris G. Demetriou
- * 
- * Permission to use, copy, modify and distribute this software and
- * its documentation is hereby granted, provided that both the copyright
- * notice and this permission notice appear in all copies of the
- * software, derivative works or modified versions, and any portions
- * thereof, and that both notices appear in supporting documentation.
- * 
- * CARNEGIE MELLON ALLOWS FREE USE OF THIS SOFTWARE IN ITS "AS IS" 
- * CONDITION.  CARNEGIE MELLON DISCLAIMS ANY LIABILITY OF ANY KIND 
- * FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE USE OF THIS SOFTWARE.
- * 
- * Carnegie Mellon requests users of this software to return to
- *
- *  Software Distribution Coordinator  or  Software.Distribution@CS.CMU.EDU
- *  School of Computer Science
- *  Carnegie Mellon University
- *  Pittsburgh PA 15213-3890
- *
- * any improvements or extensions that they make and grant Carnegie the
- * rights to redistribute these changes.
- */
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
-
-#include <sys/types.h>
-#include <machine/ieee.h>
-#include <errno.h>
-#include <math.h>
-
-/*
- * double ldexp(double val, int exp)
- * returns: val * (2**exp)
- */
-double
-ldexp(val, exp)
-	double val;
-	int exp;
-{
-	int oldexp, newexp, mulexp;
-	union doub {
-		double v;
-		struct ieee_double s;
-	} u, mul;
-
-	/*
-	 * If input is zero, or no change, just return input.
-	 * Likewise, if input is Inf or NaN, just return it.
-	 */
-	u.v = val;
-	oldexp = u.s.dbl_exp;
-	if (val == 0 || exp == 0 || oldexp == DBL_EXP_INFNAN)
-		return (val);
-
-	/*
-	 * Compute new exponent and check for over/under flow.
-	 * Underflow, unfortunately, could mean switching to denormal.
-	 * If result out of range, set ERANGE and return 0 if too small
-	 * or Inf if too big, with the same sign as the input value.
-	 */
-	newexp = oldexp + exp;
-	if (newexp >= DBL_EXP_INFNAN) {
-		/* u.s.dbl_sign = val < 0; -- already set */
-		u.s.dbl_exp = DBL_EXP_INFNAN;
-		u.s.dbl_frach = u.s.dbl_fracl = 0;
-		errno = ERANGE;
-		return (u.v);		/* Inf */
-	}
-	if (newexp <= 0) {
-		/*
-		 * The output number is either a denormal or underflows
-		 * (see comments in machine/ieee.h).
-		 */
-		if (newexp <= -DBL_FRACBITS) {
-			/* u.s.dbl_sign = val < 0; -- already set */
-			u.s.dbl_exp = 0;
-			u.s.dbl_frach = u.s.dbl_fracl = 0;
-			errno = ERANGE;
-			return (u.v);		/* zero */
-		}
-		/*
-		 * We are going to produce a denorm.  Our `exp' argument
-		 * might be as small as -2097, and we cannot compute
-		 * 2^-2097, so we may have to do this as many as three
-		 * steps (not just two, as for positive `exp's below).
-		 */
-		mul.v = 0;
-		while (exp <= -DBL_EXP_BIAS) {
-			mul.s.dbl_exp = 1;
-			val *= mul.v;
-			exp += DBL_EXP_BIAS - 1;
-		}
-		mul.s.dbl_exp = exp + DBL_EXP_BIAS;
-		val *= mul.v;
-		return (val);
-	}
-
-	/*
-	 * Newexp is positive.
-	 *
-	 * If oldexp is zero, we are starting with a denorm, and simply
-	 * adjusting the exponent will produce bogus answers.  We need
-	 * to fix that first.
-	 */
-	if (oldexp == 0) {
-		/*
-		 * Multiply by 2^mulexp to make the number normalizable.
-		 * We cannot multiply by more than 2^1023, but `exp'
-		 * argument might be as large as 2046.  A single
-		 * adjustment, however, will normalize the number even
-		 * for huge `exp's, and then we can use exponent
-		 * arithmetic just as for normal `double's.
-		 */
-		mulexp = exp <= DBL_EXP_BIAS ? exp : DBL_EXP_BIAS;
-		mul.v = 0;
-		mul.s.dbl_exp = mulexp + DBL_EXP_BIAS;
-		val *= mul.v;
-		if (mulexp == exp)
-			return (val);
-		u.v = val;
-		newexp -= mulexp;
-	}
-
-	/*
-	 * Both oldexp and newexp are positive; just replace the
-	 * old exponent with the new one.
-	 */
-	u.s.dbl_exp = newexp;
-	return (u.v);
-}
diff --git a/lib/libc/powerpc/gen/ldexp.c b/lib/libc/powerpc/gen/ldexp.c
deleted file mode 100644
index 3a85612..0000000
--- a/lib/libc/powerpc/gen/ldexp.c
+++ /dev/null
@@ -1,155 +0,0 @@
-/*-
- * Copyright (c) 1999 The NetBSD Foundation, Inc.
- * All rights reserved.
- *
- * This code is derived from software contributed to The NetBSD Foundation
- * by Charles M. Hannum.
- *
- * 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. All advertising materials mentioning features or use of this software
- *    must display the following acknowledgement:
- *        This product includes software developed by the NetBSD
- *        Foundation, Inc. and its contributors.
- * 4. Neither the name of The NetBSD Foundation 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 NETBSD FOUNDATION, INC. 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 FOUNDATION 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/cdefs.h>
-#if 0
-#if defined(LIBC_SCCS) && !defined(lint)
-__RCSID("$NetBSD: ldexp.c,v 1.8 1999/08/30 18:28:26 mycroft Exp $");
-#endif /* LIBC_SCCS and not lint */
-#endif
-__FBSDID("$FreeBSD$");
-
-#include <sys/types.h>
-#include <machine/ieee.h>
-#include <errno.h>
-#include <math.h>
-
-/*
- * Multiply the given value by 2^exp.
- */
-double
-ldexp(val, exp)
-	double val;
-	int exp;
-{
-	int oldexp, newexp;
-	union {
-		double v;
-		struct ieee_double s;
-	} u, mul;
-
-	u.v = val;
-	oldexp = u.s.dbl_exp;
-
-	/*
-	 * If input is zero, Inf or NaN, just return it.
-	 */
-	if (u.v == 0.0 || oldexp == DBL_EXP_INFNAN)
-		return (val);
-
-	if (oldexp == 0) {
-		/*
-		 * u.v is denormal.  We must adjust it so that the exponent
-		 * arithmetic below will work.
-		 */
-		if (exp <= DBL_EXP_BIAS) {
-			/*
-			 * Optimization: if the scaling can be done in a single
-			 * multiply, or underflows, just do it now.
-			 */
-			if (exp <= -DBL_FRACBITS) {
-				errno = ERANGE;
-				return (0.0);
-			}
-			mul.v = 0.0;
-			mul.s.dbl_exp = exp + DBL_EXP_BIAS;
-			u.v *= mul.v;
-			if (u.v == 0.0) {
-				errno = ERANGE;
-				return (0.0);
-			}
-			return (u.v);
-		} else {
-			/*
-			 * We know that exp is very large, and therefore the
-			 * result cannot be denormal (though it may be Inf).
-			 * Shift u.v by just enough to make it normal.
-			 */
-			mul.v = 0.0;
-			mul.s.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS;
-			u.v *= mul.v;
-			exp -= DBL_FRACBITS;
-			oldexp = u.s.dbl_exp;
-		}
-	}
-
-	/*
-	 * u.v is now normalized and oldexp has been adjusted if necessary.
-	 * Calculate the new exponent and check for underflow and overflow.
-	 */
-	newexp = oldexp + exp;
-
-	if (newexp <= 0) {
-		/*
-		 * The output number is either denormal or underflows (see
-		 * comments in machine/ieee.h).
-		 */
-		if (newexp <= -DBL_FRACBITS) {
-			errno = ERANGE;
-			return (0.0);
-		}
-		/*
-		 * Denormalize the result.  We do this with a multiply.  If exp
-		 * is very large, it won't fit in a double, so we have to
-		 * adjust the exponent first.  This is safe because we know
-		 * that u.v is normal at this point.
-		 */
-		if (exp <= -DBL_EXP_BIAS) {
-			u.s.dbl_exp = 1;
-			exp += oldexp - 1;
-		}
-		mul.v = 0.0;
-		mul.s.dbl_exp = exp + DBL_EXP_BIAS;
-		u.v *= mul.v;
-		return (u.v);
-	} else if (newexp >= DBL_EXP_INFNAN) {
-		/*
-		 * The result overflowed; return +/-Inf.
-		 */
-		u.s.dbl_exp = DBL_EXP_INFNAN;
-		u.s.dbl_frach = 0;
-		u.s.dbl_fracl = 0;
-		errno = ERANGE;
-		return (u.v);
-	} else {
-		/*
-		 * The result is normal; just replace the old exponent with the
-		 * new one.
-		 */
-		u.s.dbl_exp = newexp;
-		return (u.v);
-	}
-}
diff --git a/lib/libc/sparc64/gen/ldexp.c b/lib/libc/sparc64/gen/ldexp.c
deleted file mode 100644
index 3a85612..0000000
--- a/lib/libc/sparc64/gen/ldexp.c
+++ /dev/null
@@ -1,155 +0,0 @@
-/*-
- * Copyright (c) 1999 The NetBSD Foundation, Inc.
- * All rights reserved.
- *
- * This code is derived from software contributed to The NetBSD Foundation
- * by Charles M. Hannum.
- *
- * 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. All advertising materials mentioning features or use of this software
- *    must display the following acknowledgement:
- *        This product includes software developed by the NetBSD
- *        Foundation, Inc. and its contributors.
- * 4. Neither the name of The NetBSD Foundation 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 NETBSD FOUNDATION, INC. 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 FOUNDATION 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/cdefs.h>
-#if 0
-#if defined(LIBC_SCCS) && !defined(lint)
-__RCSID("$NetBSD: ldexp.c,v 1.8 1999/08/30 18:28:26 mycroft Exp $");
-#endif /* LIBC_SCCS and not lint */
-#endif
-__FBSDID("$FreeBSD$");
-
-#include <sys/types.h>
-#include <machine/ieee.h>
-#include <errno.h>
-#include <math.h>
-
-/*
- * Multiply the given value by 2^exp.
- */
-double
-ldexp(val, exp)
-	double val;
-	int exp;
-{
-	int oldexp, newexp;
-	union {
-		double v;
-		struct ieee_double s;
-	} u, mul;
-
-	u.v = val;
-	oldexp = u.s.dbl_exp;
-
-	/*
-	 * If input is zero, Inf or NaN, just return it.
-	 */
-	if (u.v == 0.0 || oldexp == DBL_EXP_INFNAN)
-		return (val);
-
-	if (oldexp == 0) {
-		/*
-		 * u.v is denormal.  We must adjust it so that the exponent
-		 * arithmetic below will work.
-		 */
-		if (exp <= DBL_EXP_BIAS) {
-			/*
-			 * Optimization: if the scaling can be done in a single
-			 * multiply, or underflows, just do it now.
-			 */
-			if (exp <= -DBL_FRACBITS) {
-				errno = ERANGE;
-				return (0.0);
-			}
-			mul.v = 0.0;
-			mul.s.dbl_exp = exp + DBL_EXP_BIAS;
-			u.v *= mul.v;
-			if (u.v == 0.0) {
-				errno = ERANGE;
-				return (0.0);
-			}
-			return (u.v);
-		} else {
-			/*
-			 * We know that exp is very large, and therefore the
-			 * result cannot be denormal (though it may be Inf).
-			 * Shift u.v by just enough to make it normal.
-			 */
-			mul.v = 0.0;
-			mul.s.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS;
-			u.v *= mul.v;
-			exp -= DBL_FRACBITS;
-			oldexp = u.s.dbl_exp;
-		}
-	}
-
-	/*
-	 * u.v is now normalized and oldexp has been adjusted if necessary.
-	 * Calculate the new exponent and check for underflow and overflow.
-	 */
-	newexp = oldexp + exp;
-
-	if (newexp <= 0) {
-		/*
-		 * The output number is either denormal or underflows (see
-		 * comments in machine/ieee.h).
-		 */
-		if (newexp <= -DBL_FRACBITS) {
-			errno = ERANGE;
-			return (0.0);
-		}
-		/*
-		 * Denormalize the result.  We do this with a multiply.  If exp
-		 * is very large, it won't fit in a double, so we have to
-		 * adjust the exponent first.  This is safe because we know
-		 * that u.v is normal at this point.
-		 */
-		if (exp <= -DBL_EXP_BIAS) {
-			u.s.dbl_exp = 1;
-			exp += oldexp - 1;
-		}
-		mul.v = 0.0;
-		mul.s.dbl_exp = exp + DBL_EXP_BIAS;
-		u.v *= mul.v;
-		return (u.v);
-	} else if (newexp >= DBL_EXP_INFNAN) {
-		/*
-		 * The result overflowed; return +/-Inf.
-		 */
-		u.s.dbl_exp = DBL_EXP_INFNAN;
-		u.s.dbl_frach = 0;
-		u.s.dbl_fracl = 0;
-		errno = ERANGE;
-		return (u.v);
-	} else {
-		/*
-		 * The result is normal; just replace the old exponent with the
-		 * new one.
-		 */
-		u.s.dbl_exp = newexp;
-		return (u.v);
-	}
-}