Floating point support functions

lib/ack/math/ldexp.c

@@ -14,7 +14,7 @@ ldexp(double fl, int exp)
 	int sign = 1;
 	int currexp;
 
-	if (__IsNan(fl)) {
+	if (isnan(fl)) {
 		errno = EDOM;
 		return fl;
 	}

lib/i386/Makefile.in

@@ -2,6 +2,7 @@
 
 SUBDIRS="\
 	int64 \
+	math \
 	misc \
 	rts \
 	string"

lib/i386/math/Makefile.in

@@ -0,0 +1,17 @@
+# Makefile for lib/i386/math.
+
+CFLAGS="-O -D_MINIX -D_POSIX_SOURCE"
+
+LIBRARIES=libc
+
+libc_FILES="           \
+	arch_compare.c \
+	arch_round.c   \
+	fpu_cw.s       \
+	fpu_sw.s       \
+	fpu_round.s    \
+	fegetround.c   \
+	feholdexcept.c \
+	fesetround.c"
+
+TYPE=both

lib/i386/math/arch_compare.c

@@ -0,0 +1,118 @@
+#include <assert.h>
+#include <math.h>
+
+#include "fpu_cw.h"
+#include "fpu_sw.h"
+
+#define FPUSW_FLAG_MASK	\
+	(FPUSW_CONDITION_C0 | FPUSW_CONDITION_C2 | FPUSW_CONDITION_C3)
+
+#define DOUBLE_NORMAL_MIN 2.2250738585072013830902327173324e-308 /* 2^-1022 */
+
+int fpclassify(double x)
+{
+	/* use status word returned by fpu_xam to determine type */
+	switch (fpu_xam(x) & FPUSW_FLAG_MASK)
+	{
+		case FPUSW_CONDITION_C0:
+			return FP_NAN;
+
+		case FPUSW_CONDITION_C2:
+			/* 
+			 * unfortunately, fxam always operates on long doubles
+			 * regardless of the precision setting. This means some
+			 * subnormals are incorrectly classified as normals, 
+			 * since they can be normalized using the additional
+			 * exponent bits available. However, if we already know
+			 * that the number is normal as a long double, finding
+			 * out whether it would be subnormal as a double is just
+			 * a simple comparison.
+			 */
+			if (-DOUBLE_NORMAL_MIN < x && x < DOUBLE_NORMAL_MIN)
+				return FP_SUBNORMAL;
+			else
+				return FP_NORMAL;
+
+		case FPUSW_CONDITION_C0 | FPUSW_CONDITION_C2:
+			return FP_INFINITE;
+
+		case FPUSW_CONDITION_C3:
+			return FP_ZERO;
+
+		case FPUSW_CONDITION_C3 | FPUSW_CONDITION_C2:
+			return FP_SUBNORMAL;
+	}
+
+	/* we don't expect any other case: unsupported, emtpy or reserved */
+	assert(0);
+	return -1;
+}
+
+int signbit(double x)
+{
+	u16_t sw;
+
+	/* examine and use status word to determine sign */
+	sw = fpu_xam(x);
+	return (sw & FPUSW_CONDITION_C1) ? 1 : 0;
+}
+
+#define FPUSW_GREATER	0
+#define FPUSW_LESS	FPUSW_CONDITION_C0
+#define FPUSW_EQUAL	FPUSW_CONDITION_C3
+#define FPUSW_UNORDERED	\
+	(FPUSW_CONDITION_C0 | FPUSW_CONDITION_C2 | FPUSW_CONDITION_C3)
+
+static int fpcompare(double x, double y, u16_t sw1, u16_t sw2)
+{
+	u16_t sw;
+
+	/* compare and check sanity */
+	sw = fpu_compare(x, y) & FPUSW_FLAG_MASK;
+	switch (sw)
+	{
+		case FPUSW_GREATER:
+		case FPUSW_LESS:
+		case FPUSW_EQUAL:
+		case FPUSW_UNORDERED:
+			break;
+
+		default:
+			/* other values are not possible (see IA32 Dev Man) */
+			assert(0);
+			return -1;
+	}
+
+	/* test whether FPUSW equals either sw1 or sw2 */
+	return sw == sw1 || sw == sw2;
+}
+
+int isgreater(double x, double y)
+{
+	return fpcompare(x, y, FPUSW_GREATER, -1);
+}
+
+int isgreaterequal(double x, double y)
+{
+	return fpcompare(x, y, FPUSW_GREATER, FPUSW_EQUAL);
+}
+
+int isless(double x, double y)
+{
+	return fpcompare(x, y, FPUSW_LESS, -1);
+}
+
+int islessequal(double x, double y)
+{
+	return fpcompare(x, y, FPUSW_LESS, FPUSW_EQUAL);
+}
+
+int islessgreater(double x, double y)
+{
+	return fpcompare(x, y, FPUSW_LESS, FPUSW_GREATER);
+}
+
+int isunordered(double x, double y)
+{
+	return fpcompare(x, y, FPUSW_UNORDERED, -1);
+}

lib/i386/math/arch_round.c

@@ -0,0 +1,57 @@
+#include <errno.h>
+#include <fenv.h>
+#include <math.h>
+
+#include "fpu_cw.h"
+#include "fpu_round.h"
+
+static double rndint(double x, u16_t cw_bits, u16_t cw_mask)
+{
+	u16_t cw;
+
+	/* set FPUCW to the right value */
+	cw = fpu_cw_get();
+	fpu_cw_set((cw & cw_mask) | cw_bits);
+
+	/* perform the round */
+	fpu_rndint(&x);
+	
+	/* restore FPUCW */
+	fpu_cw_set(cw);
+	return x;
+}
+
+double nearbyint(double x)
+{
+	/* round, disabling floating point precision error */
+	return rndint(x, FPUCW_EXCEPTION_MASK_PM, ~0);
+}
+
+double remainder(double x, double y)
+{
+	int xclass, yclass;
+
+	/* check arguments */
+	xclass = fpclassify(x);
+	yclass = fpclassify(y);
+	if (xclass == FP_NAN || yclass == FP_NAN)
+		return NAN;
+
+	if (xclass == FP_INFINITE || yclass == FP_ZERO)
+	{
+		errno = EDOM;
+		return NAN;
+	}
+
+	/* call the assembly implementation */
+	fpu_remainder(&x, y);
+	return x;
+}
+
+double trunc(double x)
+{
+	/* round in truncate mode, disabling floating point precision error */
+	return rndint(x, 
+		FPUCW_EXCEPTION_MASK_PM | FPUCW_ROUNDING_CONTROL_TRUNC, 
+		~FPUCW_ROUNDING_CONTROL);
+}

lib/i386/math/fegetround.c

@@ -0,0 +1,23 @@
+#include <assert.h>
+#include <fenv.h>
+
+#include "fpu_cw.h"
+
+int fegetround(void)
+{
+	u16_t cw;
+
+	/* read and categorize FPUCW */
+	cw = fpu_cw_get();
+	switch (cw & FPUCW_ROUNDING_CONTROL)
+	{
+		case FPUCW_ROUNDING_CONTROL_NEAREST: return FE_TONEAREST;
+		case FPUCW_ROUNDING_CONTROL_DOWN:    return FE_DOWNWARD;
+		case FPUCW_ROUNDING_CONTROL_UP:      return FE_UPWARD;
+		case FPUCW_ROUNDING_CONTROL_TRUNC:   return FE_TOWARDZERO;
+	}
+
+	/* each case has been handled, otherwise the constants are wrong */
+	assert(0);
+	return -1;
+}

lib/i386/math/feholdexcept.c

@@ -0,0 +1,17 @@
+#include <errno.h>
+#include <fenv.h>
+
+#include "fpu_cw.h"
+#include "fpu_sw.h"
+
+int feholdexcept(fenv_t *envp)
+{
+	/* read FPUCW and FPUSW */
+	envp->cw = fpu_cw_get();
+	envp->sw = fpu_sw_get();
+	
+	/* update FPUCW to block exceptions */
+	fpu_cw_set(envp->cw | FPUCW_EXCEPTION_MASK);
+
+	return 0;
+}

lib/i386/math/fesetround.c

@@ -0,0 +1,27 @@
+#include <errno.h>
+#include <fenv.h>
+
+#include "fpu_cw.h"
+
+int fesetround(int round)
+{
+	u16_t cw;
+
+	/* read and update FPUCW */
+	cw = fpu_cw_get() & ~FPUCW_ROUNDING_CONTROL;
+	switch (round)
+	{
+		case FE_TONEAREST: cw |= FPUCW_ROUNDING_CONTROL_NEAREST; break;
+		case FE_DOWNWARD: cw |= FPUCW_ROUNDING_CONTROL_DOWN; break;
+		case FE_UPWARD: cw |= FPUCW_ROUNDING_CONTROL_UP; break;
+		case FE_TOWARDZERO: cw |= FPUCW_ROUNDING_CONTROL_TRUNC; break;
+
+		default:
+			errno = EINVAL;
+			return -1;
+	}
+
+	/* set FPUCW to the updated value */
+	fpu_cw_set(cw);
+	return 0;
+}

lib/i386/math/fpu_cw.h

@@ -0,0 +1,33 @@
+#ifndef __FPU_CW__
+#define __FPU_CW__
+
+#include <stdint.h>
+
+/* 
+ * see section 8.1.5 "x87 FPU Control Word" in "Intel 64 and IA-32 Architectures 
+ * Software Developer's Manual Volume 1 Basic Architecture" 
+ */
+#define FPUCW_EXCEPTION_MASK    0x003f
+#define FPUCW_EXCEPTION_MASK_IM 0x0001
+#define FPUCW_EXCEPTION_MASK_DM 0x0002
+#define FPUCW_EXCEPTION_MASK_ZM 0x0004
+#define FPUCW_EXCEPTION_MASK_OM 0x0008
+#define FPUCW_EXCEPTION_MASK_UM 0x0010
+#define FPUCW_EXCEPTION_MASK_PM 0x0020
+
+#define FPUCW_PRECISION_CONTROL         0x0300
+#define FPUCW_PRECISION_CONTROL_SINGLE  0x0000
+#define FPUCW_PRECISION_CONTROL_DOUBLE  0x0200
+#define FPUCW_PRECISION_CONTROL_XDOUBLE 0x0300
+
+#define FPUCW_ROUNDING_CONTROL         0x0c00
+#define FPUCW_ROUNDING_CONTROL_NEAREST 0x0000
+#define FPUCW_ROUNDING_CONTROL_DOWN    0x0400
+#define FPUCW_ROUNDING_CONTROL_UP      0x0800
+#define FPUCW_ROUNDING_CONTROL_TRUNC   0x0c00
+
+/* get and set FPU control word */
+u16_t fpu_cw_get(void);
+void fpu_cw_set(u16_t fpu_cw);
+
+#endif /* !defined(__FPU_CW__) */

lib/i386/math/fpu_cw.s

@@ -0,0 +1,20 @@
+!	fpu_cw_get() - get FPU control word	Author: Erik van der Kouwe
+!	fpu_cw_set() - set FPU control word	9 Dec 2009
+.sect .text
+.define _fpu_cw_get
+.define _fpu_cw_set
+
+! u16_t fpu_cw_get(void)
+_fpu_cw_get:
+	! clear unused bits just to be sure
+	xor	eax,	eax
+	push	eax
+	fstcw	(esp)
+	pop	eax
+	ret
+
+! void fpu_cw_set(u16_t fpu_cw)
+_fpu_cw_set:
+	! load control word from parameter
+	fldcw	4(esp)
+	ret

lib/i386/math/fpu_round.h

@@ -0,0 +1,7 @@
+#ifndef __FPU_RNDINT__
+#define __FPU_RNDINT__
+
+void fpu_rndint(double *value);
+void fpu_remainder(double *x, double y);
+
+#endif /* !defined(__FPU_RNDINT__) */

lib/i386/math/fpu_round.s

@@ -0,0 +1,36 @@
+!	fpu_rndint() - round integer	Author: Erik van der Kouwe
+!	                            	17 Dec 2009
+.sect .text
+.define _fpu_rndint
+.define _fpu_remainder
+
+! void fpu_rndint(double *value)
+_fpu_rndint:
+	! move the value onto the floating point stack
+	mov	eax,	4(esp)
+	fldd	(eax)
+
+	! round it (beware of precision exception!)
+	frndint
+
+	! store the result
+	fstpd	(eax)
+	ret
+
+! void fpu_remainder(double *x, double y)
+_fpu_remainder:
+	! move the values onto the floating point stack
+	fldd	8(esp)
+	mov	edx,	4(esp)
+	fldd	(edx)
+
+	! compute remainder, multiple iterations may be needed
+1:	fprem1
+	.data1	0xdf, 0xe0	! fnstsw	ax	
+	sahf
+	jp	1b
+
+	! store the result and pop the divisor
+	fstpd	(edx)
+	fstp	st
+	ret

lib/i386/math/fpu_sw.h

@@ -0,0 +1,26 @@
+#ifndef __FPU_SW__
+#define __FPU_SW__
+
+#include <stdint.h>
+
+#define FPUSW_EXCEPTION_IE  0x0001
+#define FPUSW_EXCEPTION_DE  0x0002
+#define FPUSW_EXCEPTION_ZE  0x0004
+#define FPUSW_EXCEPTION_OE  0x0008
+#define FPUSW_EXCEPTION_UE  0x0010
+#define FPUSW_EXCEPTION_PE  0x0020
+#define FPUSW_STACK_FAULT   0x0040
+#define FPUSW_ERROR_SUMMARY 0x0080
+#define FPUSW_CONDITION_C0  0x0100
+#define FPUSW_CONDITION_C1  0x0200
+#define FPUSW_CONDITION_C2  0x0400
+#define FPUSW_CONDITION_C3  0x4000
+#define FPUSW_BUSY          0x8000
+
+u16_t fpu_compare(double x, double y);
+u16_t fpu_sw_get(void);
+void fpu_sw_set(u16_t value);
+u16_t fpu_xam(double value);
+
+#endif /* !defined(__FPU_SW__) */
+

lib/i386/math/fpu_sw.s

@@ -0,0 +1,38 @@
+! fpu_compare() - compare doubles	Author: Erik van der Kouwe
+! fpu_sw_get() - get FPU status      	17 Dec 2009
+! fpu_xam() - examine double
+.sect .text
+.define _fpu_compare
+.define _fpu_sw_get
+.define _fpu_xam
+
+! u16_t fpu_compare(double x, double y)
+_fpu_compare:
+	! move the values onto the floating point stack
+	fldd	12(esp)
+	fldd	4(esp)
+
+	! compare values and return status word
+	fcompp
+	jmp	_fpu_sw_get
+
+! u16_t fpu_sw_get(void)
+_fpu_sw_get:
+	! clear unused high-order word and get status word
+	xor	eax,	eax
+	.data1	0xdf, 0xe0	! fnstsw	ax
+	ret
+	
+! u16_t fpu_xam(double value)
+_fpu_xam:
+	! move the value onto the floating point stack
+	fldd	4(esp)
+
+	! examine value and get status word
+	fxam
+	call	_fpu_sw_get
+
+	! pop the value
+	fstp	st
+	ret
+

lib/math/Makefile.in

@@ -9,6 +9,7 @@ libc_FILES=" \
 	atan.c \
 	atan2.c \
 	ceil.c \
+	compare.c \
 	exp.c \
 	fabs.c \
 	floor.c \

lib/math/compare.c

@@ -0,0 +1,39 @@
+#include <assert.h>
+#include <math.h>
+
+/* functions missing here are architecture-specific and are in i386/float */
+
+int isfinite(double x)
+{
+	/* return value based on classification */
+	switch (fpclassify(x))
+	{
+		case FP_INFINITE:
+		case FP_NAN:
+			return 0;
+
+		case FP_NORMAL:
+		case FP_SUBNORMAL:
+		case FP_ZERO:
+			return 1;
+	}
+
+	/* if we get here, fpclassify is buggy */
+	assert(0);
+	return -1;
+}
+
+int isinf(double x)
+{
+	return fpclassify(x) == FP_INFINITE;
+}
+
+int isnan(double x)
+{
+	return fpclassify(x) == FP_NAN;
+}
+
+int isnormal(double x)
+{
+	return fpclassify(x) == FP_NORMAL;
+}

lib/math/hugeval.c

@@ -9,7 +9,7 @@
 #include	<math.h>
 
 double
-__huge_val(void)
+__infinity(void)
 {
 #if (CHIP == INTEL)
 	static unsigned char ieee_infinity[] = { 
@@ -21,3 +21,18 @@ __huge_val(void)
 	return 1.0e+1000;	/* This will generate a warning */
 #endif
 }
+
+double
+__qnan(void)
+{
+#if (CHIP == INTEL)
+	static unsigned char ieee_qnan[] = { 
+		0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x7f };
+
+	assert(sizeof(double) == sizeof(ieee_qnan));
+	return *(double *) ieee_qnan;
+#else
+#error QNaN not defined on this architecture
+#endif
+}
+

lib/math/localmath.h

@@ -3,6 +3,8 @@
  */
 /* $Header$ */
 
+#define __IsNan isnan
+
 /* some constants (Hart & Cheney) */
 #define	M_PI		3.14159265358979323846264338327950288
 #define	M_2PI		6.28318530717958647692528676655900576

lib/math/sqrt.c

@@ -9,6 +9,7 @@
 #include	<math.h>
 #include	<float.h>
 #include	<errno.h>
+#include	"localmath.h"
 
 #define NITER	5