about summary refs log tree commit diff
path: root/sysdeps/libm-i387/e_expf.S
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/libm-i387/e_expf.S')
-rw-r--r--sysdeps/libm-i387/e_expf.S42
1 files changed, 42 insertions, 0 deletions
diff --git a/sysdeps/libm-i387/e_expf.S b/sysdeps/libm-i387/e_expf.S
new file mode 100644
index 0000000000..c8233eba50
--- /dev/null
+++ b/sysdeps/libm-i387/e_expf.S
@@ -0,0 +1,42 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ * Adapted for float type by Ulrich Drepper <drepper@cygnus.com>.
+ */
+
+#include <machine/asm.h>
+
+RCSID("$NetBSD: $")
+
+/* e^x = 2^(x * log2(e)) */
+ENTRY(__ieee754_expf)
+	flds	4(%esp)
+/* I added the following ugly construct because exp(+-Inf) resulted
+   in NaN.  The ugliness results from the bright minds at Intel.
+   For the i686 the code can be written better.
+   -- drepper@cygnus.com.  */
+	fxam				/* Is NaN or +-Inf?  */
+	fstsw	%ax
+	sahf
+	jnc	.LnoInfNaN		/* No, jump.   */
+	jp	.LisInf			/* Is +-Inf, jump.  */
+.LnoInfNaN:
+	fldl2e
+	fmulp				/* x * log2(e) */
+	fstl	%st(1)
+	frndint				/* int(x * log2(e)) */
+	fstl	%st(2)
+	fsubrp				/* fract(x * log2(e)) */
+	f2xm1				/* 2^(fract(x * log2(e))) - 1 */
+	fld1
+	faddp				/* 2^(fract(x * log2(e))) */
+	fscale				/* e^x */
+	fstp	%st(1)
+	ret
+
+.LisInf:
+	andb	$2, %ah			/* Test sign.  */
+	jz	.LpInf			/* If positive, jump.  */
+	fldz				/* Set result to 0.  */
+.LpInf:	ret
+END (__ieee754_expf)