diff options
Diffstat (limited to 'sysdeps')
40 files changed, 1005 insertions, 166 deletions
diff --git a/sysdeps/generic/abort.c b/sysdeps/generic/abort.c index 604eadc4f5..366a543791 100644 --- a/sysdeps/generic/abort.c +++ b/sysdeps/generic/abort.c @@ -90,7 +90,7 @@ abort (void) if (stage == 4) { ++stage; - fclose (NULL); + __fcloseall (); } /* Try again. */ diff --git a/sysdeps/gnu/utmpbits.h b/sysdeps/gnu/utmpbits.h index 9fcc77de6b..de4a830bcc 100644 --- a/sysdeps/gnu/utmpbits.h +++ b/sysdeps/gnu/utmpbits.h @@ -96,7 +96,7 @@ struct utmp long ut_session; /* Session ID, used for windowing. */ struct timeval ut_tv; /* Time entry was made. */ int32_t ut_addr_v6[4]; /* Internet address of remote host. */ - char pad[18]; /* Reserved for future use. */ + char pad[20]; /* Reserved for future use. */ }; /* Backwards compatibility hacks. */ diff --git a/sysdeps/ieee754/Makefile b/sysdeps/ieee754/Makefile new file mode 100644 index 0000000000..3796e728aa --- /dev/null +++ b/sysdeps/ieee754/Makefile @@ -0,0 +1,3 @@ +ifeq ($(subdir),math) +sysdep_headers += ieee754.h +endif diff --git a/sysdeps/ieee754/ieee754.h b/sysdeps/ieee754/ieee754.h index 9e6efe4973..ba40cc7683 100644 --- a/sysdeps/ieee754/ieee754.h +++ b/sysdeps/ieee754/ieee754.h @@ -16,8 +16,15 @@ write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#ifndef _IEEE754_H + +#define _IEEE754_H 1 +#include <features.h> + #include <endian.h> +__BEGIN_DECLS + union ieee754_float { float f; @@ -153,3 +160,7 @@ union ieee854_long_double }; #define IEEE854_LONG_DOUBLE_BIAS 0x3fff + +__END_DECLS + +#endif /* ieee754.h */ diff --git a/sysdeps/libm-i387/e_exp.S b/sysdeps/libm-i387/e_exp.S index 4505ebd309..3c9d45b746 100644 --- a/sysdeps/libm-i387/e_exp.S +++ b/sysdeps/libm-i387/e_exp.S @@ -5,7 +5,7 @@ #include <machine/asm.h> -RCSID("$NetBSD: e_exp.S,v 1.4 1995/05/08 23:47:04 jtc Exp $") +RCSID("$NetBSD: e_exp.S,v 1.7 1996/07/03 17:31:28 jtc Exp $") /* e^x = 2^(x * log2(e)) */ ENTRY(__ieee754_exp) @@ -16,16 +16,16 @@ ENTRY(__ieee754_exp) -- drepper@cygnus.com. */ fxam /* Is NaN or +-Inf? */ fstsw %ax - sahf - jnc .LnoInfNaN /* No, jump. */ - jp .LisInf /* Is +-Inf, jump. */ -.LnoInfNaN: + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ fldl2e fmulp /* x * log2(e) */ fstl %st(1) frndint /* int(x * log2(e)) */ - fstl %st(2) - fsubrp /* fract(x * log2(e)) */ + fsub %st,%st(1) /* fract(x * log2(e)) */ + fxch f2xm1 /* 2^(fract(x * log2(e))) - 1 */ fld1 faddp /* 2^(fract(x * log2(e))) */ @@ -33,9 +33,8 @@ ENTRY(__ieee754_exp) fstp %st(1) ret -.LisInf: - andb $2, %ah /* Test sign. */ - jz .LpInf /* If positive, jump. */ +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ fldz /* Set result to 0. */ -.LpInf: ret +2: ret END (__ieee754_exp) diff --git a/sysdeps/libm-i387/e_expf.S b/sysdeps/libm-i387/e_expf.S index c8233eba50..cd8dc7f0b9 100644 --- a/sysdeps/libm-i387/e_expf.S +++ b/sysdeps/libm-i387/e_expf.S @@ -17,16 +17,16 @@ ENTRY(__ieee754_expf) -- drepper@cygnus.com. */ fxam /* Is NaN or +-Inf? */ fstsw %ax - sahf - jnc .LnoInfNaN /* No, jump. */ - jp .LisInf /* Is +-Inf, jump. */ -.LnoInfNaN: + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ fldl2e fmulp /* x * log2(e) */ fstl %st(1) frndint /* int(x * log2(e)) */ - fstl %st(2) - fsubrp /* fract(x * log2(e)) */ + fsub %st,%st(1) /* fract(x * log2(e)) */ + fxch f2xm1 /* 2^(fract(x * log2(e))) - 1 */ fld1 faddp /* 2^(fract(x * log2(e))) */ @@ -34,9 +34,8 @@ ENTRY(__ieee754_expf) fstp %st(1) ret -.LisInf: - andb $2, %ah /* Test sign. */ - jz .LpInf /* If positive, jump. */ +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ fldz /* Set result to 0. */ -.LpInf: ret +2: ret END (__ieee754_expf) diff --git a/sysdeps/libm-i387/e_expl.S b/sysdeps/libm-i387/e_expl.S index 9103a923aa..e398a90ef6 100644 --- a/sysdeps/libm-i387/e_expl.S +++ b/sysdeps/libm-i387/e_expl.S @@ -18,16 +18,16 @@ ENTRY(__ieee754_expl) -- drepper@cygnus.com. */ fxam /* Is NaN or +-Inf? */ fstsw %ax - sahf - jnc .LnoInfNaN /* No, jump. */ - jp .LisInf /* Is +-Inf, jump. */ -.LnoInfNaN: + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ fldl2e fmulp /* x * log2(e) */ fstl %st(1) frndint /* int(x * log2(e)) */ - fstl %st(2) - fsubrp /* fract(x * log2(e)) */ + fsub %st,%st(1) /* fract(x * log2(e)) */ + fxch f2xm1 /* 2^(fract(x * log2(e))) - 1 */ fld1 faddp /* 2^(fract(x * log2(e))) */ @@ -35,9 +35,8 @@ ENTRY(__ieee754_expl) fstp %st(1) ret -.LisInf: - andb $2, %ah /* Test sign. */ - jz .LpInf /* If positive, jump. */ +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ fldz /* Set result to 0. */ -.LpInf: ret +2: ret END (__ieee754_expl) diff --git a/sysdeps/libm-i387/e_log.S b/sysdeps/libm-i387/e_log.S index 43d430332c..e7f567d950 100644 --- a/sysdeps/libm-i387/e_log.S +++ b/sysdeps/libm-i387/e_log.S @@ -53,7 +53,7 @@ ENTRY(__ieee754_log) fyl2xp1 // log(x) ret -2: fincstp // x : log(2) +2: fstp %st(0) // x : log(2) fyl2x // log(x) ret END (__ieee754_log) diff --git a/sysdeps/libm-i387/e_log10.S b/sysdeps/libm-i387/e_log10.S index 95efc8f7a2..ecb691b905 100644 --- a/sysdeps/libm-i387/e_log10.S +++ b/sysdeps/libm-i387/e_log10.S @@ -53,7 +53,7 @@ ENTRY(__ieee754_log10) fyl2xp1 // log10(x) ret -2: fincstp // x : log10(2) +2: fstp %st(0) // x : log10(2) fyl2x // log10(x) ret END (__ieee754_log10) diff --git a/sysdeps/libm-i387/e_log10f.S b/sysdeps/libm-i387/e_log10f.S index bddb5011bc..aac58d0293 100644 --- a/sysdeps/libm-i387/e_log10f.S +++ b/sysdeps/libm-i387/e_log10f.S @@ -54,7 +54,7 @@ ENTRY(__ieee754_log10f) fyl2xp1 // log10(x) ret -2: fincstp // x : log10(2) +2: fstp %st(0) // x : log10(2) fyl2x // log10(x) ret END (__ieee754_log10f) diff --git a/sysdeps/libm-i387/e_log10l.S b/sysdeps/libm-i387/e_log10l.S index f0b8f325a7..4f51818bdd 100644 --- a/sysdeps/libm-i387/e_log10l.S +++ b/sysdeps/libm-i387/e_log10l.S @@ -55,7 +55,7 @@ ENTRY(__ieee754_log10l) fyl2xp1 // log10(x) ret -2: fincstp // x : log10(2) +2: fstp %st(0) // x : log10(2) fyl2x // log10(x) ret END(__ieee754_log10l) diff --git a/sysdeps/libm-i387/e_logf.S b/sysdeps/libm-i387/e_logf.S index b79f478cdf..4459b7fc08 100644 --- a/sysdeps/libm-i387/e_logf.S +++ b/sysdeps/libm-i387/e_logf.S @@ -54,7 +54,7 @@ ENTRY(__ieee754_logf) fyl2xp1 // log(x) ret -2: fincstp // x : log(2) +2: fstp %st(0) // x : log(2) fyl2x // log(x) ret END (__ieee754_logf) diff --git a/sysdeps/libm-i387/e_logl.S b/sysdeps/libm-i387/e_logl.S index 1657dfd5fb..08447a27e7 100644 --- a/sysdeps/libm-i387/e_logl.S +++ b/sysdeps/libm-i387/e_logl.S @@ -54,7 +54,7 @@ ENTRY(__ieee754_logl) fyl2xp1 // log(x) ret -2: fincstp // x : log(2) +2: fstp %st(0) // x : log(2) fyl2x // log(x) ret END (__ieee754_logl) diff --git a/sysdeps/libm-i387/e_pow.S b/sysdeps/libm-i387/e_pow.S new file mode 100644 index 0000000000..f6c7562d9c --- /dev/null +++ b/sysdeps/libm-i387/e_pow.S @@ -0,0 +1,120 @@ +/* ix87 specific implementation of pow function. + Copyright (C) 1996 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <machine/asm.h> + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%ecx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_pow) + fldl 4(%esp) // x + fldl 12(%esp) // y : x + +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + subl $8,%esp + + /* First see whether `y' is a natural number. In this case we + can use a more precise algorithm. */ + fld %st // y : y : x + fistpll (%esp) // y : x + fildll (%esp) // int(y) : y : x + fucomp %st(1) // y : x + fnstsw + sahf + jne 2f + + /* OK, we have an integer value for y. */ + ftst // y : x + fstp %st(0) // x + fnstsw + sahf + popl %eax + popl %edx + jnc 4f // y >= 0, jump + fdivrl MO(one) // 1/x (now referred to as x) + negl %eax + adcl $0, %edx + negl %edx +4: fldl MO(one) // 1 : x + fxch + +6: shrdl $1, %edx, %eax + jnc 5f + fxch + fmul %st(1) // x : ST*x + fxch +5: fmul %st(0), %st // x*x : ST*x + movl %eax, %ecx + orl %edx, %ecx + jnz 6b + fstp %st(0) // ST*x + ret + + .align ALIGNARG(4) +2: /* y is a real number. */ + fxch // x : y + fldl MO(one) // 1.0 : x : y + fld %st(1) // x : 1.0 : x : y + fsub %st(1) // x-1 : 1.0 : x : y + fabs // |x-1| : 1.0 : x : y + fcompl MO(limit) // 1.0 : x : y + fnstsw + fxch // x : 1.0 : y + sahf + ja 7f + fsub %st(1) // x-1 : 1.0 : y + fyl2xp1 // log2(x) : y + jmp 8f + +7: fyl2x // log2(x) : y +8: fmul %st(1) // y*log2(x) : y + fst %st(1) // y*log2(x) : y*log2(x) + frndint // int(y*log2(x)) : y*log2(x) + fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x)) + fxch // fract(y*log2(x)) : int(y*log2(x)) + f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x)) + faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x)) + fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x)) + addl $8, %esp + fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x)) + ret +END(__ieee754_pow) diff --git a/sysdeps/libm-i387/e_powf.S b/sysdeps/libm-i387/e_powf.S new file mode 100644 index 0000000000..9d6bc510b6 --- /dev/null +++ b/sysdeps/libm-i387/e_powf.S @@ -0,0 +1,120 @@ +/* ix87 specific implementation of pow function. + Copyright (C) 1996 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <machine/asm.h> + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%ecx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_powf) + flds 4(%esp) // x + flds 8(%esp) // y : x + +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + subl $8,%esp + + /* First see whether `y' is a natural number. In this case we + can use a more precise algorithm. */ + fld %st // y : y : x + fistpll (%esp) // y : x + fildll (%esp) // int(y) : y : x + fucomp %st(1) // y : x + fnstsw + sahf + jne 2f + + /* OK, we have an integer value for y. */ + ftst // y : x + fstp %st(0) // x + fnstsw + sahf + popl %eax + popl %edx + jnc 4f // y >= 0, jump + fdivrl MO(one) // 1/x (now referred to as x) + negl %eax + adcl $0, %edx + negl %edx +4: fldl MO(one) // 1 : x + fxch + +6: shrdl $1, %edx, %eax + jnc 5f + fxch + fmul %st(1) // x : ST*x + fxch +5: fmul %st(0), %st // x*x : ST*x + movl %eax, %ecx + orl %edx, %ecx + jnz 6b + fstp %st(0) // ST*x + ret + + .align ALIGNARG(4) +2: /* y is a real number. */ + fxch // x : y + fldl MO(one) // 1.0 : x : y + fld %st(1) // x : 1.0 : x : y + fsub %st(1) // x-1 : 1.0 : x : y + fabs // |x-1| : 1.0 : x : y + fcompl MO(limit) // 1.0 : x : y + fnstsw + fxch // x : 1.0 : y + sahf + ja 7f + fsub %st(1) // x-1 : 1.0 : y + fyl2xp1 // log2(x) : y + jmp 8f + +7: fyl2x // log2(x) : y +8: fmul %st(1) // y*log2(x) : y + fst %st(1) // y*log2(x) : y*log2(x) + frndint // int(y*log2(x)) : y*log2(x) + fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x)) + fxch // fract(y*log2(x)) : int(y*log2(x)) + f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x)) + faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x)) + fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x)) + addl $8, %esp + fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x)) + ret +END(__ieee754_powf) diff --git a/sysdeps/libm-i387/e_powl.S b/sysdeps/libm-i387/e_powl.S new file mode 100644 index 0000000000..ea445804f0 --- /dev/null +++ b/sysdeps/libm-i387/e_powl.S @@ -0,0 +1,120 @@ +/* ix87 specific implementation of pow function. + Copyright (C) 1996 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <machine/asm.h> + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%ecx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_powl) + fldt 4(%esp) // x + fldt 16(%esp) // y : x + +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + subl $8,%esp + + /* First see whether `y' is a natural number. In this case we + can use a more precise algorithm. */ + fld %st // y : y : x + fistpll (%esp) // y : x + fildll (%esp) // int(y) : y : x + fucomp %st(1) // y : x + fnstsw + sahf + jne 2f + + /* OK, we have an integer value for y. */ + ftst // y : x + fstp %st(0) // x + fnstsw + sahf + popl %eax + popl %edx + jnc 4f // y >= 0, jump + fdivrl MO(one) // 1/x (now referred to as x) + negl %eax + adcl $0, %edx + negl %edx +4: fldl MO(one) // 1 : x + fxch + +6: shrdl $1, %edx, %eax + jnc 5f + fxch + fmul %st(1) // x : ST*x + fxch +5: fmul %st(0), %st // x*x : ST*x + movl %eax, %ecx + orl %edx, %ecx + jnz 6b + fstp %st(0) // ST*x + ret + + .align ALIGNARG(4) +2: /* y is a real number. */ + fxch // x : y + fldl MO(one) // 1.0 : x : y + fld %st(1) // x : 1.0 : x : y + fsub %st(1) // x-1 : 1.0 : x : y + fabs // |x-1| : 1.0 : x : y + fcompl MO(limit) // 1.0 : x : y + fnstsw + fxch // x : 1.0 : y + sahf + ja 7f + fsub %st(1) // x-1 : 1.0 : y + fyl2xp1 // log2(x) : y + jmp 8f + +7: fyl2x // log2(x) : y +8: fmul %st(1) // y*log2(x) : y + fst %st(1) // y*log2(x) : y*log2(x) + frndint // int(y*log2(x)) : y*log2(x) + fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x)) + fxch // fract(y*log2(x)) : int(y*log2(x)) + f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x)) + faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x)) + fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x)) + addl $8, %esp + fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x)) + ret +END(__ieee754_powl) diff --git a/sysdeps/libm-i387/s_asinh.S b/sysdeps/libm-i387/s_asinh.S index 9b58eaa424..a71d0b66f4 100644 --- a/sysdeps/libm-i387/s_asinh.S +++ b/sysdeps/libm-i387/s_asinh.S @@ -101,7 +101,7 @@ ENTRY(__asinh) fchs // x 4: fld %st // x : x faddl MO(huge) // huge+x : x - fincstp // x + fstp %st(0) // x ret // |x| > 2^28 => y = sign(x) * (log(|x|) + log(2)) diff --git a/sysdeps/libm-i387/s_asinhf.S b/sysdeps/libm-i387/s_asinhf.S index d0a99a9461..f47bc9e0ac 100644 --- a/sysdeps/libm-i387/s_asinhf.S +++ b/sysdeps/libm-i387/s_asinhf.S @@ -101,7 +101,7 @@ ENTRY(__asinhf) fchs // x 4: fld %st // x : x faddl MO(huge) // huge+x : x - fincstp // x + fstp %st(0) // x ret // |x| > 2^14 => y = sign(x) * (log(|x|) + log(2)) diff --git a/sysdeps/libm-i387/s_asinhl.S b/sysdeps/libm-i387/s_asinhl.S index 28f76c2649..e112244908 100644 --- a/sysdeps/libm-i387/s_asinhl.S +++ b/sysdeps/libm-i387/s_asinhl.S @@ -109,7 +109,7 @@ ENTRY(__asinhl) 4: fld %st // x : x fldt MO(huge) // huge : x : x faddp // huge+x : x - fincstp // x + fstp %st(0) // x ret // |x| > 2^34 => y = sign(x) * (log(|x|) + log(2)) diff --git a/sysdeps/libm-i387/s_ceil.S b/sysdeps/libm-i387/s_ceil.S index 309fe882f8..b0159128aa 100644 --- a/sysdeps/libm-i387/s_ceil.S +++ b/sysdeps/libm-i387/s_ceil.S @@ -8,27 +8,25 @@ RCSID("$NetBSD: s_ceil.S,v 1.4 1995/05/08 23:52:13 jtc Exp $") ENTRY(__ceil) - pushl %ebp - movl %esp,%ebp + fldl 4(%esp) subl $8,%esp - fstcw -4(%ebp) /* store fpu control word */ + fstcw 4(%esp) /* store fpu control word */ /* We use here %edx although only the low 1 bits are defined. But none of the operations should care and they are faster than the 16 bit operations. */ - movl -4(%ebp),%edx - orl $0x0800,%edx /* round towards +oo */ + movl $0x0800,%edx /* round towards +oo */ + orl 4(%esp),%edx andl $0xfbff,%edx - movl %edx,-8(%ebp) - fldcw -8(%ebp) /* load modified control word */ + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ - fldl 8(%ebp); /* round */ - frndint + frndint /* round */ - fldcw -4(%ebp) /* restore original control word */ + fldcw 4(%esp) /* restore original control word */ - leave + addl $8,%esp ret END (__ceil) weak_alias (__ceil, ceil) diff --git a/sysdeps/libm-i387/s_ceilf.S b/sysdeps/libm-i387/s_ceilf.S index 768d117021..352d40d7ce 100644 --- a/sysdeps/libm-i387/s_ceilf.S +++ b/sysdeps/libm-i387/s_ceilf.S @@ -8,27 +8,25 @@ RCSID("$NetBSD: s_ceilf.S,v 1.3 1995/05/08 23:52:44 jtc Exp $") ENTRY(__ceilf) - pushl %ebp - movl %esp,%ebp + flds 4(%esp) subl $8,%esp - fstcw -4(%ebp) /* store fpu control word */ + fstcw 4(%esp) /* store fpu control word */ /* We use here %edx although only the low 1 bits are defined. But none of the operations should care and they are faster than the 16 bit operations. */ - movl -4(%ebp),%edx - orl $0x0800,%edx /* round towards +oo */ + movl $0x0800,%edx /* round towards +oo */ + orl 4(%esp),%edx andl $0xfbff,%edx - movl %edx,-8(%ebp) - fldcw -8(%ebp) /* load modified control word */ + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ - flds 8(%ebp); /* round */ - frndint + frndint /* round */ - fldcw -4(%ebp) /* restore original control word */ + fldcw 4(%esp) /* restore original control word */ - leave + addl $8,%esp ret END (__ceilf) weak_alias (__ceilf, ceilf) diff --git a/sysdeps/libm-i387/s_ceill.S b/sysdeps/libm-i387/s_ceill.S index dcf8719af3..0128966ebe 100644 --- a/sysdeps/libm-i387/s_ceill.S +++ b/sysdeps/libm-i387/s_ceill.S @@ -9,27 +9,25 @@ RCSID("$NetBSD: $") ENTRY(__ceill) - pushl %ebp - movl %esp,%ebp + fldt 4(%esp) subl $8,%esp - fstcw -4(%ebp) /* store fpu control word */ + fstcw 4(%esp) /* store fpu control word */ /* We use here %edx although only the low 1 bits are defined. But none of the operations should care and they are faster than the 16 bit operations. */ - movl -4(%ebp),%edx - orl $0x0800,%edx /* round towards +oo */ + movl $0x0800,%edx /* round towards +oo */ + orl 4(%esp),%edx andl $0xfbff,%edx - movl %edx,-8(%ebp) - fldcw -8(%ebp) /* load modified control word */ + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ - fldt 8(%ebp); /* round */ - frndint + frndint /* round */ - fldcw -4(%ebp) /* restore original control word */ + fldcw 4(%esp) /* restore original control word */ - leave + addl $8,%esp ret END (__ceill) weak_alias (__ceill, ceill) diff --git a/sysdeps/libm-i387/s_expm1.S b/sysdeps/libm-i387/s_expm1.S new file mode 100644 index 0000000000..2685d6e770 --- /dev/null +++ b/sysdeps/libm-i387/s_expm1.S @@ -0,0 +1,83 @@ +/* ix87 specific implementation of exp(x)-1. + Copyright (C) 1996 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996. + Based on code by John C. Bowman <bowman@ipp-garching.mpg.de>. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + + /* Using: e^x - 1 = 2^(x * log2(e)) - 1 */ + +#include <machine/asm.h> + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(zero,@object) +zero: .double 0.0 + ASM_SIZE_DIRECTIVE(zero) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + ASM_TYPE_DIRECTIVE(l2e,@object) +l2e: .tfloat 1.442695040888963407359924681002 + ASM_SIZE_DIRECTIVE(l2e) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__expm1) + fldl 4(%esp) // x + fxam // Is NaN or +-Inf? + fstsw %ax +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + movb $0x45, %ch + andb %ah, %ch + cmpb $0x05, %ch + je 2f // Is +-Inf, jump. + + fldt MO(l2e) // log2(e) : x + fmulp // log2(e)*x + fld %st // log2(e)*x : log2(e)*x + frndint // int(log2(e)*x) : log2(e)*x + fsub %st, %st(1) // int(log2(e)*x) : fract(log2(e)*x) + fxch // fract(log2(e)*x) : int(log2(e)*x) + f2xm1 // 2^fract(log2(e)*x)-1 : int(log2(e)*x) + fscale // 2^(log2(e)*x)-2^int(log2(e)*x) : int(log2(e)*x) + fxch // int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fldl MO(one) // 1 : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fscale // 2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + faddp // 2^(log2(e)*x) + fsubl MO(one) + ret + +2: testl $0x200, %eax // Test sign. + jz 3f // If positive, jump. + fldl MO(zero) // Set result to 0. +3: ret +END(__expm1) +weak_alias (__expm1, expm1) diff --git a/sysdeps/libm-i387/s_expm1f.S b/sysdeps/libm-i387/s_expm1f.S new file mode 100644 index 0000000000..85c3c849f1 --- /dev/null +++ b/sysdeps/libm-i387/s_expm1f.S @@ -0,0 +1,83 @@ +/* ix87 specific implementation of exp(x)-1. + Copyright (C) 1996 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996. + Based on code by John C. Bowman <bowman@ipp-garching.mpg.de>. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + + /* Using: e^x - 1 = 2^(x * log2(e)) - 1 */ + +#include <machine/asm.h> + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(zero,@object) +zero: .double 0.0 + ASM_SIZE_DIRECTIVE(zero) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + ASM_TYPE_DIRECTIVE(l2e,@object) +l2e: .tfloat 1.442695040888963407359924681002 + ASM_SIZE_DIRECTIVE(l2e) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__expm1f) + flds 4(%esp) // x + fxam // Is NaN or +-Inf? + fstsw %ax +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + movb $0x45, %ch + andb %ah, %ch + cmpb $0x05, %ch + je 2f // Is +-Inf, jump. + + fldt MO(l2e) // log2(e) : x + fmulp // log2(e)*x + fld %st // log2(e)*x : log2(e)*x + frndint // int(log2(e)*x) : log2(e)*x + fsub %st, %st(1) // int(log2(e)*x) : fract(log2(e)*x) + fxch // fract(log2(e)*x) : int(log2(e)*x) + f2xm1 // 2^fract(log2(e)*x)-1 : int(log2(e)*x) + fscale // 2^(log2(e)*x)-2^int(log2(e)*x) : int(log2(e)*x) + fxch // int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fldl MO(one) // 1 : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fscale // 2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + faddp // 2^(log2(e)*x) + fsubl MO(one) + ret + +2: testl $0x200, %eax // Test sign. + jz 3f // If positive, jump. + fldl MO(zero) // Set result to 0. +3: ret +END(__expm1f) +weak_alias (__expm1f, expm1f) diff --git a/sysdeps/libm-i387/s_expm1l.S b/sysdeps/libm-i387/s_expm1l.S new file mode 100644 index 0000000000..45aa753094 --- /dev/null +++ b/sysdeps/libm-i387/s_expm1l.S @@ -0,0 +1,83 @@ +/* ix87 specific implementation of exp(x)-1. + Copyright (C) 1996 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996. + Based on code by John C. Bowman <bowman@ipp-garching.mpg.de>. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + + /* Using: e^x - 1 = 2^(x * log2(e)) - 1 */ + +#include <machine/asm.h> + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(zero,@object) +zero: .double 0.0 + ASM_SIZE_DIRECTIVE(zero) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + ASM_TYPE_DIRECTIVE(l2e,@object) +l2e: .tfloat 1.442695040888963407359924681002 + ASM_SIZE_DIRECTIVE(l2e) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__expm1l) + fldt 4(%esp) // x + fxam // Is NaN or +-Inf? + fstsw %ax +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + movb $0x45, %ch + andb %ah, %ch + cmpb $0x05, %ch + je 2f // Is +-Inf, jump. + + fldt MO(l2e) // log2(e) : x + fmulp // log2(e)*x + fld %st // log2(e)*x : log2(e)*x + frndint // int(log2(e)*x) : log2(e)*x + fsub %st, %st(1) // int(log2(e)*x) : fract(log2(e)*x) + fxch // fract(log2(e)*x) : int(log2(e)*x) + f2xm1 // 2^fract(log2(e)*x)-1 : int(log2(e)*x) + fscale // 2^(log2(e)*x)-2^int(log2(e)*x) : int(log2(e)*x) + fxch // int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fldl MO(one) // 1 : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fscale // 2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + faddp // 2^(log2(e)*x) + fsubl MO(one) + ret + +2: testl $0x200, %eax // Test sign. + jz 3f // If positive, jump. + fldl MO(zero) // Set result to 0. +3: ret +END(__expm1l) +weak_alias (__expm1l, expm1l) diff --git a/sysdeps/libm-i387/s_floor.S b/sysdeps/libm-i387/s_floor.S index 561f83d98a..20a8660424 100644 --- a/sysdeps/libm-i387/s_floor.S +++ b/sysdeps/libm-i387/s_floor.S @@ -8,23 +8,25 @@ RCSID("$NetBSD: s_floor.S,v 1.4 1995/05/09 00:01:59 jtc Exp $") ENTRY(__floor) - pushl %ebp - movl %esp,%ebp + fldl 4(%esp) subl $8,%esp - fstcw -4(%ebp) /* store fpu control word */ - movw -4(%ebp),%dx - orw $0x0400,%dx /* round towards -oo */ - andw $0xf7ff,%dx - movw %dx,-8(%ebp) - fldcw -8(%ebp) /* load modified control word */ + fstcw 4(%esp) /* store fpu control word */ - fldl 8(%ebp); /* round */ - frndint + /* We use here %edx although only the low 1 bits are defined. + But none of the operations should care and they are faster + than the 16 bit operations. */ + movl $0x400,%edx /* round towards -oo */ + orl 4(%esp),%edx + andl $0xf7ff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ - fldcw -4(%ebp) /* restore original control word */ + frndint /* round */ - leave + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp ret END (__floor) weak_alias (__floor, floor) diff --git a/sysdeps/libm-i387/s_floorf.S b/sysdeps/libm-i387/s_floorf.S index 059f4f3a15..eca93a2aa3 100644 --- a/sysdeps/libm-i387/s_floorf.S +++ b/sysdeps/libm-i387/s_floorf.S @@ -8,23 +8,25 @@ RCSID("$NetBSD: s_floorf.S,v 1.3 1995/05/09 00:04:32 jtc Exp $") ENTRY(__floorf) - pushl %ebp - movl %esp,%ebp + flds 4(%esp) subl $8,%esp - fstcw -4(%ebp) /* store fpu control word */ - movw -4(%ebp),%dx - orw $0x0400,%dx /* round towards -oo */ - andw $0xf7ff,%dx - movw %dx,-8(%ebp) - fldcw -8(%ebp) /* load modified control word */ + fstcw 4(%esp) /* store fpu control word */ - flds 8(%ebp); /* round */ - frndint + /* We use here %edx although only the low 1 bits are defined. + But none of the operations should care and they are faster + than the 16 bit operations. */ + movl $0x400,%edx /* round towards -oo */ + orl 4(%esp),%edx + andl $0xf7ff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ - fldcw -4(%ebp) /* restore original control word */ + frndint /* round */ - leave + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp ret END (__floorf) weak_alias (__floorf, floorf) diff --git a/sysdeps/libm-i387/s_floorl.S b/sysdeps/libm-i387/s_floorl.S index 7b7cff777a..c2bf091d73 100644 --- a/sysdeps/libm-i387/s_floorl.S +++ b/sysdeps/libm-i387/s_floorl.S @@ -9,23 +9,25 @@ RCSID("$NetBSD: $") ENTRY(__floorl) - pushl %ebp - movl %esp,%ebp + fldt 4(%esp) subl $8,%esp - fstcw -4(%ebp) /* store fpu control word */ - movw -4(%ebp),%dx - orw $0x0400,%dx /* round towards -oo */ - andw $0xf7ff,%dx - movw %dx,-8(%ebp) - fldcw -8(%ebp) /* load modified control word */ + fstcw 4(%esp) /* store fpu control word */ - fldt 8(%ebp) /* round */ - frndint + /* We use here %edx although only the low 1 bits are defined. + But none of the operations should care and they are faster + than the 16 bit operations. */ + movl $0x400,%edx /* round towards -oo */ + orl 4(%esp),%edx + andl $0xf7ff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ - fldcw -4(%ebp) /* restore original control word */ + frndint /* round */ - leave + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp ret END (__floorl) weak_alias (__floorl, floorl) diff --git a/sysdeps/libm-i387/s_significandl.S b/sysdeps/libm-i387/s_significandl.S index 1d5f4a3f6a..e3a69cba73 100644 --- a/sysdeps/libm-i387/s_significandl.S +++ b/sysdeps/libm-i387/s_significandl.S @@ -11,7 +11,7 @@ RCSID("$NetBSD: $") ENTRY(__significandl) fldt 4(%esp) fxtract - fincstp + fstp %st(0) ret END (__significandl) weak_alias (__significandl, significandl) diff --git a/sysdeps/libm-ieee754/e_coshf.c b/sysdeps/libm-ieee754/e_coshf.c index e6f4ed6f37..223fbeea20 100644 --- a/sysdeps/libm-ieee754/e_coshf.c +++ b/sysdeps/libm-ieee754/e_coshf.c @@ -8,20 +8,20 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_coshf.c,v 1.5 1995/05/10 20:45:01 jtc Exp $"; +static char rcsid[] = "$NetBSD: e_coshf.c,v 1.6 1996/04/08 15:43:41 phil Exp $"; #endif #include "math.h" #include "math_private.h" #ifdef __STDC__ -static const volatile float huge = 1.0e30; +static const float huge = 1.0e30; static const float one = 1.0, half=0.5; #else static float one = 1.0, half=0.5, huge = 1.0e30; @@ -33,7 +33,7 @@ static float one = 1.0, half=0.5, huge = 1.0e30; float __ieee754_coshf(x) float x; #endif -{ +{ float t,w; int32_t ix; @@ -41,7 +41,7 @@ static float one = 1.0, half=0.5, huge = 1.0e30; ix &= 0x7fffffff; /* x is INF or NaN */ - if(ix>=0x7f800000) return x*x; + if(ix>=0x7f800000) return x*x; /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ if(ix<0x3eb17218) { diff --git a/sysdeps/libm-ieee754/e_expf.c b/sysdeps/libm-ieee754/e_expf.c index 80af17918f..fbf2691bf9 100644 --- a/sysdeps/libm-ieee754/e_expf.c +++ b/sysdeps/libm-ieee754/e_expf.c @@ -8,19 +8,19 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_expf.c,v 1.5 1995/05/10 20:45:05 jtc Exp $"; +static char rcsid[] = "$NetBSD: e_expf.c,v 1.6 1996/04/08 15:43:43 phil Exp $"; #endif #include "math.h" #include "math_private.h" -static const volatile float huge = 1.0e+30; +static const float huge = 1.0e+30; #ifdef __STDC__ static const float @@ -69,7 +69,7 @@ P5 = 4.1381369442e-08; /* 0x3331bb4c */ } /* argument reduction */ - if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ + if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb; } else { @@ -79,7 +79,7 @@ P5 = 4.1381369442e-08; /* 0x3331bb4c */ lo = t*ln2LO[0]; } x = hi - lo; - } + } else if(hx < 0x31800000) { /* when |x|<2**-28 */ if(huge+x>one) return one+x;/* trigger inexact */ } @@ -88,7 +88,7 @@ P5 = 4.1381369442e-08; /* 0x3331bb4c */ /* x is now in primary range */ t = x*x; c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); - if(k==0) return one-((x*c)/(c-(float)2.0)-x); + if(k==0) return one-((x*c)/(c-(float)2.0)-x); else y = one-((lo-(x*c)/((float)2.0-c))-hi); if(k >= -125) { u_int32_t hy; diff --git a/sysdeps/libm-ieee754/e_hypot.c b/sysdeps/libm-ieee754/e_hypot.c index 24c8ae452f..76a77ec33a 100644 --- a/sysdeps/libm-ieee754/e_hypot.c +++ b/sysdeps/libm-ieee754/e_hypot.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -16,12 +16,12 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; /* __ieee754_hypot(x,y) * - * Method : - * If (assume round-to-nearest) z=x*x+y*y - * has error less than sqrt(2)/2 ulp, than + * Method : + * If (assume round-to-nearest) z=x*x+y*y + * has error less than sqrt(2)/2 ulp, than * sqrt(z) has error less than 1 ulp (exercise). * - * So, compute sqrt(x*x+y*y) with some care as + * So, compute sqrt(x*x+y*y) with some care as * follows to get the error below 1 ulp: * * Assume x>y>0; @@ -31,10 +31,10 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; * where x1 = x with lower 32 bits cleared, x2 = x-x1; else * 2. if x <= 2y use * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y)) - * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, + * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, * y1= y with lower 32 bits chopped, y2 = y-y1. - * - * NOTE: scaling may be necessary if some argument is too + * + * NOTE: scaling may be necessary if some argument is too * large or too tiny * * Special cases: @@ -42,8 +42,8 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; * hypot(x,y) is NAN if x or y is NAN. * * Accuracy: - * hypot(x,y) returns sqrt(x^2+y^2) with error less - * than 1 ulps (units in the last place) + * hypot(x,y) returns sqrt(x^2+y^2) with error less + * than 1 ulps (units in the last place) */ #include "math.h" @@ -56,7 +56,7 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; double x, y; #endif { - double a=x,b=y,t1,t2,y1,y2,w; + double a,b,t1,t2,y1,y2,w; int32_t j,k,ha,hb; GET_HIGH_WORD(ha,x); @@ -84,7 +84,7 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; SET_HIGH_WORD(b,hb); } if(hb < 0x20b00000) { /* b < 2**-500 */ - if(hb <= 0x000fffff) { /* subnormal b or 0 */ + if(hb <= 0x000fffff) { /* subnormal b or 0 */ u_int32_t low; GET_LOW_WORD(low,b); if((hb|low)==0) return a; diff --git a/sysdeps/libm-ieee754/e_hypotf.c b/sysdeps/libm-ieee754/e_hypotf.c index 434cf94e8d..d6b1520cb8 100644 --- a/sysdeps/libm-ieee754/e_hypotf.c +++ b/sysdeps/libm-ieee754/e_hypotf.c @@ -8,7 +8,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -27,7 +27,7 @@ static char rcsid[] = "$NetBSD: e_hypotf.c,v 1.5 1995/05/12 04:57:30 jtc Exp $"; float x, y; #endif { - float a=x,b=y,t1,t2,y1,y2,w; + float a,b,t1,t2,y1,y2,w; int32_t j,k,ha,hb; GET_FLOAT_WORD(ha,x); @@ -52,7 +52,7 @@ static char rcsid[] = "$NetBSD: e_hypotf.c,v 1.5 1995/05/12 04:57:30 jtc Exp $"; SET_FLOAT_WORD(b,hb); } if(hb < 0x26800000) { /* b < 2**-50 */ - if(hb <= 0x007fffff) { /* subnormal b or 0 */ + if(hb <= 0x007fffff) { /* subnormal b or 0 */ if(hb==0) return a; SET_FLOAT_WORD(t1,0x3f000000); /* t1=2^126 */ b *= t1; diff --git a/sysdeps/libm-ieee754/e_hypotl.c b/sysdeps/libm-ieee754/e_hypotl.c new file mode 100644 index 0000000000..1a40c556dc --- /dev/null +++ b/sysdeps/libm-ieee754/e_hypotl.c @@ -0,0 +1,133 @@ +/* e_hypotl.c -- long double version of e_hypot.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: $"; +#endif + +/* __ieee754_hypotl(x,y) + * + * Method : + * If (assume round-to-nearest) z=x*x+y*y + * has error less than sqrt(2)/2 ulp, than + * sqrt(z) has error less than 1 ulp (exercise). + * + * So, compute sqrt(x*x+y*y) with some care as + * follows to get the error below 1 ulp: + * + * Assume x>y>0; + * (if possible, set rounding to round-to-nearest) + * 1. if x > 2y use + * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y + * where x1 = x with lower 32 bits cleared, x2 = x-x1; else + * 2. if x <= 2y use + * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y)) + * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, + * y1= y with lower 32 bits chopped, y2 = y-y1. + * + * NOTE: scaling may be necessary if some argument is too + * large or too tiny + * + * Special cases: + * hypot(x,y) is INF if x or y is +INF or -INF; else + * hypot(x,y) is NAN if x or y is NAN. + * + * Accuracy: + * hypot(x,y) returns sqrt(x^2+y^2) with error less + * than 1 ulps (units in the last place) + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + long double __ieee754_hypotl(long double x, long double y) +#else + long double __ieee754_hypotl(x,y) + long double x, y; +#endif +{ + long double a,b,t1,t2,y1,y2,w; + u_int32_t j,k,ea,eb; + + GET_LDOUBLE_EXP(ea,x); + ea &= 0x7fff; + GET_LDOUBLE_EXP(eb,y); + eb &= 0x7fff; + if(eb > ea) {a=y;b=x;j=ea; ea=eb;eb=j;} else {a=x;b=y;} + SET_LDOUBLE_EXP(a,ea); /* a <- |a| */ + SET_LDOUBLE_EXP(b,eb); /* b <- |b| */ + if((ea-eb)>0x46) {return a+b;} /* x/y > 2**70 */ + k=0; + if(ea > 0x5f3f) { /* a>2**8000 */ + if(ea == 0x7fff) { /* Inf or NaN */ + u_int32_t exp,high,low; + w = a+b; /* for sNaN */ + GET_LDOUBLE_WORDS(exp,high,low,a); + if(((high&0x7fffffff)|low)==0) w = a; + GET_LDOUBLE_WORDS(exp,high,low,b); + if(((eb^0x7fff)|(high&0x7fffffff)|low)==0) w = b; + return w; + } + /* scale a and b by 2**-9600 */ + ea -= 0x2580; eb -= 0x2580; k += 9600; + SET_LDOUBLE_EXP(a,ea); + SET_LDOUBLE_EXP(b,eb); + } + if(eb < 0x20bf) { /* b < 2**-8000 */ + if(eb == 0) { /* subnormal b or 0 */ + u_int32_t exp,high,low; + GET_LDOUBLE_WORDS(exp,high,low,b); + if((high|low)==0) return a; + SET_LDOUBLE_WORDS(t1, 0x7ffd, 0, 0); /* t1=2^16382 */ + b *= t1; + a *= t1; + k -= 16382; + } else { /* scale a and b by 2^9600 */ + ea += 0x2580; /* a *= 2^9600 */ + eb += 0x2580; /* b *= 2^9600 */ + k -= 9600; + SET_LDOUBLE_EXP(a,ea); + SET_LDOUBLE_EXP(b,eb); + } + } + /* medium size a and b */ + w = a-b; + if (w>b) { + u_int32_t high; + GET_LDOUBLE_MSW(high,a); + SET_LDOUBLE_WORDS(t1,ea,high,0); + t2 = a-t1; + w = __ieee754_sqrtl(t1*t1-(b*(-b)-t2*(a+t1))); + } else { + u_int32_t high; + GET_LDOUBLE_MSW(high,b); + a = a+a; + SET_LDOUBLE_WORDS(y1,eb,high,0); + y2 = b - y1; + GET_LDOUBLE_MSW(high,a); + SET_LDOUBLE_WORDS(t1,ea+1,high,0); + t2 = a - t1; + w = __ieee754_sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b))); + } + if(k!=0) { + u_int32_t exp; + t1 = 1.0; + GET_LDOUBLE_EXP(exp,t1); + SET_LDOUBLE_EXP(t1,exp+k); + return t1*w; + } else return w; +} diff --git a/sysdeps/libm-ieee754/e_powf.c b/sysdeps/libm-ieee754/e_powf.c index 1358555128..4798340c84 100644 --- a/sysdeps/libm-ieee754/e_powf.c +++ b/sysdeps/libm-ieee754/e_powf.c @@ -14,13 +14,13 @@ */ #if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_powf.c,v 1.6 1995/05/12 04:57:35 jtc Exp $"; +static char rcsid[] = "$NetBSD: e_powf.c,v 1.7 1996/04/08 15:43:44 phil Exp $"; #endif #include "math.h" #include "math_private.h" -static const volatile float huge = 1.0e+30, tiny = 1.0e-30; +static const float huge = 1.0e+30, tiny = 1.0e-30; #ifdef __STDC__ static const float diff --git a/sysdeps/libm-ieee754/e_sinhl.c b/sysdeps/libm-ieee754/e_sinhl.c new file mode 100644 index 0000000000..2d98562f94 --- /dev/null +++ b/sysdeps/libm-ieee754/e_sinhl.c @@ -0,0 +1,91 @@ +/* e_asinhl.c -- long double version of e_asinh.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: $"; +#endif + +/* __ieee754_sinhl(x) + * Method : + * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 + * 1. Replace x by |x| (sinhl(-x) = -sinhl(x)). + * 2. + * E + E/(E+1) + * 0 <= x <= 25 : sinhl(x) := --------------, E=expm1l(x) + * 2 + * + * 25 <= x <= lnovft : sinhl(x) := expl(x)/2 + * lnovft <= x <= ln2ovft: sinhl(x) := expl(x/2)/2 * expl(x/2) + * ln2ovft < x : sinhl(x) := x*shuge (overflow) + * + * Special cases: + * sinhl(x) is |x| if x is +INF, -INF, or NaN. + * only sinhl(0)=0 is exact for finite x. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double one = 1.0, shuge = 1.0e4931L; +#else +static long double one = 1.0, shuge = 1.0e4931L; +#endif + +#ifdef __STDC__ + long double __ieee754_sinhl(long double x) +#else + long double __ieee754_sinhl(x) + long double x; +#endif +{ + long double t,w,h; + u_int32_t jx,ix,i0,i1; + + /* Words of |x|. */ + GET_LDOUBLE_WORDS(jx,i0,i1,x); + ix = jx&0x7fff; + + /* x is INF or NaN */ + if(ix==0x7fff) return x+x; + + h = 0.5; + if (jx<0) h = -h; + /* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */ + if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x|<25 */ + if (ix<0x3fe3) /* |x|<2**-28 */ + if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */ + t = __expm1l(fabsl(x)); + if(ix<0x3fff) return h*(2.0*t-t*t/(t+one)); + return h*(t+t/(t+one)); + } + + /* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */ + if (ix < 0x400c || (ix == 0x400c && i0 < 0xb17217f7)) + return h*__ieee754_expl(fabsl(x)); + + /* |x| in [log(maxdouble), overflowthresold] */ + if (ix<0x400c || (ix == 0x400c && (i0 < 0xb174ddc0 + || (i0 == 0xb174ddc0 + && i1 <= 0x31aec0ea)))) { + w = __ieee754_expl(0.5*fabsl(x)); + t = h*w; + return t*w; + } + + /* |x| > overflowthresold, sinhl(x) overflow */ + return x*shuge; +} diff --git a/sysdeps/m68k/Makefile b/sysdeps/m68k/Makefile index a5e9064767..dc502912bf 100644 --- a/sysdeps/m68k/Makefile +++ b/sysdeps/m68k/Makefile @@ -12,9 +12,9 @@ # Library General Public License for more details. # You should have received a copy of the GNU Library General Public -# License along with the GNU C Library; see the file COPYING.LIB. If -# not, write to the Free Software Foundation, Inc., 675 Mass Ave, -# Cambridge, MA 02139, USA. +# License along with the GNU C Library; see the file COPYING.LIB. If not, +# write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, +# Boston, MA 02111-1307, USA. # The mpn functions need this. All existing 68k ports use MIT syntax. If # a new port wants to use Motorola or Sony syntax, it can redefine this @@ -30,3 +30,7 @@ CFLAGS-setjmp.c := -fno-omit-frame-pointer # The 68k `long double' is a distinct type we support. long-double-fcts = yes + +ifeq ($(subdir),elf) +CFLAGS-rtld.c += -Wno-uninitialized +endif diff --git a/sysdeps/posix/sigpause.c b/sysdeps/posix/sigpause.c index b3eec8f39a..8b820ecad9 100644 --- a/sysdeps/posix/sigpause.c +++ b/sysdeps/posix/sigpause.c @@ -34,7 +34,7 @@ __sigpause (sig_or_mask, is_sig) /* The modern X/Open implementation is requested. */ if (sigprocmask (0, NULL, &set) < 0 /* Yes, we call `sigaddset' and not `__sigaddset'. */ - || sigaddset (&set, sig_or_mask) < 0) + || sigdelset (&set, sig_or_mask) < 0) return -1; } else @@ -54,3 +54,16 @@ __sigpause (sig_or_mask, is_sig) return sigsuspend (&set); } + + +/* We have to provide a default version of this function since the + standards demand it. The version which is a bit more reasonable is + the BSD version. So make this the default. */ +int __default_sigpause __P ((int mask)); +int +__default_sigpause (mask) + int mask; +{ + return __sigpause (mask, 0); +} +weak_alias (__default_sigpause, sigpause) diff --git a/sysdeps/stub/e_hypotl.c b/sysdeps/stub/e_hypotl.c deleted file mode 100644 index a7a6e124d3..0000000000 --- a/sysdeps/stub/e_hypotl.c +++ /dev/null @@ -1,11 +0,0 @@ -#include <math.h> -#include <stdio.h> - -long double -__ieee754_hypotl (long double x, long double y) -{ - fputs ("__ieee754_hypotl not implemented\n", stderr); - return 0.0; -} - -stub_warning (__ieee754_hypotl) diff --git a/sysdeps/stub/e_sinhl.c b/sysdeps/stub/e_sinhl.c deleted file mode 100644 index 8fe11d4ea9..0000000000 --- a/sysdeps/stub/e_sinhl.c +++ /dev/null @@ -1,11 +0,0 @@ -#include <math.h> -#include <stdio.h> - -long double -__ieee754_sinhl (long double x) -{ - fputs ("__ieee754_sinhl not implemented\n", stderr); - return 0.0; -} - -stub_warning (__ieee754_sinhl) |