- ldexp implementation has issues. switch to the one from libm

- remove frexp in hppa64, cloned from hppa
- move generic ieee754 implementations of modf and ldexp to gen
ok kettenis@, "looks good" millert@

3 files changed, 241 insertions, 4 deletions

diff --git a/lib/libc/gen/Makefile.inc b/lib/libc/gen/Makefile.inc
index 6a55a7d1f2c..b92d8d515c7 100644
--- a/lib/libc/gen/Makefile.inc
+++ b/lib/libc/gen/Makefile.inc
@@ -1,4 +1,4 @@
-#	$OpenBSD: Makefile.inc,v 1.44 2008/12/12 00:12:47 martynas Exp $
+#	$OpenBSD: Makefile.inc,v 1.45 2009/04/19 16:42:06 martynas Exp $
 
 # gen sources
 .PATH: ${LIBCSRCDIR}/arch/${MACHINE_ARCH}/gen ${LIBCSRCDIR}/gen
@@ -12,9 +12,9 @@ SRCS+=  alarm.c assert.c auth_subr.c authenticate.c \
         getloadavg.c getlogin.c getmntinfo.c getnetgrent.c getpagesize.c \
         getpwent.c getttyent.c getusershell.c glob.c initgroups.c \
         isatty.c isctype.c isfdtype.c isfinite.c isinf.c isnan.c \
-        isnormal.c signbit.c lockf.c login_cap.c nice.c nlist.c \
-        nftw.c opendir.c pause.c popen.c psignal.c pw_dup.c pwcache.c \
-        raise.c readdir.c readpassphrase.c rewinddir.c scandir.c \
+        isnormal.c signbit.c ldexp.c lockf.c login_cap.c modf.c nice.c \
+        nlist.c nftw.c opendir.c pause.c popen.c psignal.c pw_dup.c \
+        pwcache.c raise.c readdir.c readpassphrase.c rewinddir.c scandir.c \
         seekdir.c setdomainname.c sethostname.c setjmperr.c setmode.c \
         setproctitle.c siginterrupt.c siglist.c signal.c signame.c \
         sigsetops.c sleep.c statvfs.c strtofflags.c sysconf.c sysctl.c \
diff --git a/lib/libc/gen/ldexp.c b/lib/libc/gen/ldexp.c
new file mode 100644
index 00000000000..451950889b2
--- /dev/null
+++ b/lib/libc/gen/ldexp.c
@@ -0,0 +1,133 @@
+/*	$OpenBSD: ldexp.c,v 1.1 2009/04/19 16:42:06 martynas Exp $	*/
+/* @(#)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.
+ * ====================================================
+ */
+
+#if 0
+__FBSDID("$FreeBSD: src/lib/libc/gen/ldexp.c,v 1.1 2005/01/22 06:03:40 das Exp $");
+#endif
+
+#include <sys/types.h>
+#include <sys/cdefs.h>
+#include <machine/endian.h>
+#include <float.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;
+}
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(ldexpl, ldexp);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libc/gen/modf.c b/lib/libc/gen/modf.c
new file mode 100644
index 00000000000..022cf94d34c
--- /dev/null
+++ b/lib/libc/gen/modf.c
@@ -0,0 +1,104 @@
+/*	$OpenBSD: modf.c,v 1.1 2009/04/19 16:42:06 martynas Exp $	*/
+/*	$NetBSD: modf.c,v 1.1 1995/02/10 17:50:25 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/types.h>
+#include <machine/ieee.h>
+#include <errno.h>
+#include <math.h>
+
+/*
+ * double modf(double val, double *iptr)
+ * returns: f and i such that |f| < 1.0, (f + i) = val, and
+ *	sign(f) == sign(i) == sign(val).
+ *
+ * Beware signedness when doing subtraction, and also operand size!
+ */
+double
+modf(double val, double *iptr)
+{
+	union doub {
+		double v;
+		struct ieee_double s;
+	} u, v;
+	u_int64_t frac;
+
+	/*
+	 * If input is Inf or NaN, return it and leave i alone.
+	 */
+	u.v = val;
+	if (u.s.dbl_exp == DBL_EXP_INFNAN)
+		return (u.v);
+
+	/*
+	 * If input can't have a fractional part, return
+	 * (appropriately signed) zero, and make i be the input.
+	 */
+	if ((int)u.s.dbl_exp - DBL_EXP_BIAS > DBL_FRACBITS - 1) {
+		*iptr = u.v;
+		v.v = 0.0;
+		v.s.dbl_sign = u.s.dbl_sign;
+		return (v.v);
+	}
+
+	/*
+	 * If |input| < 1.0, return it, and set i to the appropriately
+	 * signed zero.
+	 */
+	if (u.s.dbl_exp < DBL_EXP_BIAS) {
+		v.v = 0.0;
+		v.s.dbl_sign = u.s.dbl_sign;
+		*iptr = v.v;
+		return (u.v);
+	}
+
+	/*
+	 * There can be a fractional part of the input.
+	 * If you look at the math involved for a few seconds, it's
+	 * plain to see that the integral part is the input, with the
+	 * low (DBL_FRACBITS - (exponent - DBL_EXP_BIAS)) bits zeroed,
+	 * the fractional part is the part with the rest of the
+	 * bits zeroed.  Just zeroing the high bits to get the
+	 * fractional part would yield a fraction in need of
+	 * normalization.  Therefore, we take the easy way out, and
+	 * just use subtraction to get the fractional part.
+	 */
+	v.v = u.v;
+	/* Zero the low bits of the fraction, the sleazy way. */
+	frac = ((u_int64_t)v.s.dbl_frach << 32) + v.s.dbl_fracl;
+	frac >>= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS);
+	frac <<= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS);
+	v.s.dbl_fracl = frac & 0xffffffff;
+	v.s.dbl_frach = frac >> 32;
+	*iptr = v.v;
+
+	u.v -= v.v;
+	u.s.dbl_sign = v.s.dbl_sign;
+	return (u.v);
+}