diff options
Diffstat (limited to 'sysdeps')
| -rw-r--r-- | sysdeps/i386/fpu/s_log1p.S | 27 | ||||
| -rw-r--r-- | sysdeps/i386/fpu/s_log1pf.S | 27 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_log1p.c | 10 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_log1pf.c | 8 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128/s_log1pl.c | 6 |
5 files changed, 75 insertions, 3 deletions
diff --git a/sysdeps/i386/fpu/s_log1p.S b/sysdeps/i386/fpu/s_log1p.S index 8624249dba..c2559a3f18 100644 --- a/sysdeps/i386/fpu/s_log1p.S +++ b/sysdeps/i386/fpu/s_log1p.S @@ -17,6 +17,13 @@ RCSID("$NetBSD: s_log1p.S,v 1.7 1995/05/09 00:10:58 jtc Exp $") limit: .double 0.29 one: .double 1.0 + .section .rodata.cst8,"aM",@progbits,8 + + .p2align 3 + .type dbl_min,@object +dbl_min: .byte 0, 0, 0, 0, 0, 0, 0x10, 0 + ASM_SIZE_DIRECTIVE(dbl_min) + /* * Use the fyl2xp1 function when the argument is in the range -0.29 to 0.29, * otherwise fyl2x with the needed extra computation. @@ -55,7 +62,25 @@ ENTRY(__log1p) ret 2: fyl2xp1 - ret +#ifdef PIC + fldl dbl_min@GOTOFF(%edx) +#else + fldl dbl_min +#endif + fld %st(1) + fabs + fucompp + fnstsw + sahf + jnc 1f + subl $8, %esp + cfi_adjust_cfa_offset (8) + fld %st(0) + fmul %st(0) + fstpl (%esp) + addl $8, %esp + cfi_adjust_cfa_offset (-8) +1: ret 3: jp 4b // in case x is ħInf fstp %st(1) diff --git a/sysdeps/i386/fpu/s_log1pf.S b/sysdeps/i386/fpu/s_log1pf.S index b071e7372d..8fca22e4ff 100644 --- a/sysdeps/i386/fpu/s_log1pf.S +++ b/sysdeps/i386/fpu/s_log1pf.S @@ -17,6 +17,13 @@ RCSID("$NetBSD: s_log1pf.S,v 1.4 1995/05/09 00:13:05 jtc Exp $") limit: .float 0.29 one: .float 1.0 + .section .rodata.cst4,"aM",@progbits,4 + + .p2align 2 + .type flt_min,@object +flt_min: .byte 0, 0, 0x80, 0 + ASM_SIZE_DIRECTIVE(flt_min) + /* * Use the fyl2xp1 function when the argument is in the range -0.29 to 0.29, * otherwise fyl2x with the needed extra computation. @@ -55,7 +62,25 @@ ENTRY(__log1pf) ret 2: fyl2xp1 - ret +#ifdef PIC + flds flt_min@GOTOFF(%edx) +#else + flds flt_min +#endif + fld %st(1) + fabs + fucompp + fnstsw + sahf + jnc 1f + subl $4, %esp + cfi_adjust_cfa_offset (4) + fld %st(0) + fmul %st(0) + fstps (%esp) + addl $4, %esp + cfi_adjust_cfa_offset (-4) +1: ret 3: jp 4b // in case x is ħInf fstp %st(1) diff --git a/sysdeps/ieee754/dbl-64/s_log1p.c b/sysdeps/ieee754/dbl-64/s_log1p.c index 86bbfbacaf..cff555b0aa 100644 --- a/sysdeps/ieee754/dbl-64/s_log1p.c +++ b/sysdeps/ieee754/dbl-64/s_log1p.c @@ -78,6 +78,7 @@ * See HP-15C Advanced Functions Handbook, p.193. */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -118,7 +119,14 @@ __log1p (double x) { math_force_eval (two54 + x); /* raise inexact */ if (ax < 0x3c900000) /* |x| < 2**-54 */ - return x; + { + if (fabs (x) < DBL_MIN) + { + double force_underflow = x * x; + math_force_eval (force_underflow); + } + return x; + } else return x - x * x * 0.5; } diff --git a/sysdeps/ieee754/flt-32/s_log1pf.c b/sysdeps/ieee754/flt-32/s_log1pf.c index 94c33fca16..83a09f1414 100644 --- a/sysdeps/ieee754/flt-32/s_log1pf.c +++ b/sysdeps/ieee754/flt-32/s_log1pf.c @@ -13,6 +13,7 @@ * ==================================================== */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -48,7 +49,14 @@ __log1pf(float x) if(ax<0x31000000) { /* |x| < 2**-29 */ math_force_eval(two25+x); /* raise inexact */ if (ax<0x24800000) /* |x| < 2**-54 */ + { + if (fabsf (x) < FLT_MIN) + { + float force_underflow = x * x; + math_force_eval (force_underflow); + } return x; + } else return x - x*x*(float)0.5; } diff --git a/sysdeps/ieee754/ldbl-128/s_log1pl.c b/sysdeps/ieee754/ldbl-128/s_log1pl.c index b70a55b758..ff759bc000 100644 --- a/sysdeps/ieee754/ldbl-128/s_log1pl.c +++ b/sysdeps/ieee754/ldbl-128/s_log1pl.c @@ -53,6 +53,7 @@ <http://www.gnu.org/licenses/>. */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -140,6 +141,11 @@ __log1pl (long double xm1) if ((hx & 0x7fffffff) < 0x3f8e0000) { + if (fabsl (xm1) < LDBL_MIN) + { + long double force_underflow = xm1 * xm1; + math_force_eval (force_underflow); + } if ((int) xm1 == 0) return xm1; } |