From 60fe3744a17a2b618cfea7c4f28c8899518f5a3b Mon Sep 17 00:00:00 2001
From: das <das@FreeBSD.org>
Date: Mon, 7 Mar 2005 04:56:46 +0000
Subject: [PATCH] Implement nexttoward and nextafterl; the latter is also known
 as nexttowardl.  These are not needed on machines where long doubles look
 like IEEE-754 doubles, so the implementation only supports the usual long
 double formats with 15-bit exponents.

Anything bizarre, such as machines where floating-point and integer
data have different endianness, will cause problems.  This is the case
with big endian ia64 according to libc/ia64/_fpmath.h.  Please contact
me if you managed to get a machine running this way.
---
 lib/msun/src/s_nextafterl.c | 82 +++++++++++++++++++++++++++++++++++++
 lib/msun/src/s_nexttoward.c | 73 +++++++++++++++++++++++++++++++++
 2 files changed, 155 insertions(+)
 create mode 100644 lib/msun/src/s_nextafterl.c
 create mode 100644 lib/msun/src/s_nexttoward.c

diff --git a/lib/msun/src/s_nextafterl.c b/lib/msun/src/s_nextafterl.c
new file mode 100644
index 000000000000..a74dd79a9090
--- /dev/null
+++ b/lib/msun/src/s_nextafterl.c
@@ -0,0 +1,82 @@
+/* @(#)s_nextafter.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.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD$";
+#endif
+
+/* IEEE functions
+ *	nextafter(x,y)
+ *	return the next machine floating-point number of x in the
+ *	direction toward y.
+ *   Special cases:
+ */
+
+#include <sys/cdefs.h>
+#include <float.h>
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+#if LDBL_MAX_EXP != 0x4000
+#error "Unsupported long double format"
+#endif
+
+long double
+nextafterl(long double x, long double y)
+{
+	volatile long double t;
+	union IEEEl2bits ux, uy;
+
+	ux.e = x;
+	uy.e = y;
+
+	if ((ux.bits.exp == 0x7fff &&
+	     ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl) != 0) ||
+	    (uy.bits.exp == 0x7fff &&
+	     ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0))
+	   return x+y;	/* x or y is nan */
+	if(x==y) return y;		/* x=y, return y */
+	if(x==0.0) {
+	    ux.bits.manh = 0;			/* return +-minsubnormal */
+	    ux.bits.manl = 1;
+	    ux.bits.sign = uy.bits.sign;
+	    t = ux.e*ux.e;
+	    if(t==ux.e) return t; else return ux.e; /* raise underflow flag */
+	}
+	if(x>0.0 ^ x<y) {			/* x -= ulp */
+	    if(ux.bits.manl==0) {
+		if ((ux.bits.manh&~LDBL_NBIT)==0)
+		    ux.bits.exp -= 1;
+		ux.bits.manh = (ux.bits.manh - 1) | (ux.bits.manh & LDBL_NBIT);
+	    }
+	    ux.bits.manl -= 1;
+	} else {				/* x += ulp */
+	    ux.bits.manl += 1;
+	    if(ux.bits.manl==0) {
+		ux.bits.manh = (ux.bits.manh + 1) | (ux.bits.manh & LDBL_NBIT);
+		if ((ux.bits.manh&~LDBL_NBIT)==0)
+		    ux.bits.exp += 1;
+	    }
+	}
+	if(ux.bits.exp==0x7fff) return x+x;	/* overflow  */
+	if(ux.bits.exp==0) {			/* underflow */
+	    mask_nbit_l(ux);
+	    t = ux.e * ux.e;
+	    if(t!=ux.e)			/* raise underflow flag */
+		return ux.e;
+	}
+	return ux.e;
+}
+
+__strong_reference(nextafterl, nexttowardl);
diff --git a/lib/msun/src/s_nexttoward.c b/lib/msun/src/s_nexttoward.c
new file mode 100644
index 000000000000..a493c83ff317
--- /dev/null
+++ b/lib/msun/src/s_nexttoward.c
@@ -0,0 +1,73 @@
+/* @(#)s_nextafter.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.
+ * ====================================================
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD$";
+#endif
+
+/*
+ * We assume that a long double has a 15-bit exponent.  On systems
+ * where long double is the same as double, nexttoward() is an alias
+ * for nextafter(), so we don't use this routine.
+ */
+
+#include <float.h>
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+#if LDBL_MAX_EXP != 0x4000
+#error "Unsupported long double format"
+#endif
+
+double
+nexttoward(double x, long double y)
+{
+	union IEEEl2bits uy;
+	volatile double t;
+	int32_t hx,ix;
+	u_int32_t lx;
+
+	EXTRACT_WORDS(hx,lx,x);
+	ix = hx&0x7fffffff;		/* |x| */
+	uy.e = y;
+
+	if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) ||
+	    (uy.bits.exp == 0x7fff &&
+	     ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0))
+	   return x+y;	/* x or y is nan */
+	if(x==y) return (double)y;		/* x=y, return y */
+	if(x==0.0) {
+	    INSERT_WORDS(x,uy.bits.sign<<31,1);	/* return +-minsubnormal */
+	    t = x*x;
+	    if(t==x) return t; else return x;	/* raise underflow flag */
+	}
+	if(hx>0.0 ^ x < y) {			/* x -= ulp */
+	    if(lx==0) hx -= 1;
+	    lx -= 1;
+	} else {				/* x += ulp */
+	    lx += 1;
+	    if(lx==0) hx += 1;
+	}
+	ix = hx&0x7ff00000;
+	if(ix>=0x7ff00000) return x+x;	/* overflow  */
+	if(ix<0x00100000) {		/* underflow */
+	    t = x*x;
+	    if(t!=x) {		/* raise underflow flag */
+	        INSERT_WORDS(y,hx,lx);
+		return y;
+	    }
+	}
+	INSERT_WORDS(x,hx,lx);
+	return x;
+}