summaryrefslogtreecommitdiff
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)