math: fix x86 asin, atan, exp, log1p to raise underflow

underflow is raised by an inexact subnormal float store,
since subnormal operations are slow, check the underflow
flag and skip the store if it's already raised

diff --git a/src/math/i386/asin.s b/src/math/i386/asin.s
index 932c7542442787f30361489b3f4e9d57fed7ffec..a9f691bff3b16c21e0ab10a75b2243707117b92a 100644 (file)
--- a/src/math/i386/asin.s
+++ b/src/math/i386/asin.s
@@ -2,7 +2,18 @@
 .type asinf,@function
 asinf:
        flds 4(%esp)
-       jmp 1f
+       mov 4(%esp),%eax
+       add %eax,%eax
+       cmp $0x01000000,%eax
+       jae 1f
+               # subnormal x, return x with underflow
+       fnstsw %ax
+       and $16,%ax
+       jnz 2f
+       fld %st(0)
+       fmul %st(1)
+       fstps 4(%esp)
+2:     ret
 
 .global asinl
 .type asinl,@function
@@ -14,6 +25,16 @@ asinl:
 .type asin,@function
 asin:
        fldl 4(%esp)
+       mov 8(%esp),%eax
+       add %eax,%eax
+       cmp $0x00200000,%eax
+       jae 1f
+               # subnormal x, return x with underflow
+       fnstsw %ax
+       and $16,%ax
+       jnz 2f
+       fsts 4(%esp)
+2:     ret
 1:     fld %st(0)
        fld1
        fsub %st(0),%st(1)

diff --git a/src/math/i386/atan.s b/src/math/i386/atan.s
index 7e28b39510140703ade798cb5a034f12efd54c23..d73137b28eb545c543338383b0a5bbc8947345a7 100644 (file)
--- a/src/math/i386/atan.s
+++ b/src/math/i386/atan.s
@@ -2,6 +2,16 @@
 .type atan,@function
 atan:
        fldl 4(%esp)
+       mov 8(%esp),%eax
+       add %eax,%eax
+       cmp $0x00200000,%eax
+       jb 1f
        fld1
        fpatan
        ret
+               # subnormal x, return x with underflow
+1:     fnstsw %ax
+       and $16,%ax
+       jnz 2f
+       fsts 4(%esp)
+2:     ret

diff --git a/src/math/i386/atanf.s b/src/math/i386/atanf.s
index 3cd402332e75f61a3a29eeaa09cebcc62f58cbea..8caddefa26efc0256e756b4b164def3b317dd148 100644 (file)
--- a/src/math/i386/atanf.s
+++ b/src/math/i386/atanf.s
@@ -2,6 +2,18 @@
 .type atanf,@function
 atanf:
        flds 4(%esp)
+       mov 4(%esp),%eax
+       add %eax,%eax
+       cmp $0x01000000,%eax
+       jb 1f
        fld1
        fpatan
        ret
+               # subnormal x, return x with underflow
+1:     fnstsw %ax
+       and $16,%ax
+       jnz 2f
+       fld %st(0)
+       fmul %st(1)
+       fstps 4(%esp)
+2:     ret

diff --git a/src/math/i386/exp.s b/src/math/i386/exp.s
index e3b42af5ae244fde2e2e30b048d9877d9d1a1177..e5f545888e35e9391950cf6ca52c29abd5460b7e 100644 (file)
--- a/src/math/i386/exp.s
+++ b/src/math/i386/exp.s
@@ -2,7 +2,18 @@
 .type expm1f,@function
 expm1f:
        flds 4(%esp)
-       jmp 1f
+       mov 4(%esp),%eax
+       add %eax,%eax
+       cmp $0x01000000,%eax
+       jae 1f
+               # subnormal x, return x with underflow
+       fnstsw %ax
+       and $16,%ax
+       jnz 2f
+       fld %st(0)
+       fmul %st(1)
+       fstps 4(%esp)
+2:     ret
 
 .global expm1l
 .type expm1l,@function
@@ -14,10 +25,32 @@ expm1l:
 .type expm1,@function
 expm1:
        fldl 4(%esp)
+       mov 8(%esp),%eax
+       add %eax,%eax
+       cmp $0x00200000,%eax
+       jae 1f
+               # subnormal x, return x with underflow
+       fnstsw %ax
+       and $16,%ax
+       jnz 2f
+       fsts 4(%esp)
+2:     ret
 1:     fldl2e
        fmulp
+       mov $0xc2820000,%eax
+       push %eax
+       flds (%esp)
+       pop %eax
+       fucomp %st(1)
+       fnstsw %ax
+       sahf
        fld1
-       fld %st(1)
+       jb 1f
+               # x*log2e < -65, return -1 without underflow
+       fstp %st(1)
+       fchs
+       ret
+1:     fld %st(1)
        fabs
        fucom %st(1)
        fnstsw %ax

diff --git a/src/math/i386/log1p.s b/src/math/i386/log1p.s
index 9971e53c5446f8ace46fd9d1b591f9f927d8f023..6b6929c72800ab35c07bc331995fe1ed2de1ba61 100644 (file)
--- a/src/math/i386/log1p.s
+++ b/src/math/i386/log1p.s
@@ -7,9 +7,18 @@ log1p:
        fldl 4(%esp)
        cmp $0x3fd28f00,%eax
        ja 1f
+       cmp $0x00100000,%eax
+       jb 2f
        fyl2xp1
        ret
 1:     fld1
        faddp
        fyl2x
        ret
+               # subnormal x, return x with underflow
+2:     fnstsw %ax
+       and $16,%ax
+       jnz 1f
+       fsts 4(%esp)
+       fstp %st(1)
+1:     ret

diff --git a/src/math/i386/log1pf.s b/src/math/i386/log1pf.s
index 2680a8a624ed9f4e35555c97250f8f1fbe63f749..c0bcd30f0e6efd46b0b1223d5439097f7bc7febd 100644 (file)
--- a/src/math/i386/log1pf.s
+++ b/src/math/i386/log1pf.s
@@ -7,9 +7,19 @@ log1pf:
        flds 4(%esp)
        cmp $0x3e940000,%eax
        ja 1f
+       cmp $0x00800000,%eax
+       jb 2f
        fyl2xp1
        ret
 1:     fld1
        faddp
        fyl2x
        ret
+               # subnormal x, return x with underflow
+2:     fnstsw %ax
+       and $16,%ax
+       jnz 1f
+       fxch
+       fmul %st(1)
+       fstps 4(%esp)
+1:     ret