fix broken exponential asm

infinities were getting converted into nans. the new code simply tests
for infinity and replaces it with a large magnitude value of the same
sign.

also, the fcomi instruction is apparently not part of the i387
instruction set, so avoid using it.

diff --git a/src/math/i386/exp.s b/src/math/i386/exp.s
index 18f6cd6..f4769d5 100644 (file)
--- a/src/math/i386/exp.s
+++ b/src/math/i386/exp.s
@@ -34,6 +34,15 @@ exp:
 .type exp2,@function
 exp2:
        fldl 4(%esp)
 .type exp2,@function
 exp2:
        fldl 4(%esp)

+1:     fxam
+       fnstsw %ax
+       sahf
+       jnp 1f
+       jnc 1f
+       fstps 4(%esp)
+       mov $0xfe,%al
+       and %al,7(%esp)
+       flds 4(%esp)

 1:     fld %st(0)
        frndint
        fxch %st(1)
 1:     fld %st(0)
        frndint
        fxch %st(1)

diff --git a/src/math/i386/expm1.s b/src/math/i386/expm1.s
index d6d511e..bbb5d12 100644 (file)
--- a/src/math/i386/expm1.s
+++ b/src/math/i386/expm1.s
@@ -14,12 +14,23 @@ expm1l:
 .type expm1,@function
 expm1:
        fldl 4(%esp)
 .type expm1,@function
 expm1:
        fldl 4(%esp)

+1:     fxam
+       fnstsw %ax
+       sahf
+       jnp 1f
+       jnc 1f
+       fstps 4(%esp)
+       mov $0xfe,%al
+       and %al,7(%esp)
+       flds 4(%esp)

 1:     fldl2e
        fmulp
        fld %st(0)
        frndint
        fldz
 1:     fldl2e
        fmulp
        fld %st(0)
        frndint
        fldz

-       fcompi
+       fcomp
+       fnstsw %ax
+       sahf

        jnz 1f
        fstp %st(0)
        f2xm1
        jnz 1f
        fstp %st(0)
        f2xm1