diff options
Diffstat (limited to 'sysdeps/i386/fpu')
151 files changed, 7920 insertions, 2 deletions
diff --git a/sysdeps/i386/fpu/Implies b/sysdeps/i386/fpu/Implies deleted file mode 100644 index de9b0b27aa..0000000000 --- a/sysdeps/i386/fpu/Implies +++ /dev/null @@ -1,2 +0,0 @@ -# For x86 machines with FPU, use the i387 port of libm by JT Conklin. -libm-i387 diff --git a/sysdeps/i386/fpu/e_acos.S b/sysdeps/i386/fpu/e_acos.S new file mode 100644 index 0000000000..b9d07b1091 --- /dev/null +++ b/sysdeps/i386/fpu/e_acos.S @@ -0,0 +1,21 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_acos.S,v 1.4 1995/05/08 23:44:37 jtc Exp $") + +/* acos = atan (sqrt(1 - x^2) / x) */ +ENTRY(__ieee754_acos) + fldl 4(%esp) /* x */ + fld %st /* x : x */ + fmul %st(0) /* x^2 : x */ + fld1 /* 1 : x^2 : x */ + fsubp /* 1 - x^2 : x */ + fsqrt /* sqrt (1 - x^2) : x */ + fxch %st(1) /* x : sqrt (1 - x^2) */ + fpatan /* atan (sqrt(1 - x^2) / x) */ + ret +END (__ieee754_acos) diff --git a/sysdeps/i386/fpu/e_acosf.S b/sysdeps/i386/fpu/e_acosf.S new file mode 100644 index 0000000000..50b13fd1bd --- /dev/null +++ b/sysdeps/i386/fpu/e_acosf.S @@ -0,0 +1,22 @@ +/* + * 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: $") + +/* acos = atan (sqrt(1 - x^2) / x) */ +ENTRY(__ieee754_acosf) + flds 4(%esp) /* x */ + fld %st + fmul %st(0) /* x^2 */ + fld1 + fsubp /* 1 - x^2 */ + fsqrt /* sqrt (1 - x^2) */ + fxch %st(1) + fpatan + ret +END (__ieee754_acosf) diff --git a/sysdeps/i386/fpu/e_acosh.S b/sysdeps/i386/fpu/e_acosh.S new file mode 100644 index 0000000000..a3397b365c --- /dev/null +++ b/sysdeps/i386/fpu/e_acosh.S @@ -0,0 +1,105 @@ +/* ix87 specific implementation of arcsinh. + 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(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_acosh) + movl 8(%esp), %ecx + cmpl $0x3ff00000, %ecx + jl 5f // < 1 => invalid + fldln2 // log(2) + fldl 4(%esp) // x : log(2) + cmpl $0x41b00000, %ecx + ja 3f // x > 2^28 +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + cmpl $0x40000000, %ecx + ja 4f // x > 2 + + // 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2)) + fsubl MO(one) // x-1 : log(2) + fld %st // x-1 : x-1 : log(2) + fmul %st(1) // (x-1)^2 : x-1 : log(2) + fadd %st(1) // x-1+(x-1)^2 : x-1 : log(2) + fadd %st(1) // 2*(x-1)+(x-1)^2 : x-1 : log(2) + fsqrt // sqrt(2*(x-1)+(x-1)^2) : x-1 : log(2) + faddp // x-1+sqrt(2*(x-1)+(x-1)^2) : log(2) + fcoml MO(limit) + fnstsw + sahf + ja 2f + fyl2xp1 // log1p(x-1+sqrt(2*(x-1)+(x-1)^2)) + ret + +2: faddl MO(one) // x+sqrt(2*(x-1)+(x-1)^2) : log(2) + fyl2x // log(x+sqrt(2*(x-1)+(x-1)^2)) + ret + + // x > 2^28 => y = log(x) + log(2) + .align ALIGNARG(4) +3: fyl2x // log(x) + fldln2 // log(2) : log(x) + faddp // log(x)+log(2) + ret + + // 2^28 > x > 2 => y = log(2*x - 1/(x+sqrt(x*x-1))) + .align ALIGNARG(4) +4: fld %st // x : x : log(2) + fadd %st, %st(1) // x : 2*x : log(2) + fld %st // x : x : 2*x : log(2) + fmul %st(1) // x^2 : x : 2*x : log(2) + fsubl MO(one) // x^2-1 : x : 2*x : log(2) + fsqrt // sqrt(x^2-1) : x : 2*x : log(2) + faddp // x+sqrt(x^2-1) : 2*x : log(2) + fdivrl MO(one) // 1/(x+sqrt(x^2-1)) : 2*x : log(2) + fsubrp // 2*x+1/(x+sqrt(x^2)-1) : log(2) + fyl2x // log(2*x+1/(x+sqrt(x^2-1))) + ret + + // x < 1 => NaN + .align ALIGNARG(4) +5: fldz + fdiv %st, %st(0) + ret +END(__ieee754_acosh) diff --git a/sysdeps/i386/fpu/e_acoshf.S b/sysdeps/i386/fpu/e_acoshf.S new file mode 100644 index 0000000000..8aa78957e2 --- /dev/null +++ b/sysdeps/i386/fpu/e_acoshf.S @@ -0,0 +1,105 @@ +/* ix87 specific implementation of arcsinh. + Copyright (C) 1996, 1997 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(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_acoshf) + movl 4(%esp), %ecx + cmpl $0x3f800000, %ecx + jl 5f // < 1 => invalid + fldln2 // log(2) + flds 4(%esp) // x : log(2) + cmpl $0x47000000, %ecx + ja 3f // x > 2^14 +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + cmpl $0x40000000, %ecx + ja 4f // x > 2 + + // 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2)) + fsubl MO(one) // x-1 : log(2) + fld %st // x-1 : x-1 : log(2) + fmul %st(1) // (x-1)^2 : x-1 : log(2) + fadd %st(1) // x-1+(x-1)^2 : x-1 : log(2) + fadd %st(1) // 2*(x-1)+(x-1)^2 : x-1 : log(2) + fsqrt // sqrt(2*(x-1)+(x-1)^2) : x-1 : log(2) + faddp // x-1+sqrt(2*(x-1)+(x-1)^2) : log(2) + fcoml MO(limit) + fnstsw + sahf + ja 2f + fyl2xp1 // log1p(x-1+sqrt(2*(x-1)+(x-1)^2)) + ret + +2: faddl MO(one) // x+sqrt(2*(x-1)+(x-1)^2) : log(2) + fyl2x // log(x+sqrt(2*(x-1)+(x-1)^2)) + ret + + // x > 2^14 => y = log(x) + log(2) + .align ALIGNARG(4) +3: fyl2x // log(x) + fldln2 // log(2) : log(x) + faddp // log(x)+log(2) + ret + + // 2^28 > x > 2 => y = log(2*x - 1/(x+sqrt(x*x-1))) + .align ALIGNARG(4) +4: fld %st // x : x : log(2) + fadd %st, %st(1) // x : 2*x : log(2) + fld %st // x : x : 2*x : log(2) + fmul %st(1) // x^2 : x : 2*x : log(2) + fsubl MO(one) // x^2-1 : x : 2*x : log(2) + fsqrt // sqrt(x^2-1) : x : 2*x : log(2) + faddp // x+sqrt(x^2-1) : 2*x : log(2) + fdivrl MO(one) // 1/(x+sqrt(x^2-1)) : 2*x : log(2) + fsubrp // 2*x+1/(x+sqrt(x^2)-1) : log(2) + fyl2x // log(2*x+1/(x+sqrt(x^2-1))) + ret + + // x < 1 => NaN + .align ALIGNARG(4) +5: fldz + fdiv %st, %st(0) + ret +END(__ieee754_acoshf) diff --git a/sysdeps/i386/fpu/e_acoshl.S b/sysdeps/i386/fpu/e_acoshl.S new file mode 100644 index 0000000000..0c81daaebe --- /dev/null +++ b/sysdeps/i386/fpu/e_acoshl.S @@ -0,0 +1,112 @@ +/* ix87 specific implementation of arcsinh. + Copyright (C) 1996, 1997 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) + /* Please note that we use double value for 1.0. This number + has an exact representation and so we don't get accuracy + problems. The advantage is that the code is simpler. */ + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_acoshl) + movl 12(%esp), %ecx + andl $0xffff, %ecx + cmpl $0x3fff, %ecx + jl 5f // < 1 => invalid + fldln2 // log(2) + fldt 4(%esp) // x : log(2) + cmpl $0x4020, %ecx + ja 3f // x > 2^34 +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + cmpl $0x4000, %ecx + ja 4f // x > 2 + + // 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2)) + fsubl MO(one) // x-1 : log(2) + fld %st // x-1 : x-1 : log(2) + fmul %st(1) // (x-1)^2 : x-1 : log(2) + fadd %st(1) // x-1+(x-1)^2 : x-1 : log(2) + fadd %st(1) // 2*(x-1)+(x-1)^2 : x-1 : log(2) + fsqrt // sqrt(2*(x-1)+(x-1)^2) : x-1 : log(2) + faddp // x-1+sqrt(2*(x-1)+(x-1)^2) : log(2) + fcoml MO(limit) + fnstsw + sahf + ja 2f + fyl2xp1 // log1p(x-1+sqrt(2*(x-1)+(x-1)^2)) + ret + +2: faddl MO(one) // x+sqrt(2*(x-1)+(x-1)^2) : log(2) + fyl2x // log(x+sqrt(2*(x-1)+(x-1)^2)) + ret + + // x > 2^34 => y = log(x) + log(2) + .align ALIGNARG(4) +3: fyl2x // log(x) + fldln2 // log(2) : log(x) + faddp // log(x)+log(2) + ret + + // 2^34 > x > 2 => y = log(2*x - 1/(x+sqrt(x*x-1))) + .align ALIGNARG(4) +4: fld %st // x : x : log(2) + fadd %st, %st(1) // x : 2*x : log(2) + fld %st // x : x : 2*x : log(2) + fmul %st(1) // x^2 : x : 2*x : log(2) + fsubl MO(one) // x^2-1 : x : 2*x : log(2) + fsqrt // sqrt(x^2-1) : x : 2*x : log(2) + faddp // x+sqrt(x^2-1) : 2*x : log(2) + fdivrl MO(one) // 1/(x+sqrt(x^2-1)) : 2*x : log(2) + fsubrp // 2*x+1/(x+sqrt(x^2)-1) : log(2) + fyl2x // log(2*x+1/(x+sqrt(x^2-1))) + ret + + // x < 1 => NaN + .align ALIGNARG(4) +5: fldz + fdiv %st, %st(0) + ret +END(__ieee754_acoshl) diff --git a/sysdeps/i386/fpu/e_acosl.S b/sysdeps/i386/fpu/e_acosl.S new file mode 100644 index 0000000000..d69f056556 --- /dev/null +++ b/sysdeps/i386/fpu/e_acosl.S @@ -0,0 +1,22 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + + +/* acosl = atanl (sqrtl(1 - x^2) / x) */ +ENTRY(__ieee754_acosl) + fldt 4(%esp) /* x */ + fld %st + fmul %st(0) /* x^2 */ + fld1 + fsubp /* 1 - x^2 */ + fsqrt /* sqrtl (1 - x^2) */ + fxch %st(1) + fpatan + ret +END (__ieee754_acosl) diff --git a/sysdeps/i386/fpu/e_asin.S b/sysdeps/i386/fpu/e_asin.S new file mode 100644 index 0000000000..945e308245 --- /dev/null +++ b/sysdeps/i386/fpu/e_asin.S @@ -0,0 +1,20 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_asin.S,v 1.4 1995/05/08 23:45:40 jtc Exp $") + +/* asin = atan (x / sqrt(1 - x^2)) */ +ENTRY(__ieee754_asin) + fldl 4(%esp) /* x */ + fld %st + fmul %st(0) /* x^2 */ + fld1 + fsubp /* 1 - x^2 */ + fsqrt /* sqrt (1 - x^2) */ + fpatan + ret +END (__ieee754_asin) diff --git a/sysdeps/i386/fpu/e_asinf.S b/sysdeps/i386/fpu/e_asinf.S new file mode 100644 index 0000000000..d450e9a740 --- /dev/null +++ b/sysdeps/i386/fpu/e_asinf.S @@ -0,0 +1,21 @@ +/* + * 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: $") + +/* asin = atan (x / sqrt(1 - x^2)) */ +ENTRY(__ieee754_asinf) + flds 4(%esp) /* x */ + fld %st + fmul %st(0) /* x^2 */ + fld1 + fsubp /* 1 - x^2 */ + fsqrt /* sqrt (1 - x^2) */ + fpatan + ret +END (__ieee754_asinf) diff --git a/sysdeps/i386/fpu/e_asinl.S b/sysdeps/i386/fpu/e_asinl.S new file mode 100644 index 0000000000..3919fbcf58 --- /dev/null +++ b/sysdeps/i386/fpu/e_asinl.S @@ -0,0 +1,22 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +/* asinl = atanl (x / sqrtl(1 - x^2)) */ +ENTRY(__ieee754_asinl) + fldt 4(%esp) /* x */ + fld %st + fmul %st(0) /* x^2 */ + fld1 + fsubp /* 1 - x^2 */ + fsqrt /* sqrt (1 - x^2) */ + fpatan + ret +END (__ieee754_asinl) diff --git a/sysdeps/i386/fpu/e_atan2.S b/sysdeps/i386/fpu/e_atan2.S new file mode 100644 index 0000000000..8df04e485e --- /dev/null +++ b/sysdeps/i386/fpu/e_atan2.S @@ -0,0 +1,15 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_atan2.S,v 1.4 1995/05/08 23:46:28 jtc Exp $") + +ENTRY(__ieee754_atan2) + fldl 4(%esp) + fldl 12(%esp) + fpatan + ret +END (__ieee754_atan2) diff --git a/sysdeps/i386/fpu/e_atan2f.S b/sysdeps/i386/fpu/e_atan2f.S new file mode 100644 index 0000000000..fc6621f183 --- /dev/null +++ b/sysdeps/i386/fpu/e_atan2f.S @@ -0,0 +1,15 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_atan2f.S,v 1.1 1995/05/08 23:35:10 jtc Exp $") + +ENTRY(__ieee754_atan2f) + flds 4(%esp) + flds 8(%esp) + fpatan + ret +END (__ieee754_atan2f) diff --git a/sysdeps/i386/fpu/e_atan2l.S b/sysdeps/i386/fpu/e_atan2l.S new file mode 100644 index 0000000000..f58eaa94a9 --- /dev/null +++ b/sysdeps/i386/fpu/e_atan2l.S @@ -0,0 +1,17 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__ieee754_atan2l) + fldt 4(%esp) + fldt 16(%esp) + fpatan + ret +END (__ieee754_atan2l) diff --git a/sysdeps/i386/fpu/e_atanh.S b/sysdeps/i386/fpu/e_atanh.S new file mode 100644 index 0000000000..231e96f57f --- /dev/null +++ b/sysdeps/i386/fpu/e_atanh.S @@ -0,0 +1,101 @@ +/* ix87 specific implementation of arctanh 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(half,@object) +half: .double 0.5 + ASM_SIZE_DIRECTIVE(half) + 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) + ASM_TYPE_DIRECTIVE(ln2_2,@object) +ln2_2: .tfloat 0.3465735902799726547086160 + ASM_SIZE_DIRECTIVE(ln2_2) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_atanh) + movl 8(%esp), %ecx + +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + + andl $0x80000000, %ecx // ECX == 0 iff X >= 0 + + fldt MO(ln2_2) // 0.5*ln2 + xorl %ecx, 8(%esp) + fldl 4(%esp) // |x| : 0.5*ln2 + fcoml MO(half) // |x| : 0.5*ln2 + fld %st // |x| : |x| : 0.5*ln2 + fnstsw // |x| : |x| : 0.5*ln2 + sahf + jae 2f + fadd %st, %st(1) // |x| : 2*|x| : 0.5*ln2 + fld %st // |x| : |x| : 2*|x| : 0.5*ln2 + fsubrl MO(one) // 1-|x| : |x| : 2*|x| : 0.5*ln2 + fxch // |x| : 1-|x| : 2*|x| : 0.5*ln2 + fmul %st(2) // 2*|x|^2 : 1-|x| : 2*|x| : 0.5*ln2 + fdivp // (2*|x|^2)/(1-|x|) : 2*|x| : 0.5*ln2 + faddp // 2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2 + fcoml MO(limit) // 2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2 + fnstsw // 2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2 + sahf + jae 4f + fyl2xp1 // 0.5*ln2*ld(1+2*|x|+(2*|x|^2)/(1-|x|)) + jecxz 3f + fchs // 0.5*ln2*ld(1+2*x+(2*x^2)/(1-x)) +3: ret + + .align ALIGNARG(4) +4: faddl MO(one) // 1+2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2 + fyl2x // 0.5*ln2*ld(1+2*|x|+(2*|x|^2)/(1-|x|)) + jecxz 3f + fchs // 0.5*ln2*ld(1+2*x+(2*x^2)/(1-x)) +3: ret + + .align ALIGNARG(4) +2: faddl MO(one) // 1+|x| : |x| : 0.5*ln2 + fxch // |x| : 1+|x| : 0.5*ln2 + fsubrl MO(one) // 1-|x| : 1+|x| : 0.5*ln2 + fdivrp // (1+|x|)/(1-|x|) : 0.5*ln2 + fyl2x // 0.5*ln2*ld((1+|x|)/(1-|x|)) + jecxz 3f + fchs // 0.5*ln2*ld((1+x)/(1-x)) +3: ret +END(__ieee754_atanh) diff --git a/sysdeps/i386/fpu/e_atanhf.S b/sysdeps/i386/fpu/e_atanhf.S new file mode 100644 index 0000000000..687d4c97fb --- /dev/null +++ b/sysdeps/i386/fpu/e_atanhf.S @@ -0,0 +1,102 @@ +/* ix87 specific implementation of arctanh 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(half,@object) +half: .double 0.5 + ASM_SIZE_DIRECTIVE(half) + 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) + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(ln2_2,@object) +ln2_2: .tfloat 0.3465735902799726547086160 + ASM_SIZE_DIRECTIVE(ln2_2) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_atanhf) + movl 4(%esp), %ecx + +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + + andl $0x80000000, %ecx // ECX == 0 iff X >= 0 + + fldt MO(ln2_2) // 0.5*ln2 + xorl %ecx, 4(%esp) + flds 4(%esp) // |x| : 0.5*ln2 + fcoml MO(half) // |x| : 0.5*ln2 + fld %st(0) // |x| : |x| : 0.5*ln2 + fnstsw // |x| : |x| : 0.5*ln2 + sahf + jae 2f + fadd %st, %st(1) // |x| : 2*|x| : 0.5*ln2 + fld %st // |x| : |x| : 2*|x| : 0.5*ln2 + fsubrl MO(one) // 1-|x| : |x| : 2*|x| : 0.5*ln2 + fxch // |x| : 1-|x| : 2*|x| : 0.5*ln2 + fmul %st(2) // 2*|x|^2 : 1-|x| : 2*|x| : 0.5*ln2 + fdivp // (2*|x|^2)/(1-|x|) : 2*|x| : 0.5*ln2 + faddp // 2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2 + fcoml MO(limit) // 2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2 + fnstsw // 2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2 + sahf + jae 4f + fyl2xp1 // 0.5*ln2*ld(1+2*|x|+(2*|x|^2)/(1-|x|)) + jecxz 3f + fchs // 0.5*ln2*ld(1+2*x+(2*x^2)/(1-x)) +3: ret + + .align ALIGNARG(4) +4: faddl MO(one) // 1+2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2 + fyl2x // 0.5*ln2*ld(1+2*|x|+(2*|x|^2)/(1-|x|)) + jecxz 3f + fchs // 0.5*ln2*ld(1+2*x+(2*x^2)/(1-x)) +3: ret + + .align ALIGNARG(4) +2: faddl MO(one) // 1+|x| : |x| : 0.5*ln2 + fxch // |x| : 1+|x| : 0.5*ln2 + fsubrl MO(one) // 1-|x| : 1+|x| : 0.5*ln2 + fdivrp // (1+|x|)/(1-|x|) : 0.5*ln2 + fyl2x // 0.5*ln2*ld((1+|x|)/(1-|x|)) + jecxz 3f + fchs // 0.5*ln2*ld((1+x)/(1-x)) +3: ret +END(__ieee754_atanhf) diff --git a/sysdeps/i386/fpu/e_atanhl.S b/sysdeps/i386/fpu/e_atanhl.S new file mode 100644 index 0000000000..8a2bd11ce4 --- /dev/null +++ b/sysdeps/i386/fpu/e_atanhl.S @@ -0,0 +1,108 @@ +/* ix87 specific implementation of arctanh 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) + /* Please note that we use double values for 0.5 and 1.0. These + numbers have exact representations and so we don't get accuracy + problems. The advantage is that the code is simpler. */ + ASM_TYPE_DIRECTIVE(half,@object) +half: .double 0.5 + ASM_SIZE_DIRECTIVE(half) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(ln2_2,@object) +ln2_2: .tfloat 0.3465735902799726547086160 + ASM_SIZE_DIRECTIVE(ln2_2) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_atanhl) + movl 12(%esp), %ecx + +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + + andl $0x8000, %ecx // ECX == 0 iff X >= 0 + + fldt MO(ln2_2) // 0.5*ln2 + xorl %ecx, 12(%esp) + fldt 4(%esp) // |x| : 0.5*ln2 + fcoml MO(half) // |x| : 0.5*ln2 + fld %st(0) // |x| : |x| : 0.5*ln2 + fnstsw // |x| : |x| : 0.5*ln2 + sahf + jae 2f + fadd %st, %st(1) // |x| : 2*|x| : 0.5*ln2 + fld %st // |x| : |x| : 2*|x| : 0.5*ln2 + fsubrl MO(one) // 1-|x| : |x| : 2*|x| : 0.5*ln2 + fxch // |x| : 1-|x| : 2*|x| : 0.5*ln2 + fmul %st(2) // 2*|x|^2 : 1-|x| : 2*|x| : 0.5*ln2 + fdivp // (2*|x|^2)/(1-|x|) : 2*|x| : 0.5*ln2 + faddp // 2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2 + fcoml MO(limit) // 2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2 + fnstsw // 2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2 + sahf + jae 4f + fyl2xp1 // 0.5*ln2*ld(1+2*|x|+(2*|x|^2)/(1-|x|)) + jecxz 3f + fchs // 0.5*ln2*ld(1+2*x+(2*x^2)/(1-x)) +3: ret + + .align ALIGNARG(4) +4: faddl MO(one) // 1+2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2 + fyl2x // 0.5*ln2*ld(1+2*|x|+(2*|x|^2)/(1-|x|)) + jecxz 3f + fchs // 0.5*ln2*ld(1+2*x+(2*x^2)/(1-x)) +3: ret + + .align ALIGNARG(4) +2: faddl MO(one) // 1+|x| : |x| : 0.5*ln2 + fxch // |x| : 1+|x| : 0.5*ln2 + fsubrl MO(one) // 1-|x| : 1+|x| : 0.5*ln2 + fdivrp // (1+|x|)/(1-|x|) : 0.5*ln2 + fyl2x // 0.5*ln2*ld((1+|x|)/(1-|x|)) + jecxz 3f + fchs // 0.5*ln2*ld((1+x)/(1-x)) +3: ret +END(__ieee754_atanhl) diff --git a/sysdeps/i386/fpu/e_exp.S b/sysdeps/i386/fpu/e_exp.S new file mode 100644 index 0000000000..4a75fa1d1c --- /dev/null +++ b/sysdeps/i386/fpu/e_exp.S @@ -0,0 +1,41 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +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) + fldl 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 + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + fldl2e + fmulp /* x * log2(e) */ + fld %st + frndint /* int(x * log2(e)) */ + fsubr %st,%st(1) /* fract(x * log2(e)) */ + fxch + f2xm1 /* 2^(fract(x * log2(e))) - 1 */ + fld1 + faddp /* 2^(fract(x * log2(e))) */ + fscale /* e^x */ + fstp %st(1) + ret + +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +2: ret +END (__ieee754_exp) diff --git a/sysdeps/i386/fpu/e_exp10.S b/sysdeps/i386/fpu/e_exp10.S new file mode 100644 index 0000000000..6bfcdbb723 --- /dev/null +++ b/sysdeps/i386/fpu/e_exp10.S @@ -0,0 +1,38 @@ +/* + * Written by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +/* 10^x = 2^(x * log2(10)) */ +ENTRY(__ieee754_exp10) + fldl 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 + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + fldl2t + fmulp /* x * log2(10) */ + fld %st + frndint /* int(x * log2(10)) */ + fsubr %st,%st(1) /* fract(x * log2(10)) */ + fxch + f2xm1 /* 2^(fract(x * log2(10))) - 1 */ + fld1 + faddp /* 2^(fract(x * log2(10))) */ + fscale /* e^x */ + fstp %st(1) + ret + +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +2: ret +END (__ieee754_exp10) diff --git a/sysdeps/i386/fpu/e_exp10f.S b/sysdeps/i386/fpu/e_exp10f.S new file mode 100644 index 0000000000..4791b99afa --- /dev/null +++ b/sysdeps/i386/fpu/e_exp10f.S @@ -0,0 +1,38 @@ +/* + * Written by Ulrich Drepper. + */ + +#include <machine/asm.h> + +/* e^x = 2^(x * log2(10)) */ +ENTRY(__ieee754_exp10f) + 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 + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + fldl2t + fmulp /* x * log2(10) */ + fld %st + frndint /* int(x * log2(10)) */ + fsubr %st,%st(1) /* fract(x * log2(10)) */ + fxch + f2xm1 /* 2^(fract(x * log2(10))) - 1 */ + fld1 + faddp /* 2^(fract(x * log2(10))) */ + fscale /* e^x */ + fstp %st(1) + ret + +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +2: ret +END (__ieee754_exp10f) diff --git a/sysdeps/i386/fpu/e_exp10l.S b/sysdeps/i386/fpu/e_exp10l.S new file mode 100644 index 0000000000..71f0da792d --- /dev/null +++ b/sysdeps/i386/fpu/e_exp10l.S @@ -0,0 +1,38 @@ +/* + * Written by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +/* e^x = 2^(x * log2l(10)) */ +ENTRY(__ieee754_exp10l) + fldt 4(%esp) +/* I added the following ugly construct because expl(+-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 + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + fldl2t + fmulp /* x * log2(10) */ + fld %st + frndint /* int(x * log2(10)) */ + fsubr %st,%st(1) /* fract(x * log2(10)) */ + fxch + f2xm1 /* 2^(fract(x * log2(10))) - 1 */ + fld1 + faddp /* 2^(fract(x * log2(10))) */ + fscale /* e^x */ + fstp %st(1) + ret + +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +2: ret +END (__ieee754_exp10l) diff --git a/sysdeps/i386/fpu/e_expf.S b/sysdeps/i386/fpu/e_expf.S new file mode 100644 index 0000000000..5fd49b89fd --- /dev/null +++ b/sysdeps/i386/fpu/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 + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + fldl2e + fmulp /* x * log2(e) */ + fld %st + frndint /* int(x * log2(e)) */ + fsubr %st,%st(1) /* fract(x * log2(e)) */ + fxch + f2xm1 /* 2^(fract(x * log2(e))) - 1 */ + fld1 + faddp /* 2^(fract(x * log2(e))) */ + fscale /* e^x */ + fstp %st(1) + ret + +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +2: ret +END (__ieee754_expf) diff --git a/sysdeps/i386/fpu/e_expl.S b/sysdeps/i386/fpu/e_expl.S new file mode 100644 index 0000000000..2bcdf58c58 --- /dev/null +++ b/sysdeps/i386/fpu/e_expl.S @@ -0,0 +1,43 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +/* e^x = 2^(x * log2l(e)) */ +ENTRY(__ieee754_expl) + fldt 4(%esp) +/* I added the following ugly construct because expl(+-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 + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + fldl2e + fmulp /* x * log2(e) */ + fld %st + frndint /* int(x * log2(e)) */ + fsubr %st,%st(1) /* fract(x * log2(e)) */ + fxch + f2xm1 /* 2^(fract(x * log2(e))) - 1 */ + fld1 + faddp /* 2^(fract(x * log2(e))) */ + fscale /* e^x */ + fstp %st(1) + ret + +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +2: ret +END (__ieee754_expl) diff --git a/sysdeps/i386/fpu/e_fmod.S b/sysdeps/i386/fpu/e_fmod.S new file mode 100644 index 0000000000..4cf6e92054 --- /dev/null +++ b/sysdeps/i386/fpu/e_fmod.S @@ -0,0 +1,19 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_fmod.S,v 1.4 1995/05/08 23:47:56 jtc Exp $") + +ENTRY(__ieee754_fmod) + fldl 12(%esp) + fldl 4(%esp) +1: fprem + fstsw %ax + sahf + jp 1b + fstp %st(1) + ret +END (__ieee754_fmod) diff --git a/sysdeps/i386/fpu/e_fmodf.S b/sysdeps/i386/fpu/e_fmodf.S new file mode 100644 index 0000000000..bbce40976d --- /dev/null +++ b/sysdeps/i386/fpu/e_fmodf.S @@ -0,0 +1,20 @@ +/* + * 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: $") + +ENTRY(__ieee754_fmodf) + flds 8(%esp) + flds 4(%esp) +1: fprem + fstsw %ax + sahf + jp 1b + fstp %st(1) + ret +END(__ieee754_fmodf) diff --git a/sysdeps/i386/fpu/e_fmodl.S b/sysdeps/i386/fpu/e_fmodl.S new file mode 100644 index 0000000000..7ae63a40ab --- /dev/null +++ b/sysdeps/i386/fpu/e_fmodl.S @@ -0,0 +1,21 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__ieee754_fmodl) + fldt 16(%esp) + fldt 4(%esp) +1: fprem + fstsw %ax + sahf + jp 1b + fstp %st(1) + ret +END (__ieee754_fmodl) diff --git a/sysdeps/i386/fpu/e_hypot.S b/sysdeps/i386/fpu/e_hypot.S new file mode 100644 index 0000000000..07a32878a0 --- /dev/null +++ b/sysdeps/i386/fpu/e_hypot.S @@ -0,0 +1,62 @@ +/* Compute the hypothenuse of X and Y. + Copyright (C) 1998 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1998. + + 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 <sysdep.h> + + .text +ENTRY(__ieee754_hypot) + fldl 4(%esp) // x + fxam + fnstsw + fldl 12(%esp) // y : x + movb %ah, %ch + fxam + fnstsw + movb %ah, %al + orb %ch, %ah + sahf + jc 1f + fmul %st(0) // y * y : x + fxch // x : y * y + fmul %st(0) // x * x : y * y + faddp // x * x + y * y + fsqrt +2: ret + + // We have to test whether any of the parameters is Inf. + // In this case the result is infinity. +1: andb $0x45, %al + cmpb $5, %al + je 3f // jump if y is Inf + andb $0x45, %ch + cmpb $5, %ch + jne 4f // jump if x is not Inf + fxch +3: fstp %st(1) + fabs + jmp 2b + +4: testb $1, %al + jnz 5f // y is NaN + fxch +5: fstp %st(1) + jmp 2b + +END(__ieee754_hypot) diff --git a/sysdeps/i386/fpu/e_hypotf.S b/sysdeps/i386/fpu/e_hypotf.S new file mode 100644 index 0000000000..bf5416b664 --- /dev/null +++ b/sysdeps/i386/fpu/e_hypotf.S @@ -0,0 +1,62 @@ +/* Compute the hypothenuse of X and Y. + Copyright (C) 1998 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1998. + + 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 <sysdep.h> + + .text +ENTRY(__ieee754_hypotf) + flds 4(%esp) // x + fxam + fnstsw + flds 8(%esp) // y : x + movb %ah, %ch + fxam + fnstsw + movb %ah, %al + orb %ch, %ah + sahf + jc 1f + fmul %st(0) // y * y : x + fxch // x : y * y + fmul %st(0) // x * x : y * y + faddp // x * x + y * y + fsqrt +2: ret + + // We have to test whether any of the parameters is Inf. + // In this case the result is infinity. +1: andb $0x45, %al + cmpb $5, %al + je 3f // jump if y is Inf + andb $0x45, %ch + cmpb $5, %ch + jne 4f // jump if x is not Inf + fxch +3: fstp %st(1) + fabs + jmp 2b + +4: testb $1, %al + jnz 5f // y is NaN + fxch +5: fstp %st(1) + jmp 2b + +END(__ieee754_hypotf) diff --git a/sysdeps/i386/fpu/e_log.S b/sysdeps/i386/fpu/e_log.S new file mode 100644 index 0000000000..c7cacdfb0a --- /dev/null +++ b/sysdeps/i386/fpu/e_log.S @@ -0,0 +1,59 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_log.S,v 1.4 1995/05/08 23:48:39 jtc Exp $") + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_log) + fldln2 // log(2) + fldl 4(%esp) // x : log(2) +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + fld %st // x : x : log(2) + fsubl MO(one) // x-1 : x : log(2) + fld %st // x-1 : x-1 : x : log(2) + fabs // |x-1| : x-1 : x : log(2) + fcompl MO(limit) // x-1 : x : log(2) + fnstsw // x-1 : x : log(2) + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : log(2) + fyl2xp1 // log(x) + ret + +2: fstp %st(0) // x : log(2) + fyl2x // log(x) + ret +END (__ieee754_log) diff --git a/sysdeps/i386/fpu/e_log10.S b/sysdeps/i386/fpu/e_log10.S new file mode 100644 index 0000000000..2c8488c3a9 --- /dev/null +++ b/sysdeps/i386/fpu/e_log10.S @@ -0,0 +1,59 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_log10.S,v 1.4 1995/05/08 23:49:24 jtc Exp $") + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_log10) + fldlg2 // log10(2) + fldl 4(%esp) // x : log10(2) +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + fld %st // x : x : log10(2) + fsubl MO(one) // x-1 : x : log10(2) + fld %st // x-1 : x-1 : x : log10(2) + fabs // |x-1| : x-1 : x : log10(2) + fcompl MO(limit) // x-1 : x : log10(2) + fnstsw // x-1 : x : log10(2) + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : log10(2) + fyl2xp1 // log10(x) + ret + +2: fstp %st(0) // x : log10(2) + fyl2x // log10(x) + ret +END (__ieee754_log10) diff --git a/sysdeps/i386/fpu/e_log10f.S b/sysdeps/i386/fpu/e_log10f.S new file mode 100644 index 0000000000..2c07161085 --- /dev/null +++ b/sysdeps/i386/fpu/e_log10f.S @@ -0,0 +1,60 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * Adapted for float type by Ulrich Drepper <drepper@cygnus.com>. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_log10f) + fldlg2 // log10(2) + flds 4(%esp) // x : log10(2) +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + fld %st // x : x : log10(2) + fsubl MO(one) // x-1 : x : log10(2) + fld %st // x-1 : x-1 : x : log10(2) + fabs // |x-1| : x-1 : x : log10(2) + fcompl MO(limit) // x-1 : x : log10(2) + fnstsw // x-1 : x : log10(2) + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : log10(2) + fyl2xp1 // log10(x) + ret + +2: fstp %st(0) // x : log10(2) + fyl2x // log10(x) + ret +END (__ieee754_log10f) diff --git a/sysdeps/i386/fpu/e_log10l.S b/sysdeps/i386/fpu/e_log10l.S new file mode 100644 index 0000000000..6fe7c5a6f7 --- /dev/null +++ b/sysdeps/i386/fpu/e_log10l.S @@ -0,0 +1,61 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_log10l) + fldlg2 // log10(2) + fldt 4(%esp) // x : log10(2) +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + fld %st // x : x : log10(2) + fsubl MO(one) // x-1 : x : log10(2) + fld %st // x-1 : x-1 : x : log10(2) + fabs // |x-1| : x-1 : x : log10(2) + fcompl MO(limit) // x-1 : x : log10(2) + fnstsw // x-1 : x : log10(2) + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : log10(2) + fyl2xp1 // log10(x) + ret + +2: fstp %st(0) // x : log10(2) + fyl2x // log10(x) + ret +END(__ieee754_log10l) diff --git a/sysdeps/i386/fpu/e_logf.S b/sysdeps/i386/fpu/e_logf.S new file mode 100644 index 0000000000..bdba1d3225 --- /dev/null +++ b/sysdeps/i386/fpu/e_logf.S @@ -0,0 +1,60 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * Adapted for float by Ulrich Drepper <drepper@cygnus.com>. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_log.S,v 1.4 1995/05/08 23:48:39 jtc Exp $") + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_logf) + fldln2 // log(2) + flds 4(%esp) // x : log(2) +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + fld %st // x : x : log(2) + fsubl MO(one) // x-1 : x : log(2) + fld %st // x-1 : x-1 : x : log(2) + fabs // |x-1| : x-1 : x : log(2) + fcompl MO(limit) // x-1 : x : log(2) + fnstsw // x-1 : x : log(2) + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : log(2) + fyl2xp1 // log(x) + ret + +2: fstp %st(0) // x : log(2) + fyl2x // log(x) + ret +END (__ieee754_logf) diff --git a/sysdeps/i386/fpu/e_logl.S b/sysdeps/i386/fpu/e_logl.S new file mode 100644 index 0000000000..bda3ea508e --- /dev/null +++ b/sysdeps/i386/fpu/e_logl.S @@ -0,0 +1,60 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__ieee754_logl) + fldln2 // log(2) + fldt 4(%esp) // x : log(2) +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + fld %st // x : x : log(2) + fsubl MO(one) // x-1 : x : log(2) + fld %st // x-1 : x-1 : x : log(2) + fabs // |x-1| : x-1 : x : log(2) + fcompl MO(limit) // x-1 : x : log(2) + fnstsw // x-1 : x : log(2) + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : log(2) + fyl2xp1 // log(x) + ret + +2: fstp %st(0) // x : log(2) + fyl2x // log(x) + ret +END (__ieee754_logl) diff --git a/sysdeps/i386/fpu/e_pow.S b/sysdeps/i386/fpu/e_pow.S new file mode 100644 index 0000000000..75ad211872 --- /dev/null +++ b/sysdeps/i386/fpu/e_pow.S @@ -0,0 +1,320 @@ +/* ix87 specific implementation of pow function. + Copyright (C) 1996, 1997, 1998 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(infinity,@object) +inf_zero: +infinity: + .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f + ASM_SIZE_DIRECTIVE(infinity) + ASM_TYPE_DIRECTIVE(zero,@object) +zero: .double 0.0 + ASM_SIZE_DIRECTIVE(zero) + ASM_TYPE_DIRECTIVE(minf_mzero,@object) +minf_mzero: +minfinity: + .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff +mzero: + .byte 0, 0, 0, 0, 0, 0, 0, 0x80 + ASM_SIZE_DIRECTIVE(minf_mzero) + 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) +#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) +#else +#define MO(op) op +#define MOX(op,x,f) op(,x,f) +#endif + + .text +ENTRY(__ieee754_pow) + fldl 12(%esp) // y + fxam + +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + + fnstsw + movb %ah, %dl + andb $0x45, %ah + cmpb $0x40, %ah // is y == 0 ? + je 11f + + cmpb $0x05, %ah // is y == ±inf ? + je 12f + + cmpb $0x01, %ah // is y == NaN ? + je 30f + + fldl 4(%esp) // x : y + + subl $8,%esp + + fxam + fnstsw + movb %ah, %dh + andb $0x45, %ah + cmpb $0x40, %ah + je 20f // x is ±0 + + cmpb $0x05, %ah + je 15f // x is ±inf + + fxch // y : x + + /* 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. */ + popl %eax + popl %edx + orl $0, %edx + fstp %st(0) // x + jns 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 + shrl $1, %edx + movl %eax, %ecx + orl %edx, %ecx + jnz 6b + fstp %st(0) // ST*x +30: 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 + + + // pow(x,±0) = 1 + .align ALIGNARG(4) +11: fstp %st(0) // pop y + fldl MO(one) + ret + + // y == ±inf + .align ALIGNARG(4) +12: fstp %st(0) // pop y + fldl 4(%esp) // x + fabs + fcompl MO(one) // < 1, == 1, or > 1 + fnstsw + andb $0x45, %ah + cmpb $0x45, %ah + je 13f // jump if x is NaN + + cmpb $0x40, %ah + je 14f // jump if |x| == 1 + + shlb $1, %ah + xorb %ah, %dl + andl $2, %edx + fldl MOX(inf_zero, %edx, 4) + ret + + .align ALIGNARG(4) +14: fldl MO(infinity) + fmull MO(zero) // raise invalid exception + ret + + .align ALIGNARG(4) +13: fldl 4(%esp) // load x == NaN + ret + + .align ALIGNARG(4) + // x is ±inf +15: fstp %st(0) // y + testb $2, %dh + jz 16f // jump if x == +inf + + // We must find out whether y is an odd integer. + fld %st // y : y + fistpll (%esp) // y + fildll (%esp) // int(y) : y + fucompp // <empty> + fnstsw + sahf + jne 17f + + // OK, the value is an integer, but is the number of bits small + // enough so that all are coming from the mantissa? + popl %eax + popl %edx + andb $1, %al + jz 18f // jump if not odd + movl %edx, %eax + orl %edx, %edx + jns 155f + negl %eax +155: cmpl $0x00200000, %eax + ja 18f // does not fit in mantissa bits + // It's an odd integer. + shrl $31, %edx + fldl MOX(minf_mzero, %edx, 8) + ret + + .align ALIGNARG(4) +16: fcompl MO(zero) + addl $8, %esp + fnstsw + shrl $5, %eax + andl $8, %eax + fldl MOX(inf_zero, %eax, 1) + ret + + .align ALIGNARG(4) +17: shll $30, %edx // sign bit for y in right position + addl $8, %esp +18: shrl $31, %edx + fldl MOX(inf_zero, %edx, 8) + ret + + .align ALIGNARG(4) + // x is ±0 +20: fstp %st(0) // y + testb $2, %dl + jz 21f // y > 0 + + // x is ±0 and y is < 0. We must find out whether y is an odd integer. + testb $2, %dh + jz 25f + + fld %st // y : y + fistpll (%esp) // y + fildll (%esp) // int(y) : y + fucompp // <empty> + fnstsw + sahf + jne 26f + + // OK, the value is an integer, but is the number of bits small + // enough so that all are coming from the mantissa? + popl %eax + popl %edx + andb $1, %al + jz 27f // jump if not odd + cmpl $0xffe00000, %edx + jbe 27f // does not fit in mantissa bits + // It's an odd integer. + // Raise divide-by-zero exception and get minus infinity value. + fldl MO(one) + fdivl MO(zero) + fchs + ret + +25: fstp %st(0) +26: popl %eax + popl %edx +27: // Raise divide-by-zero exception and get infinity value. + fldl MO(one) + fdivl MO(zero) + ret + + .align ALIGNARG(4) + // x is ±0 and y is > 0. We must find out whether y is an odd integer. +21: testb $2, %dh + jz 22f + + fld %st // y : y + fistpll (%esp) // y + fildll (%esp) // int(y) : y + fucompp // <empty> + fnstsw + sahf + jne 23f + + // OK, the value is an integer, but is the number of bits small + // enough so that all are coming from the mantissa? + popl %eax + popl %edx + andb $1, %al + jz 24f // jump if not odd + cmpl $0xffe00000, %edx + jae 24f // does not fit in mantissa bits + // It's an odd integer. + fldl MO(mzero) + ret + +22: fstp %st(0) +23: popl %eax + popl %edx +24: fldl MO(zero) + ret + +END(__ieee754_pow) diff --git a/sysdeps/i386/fpu/e_powf.S b/sysdeps/i386/fpu/e_powf.S new file mode 100644 index 0000000000..d7342bf56f --- /dev/null +++ b/sysdeps/i386/fpu/e_powf.S @@ -0,0 +1,310 @@ +/* ix87 specific implementation of pow function. + Copyright (C) 1996, 1997 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(infinity,@object) +inf_zero: +infinity: + .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f + ASM_SIZE_DIRECTIVE(infinity) + ASM_TYPE_DIRECTIVE(zero,@object) +zero: .double 0.0 + ASM_SIZE_DIRECTIVE(zero) + ASM_TYPE_DIRECTIVE(minf_mzero,@object) +minf_mzero: +minfinity: + .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff +mzero: + .byte 0, 0, 0, 0, 0, 0, 0, 0x80 + ASM_SIZE_DIRECTIVE(minf_mzero) + 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) +#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) +#else +#define MO(op) op +#define MOX(op,x,f) op(,x,f) +#endif + + .text +ENTRY(__ieee754_powf) + flds 8(%esp) // y + fxam + +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + + fnstsw + movb %ah, %dl + andb $0x45, %ah + cmpb $0x40, %ah // is y == 0 ? + je 11f + + cmpb $0x05, %ah // is y == ±inf ? + je 12f + + cmpb $0x01, %ah // is y == NaN ? + je 30f + + flds 4(%esp) // x : y + + subl $4, %esp + + fxam + fnstsw + movb %ah, %dh + andb $0x45, %ah + cmpb $0x40, %ah + je 20f // x is ±0 + + cmpb $0x05, %ah + je 15f // x is ±inf + + fxch // y : x + + /* First see whether `y' is a natural number. In this case we + can use a more precise algorithm. */ + fld %st // y : y : x + fistpl (%esp) // y : x + fildl (%esp) // int(y) : y : x + fucomp %st(1) // y : x + fnstsw + sahf + jne 2f + + /* OK, we have an integer value for y. */ + popl %edx + orl $0, %edx + fstp %st(0) // x + jns 4f // y >= 0, jump + fdivrl MO(one) // 1/x (now referred to as x) + negl %edx +4: fldl MO(one) // 1 : x + fxch + +6: shrl $1, %edx + jnc 5f + fxch + fmul %st(1) // x : ST*x + fxch +5: fmul %st(0), %st // x*x : ST*x + testl %edx, %edx + jnz 6b + fstp %st(0) // ST*x +30: 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 $4, %esp + fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x)) + ret + + + // pow(x,±0) = 1 + .align ALIGNARG(4) +11: fstp %st(0) // pop y + fldl MO(one) + ret + + // y == ±inf + .align ALIGNARG(4) +12: fstp %st(0) // pop y + flds 4(%esp) // x + fabs + fcompl MO(one) // < 1, == 1, or > 1 + fnstsw + andb $0x45, %ah + cmpb $0x45, %ah + je 13f // jump if x is NaN + + cmpb $0x40, %ah + je 14f // jump if |x| == 1 + + shlb $1, %ah + xorb %ah, %dl + andl $2, %edx + fldl MOX(inf_zero, %edx, 4) + ret + + .align ALIGNARG(4) +14: fldl MO(infinity) + fmull MO(zero) // raise invalid exception + ret + + .align ALIGNARG(4) +13: flds 4(%esp) // load x == NaN + ret + + .align ALIGNARG(4) + // x is ±inf +15: fstp %st(0) // y + testb $2, %dh + jz 16f // jump if x == +inf + + // We must find out whether y is an odd integer. + fld %st // y : y + fistpl (%esp) // y + fildl (%esp) // int(y) : y + fucompp // <empty> + fnstsw + sahf + jne 17f + + // OK, the value is an integer, but is the number of bits small + // enough so that all are coming from the mantissa? + popl %edx + testb $1, %dl + jz 18f // jump if not odd + movl %edx, %eax + orl %edx, %edx + jns 155f + negl %eax +155: cmpl $0x01000000, %eax + ja 18f // does not fit in mantissa bits + // It's an odd integer. + shrl $31, %edx + fldl MOX(minf_mzero, %edx, 8) + ret + + .align ALIGNARG(4) +16: fcompl MO(zero) + addl $4, %esp + fnstsw + shrl $5, %eax + andl $8, %eax + fldl MOX(inf_zero, %eax, 1) + ret + + .align ALIGNARG(4) +17: shll $30, %edx // sign bit for y in right position + addl $4, %esp +18: shrl $31, %edx + fldl MOX(inf_zero, %edx, 8) + ret + + .align ALIGNARG(4) + // x is ±0 +20: fstp %st(0) // y + testb $2, %dl + jz 21f // y > 0 + + // x is ±0 and y is < 0. We must find out whether y is an odd integer. + testb $2, %dh + jz 25f + + fld %st // y : y + fistpl (%esp) // y + fildl (%esp) // int(y) : y + fucompp // <empty> + fnstsw + sahf + jne 26f + + // OK, the value is an integer, but is the number of bits small + // enough so that all are coming from the mantissa? + popl %edx + testb $1, %dl + jz 27f // jump if not odd + cmpl $0xff000000, %edx + jbe 27f // does not fit in mantissa bits + // It's an odd integer. + // Raise divide-by-zero exception and get minus infinity value. + fldl MO(one) + fdivl MO(zero) + fchs + ret + +25: fstp %st(0) +26: popl %eax +27: // Raise divide-by-zero exception and get infinity value. + fldl MO(one) + fdivl MO(zero) + ret + + .align ALIGNARG(4) + // x is ±0 and y is > 0. We must find out whether y is an odd integer. +21: testb $2, %dh + jz 22f + + fld %st // y : y + fistpl (%esp) // y + fildl (%esp) // int(y) : y + fucompp // <empty> + fnstsw + sahf + jne 23f + + // OK, the value is an integer, but is the number of bits small + // enough so that all are coming from the mantissa? + popl %edx + testb $1, %dl + jz 24f // jump if not odd + cmpl $0xff000000, %edx + jae 24f // does not fit in mantissa bits + // It's an odd integer. + fldl MO(mzero) + ret + +22: fstp %st(0) +23: popl %eax +24: fldl MO(zero) + ret + +END(__ieee754_powf) diff --git a/sysdeps/i386/fpu/e_powl.S b/sysdeps/i386/fpu/e_powl.S new file mode 100644 index 0000000000..2e09dcc820 --- /dev/null +++ b/sysdeps/i386/fpu/e_powl.S @@ -0,0 +1,307 @@ +/* ix87 specific implementation of pow function. + Copyright (C) 1996, 1997, 1998 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(infinity,@object) +inf_zero: +infinity: + .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f + ASM_SIZE_DIRECTIVE(infinity) + ASM_TYPE_DIRECTIVE(zero,@object) +zero: .double 0.0 + ASM_SIZE_DIRECTIVE(zero) + ASM_TYPE_DIRECTIVE(minf_mzero,@object) +minf_mzero: +minfinity: + .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff +mzero: + .byte 0, 0, 0, 0, 0, 0, 0, 0x80 + ASM_SIZE_DIRECTIVE(minf_mzero) + 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) +#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) +#else +#define MO(op) op +#define MOX(op,x,f) op(,x,f) +#endif + + .text +ENTRY(__ieee754_powl) + fldt 16(%esp) // y + fxam + +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + + fnstsw + movb %ah, %dl + andb $0x45, %ah + cmpb $0x40, %ah // is y == 0 ? + je 11f + + cmpb $0x05, %ah // is y == ±inf ? + je 12f + + cmpb $0x01, %ah // is y == NaN ? + je 30f + + fldt 4(%esp) // x : y + + subl $8,%esp + + fxam + fnstsw + movb %ah, %dh + andb $0x45, %ah + cmpb $0x40, %ah + je 20f // x is ±0 + + cmpb $0x05, %ah + je 15f // x is ±inf + + fxch // y : x + + /* 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. */ + popl %eax + popl %edx + orl $0, %edx + fstp %st(0) // x + jns 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 + shrl $1, %edx + movl %eax, %ecx + orl %edx, %ecx + jnz 6b + fstp %st(0) // ST*x +30: 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 + + + // pow(x,±0) = 1 + .align ALIGNARG(4) +11: fstp %st(0) // pop y + fldl MO(one) + ret + + // y == ±inf + .align ALIGNARG(4) +12: fstp %st(0) // pop y + fldt 4(%esp) // x + fabs + fcompl MO(one) // < 1, == 1, or > 1 + fnstsw + andb $0x45, %ah + cmpb $0x45, %ah + je 13f // jump if x is NaN + + cmpb $0x40, %ah + je 14f // jump if |x| == 1 + + shlb $1, %ah + xorb %ah, %dl + andl $2, %edx + fldl MOX(inf_zero, %edx, 4) + ret + + .align ALIGNARG(4) +14: fldl MO(infinity) + fmull MO(zero) // raise invalid exception + ret + + .align ALIGNARG(4) +13: fldt 4(%esp) // load x == NaN + ret + + .align ALIGNARG(4) + // x is ±inf +15: fstp %st(0) // y + testb $2, %dh + jz 16f // jump if x == +inf + + // We must find out whether y is an odd integer. + fld %st // y : y + fistpll (%esp) // y + fildll (%esp) // int(y) : y + fucompp // <empty> + fnstsw + sahf + jne 17f + + // OK, the value is an integer, but is it odd? + popl %eax + popl %edx + andb $1, %al + jz 18f // jump if not odd + // It's an odd integer. + shrl $31, %edx + fldl MOX(minf_mzero, %edx, 8) + ret + + .align ALIGNARG(4) +16: fcompl MO(zero) + addl $8, %esp + fnstsw + shrl $5, %eax + andl $8, %eax + fldl MOX(inf_zero, %eax, 1) + ret + + .align ALIGNARG(4) +17: shll $30, %edx // sign bit for y in right position + addl $8, %esp +18: shrl $31, %edx + fldl MOX(inf_zero, %edx, 8) + ret + + .align ALIGNARG(4) + // x is ±0 +20: fstp %st(0) // y + testb $2, %dl + jz 21f // y > 0 + + // x is ±0 and y is < 0. We must find out whether y is an odd integer. + testb $2, %dh + jz 25f + + fld %st // y : y + fistpll (%esp) // y + fildll (%esp) // int(y) : y + fucompp // <empty> + fnstsw + sahf + jne 26f + + // OK, the value is an integer, but is it odd? + popl %eax + popl %edx + andb $1, %al + jz 27f // jump if not odd + // It's an odd integer. + // Raise divide-by-zero exception and get minus infinity value. + fldl MO(one) + fdivl MO(zero) + fchs + ret + +25: fstp %st(0) +26: popl %eax + popl %edx +27: // Raise divide-by-zero exception and get infinity value. + fldl MO(one) + fdivl MO(zero) + ret + + .align ALIGNARG(4) + // x is ±0 and y is > 0. We must find out whether y is an odd integer. +21: testb $2, %dh + jz 22f + + fld %st // y : y + fistpll (%esp) // y + fildll (%esp) // int(y) : y + fucompp // <empty> + fnstsw + sahf + jne 23f + + // OK, the value is an integer, but is it odd? + popl %eax + popl %edx + andb $1, %al + jz 24f // jump if not odd + // It's an odd integer. + fldl MO(mzero) + ret + +22: fstp %st(0) +23: popl %eax + popl %edx +24: fldl MO(zero) + ret + +END(__ieee754_powl) diff --git a/sysdeps/i386/fpu/e_rem_pio2.c b/sysdeps/i386/fpu/e_rem_pio2.c new file mode 100644 index 0000000000..1347b0468c --- /dev/null +++ b/sysdeps/i386/fpu/e_rem_pio2.c @@ -0,0 +1,3 @@ +/* Empty. This file is only meant to avoid compiling the file with the + same name in the libm-ieee754 directory. The code is not used since + there is an assembler version for all users of this file. */ diff --git a/sysdeps/i386/fpu/e_rem_pio2f.c b/sysdeps/i386/fpu/e_rem_pio2f.c new file mode 100644 index 0000000000..1347b0468c --- /dev/null +++ b/sysdeps/i386/fpu/e_rem_pio2f.c @@ -0,0 +1,3 @@ +/* Empty. This file is only meant to avoid compiling the file with the + same name in the libm-ieee754 directory. The code is not used since + there is an assembler version for all users of this file. */ diff --git a/sysdeps/i386/fpu/e_rem_pio2l.c b/sysdeps/i386/fpu/e_rem_pio2l.c new file mode 100644 index 0000000000..1347b0468c --- /dev/null +++ b/sysdeps/i386/fpu/e_rem_pio2l.c @@ -0,0 +1,3 @@ +/* Empty. This file is only meant to avoid compiling the file with the + same name in the libm-ieee754 directory. The code is not used since + there is an assembler version for all users of this file. */ diff --git a/sysdeps/i386/fpu/e_remainder.S b/sysdeps/i386/fpu/e_remainder.S new file mode 100644 index 0000000000..2f43cb894c --- /dev/null +++ b/sysdeps/i386/fpu/e_remainder.S @@ -0,0 +1,19 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_remainder.S,v 1.4 1995/05/08 23:49:37 jtc Exp $") + +ENTRY(__ieee754_remainder) + fldl 12(%esp) + fldl 4(%esp) +1: fprem1 + fstsw %ax + sahf + jp 1b + fstp %st(1) + ret +END (__ieee754_remainder) diff --git a/sysdeps/i386/fpu/e_remainderf.S b/sysdeps/i386/fpu/e_remainderf.S new file mode 100644 index 0000000000..79f821993b --- /dev/null +++ b/sysdeps/i386/fpu/e_remainderf.S @@ -0,0 +1,19 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_remainderf.S,v 1.2 1995/05/08 23:49:47 jtc Exp $") + +ENTRY(__ieee754_remainderf) + flds 8(%esp) + flds 4(%esp) +1: fprem1 + fstsw %ax + sahf + jp 1b + fstp %st(1) + ret +END (__ieee754_remainderf) diff --git a/sysdeps/i386/fpu/e_remainderl.S b/sysdeps/i386/fpu/e_remainderl.S new file mode 100644 index 0000000000..5f50b626a2 --- /dev/null +++ b/sysdeps/i386/fpu/e_remainderl.S @@ -0,0 +1,21 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__ieee754_remainderl) + fldt 16(%esp) + fldt 4(%esp) +1: fprem1 + fstsw %ax + sahf + jp 1b + fstp %st(1) + ret +END (__ieee754_remainderl) diff --git a/sysdeps/i386/fpu/e_scalb.S b/sysdeps/i386/fpu/e_scalb.S new file mode 100644 index 0000000000..7ff5541e2f --- /dev/null +++ b/sysdeps/i386/fpu/e_scalb.S @@ -0,0 +1,94 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Correct handling of y==-inf <drepper@gnu> + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_scalb.S,v 1.4 1995/05/08 23:49:52 jtc Exp $") + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(zero_nan,@object) +zero_nan: + .double 0.0 +nan: .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f +minus_zero: + .byte 0, 0, 0, 0, 0, 0, 0, 0x80 + .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f + ASM_SIZE_DIRECTIVE(zero_nan) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%ecx) +#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) +#else +#define MO(op) op +#define MOX(op,x,f) op(,x,f) +#endif + + .text +ENTRY(__ieee754_scalb) + fldl 12(%esp) + fxam + fnstsw + fldl 4(%esp) + andl $0x4700, %eax + cmpl $0x0700, %eax + je 1f + andl $0x4500, %eax + cmpl $0x0100, %eax + je 2f + fxam + fnstsw + andl $0x4500, %eax + cmpl $0x0100, %eax + je 2f + fld %st(1) + frndint + fcomp %st(2) + fnstsw + sahf + jne 2f + fscale + fstp %st(1) + ret + + /* y is -inf */ +1: fxam +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + fnstsw + movl 8(%esp), %edx + shrl $5, %eax + fstp %st + fstp %st + andl $0x80000000, %edx + andl $8, %eax + shrl $27, %edx + addl %edx, %eax + fldl MOX(zero_nan, %eax, 1) + ret + + /* The result is NaN, but we must not raise an exception. + So use a variable. */ +2: fstp %st + fstp %st +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + fldl MO(nan) + ret +END(__ieee754_scalb) diff --git a/sysdeps/i386/fpu/e_scalbf.S b/sysdeps/i386/fpu/e_scalbf.S new file mode 100644 index 0000000000..4222eecc97 --- /dev/null +++ b/sysdeps/i386/fpu/e_scalbf.S @@ -0,0 +1,96 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * Adapted for float type by Ulrich Drepper <drepper@cygnus.com>. + * + * Correct handling of y==-inf <drepper@gnu> + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(zero_nan,@object) +zero_nan: + .double 0.0 +nan: .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f +minus_zero: + .byte 0, 0, 0, 0, 0, 0, 0, 0x80 + .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f + ASM_SIZE_DIRECTIVE(zero_nan) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%ecx) +#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) +#else +#define MO(op) op +#define MOX(op,x,f) op(,x,f) +#endif + + + .text +ENTRY(__ieee754_scalbf) + flds 8(%esp) + fxam + fnstsw + flds 4(%esp) + andl $0x4700, %eax + cmpl $0x0700, %eax + je 1f + andl $0x4500, %eax + cmpl $0x0100, %eax + je 2f + fxam + fnstsw + andl $0x4500, %eax + cmpl $0x0100, %eax + je 2f + fld %st(1) + frndint + fcomp %st(2) + fnstsw + sahf + jne 2f + fscale + fstp %st(1) + ret + + /* y is -inf */ +1: fxam +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + fnstsw + movl 4(%esp), %edx + shrl $5, %eax + fstp %st + fstp %st + andl $0x80000000, %edx + andl $8, %eax + shrl $27, %edx + addl %edx, %eax + fldl MOX(zero_nan, %eax, 1) + ret + + /* The result is NaN, but we must not raise an exception. + So use a variable. */ +2: fstp %st + fstp %st +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + fldl MO(nan) + ret +END(__ieee754_scalbf) diff --git a/sysdeps/i386/fpu/e_scalbl.S b/sysdeps/i386/fpu/e_scalbl.S new file mode 100644 index 0000000000..56cc833a56 --- /dev/null +++ b/sysdeps/i386/fpu/e_scalbl.S @@ -0,0 +1,96 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + * + * Correct handling of y==-inf <drepper@gnu> + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(zero_nan,@object) +zero_nan: + .double 0.0 +nan: .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f +minus_zero: + .byte 0, 0, 0, 0, 0, 0, 0, 0x80 + .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f + ASM_SIZE_DIRECTIVE(zero_nan) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%ecx) +#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) +#else +#define MO(op) op +#define MOX(op,x,f) op(,x,f) +#endif + + .text +ENTRY(__ieee754_scalbl) + fldt 16(%esp) + fxam + fnstsw + fldt 4(%esp) + andl $0x4700, %eax + cmpl $0x0700, %eax + je 1f + andl $0x4500, %eax + cmpl $0x0100, %eax + je 2f + fxam + fnstsw + andl $0x4500, %eax + cmpl $0x0100, %eax + je 2f + fld %st(1) + frndint + fcomp %st(2) + fnstsw + sahf + jne 2f + fscale + fstp %st(1) + ret + + /* y is -inf */ +1: fxam +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + fnstsw + movl 12(%esp), %edx + shrl $5, %eax + fstp %st + fstp %st + andl $0x8000, %edx + andl $8, %eax + shrl $11, %edx + addl %edx, %eax + fldl MOX(zero_nan, %eax, 1) + ret + + /* The result is NaN, but we must not raise an exception. + So use a variable. */ +2: fstp %st + fstp %st +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + fldl MO(nan) + ret +END(__ieee754_scalbl) diff --git a/sysdeps/i386/fpu/e_sqrt.S b/sysdeps/i386/fpu/e_sqrt.S new file mode 100644 index 0000000000..6f253d51aa --- /dev/null +++ b/sysdeps/i386/fpu/e_sqrt.S @@ -0,0 +1,14 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_sqrt.S,v 1.4 1995/05/08 23:49:57 jtc Exp $") + +ENTRY(__ieee754_sqrt) + fldl 4(%esp) + fsqrt + ret +END (__ieee754_sqrt) diff --git a/sysdeps/i386/fpu/e_sqrtf.S b/sysdeps/i386/fpu/e_sqrtf.S new file mode 100644 index 0000000000..5ce1ad0544 --- /dev/null +++ b/sysdeps/i386/fpu/e_sqrtf.S @@ -0,0 +1,14 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: e_sqrtf.S,v 1.2 1995/05/08 23:50:14 jtc Exp $") + +ENTRY(__ieee754_sqrtf) + flds 4(%esp) + fsqrt + ret +END (__ieee754_sqrtf) diff --git a/sysdeps/i386/fpu/e_sqrtl.S b/sysdeps/i386/fpu/e_sqrtl.S new file mode 100644 index 0000000000..d47aae5cb5 --- /dev/null +++ b/sysdeps/i386/fpu/e_sqrtl.S @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__ieee754_sqrtl) + fldt 4(%esp) + fsqrt + ret +END (__ieee754_sqrtl) diff --git a/sysdeps/i386/fpu/k_rem_pio2.c b/sysdeps/i386/fpu/k_rem_pio2.c new file mode 100644 index 0000000000..1347b0468c --- /dev/null +++ b/sysdeps/i386/fpu/k_rem_pio2.c @@ -0,0 +1,3 @@ +/* Empty. This file is only meant to avoid compiling the file with the + same name in the libm-ieee754 directory. The code is not used since + there is an assembler version for all users of this file. */ diff --git a/sysdeps/i386/fpu/k_rem_pio2f.c b/sysdeps/i386/fpu/k_rem_pio2f.c new file mode 100644 index 0000000000..1347b0468c --- /dev/null +++ b/sysdeps/i386/fpu/k_rem_pio2f.c @@ -0,0 +1,3 @@ +/* Empty. This file is only meant to avoid compiling the file with the + same name in the libm-ieee754 directory. The code is not used since + there is an assembler version for all users of this file. */ diff --git a/sysdeps/i386/fpu/k_rem_pio2l.c b/sysdeps/i386/fpu/k_rem_pio2l.c new file mode 100644 index 0000000000..1347b0468c --- /dev/null +++ b/sysdeps/i386/fpu/k_rem_pio2l.c @@ -0,0 +1,3 @@ +/* Empty. This file is only meant to avoid compiling the file with the + same name in the libm-ieee754 directory. The code is not used since + there is an assembler version for all users of this file. */ diff --git a/sysdeps/i386/fpu/s_asinh.S b/sysdeps/i386/fpu/s_asinh.S new file mode 100644 index 0000000000..a4c52cb67c --- /dev/null +++ b/sysdeps/i386/fpu/s_asinh.S @@ -0,0 +1,139 @@ +/* ix87 specific implementation of arcsinh. + Copyright (C) 1996, 1997 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(huge,@object) +huge: .double 1e+300 + ASM_SIZE_DIRECTIVE(huge) + 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(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__asinh) + movl 8(%esp), %ecx + movl $0x7fffffff, %eax + andl %ecx, %eax + andl $0x80000000, %ecx + movl %eax, %edx + orl $0x800fffff, %edx + incl %edx + jz 7f // x in ±Inf or NaN + xorl %ecx, 8(%esp) + fldl 4(%esp) // |x| + cmpl $0x3e300000, %eax + jb 2f // |x| < 2^-28 + fldln2 // log(2) : |x| + cmpl $0x41b00000, %eax + fxch // |x| : log(2) + ja 3f // |x| > 2^28 +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + cmpl $0x40000000, %eax + ja 5f // |x| > 2 + + // 2^-28 <= |x| <= 2 => y = sign(x)*log1p(|x|+|x|^2/(1+sqrt(1+|x|^2))) + fld %st // |x| : |x| : log(2) + fmul %st(1) // |x|^2 : |x| : log(2) + fld %st // |x|^2 : |x|^2 : |x| : log(2) + faddl MO(one) // 1+|x|^2 : |x|^2 : |x| : log(2) + fsqrt // sqrt(1+|x|^2) : |x|^2 : |x| : log(2) + faddl MO(one) // 1+sqrt(1+|x|^2) : |x|^2 : |x| : log(2) + fdivrp // |x|^2/(1+sqrt(1+|x|^2)) : |x| : log(2) + faddp // |x|+|x|^2/(1+sqrt(1+|x|^2)) : log(2) + fcoml MO(limit) + fnstsw + sahf + ja 6f + fyl2xp1 + jecxz 4f + fchs +4: ret + +7: fldl 4(%esp) + ret + +6: faddl MO(one) + fyl2x + jecxz 4f + fchs +4: ret + + // |x| < 2^-28 => y = x (inexact iff |x| != 0.0) + .align ALIGNARG(4) +2: +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + jecxz 4f + fchs // x +4: fld %st // x : x + faddl MO(huge) // huge+x : x + fstp %st(0) // x + ret + + // |x| > 2^28 => y = sign(x) * (log(|x|) + log(2)) + .align ALIGNARG(4) +3: fyl2x // log(|x|) + fldln2 // log(2) : log(|x|) + faddp // log(|x|)+log(2) + jecxz 4f + fchs +4: ret + + // |x| > 2 => y = sign(x) * log(2*|x| + 1/(|x|+sqrt(x*x+1))) + .align ALIGNARG(4) +5: fld %st // |x| : |x| : log(2) + fadd %st, %st(1) // |x| : 2*|x| : log(2) + fld %st // |x| : |x| : 2*|x| : log(2) + fmul %st(1) // |x|^2 : |x| : 2*|x| : log(2) + faddl MO(one) // 1+|x|^2 : |x| : 2*|x| : log(2) + fsqrt // sqrt(1+|x|^2) : |x| : 2*|x| : log(2) + faddp // |x|+sqrt(1+|x|^2) : 2*|x| : log(2) + fdivrl MO(one) // 1/(|x|+sqrt(1+|x|^2)) : 2*|x| : log(2) + faddp // 2*|x|+1/(|x|+sqrt(1+|x|^2)) : log(2) + fyl2x // log(2*|x|+1/(|x|+sqrt(1+|x|^2))) + jecxz 4f + fchs +4: ret +END(__asinh) +weak_alias (__asinh, asinh) diff --git a/sysdeps/i386/fpu/s_asinhf.S b/sysdeps/i386/fpu/s_asinhf.S new file mode 100644 index 0000000000..a6925c7b4c --- /dev/null +++ b/sysdeps/i386/fpu/s_asinhf.S @@ -0,0 +1,139 @@ +/* ix87 specific implementation of arcsinh. + Copyright (C) 1996, 1997, 1999 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(huge,@object) +huge: .double 1e+36 + ASM_SIZE_DIRECTIVE(huge) + 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(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__asinhf) + movl 4(%esp), %ecx + movl $0x7fffffff, %eax + andl %ecx, %eax + andl $0x80000000, %ecx + movl %eax, %edx + orl $0x807fffff, %edx + incl %edx + jz 7f // x in ±Inf or NaN + xorl %ecx, 4(%esp) + flds 4(%esp) // |x| + cmpl $0x38000000, %eax + jb 2f // |x| < 2^-14 + fldln2 // log(2) : |x| + cmpl $0x47000000, %eax + fxch // |x| : log(2) + ja 3f // |x| > 2^14 +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + cmpl $0x40000000, %eax + ja 5f // |x| > 2 + + // 2^-14 <= |x| <= 2 => y = sign(x)*log1p(|x|+|x|^2/(1+sqrt(1+|x|^2))) + fld %st // |x| : |x| : log(2) + fmul %st(1) // |x|^2 : |x| : log(2) + fld %st // |x|^2 : |x|^2 : |x| : log(2) + faddl MO(one) // 1+|x|^2 : |x|^2 : |x| : log(2) + fsqrt // sqrt(1+|x|^2) : |x|^2 : |x| : log(2) + faddl MO(one) // 1+sqrt(1+|x|^2) : |x|^2 : |x| : log(2) + fdivrp // |x|^2/(1+sqrt(1+|x|^2)) : |x| : log(2) + faddp // |x|+|x|^2/(1+sqrt(1+|x|^2)) : log(2) + fcoml MO(limit) + fnstsw + sahf + ja 6f + fyl2xp1 + jecxz 4f + fchs +4: ret + +7: flds 4(%esp) + ret + +6: faddl MO(one) + fyl2x + jecxz 4f + fchs +4: ret + + // |x| < 2^-14 => y = x (inexact iff |x| != 0.0) + .align ALIGNARG(4) +2: +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + jecxz 4f + fchs // x +4: fld %st // x : x + faddl MO(huge) // huge+x : x + fstp %st(0) // x + ret + + // |x| > 2^14 => y = sign(x) * (log(|x|) + log(2)) + .align ALIGNARG(4) +3: fyl2x // log(|x|) + fldln2 // log(2) : log(|x|) + faddp // log(|x|)+log(2) + jecxz 4f + fchs +4: ret + + // |x| > 2 => y = sign(x) * log(2*|x| + 1/(|x|+sqrt(x*x+1))) + .align ALIGNARG(4) +5: fld %st // |x| : |x| : log(2) + fadd %st, %st(1) // |x| : 2*|x| : log(2) + fld %st // |x| : |x| : 2*|x| : log(2) + fmul %st(1) // |x|^2 : |x| : 2*|x| : log(2) + faddl MO(one) // 1+|x|^2 : |x| : 2*|x| : log(2) + fsqrt // sqrt(1+|x|^2) : |x| : 2*|x| : log(2) + faddp // |x|+sqrt(1+|x|^2) : 2*|x| : log(2) + fdivrl MO(one) // 1/(|x|+sqrt(1+|x|^2)) : 2*|x| : log(2) + faddp // 2*|x|+1/(|x|+sqrt(1+|x|^2)) : log(2) + fyl2x // log(2*|x|+1/(|x|+sqrt(1+|x|^2))) + jecxz 4f + fchs +4: ret +END(__asinhf) +weak_alias (__asinhf, asinhf) diff --git a/sysdeps/i386/fpu/s_asinhl.S b/sysdeps/i386/fpu/s_asinhl.S new file mode 100644 index 0000000000..62e29bc58e --- /dev/null +++ b/sysdeps/i386/fpu/s_asinhl.S @@ -0,0 +1,147 @@ +/* ix87 specific implementation of arcsinh. + Copyright (C) 1996, 1997 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(huge,@object) +huge: .tfloat 1e+4930 + ASM_SIZE_DIRECTIVE(huge) + .align ALIGNARG(4) + /* Please note that we use double value for 1.0. This number + has an exact representation and so we don't get accuracy + problems. The advantage is that the code is simpler. */ + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__asinhl) + movl 12(%esp), %ecx + movl $0x7fff, %eax + andl %ecx, %eax + andl $0x8000, %ecx + movl %eax, %edx + orl $0xffff8000, %edx + incl %edx + jz 7f // x in ±Inf or NaN + xorl %ecx, 12(%esp) + fldt 4(%esp) // |x| + cmpl $0x3fde, %eax + jb 2f // |x| < 2^-34 + fldln2 // log(2) : |x| + cmpl $0x4020, %eax + fxch // |x| : log(2) + ja 3f // |x| > 2^34 +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + cmpl $0x4000, %eax + ja 5f // |x| > 2 + + // 2^-34 <= |x| <= 2 => y = sign(x)*log1p(|x|+|x|^2/(1+sqrt(1+|x|^2))) + fld %st // |x| : |x| : log(2) + fmul %st(1) // |x|^2 : |x| : log(2) + fld %st // |x|^2 : |x|^2 : |x| : log(2) + faddl MO(one) // 1+|x|^2 : |x|^2 : |x| : log(2) + fsqrt // sqrt(1+|x|^2) : |x|^2 : |x| : log(2) + faddl MO(one) // 1+sqrt(1+|x|^2) : |x|^2 : |x| : log(2) + fdivrp // |x|^2/(1+sqrt(1+|x|^2)) : |x| : log(2) + faddp // |x|+|x|^2/(1+sqrt(1+|x|^2)) : log(2) + fcoml MO(limit) + fnstsw + sahf + ja 6f + fyl2xp1 + jecxz 4f + fchs +4: ret + +7: fldt 4(%esp) + ret + +6: faddl MO(one) + fyl2x + jecxz 4f + fchs +4: ret + + // |x| < 2^-34 => y = x (inexact iff |x| != 0.0) + .align ALIGNARG(4) +2: +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + jecxz 4f + fchs // x +4: fld %st // x : x + fldt MO(huge) // huge : x : x + faddp // huge+x : x + fstp %st(0) // x + ret + + // |x| > 2^34 => y = sign(x) * (log(|x|) + log(2)) + .align ALIGNARG(4) +3: fyl2x // log(|x|) + fldln2 // log(2) : log(|x|) + faddp // log(|x|)+log(2) + jecxz 4f + fchs +4: ret + + // |x| > 2 => y = sign(x) * log(2*|x| + 1/(|x|+sqrt(x*x+1))) + .align ALIGNARG(4) +5: fld %st // |x| : |x| : log(2) + fadd %st, %st(1) // |x| : 2*|x| : log(2) + fld %st // |x| : |x| : 2*|x| : log(2) + fmul %st(1) // |x|^2 : |x| : 2*|x| : log(2) + faddl MO(one) // 1+|x|^2 : |x| : 2*|x| : log(2) + fsqrt // sqrt(1+|x|^2) : |x| : 2*|x| : log(2) + faddp // |x|+sqrt(1+|x|^2) : 2*|x| : log(2) + fdivrl MO(one) // 1/(|x|+sqrt(1+|x|^2)) : 2*|x| : log(2) + faddp // 2*|x|+1/(|x|+sqrt(1+|x|^2)) : log(2) + fyl2x // log(2*|x|+1/(|x|+sqrt(1+|x|^2))) + jecxz 4f + fchs +4: ret +END(__asinhl) +weak_alias (__asinhl, asinhl) diff --git a/sysdeps/i386/fpu/s_atan.S b/sysdeps/i386/fpu/s_atan.S new file mode 100644 index 0000000000..7502f6d828 --- /dev/null +++ b/sysdeps/i386/fpu/s_atan.S @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_atan.S,v 1.4 1995/05/08 23:50:41 jtc Exp $") + +ENTRY(__atan) + fldl 4(%esp) + fld1 + fpatan + ret +END (__atan) +weak_alias (__atan, atan) diff --git a/sysdeps/i386/fpu/s_atanf.S b/sysdeps/i386/fpu/s_atanf.S new file mode 100644 index 0000000000..70232c8240 --- /dev/null +++ b/sysdeps/i386/fpu/s_atanf.S @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_atanf.S,v 1.3 1995/05/08 23:51:33 jtc Exp $") + +ENTRY(__atanf) + flds 4(%esp) + fld1 + fpatan + ret +END (__atanf) +weak_alias (__atanf, atanf) diff --git a/sysdeps/i386/fpu/s_atanl.S b/sysdeps/i386/fpu/s_atanl.S new file mode 100644 index 0000000000..8b07272764 --- /dev/null +++ b/sysdeps/i386/fpu/s_atanl.S @@ -0,0 +1,18 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__atanl) + fldt 4(%esp) + fld1 + fpatan + ret +END (__atanl) +weak_alias (__atanl, atanl) diff --git a/sysdeps/i386/fpu/s_cbrt.S b/sysdeps/i386/fpu/s_cbrt.S new file mode 100644 index 0000000000..3f6a0174f2 --- /dev/null +++ b/sysdeps/i386/fpu/s_cbrt.S @@ -0,0 +1,201 @@ +/* Compute cubic root of double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Dirk Alboth <dirka@uni-paderborn.de> and + Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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(f7,@object) +f7: .double -0.145263899385486377 + ASM_SIZE_DIRECTIVE(f7) + ASM_TYPE_DIRECTIVE(f6,@object) +f6: .double 0.784932344976639262 + ASM_SIZE_DIRECTIVE(f6) + ASM_TYPE_DIRECTIVE(f5,@object) +f5: .double -1.83469277483613086 + ASM_SIZE_DIRECTIVE(f5) + ASM_TYPE_DIRECTIVE(f4,@object) +f4: .double 2.44693122563534430 + ASM_SIZE_DIRECTIVE(f4) + ASM_TYPE_DIRECTIVE(f3,@object) +f3: .double -2.11499494167371287 + ASM_SIZE_DIRECTIVE(f3) + ASM_TYPE_DIRECTIVE(f2,@object) +f2: .double 1.50819193781584896 + ASM_SIZE_DIRECTIVE(f2) + ASM_TYPE_DIRECTIVE(f1,@object) +f1: .double 0.354895765043919860 + ASM_SIZE_DIRECTIVE(f1) + +#define CBRT2 1.2599210498948731648 +#define ONE_CBRT2 0.793700525984099737355196796584 +#define SQR_CBRT2 1.5874010519681994748 +#define ONE_SQR_CBRT2 0.629960524947436582364439673883 + + ASM_TYPE_DIRECTIVE(factor,@object) +factor: .double ONE_SQR_CBRT2 + .double ONE_CBRT2 + .double 1.0 + .double CBRT2 + .double SQR_CBRT2 + ASM_SIZE_DIRECTIVE(factor) + + ASM_TYPE_DIRECTIVE(two54,@object) +two54: .byte 0, 0, 0, 0, 0, 0, 0x50, 0x43 + ASM_SIZE_DIRECTIVE(two54) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%ebx) +#define MOX(op,x) op##@GOTOFF(%ebx,x,1) +#else +#define MO(op) op +#define MOX(op,x) op(x) +#endif + + .text +ENTRY(__cbrt) + movl 4(%esp), %ecx + movl 8(%esp), %eax + movl %eax, %edx + andl $0x7fffffff, %eax + orl %eax, %ecx + jz 1f + xorl %ecx, %ecx + cmpl $0x7ff00000, %eax + jae 1f + +#ifdef PIC + pushl %ebx + call 3f +3: popl %ebx + addl $_GLOBAL_OFFSET_TABLE_+[.-3b], %ebx +#endif + + cmpl $0x00100000, %eax + jae 2f + +#ifdef PIC + fldl 8(%esp) +#else + fldl 4(%esp) +#endif + fmull MO(two54) + movl $-54, %ecx +#ifdef PIC + fstpl 8(%esp) + movl 12(%esp), %eax +#else + fstpl 4(%esp) + movl 8(%esp), %eax +#endif + movl %eax, %edx + andl $0x7fffffff, %eax + +2: shrl $20, %eax + andl $0x800fffff, %edx + subl $1022, %eax + orl $0x3fe00000, %edx + addl %eax, %ecx +#ifdef PIC + movl %edx, 12(%esp) + + fldl 8(%esp) /* xm */ +#else + movl %edx, 8(%esp) + + fldl 4(%esp) /* xm */ +#endif + fabs + + /* The following code has two tracks: + a) compute the normalized cbrt value + b) compute xe/3 and xe%3 + The right track computes the value for b) and this is done + in an optimized way by avoiding division. + + But why two tracks at all? Very easy: efficiency. Some FP + instruction can overlap with a certain amount of integer (and + FP) instructions. So we get (except for the imull) all + instructions for free. */ + + fld %st(0) /* xm : xm */ + + fmull MO(f7) /* f7*xm : xm */ + movl $1431655766, %eax + faddl MO(f6) /* f6+f7*xm : xm */ + imull %ecx + fmul %st(1) /* (f6+f7*xm)*xm : xm */ + movl %ecx, %eax + faddl MO(f5) /* f5+(f6+f7*xm)*xm : xm */ + sarl $31, %eax + fmul %st(1) /* (f5+(f6+f7*xm)*xm)*xm : xm */ + subl %eax, %edx + faddl MO(f4) /* f4+(f5+(f6+f7*xm)*xm)*xm : xm */ + fmul %st(1) /* (f4+(f5+(f6+f7*xm)*xm)*xm)*xm : xm */ + faddl MO(f3) /* f3+(f4+(f5+(f6+f7*xm)*xm)*xm)*xm : xm */ + fmul %st(1) /* (f3+(f4+(f5+(f6+f7*xm)*xm)*xm)*xm)*xm : xm */ + faddl MO(f2) /* f2+(f3+(f4+(f5+(f6+f7*xm)*xm)*xm)*xm)*xm : xm */ + fmul %st(1) /* (f2+(f3+(f4+(f5+(f6+f7*xm)*xm)*xm)*xm)*xm)*xm : xm */ + faddl MO(f1) /* u:=f1+(f2+(f3+(f4+(f5+(f6+f7*xm)*xm)*xm)*xm)*xm)*xm : xm */ + + fld %st /* u : u : xm */ + fmul %st(1) /* u*u : u : xm */ + fld %st(2) /* xm : u*u : u : xm */ + fadd %st /* 2*xm : u*u : u : xm */ + fxch %st(1) /* u*u : 2*xm : u : xm */ + fmul %st(2) /* t2:=u*u*u : 2*xm : u : xm */ + movl %edx, %eax + fadd %st, %st(1) /* t2 : t2+2*xm : u : xm */ + leal (%edx,%edx,2),%edx + fadd %st(0) /* 2*t2 : t2+2*xm : u : xm */ + subl %edx, %ecx + faddp %st, %st(3) /* t2+2*xm : u : 2*t2+xm */ + shll $3, %ecx + fmulp /* u*(t2+2*xm) : 2*t2+xm */ + fdivp %st, %st(1) /* u*(t2+2*xm)/(2*t2+xm) */ + fmull MOX(16+factor,%ecx) /* u*(t2+2*xm)/(2*t2+xm)*FACT */ + pushl %eax + fildl (%esp) /* xe/3 : u*(t2+2*xm)/(2*t2+xm)*FACT */ + fxch /* u*(t2+2*xm)/(2*t2+xm)*FACT : xe/3 */ + fscale /* u*(t2+2*xm)/(2*t2+xm)*FACT*2^xe/3 */ + popl %edx +#ifdef PIC + movl 12(%esp), %eax + popl %ebx +#else + movl 8(%esp), %eax +#endif + testl %eax, %eax + fstp %st(1) + jns 4f + fchs +4: ret + + /* Return the argument. */ +1: fldl 4(%esp) + ret +END(__cbrt) +weak_alias (__cbrt, cbrt) diff --git a/sysdeps/i386/fpu/s_cbrtf.S b/sysdeps/i386/fpu/s_cbrtf.S new file mode 100644 index 0000000000..a14e04ed2f --- /dev/null +++ b/sysdeps/i386/fpu/s_cbrtf.S @@ -0,0 +1,178 @@ +/* Compute cubic root of float value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Dirk Alboth <dirka@uni-paderborn.de> and + Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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(f3,@object) +f3: .double 0.191502161678719066 + ASM_SIZE_DIRECTIVE(f3) + ASM_TYPE_DIRECTIVE(f2,@object) +f2: .double 0.697570460207922770 + ASM_SIZE_DIRECTIVE(f2) + ASM_TYPE_DIRECTIVE(f1,@object) +f1: .double 0.492659620528969547 + ASM_SIZE_DIRECTIVE(f1) + +#define CBRT2 1.2599210498948731648 +#define ONE_CBRT2 0.793700525984099737355196796584 +#define SQR_CBRT2 1.5874010519681994748 +#define ONE_SQR_CBRT2 0.629960524947436582364439673883 + + ASM_TYPE_DIRECTIVE(factor,@object) + .align ALIGNARG(4) +factor: .double ONE_SQR_CBRT2 + .double ONE_CBRT2 + .double 1.0 + .double CBRT2 + .double SQR_CBRT2 + ASM_SIZE_DIRECTIVE(factor) + + ASM_TYPE_DIRECTIVE(two25,@object) +two25: .byte 0, 0, 0, 0x4c + ASM_SIZE_DIRECTIVE(two25) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%ebx) +#define MOX(op,x) op##@GOTOFF(%ebx,x,1) +#else +#define MO(op) op +#define MOX(op,x) op(x) +#endif + + .text +ENTRY(__cbrtf) + movl 4(%esp), %eax + xorl %ecx, %ecx + movl %eax, %edx + andl $0x7fffffff, %eax + jz 1f + cmpl $0x7f800000, %eax + jae 1f + +#ifdef PIC + pushl %ebx + call 3f +3: popl %ebx + addl $_GLOBAL_OFFSET_TABLE_+[.-3b], %ebx +#endif + + cmpl $0x00800000, %eax + jae 2f + +#ifdef PIC + flds 8(%esp) +#else + flds 4(%esp) +#endif + fmuls MO(two25) + movl $-25, %ecx +#ifdef PIC + fstps 8(%esp) + movl 8(%esp), %eax +#else + fstps 4(%esp) + movl 4(%esp), %eax +#endif + movl %eax, %edx + andl $0x7fffffff, %eax + +2: shrl $23, %eax + andl $0x807fffff, %edx + subl $126, %eax + orl $0x3f000000, %edx + addl %eax, %ecx +#ifdef PIC + movl %edx, 8(%esp) + + flds 8(%esp) /* xm */ +#else + movl %edx, 4(%esp) + + flds 4(%esp) /* xm */ +#endif + fabs + + /* The following code has two tracks: + a) compute the normalized cbrt value + b) compute xe/3 and xe%3 + The right track computes the value for b) and this is done + in an optimized way by avoiding division. + + But why two tracks at all? Very easy: efficiency. Some FP + instruction can overlap with a certain amount of integer (and + FP) instructions. So we get (except for the imull) all + instructions for free. */ + + fld %st(0) /* xm : xm */ + fmull MO(f3) /* f3*xm : xm */ + movl $1431655766, %eax + fsubrl MO(f2) /* f2-f3*xm : xm */ + imull %ecx + fmul %st(1) /* (f2-f3*xm)*xm : xm */ + movl %ecx, %eax + faddl MO(f1) /* u:=f1+(f2-f3*xm)*xm : xm */ + sarl $31, %eax + fld %st /* u : u : xm */ + subl %eax, %edx + fmul %st(1) /* u*u : u : xm */ + fld %st(2) /* xm : u*u : u : xm */ + fadd %st /* 2*xm : u*u : u : xm */ + fxch %st(1) /* u*u : 2*xm : u : xm */ + fmul %st(2) /* t2:=u*u*u : 2*xm : u : xm */ + movl %edx, %eax + fadd %st, %st(1) /* t2 : t2+2*xm : u : xm */ + leal (%edx,%edx,2),%edx + fadd %st(0) /* 2*t2 : t2+2*xm : u : xm */ + subl %edx, %ecx + faddp %st, %st(3) /* t2+2*xm : u : 2*t2+xm */ + shll $3, %ecx + fmulp /* u*(t2+2*xm) : 2*t2+xm */ + fdivp %st, %st(1) /* u*(t2+2*xm)/(2*t2+xm) */ + fmull MOX(16+factor,%ecx) /* u*(t2+2*xm)/(2*t2+xm)*FACT */ + pushl %eax + fildl (%esp) /* xe/3 : u*(t2+2*xm)/(2*t2+xm)*FACT */ + fxch /* u*(t2+2*xm)/(2*t2+xm)*FACT : xe/3 */ + fscale /* u*(t2+2*xm)/(2*t2+xm)*FACT*2^xe/3 */ + popl %edx +#ifdef PIC + movl 8(%esp), %eax + popl %ebx +#else + movl 4(%esp), %eax +#endif + testl %eax, %eax + fstp %st(1) + jns 4f + fchs +4: ret + + /* Return the argument. */ +1: flds 4(%esp) + ret +END(__cbrtf) +weak_alias (__cbrtf, cbrtf) diff --git a/sysdeps/i386/fpu/s_cbrtl.S b/sysdeps/i386/fpu/s_cbrtl.S new file mode 100644 index 0000000000..6a3b9a8dc5 --- /dev/null +++ b/sysdeps/i386/fpu/s_cbrtl.S @@ -0,0 +1,229 @@ +/* Compute cubic root of long double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Dirk Alboth <dirka@uni-paderborn.de> and + Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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(f8,@object) +f8: .tfloat 0.161617097923756032 + ASM_SIZE_DIRECTIVE(f8) + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(f7,@object) +f7: .tfloat -0.988553671195413709 + ASM_SIZE_DIRECTIVE(f7) + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(f6,@object) +f6: .tfloat 2.65298938441952296 + ASM_SIZE_DIRECTIVE(f6) + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(f5,@object) +f5: .tfloat -4.11151425200350531 + ASM_SIZE_DIRECTIVE(f5) + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(f4,@object) +f4: .tfloat 4.09559907378707839 + ASM_SIZE_DIRECTIVE(f4) + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(f3,@object) +f3: .tfloat -2.82414939754975962 + ASM_SIZE_DIRECTIVE(f3) + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(f2,@object) +f2: .tfloat 1.67595307700780102 + ASM_SIZE_DIRECTIVE(f2) + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(f1,@object) +f1: .tfloat 0.338058687610520237 + ASM_SIZE_DIRECTIVE(f1) + +#define CBRT2 1.2599210498948731648 +#define ONE_CBRT2 0.793700525984099737355196796584 +#define SQR_CBRT2 1.5874010519681994748 +#define ONE_SQR_CBRT2 0.629960524947436582364439673883 + + /* We make the entries in the following table all 16 bytes + wide to avoid having to implement a multiplication by 10. */ + ASM_TYPE_DIRECTIVE(factor,@object) + .align ALIGNARG(4) +factor: .tfloat ONE_SQR_CBRT2 + .byte 0, 0, 0, 0, 0, 0 + .tfloat ONE_CBRT2 + .byte 0, 0, 0, 0, 0, 0 + .tfloat 1.0 + .byte 0, 0, 0, 0, 0, 0 + .tfloat CBRT2 + .byte 0, 0, 0, 0, 0, 0 + .tfloat SQR_CBRT2 + ASM_SIZE_DIRECTIVE(factor) + + ASM_TYPE_DIRECTIVE(two64,@object) + .align ALIGNARG(4) +two64: .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43 + ASM_SIZE_DIRECTIVE(two64) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%ebx) +#define MOX(op,x) op##@GOTOFF(%ebx,x,1) +#else +#define MO(op) op +#define MOX(op,x) op(x) +#endif + + .text +ENTRY(__cbrtl) + movl 4(%esp), %ecx + movl 12(%esp), %eax + orl 8(%esp), %ecx + movl %eax, %edx + andl $0x7fff, %eax + orl %eax, %ecx + jz 1f + xorl %ecx, %ecx + cmpl $0x7fff, %eax + je 1f + +#ifdef PIC + pushl %ebx + call 3f +3: popl %ebx + addl $_GLOBAL_OFFSET_TABLE_+[.-3b], %ebx +#endif + + cmpl $0, %eax + jne 2f + +#ifdef PIC + fldt 8(%esp) +#else + fldt 4(%esp) +#endif + fmull MO(two64) + movl $-64, %ecx +#ifdef PIC + fstpt 8(%esp) + movl 16(%esp), %eax +#else + fstpt 4(%esp) + movl 12(%esp), %eax +#endif + movl %eax, %edx + andl $0x7fff, %eax + +2: andl $0x8000, %edx + subl $16382, %eax + orl $0x3ffe, %edx + addl %eax, %ecx +#ifdef PIC + movl %edx, 16(%esp) + + fldt 8(%esp) /* xm */ +#else + movl %edx, 12(%esp) + + fldt 4(%esp) /* xm */ +#endif + fabs + + /* The following code has two tracks: + a) compute the normalized cbrt value + b) compute xe/3 and xe%3 + The right track computes the value for b) and this is done + in an optimized way by avoiding division. + + But why two tracks at all? Very easy: efficiency. Some FP + instruction can overlap with a certain amount of integer (and + FP) instructions. So we get (except for the imull) all + instructions for free. */ + + fldt MO(f8) /* f8 : xm */ + fmul %st(1) /* f8*xm : xm */ + + fldt MO(f7) + faddp /* f7+f8*xm : xm */ + fmul %st(1) /* (f7+f8*xm)*xm : xm */ + movl $1431655766, %eax + fldt MO(f6) + faddp /* f6+(f7+f8*xm)*xm : xm */ + imull %ecx + fmul %st(1) /* (f6+(f7+f8*xm)*xm)*xm : xm */ + movl %ecx, %eax + fldt MO(f5) + faddp /* f5+(f6+(f7+f8*xm)*xm)*xm : xm */ + sarl $31, %eax + fmul %st(1) /* (f5+(f6+(f7+f8*xm)*xm)*xm)*xm : xm */ + subl %eax, %edx + fldt MO(f4) + faddp /* f4+(f5+(f6+(f7+f8*xm)*xm)*xm)*xm : xm */ + fmul %st(1) /* (f4+(f5+(f6+(f7+f8*xm)*xm)*xm)*xm)*xm : xm */ + fldt MO(f3) + faddp /* f3+(f4+(f5+(f6+(f7+f8*xm)*xm)*xm)*xm)*xm : xm */ + fmul %st(1) /* (f3+(f4+(f5+(f6+(f7+f8*xm)*xm)*xm)*xm)*xm)*xm : xm */ + fldt MO(f2) + faddp /* f2+(f3+(f4+(f5+(f6+(f7+f8*xm)*xm)*xm)*xm)*xm)*xm : xm */ + fmul %st(1) /* (f2+(f3+(f4+(f5+(f6+(f7+f8*xm)*xm)*xm)*xm)*xm)*xm)*xm : xm */ + fldt MO(f1) + faddp /* u:=f1+(f2+(f3+(f4+(f5+(f6+(f7+f8*xm)*xm)*xm)*xm)*xm)*xm)*xm : xm */ + + fld %st /* u : u : xm */ + fmul %st(1) /* u*u : u : xm */ + fld %st(2) /* xm : u*u : u : xm */ + fadd %st /* 2*xm : u*u : u : xm */ + fxch %st(1) /* u*u : 2*xm : u : xm */ + fmul %st(2) /* t2:=u*u*u : 2*xm : u : xm */ + movl %edx, %eax + fadd %st, %st(1) /* t2 : t2+2*xm : u : xm */ + leal (%edx,%edx,2),%edx + fadd %st(0) /* 2*t2 : t2+2*xm : u : xm */ + subl %edx, %ecx + faddp %st, %st(3) /* t2+2*xm : u : 2*t2+xm */ + shll $4, %ecx + fmulp /* u*(t2+2*xm) : 2*t2+xm */ + fdivp %st, %st(1) /* u*(t2+2*xm)/(2*t2+xm) */ + fldt MOX(32+factor,%ecx) + fmulp /* u*(t2+2*xm)/(2*t2+xm)*FACT */ + pushl %eax + fildl (%esp) /* xe/3 : u*(t2+2*xm)/(2*t2+xm)*FACT */ + fxch /* u*(t2+2*xm)/(2*t2+xm)*FACT : xe/3 */ + fscale /* u*(t2+2*xm)/(2*t2+xm)*FACT*2^xe/3 */ + popl %edx +#ifdef PIC + movl 16(%esp), %eax + popl %ebx +#else + movl 12(%esp), %eax +#endif + testl $0x8000, %eax + fstp %st(1) + jz 4f + fchs +4: ret + + /* Return the argument. */ +1: fldt 4(%esp) + ret +END(__cbrtl) +weak_alias (__cbrtl, cbrtl) diff --git a/sysdeps/i386/fpu/s_ceil.S b/sysdeps/i386/fpu/s_ceil.S new file mode 100644 index 0000000000..b0159128aa --- /dev/null +++ b/sysdeps/i386/fpu/s_ceil.S @@ -0,0 +1,32 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_ceil.S,v 1.4 1995/05/08 23:52:13 jtc Exp $") + +ENTRY(__ceil) + fldl 4(%esp) + subl $8,%esp + + 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 $0x0800,%edx /* round towards +oo */ + orl 4(%esp),%edx + andl $0xfbff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ + + frndint /* round */ + + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp + ret +END (__ceil) +weak_alias (__ceil, ceil) diff --git a/sysdeps/i386/fpu/s_ceilf.S b/sysdeps/i386/fpu/s_ceilf.S new file mode 100644 index 0000000000..352d40d7ce --- /dev/null +++ b/sysdeps/i386/fpu/s_ceilf.S @@ -0,0 +1,32 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_ceilf.S,v 1.3 1995/05/08 23:52:44 jtc Exp $") + +ENTRY(__ceilf) + flds 4(%esp) + subl $8,%esp + + 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 $0x0800,%edx /* round towards +oo */ + orl 4(%esp),%edx + andl $0xfbff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ + + frndint /* round */ + + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp + ret +END (__ceilf) +weak_alias (__ceilf, ceilf) diff --git a/sysdeps/i386/fpu/s_ceill.S b/sysdeps/i386/fpu/s_ceill.S new file mode 100644 index 0000000000..0128966ebe --- /dev/null +++ b/sysdeps/i386/fpu/s_ceill.S @@ -0,0 +1,33 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__ceill) + fldt 4(%esp) + subl $8,%esp + + 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 $0x0800,%edx /* round towards +oo */ + orl 4(%esp),%edx + andl $0xfbff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ + + frndint /* round */ + + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp + ret +END (__ceill) +weak_alias (__ceill, ceill) diff --git a/sysdeps/i386/fpu/s_cexp.S b/sysdeps/i386/fpu/s_cexp.S new file mode 100644 index 0000000000..61158d9540 --- /dev/null +++ b/sysdeps/i386/fpu/s_cexp.S @@ -0,0 +1,259 @@ +/* ix87 specific implementation of complex exponential function for double. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object) +huge_nan_null_null: + .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f + .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f + .double 0.0 +zero: .double 0.0 +infinity: + .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f + .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f + .double 0.0 + .byte 0, 0, 0, 0, 0, 0, 0, 0x80 + ASM_SIZE_DIRECTIVE(huge_nan_null_null) + + ASM_TYPE_DIRECTIVE(twopi,@object) +twopi: + .byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40 + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(twopi) + + ASM_TYPE_DIRECTIVE(l2e,@object) +l2e: + .byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(l2e) + + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%ecx) +#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) +#else +#define MO(op) op +#define MOX(op,x,f) op(,x,f) +#endif + + .text +ENTRY(__cexp) + fldl 8(%esp) /* x */ + fxam + fnstsw + fldl 16(%esp) /* y : x */ +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + movb %ah, %dh + andb $0x45, %ah + cmpb $0x05, %ah + je 1f /* Jump if real part is +-Inf */ + cmpb $0x01, %ah + je 2f /* Jump if real part is NaN */ + + fxam /* y : x */ + fnstsw + /* If the imaginary part is not finite we return NaN+i NaN, as + for the case when the real part is NaN. A test for +-Inf and + NaN would be necessary. But since we know the stack register + we applied `fxam' to is not empty we can simply use one test. + Check your FPU manual for more information. */ + andb $0x01, %ah + cmpb $0x01, %ah + je 20f + + /* We have finite numbers in the real and imaginary part. Do + the real work now. */ + fxch /* x : y */ + fldt MO(l2e) /* log2(e) : x : y */ + fmulp /* x * log2(e) : y */ + fld %st /* x * log2(e) : x * log2(e) : y */ + frndint /* int(x * log2(e)) : x * log2(e) : y */ + fsubr %st, %st(1) /* int(x * log2(e)) : frac(x * log2(e)) : y */ + fxch /* frac(x * log2(e)) : int(x * log2(e)) : y */ + f2xm1 /* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */ + faddl MO(one) /* 2^frac(x * log2(e)) : int(x * log2(e)) : y */ + fscale /* e^x : int(x * log2(e)) : y */ + fst %st(1) /* e^x : e^x : y */ + fxch %st(2) /* y : e^x : e^x */ + fsincos /* cos(y) : sin(y) : e^x : e^x */ + fnstsw + testl $0x400, %eax + jnz 7f + fmulp %st, %st(3) /* sin(y) : e^x : e^x * cos(y) */ + fmulp %st, %st(1) /* e^x * sin(y) : e^x * cos(y) */ + movl 4(%esp), %eax /* Pointer to memory for result. */ + fstpl 8(%eax) + fstpl (%eax) + ret $4 + + /* We have to reduce the argument to fsincos. */ + .align ALIGNARG(4) +7: fldt MO(twopi) /* 2*pi : y : e^x : e^x */ + fxch /* y : 2*pi : e^x : e^x */ +8: fprem1 /* y%(2*pi) : 2*pi : e^x : e^x */ + fnstsw + testl $0x400, %eax + jnz 8b + fstp %st(1) /* y%(2*pi) : e^x : e^x */ + fsincos /* cos(y) : sin(y) : e^x : e^x */ + fmulp %st, %st(3) + fmulp %st, %st(1) + movl 4(%esp), %eax /* Pointer to memory for result. */ + fstpl 8(%eax) + fstpl (%eax) + ret $4 + + /* The real part is +-inf. We must make further differences. */ + .align ALIGNARG(4) +1: fxam /* y : x */ + fnstsw + movb %ah, %dl + testb $0x01, %ah /* See above why 0x01 is usable here. */ + jne 3f + + + /* The real part is +-Inf and the imaginary part is finite. */ + andl $0x245, %edx + cmpb $0x40, %dl /* Imaginary part == 0? */ + je 4f /* Yes -> */ + + fxch /* x : y */ + shrl $5, %edx + fstp %st(0) /* y */ /* Drop the real part. */ + andl $16, %edx /* This puts the sign bit of the real part + in bit 4. So we can use it to index a + small array to select 0 or Inf. */ + fsincos /* cos(y) : sin(y) */ + fnstsw + testl $0x0400, %eax + jnz 5f + fldl MOX(huge_nan_null_null,%edx,1) + movl 4(%esp), %edx /* Pointer to memory for result. */ + fstl 8(%edx) + fstpl (%edx) + ftst + fnstsw + shll $23, %eax + andl $0x80000000, %eax + orl %eax, 4(%edx) + fstp %st(0) + ftst + fnstsw + shll $23, %eax + andl $0x80000000, %eax + orl %eax, 12(%edx) + fstp %st(0) + ret $4 + /* We must reduce the argument to fsincos. */ + .align ALIGNARG(4) +5: fldt MO(twopi) + fxch +6: fprem1 + fnstsw + testl $0x400, %eax + jnz 6b + fstp %st(1) + fsincos + fldl MOX(huge_nan_null_null,%edx,1) + movl 4(%esp), %edx /* Pointer to memory for result. */ + fstl 8(%edx) + fstpl (%edx) + ftst + fnstsw + shll $23, %eax + andl $0x80000000, %eax + orl %eax, 4(%edx) + fstp %st(0) + ftst + fnstsw + shll $23, %eax + andl $0x80000000, %eax + orl %eax, 12(%edx) + fstp %st(0) + ret $4 + + /* The real part is +-Inf and the imaginary part is +-0. So return + +-Inf+-0i. */ + .align ALIGNARG(4) +4: movl 4(%esp), %eax /* Pointer to memory for result. */ + fstpl 8(%eax) + shrl $5, %edx + fstp %st(0) + andl $16, %edx + fldl MOX(huge_nan_null_null,%edx,1) + fstpl (%eax) + ret $4 + + /* The real part is +-Inf, the imaginary is also is not finite. */ + .align ALIGNARG(4) +3: fstp %st(0) + fstp %st(0) /* <empty> */ + andb $0x45, %ah + andb $0x47, %dh + xorb %dh, %ah + jnz 30f + fldl MO(infinity) /* Raise invalid exception. */ + fmull MO(zero) + fstp %st(0) +30: movl %edx, %eax + shrl $5, %edx + shll $4, %eax + andl $16, %edx + andl $32, %eax + orl %eax, %edx + movl 4(%esp), %eax /* Pointer to memory for result. */ + + fldl MOX(huge_nan_null_null,%edx,1) + fldl MOX(huge_nan_null_null+8,%edx,1) + fxch + fstpl (%eax) + fstpl 8(%eax) + ret $4 + + /* The real part is NaN. */ + .align ALIGNARG(4) +20: fldl MO(infinity) /* Raise invalid exception. */ + fmull MO(zero) + fstp %st(0) +2: fstp %st(0) + fstp %st(0) + movl 4(%esp), %eax /* Pointer to memory for result. */ + fldl MO(huge_nan_null_null+8) + fstl (%eax) + fstpl 8(%eax) + ret $4 + +END(__cexp) +weak_alias (__cexp, cexp) diff --git a/sysdeps/i386/fpu/s_cexpf.S b/sysdeps/i386/fpu/s_cexpf.S new file mode 100644 index 0000000000..d6dcebcb23 --- /dev/null +++ b/sysdeps/i386/fpu/s_cexpf.S @@ -0,0 +1,255 @@ +/* ix87 specific implementation of complex exponential function for double. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object) +huge_nan_null_null: + .byte 0, 0, 0x80, 0x7f + .byte 0, 0, 0xc0, 0x7f + .float 0.0 +zero: .float 0.0 +infinity: + .byte 0, 0, 0x80, 0x7f + .byte 0, 0, 0xc0, 0x7f + .float 0.0 + .byte 0, 0, 0, 0x80 + ASM_SIZE_DIRECTIVE(huge_nan_null_null) + + ASM_TYPE_DIRECTIVE(twopi,@object) +twopi: + .byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40 + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(twopi) + + ASM_TYPE_DIRECTIVE(l2e,@object) +l2e: + .byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(l2e) + + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%ecx) +#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) +#else +#define MO(op) op +#define MOX(op,x,f) op(,x,f) +#endif + + .text +ENTRY(__cexpf) + flds 4(%esp) /* x */ + fxam + fnstsw + flds 8(%esp) /* y : x */ +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + movb %ah, %dh + andb $0x45, %ah + cmpb $0x05, %ah + je 1f /* Jump if real part is +-Inf */ + cmpb $0x01, %ah + je 2f /* Jump if real part is NaN */ + + fxam /* y : x */ + fnstsw + /* If the imaginary part is not finite we return NaN+i NaN, as + for the case when the real part is NaN. A test for +-Inf and + NaN would be necessary. But since we know the stack register + we applied `fxam' to is not empty we can simply use one test. + Check your FPU manual for more information. */ + andb $0x01, %ah + cmpb $0x01, %ah + je 20f + + /* We have finite numbers in the real and imaginary part. Do + the real work now. */ + fxch /* x : y */ + fldt MO(l2e) /* log2(e) : x : y */ + fmulp /* x * log2(e) : y */ + fld %st /* x * log2(e) : x * log2(e) : y */ + frndint /* int(x * log2(e)) : x * log2(e) : y */ + fsubr %st, %st(1) /* int(x * log2(e)) : frac(x * log2(e)) : y */ + fxch /* frac(x * log2(e)) : int(x * log2(e)) : y */ + f2xm1 /* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */ + faddl MO(one) /* 2^frac(x * log2(e)) : int(x * log2(e)) : y */ + fscale /* e^x : int(x * log2(e)) : y */ + fst %st(1) /* e^x : e^x : y */ + fxch %st(2) /* y : e^x : e^x */ + fsincos /* cos(y) : sin(y) : e^x : e^x */ + fnstsw + testl $0x400, %eax + jnz 7f + fmulp %st, %st(3) /* sin(y) : e^x : e^x * cos(y) */ + fmulp %st, %st(1) /* e^x * sin(y) : e^x * cos(y) */ + subl $8, %esp + fstps 4(%esp) + fstps (%esp) + popl %eax + popl %edx + ret + + /* We have to reduce the argument to fsincos. */ + .align ALIGNARG(4) +7: fldt MO(twopi) /* 2*pi : y : e^x : e^x */ + fxch /* y : 2*pi : e^x : e^x */ +8: fprem1 /* y%(2*pi) : 2*pi : e^x : e^x */ + fnstsw + testl $0x400, %eax + jnz 8b + fstp %st(1) /* y%(2*pi) : e^x : e^x */ + fsincos /* cos(y) : sin(y) : e^x : e^x */ + fmulp %st, %st(3) + fmulp %st, %st(1) + subl $8, %esp + fstps 4(%esp) + fstps (%esp) + popl %eax + popl %edx + ret + + /* The real part is +-inf. We must make further differences. */ + .align ALIGNARG(4) +1: fxam /* y : x */ + fnstsw + movb %ah, %dl + testb $0x01, %ah /* See above why 0x01 is usable here. */ + jne 3f + + + /* The real part is +-Inf and the imaginary part is finite. */ + andl $0x245, %edx + cmpb $0x40, %dl /* Imaginary part == 0? */ + je 4f /* Yes -> */ + + fxch /* x : y */ + shrl $6, %edx + fstp %st(0) /* y */ /* Drop the real part. */ + andl $8, %edx /* This puts the sign bit of the real part + in bit 3. So we can use it to index a + small array to select 0 or Inf. */ + fsincos /* cos(y) : sin(y) */ + fnstsw + testl $0x0400, %eax + jnz 5f + fxch + ftst + fnstsw + fstp %st(0) + shll $23, %eax + andl $0x80000000, %eax + orl MOX(huge_nan_null_null,%edx,1), %eax + movl MOX(huge_nan_null_null,%edx,1), %ecx + movl %eax, %edx + ftst + fnstsw + fstp %st(0) + shll $23, %eax + andl $0x80000000, %eax + orl %ecx, %eax + ret + /* We must reduce the argument to fsincos. */ + .align ALIGNARG(4) +5: fldt MO(twopi) + fxch +6: fprem1 + fnstsw + testl $0x400, %eax + jnz 6b + fstp %st(1) + fsincos + fxch + ftst + fnstsw + fstp %st(0) + shll $23, %eax + andl $0x80000000, %eax + orl MOX(huge_nan_null_null,%edx,1), %eax + movl MOX(huge_nan_null_null,%edx,1), %ecx + movl %eax, %edx + ftst + fnstsw + fstp %st(0) + shll $23, %eax + andl $0x80000000, %eax + orl %ecx, %eax + ret + + /* The real part is +-Inf and the imaginary part is +-0. So return + +-Inf+-0i. */ + .align ALIGNARG(4) +4: subl $4, %esp + fstps (%esp) + shrl $6, %edx + fstp %st(0) + andl $8, %edx + movl MOX(huge_nan_null_null,%edx,1), %eax + popl %edx + ret + + /* The real part is +-Inf, the imaginary is also is not finite. */ + .align ALIGNARG(4) +3: fstp %st(0) + fstp %st(0) /* <empty> */ + andb $0x45, %ah + andb $0x47, %dh + xorb %dh, %ah + jnz 30f + flds MO(infinity) /* Raise invalid exception. */ + fmuls MO(zero) + fstp %st(0) +30: movl %edx, %eax + shrl $6, %edx + shll $3, %eax + andl $8, %edx + andl $16, %eax + orl %eax, %edx + + movl MOX(huge_nan_null_null,%edx,1), %eax + movl MOX(huge_nan_null_null+4,%edx,1), %edx + ret + + /* The real part is NaN. */ + .align ALIGNARG(4) +20: flds MO(infinity) /* Raise invalid exception. */ + fmuls MO(zero) + fstp %st(0) +2: fstp %st(0) + fstp %st(0) + movl MO(huge_nan_null_null+4), %eax + movl %eax, %edx + ret + +END(__cexpf) +weak_alias (__cexpf, cexpf) diff --git a/sysdeps/i386/fpu/s_cexpl.S b/sysdeps/i386/fpu/s_cexpl.S new file mode 100644 index 0000000000..203a7ee779 --- /dev/null +++ b/sysdeps/i386/fpu/s_cexpl.S @@ -0,0 +1,262 @@ +/* ix87 specific implementation of complex exponential function for double. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object) +huge_nan_null_null: + .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f + .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f + .double 0.0 +zero: .double 0.0 +infinity: + .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f + .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f + .double 0.0 + .byte 0, 0, 0, 0, 0, 0, 0, 0x80 + ASM_SIZE_DIRECTIVE(huge_nan_null_null) + + ASM_TYPE_DIRECTIVE(twopi,@object) +twopi: + .byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40 + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(twopi) + + ASM_TYPE_DIRECTIVE(l2e,@object) +l2e: + .byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(l2e) + + ASM_TYPE_DIRECTIVE(one,@object) +one: .double 1.0 + ASM_SIZE_DIRECTIVE(one) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%ecx) +#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f) +#else +#define MO(op) op +#define MOX(op,x,f) op(,x,f) +#endif + + .text +ENTRY(__cexpl) + fldt 8(%esp) /* x */ + fxam + fnstsw + fldt 20(%esp) /* y : x */ +#ifdef PIC + call 1f +1: popl %ecx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx +#endif + movb %ah, %dh + andb $0x45, %ah + cmpb $0x05, %ah + je 1f /* Jump if real part is +-Inf */ + cmpb $0x01, %ah + je 2f /* Jump if real part is NaN */ + + fxam /* y : x */ + fnstsw + /* If the imaginary part is not finite we return NaN+i NaN, as + for the case when the real part is NaN. A test for +-Inf and + NaN would be necessary. But since we know the stack register + we applied `fxam' to is not empty we can simply use one test. + Check your FPU manual for more information. */ + andb $0x01, %ah + cmpb $0x01, %ah + je 20f + + /* We have finite numbers in the real and imaginary part. Do + the real work now. */ + fxch /* x : y */ + fldt MO(l2e) /* log2(e) : x : y */ + fmulp /* x * log2(e) : y */ + fld %st /* x * log2(e) : x * log2(e) : y */ + frndint /* int(x * log2(e)) : x * log2(e) : y */ + fsubr %st, %st(1) /* int(x * log2(e)) : frac(x * log2(e)) : y */ + fxch /* frac(x * log2(e)) : int(x * log2(e)) : y */ + f2xm1 /* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */ + faddl MO(one) /* 2^frac(x * log2(e)) : int(x * log2(e)) : y */ + fscale /* e^x : int(x * log2(e)) : y */ + fst %st(1) /* e^x : e^x : y */ + fxch %st(2) /* y : e^x : e^x */ + fsincos /* cos(y) : sin(y) : e^x : e^x */ + fnstsw + testl $0x400, %eax + jnz 7f + fmulp %st, %st(3) /* sin(y) : e^x : e^x * cos(y) */ + fmulp %st, %st(1) /* e^x * sin(y) : e^x * cos(y) */ + movl 4(%esp), %eax /* Pointer to memory for result. */ + fstpt 12(%eax) + fstpt (%eax) + ret $4 + + /* We have to reduce the argument to fsincos. */ + .align ALIGNARG(4) +7: fldt MO(twopi) /* 2*pi : y : e^x : e^x */ + fxch /* y : 2*pi : e^x : e^x */ +8: fprem1 /* y%(2*pi) : 2*pi : e^x : e^x */ + fnstsw + testl $0x400, %eax + jnz 8b + fstp %st(1) /* y%(2*pi) : e^x : e^x */ + fsincos /* cos(y) : sin(y) : e^x : e^x */ + fmulp %st, %st(3) + fmulp %st, %st(1) + movl 4(%esp), %eax /* Pointer to memory for result. */ + fstpt 12(%eax) + fstpt (%eax) + ret $4 + + /* The real part is +-inf. We must make further differences. */ + .align ALIGNARG(4) +1: fxam /* y : x */ + fnstsw + movb %ah, %dl + testb $0x01, %ah /* See above why 0x01 is usable here. */ + jne 3f + + + /* The real part is +-Inf and the imaginary part is finite. */ + andl $0x245, %edx + cmpb $0x40, %dl /* Imaginary part == 0? */ + je 4f /* Yes -> */ + + fxch /* x : y */ + shrl $5, %edx + fstp %st(0) /* y */ /* Drop the real part. */ + andl $16, %edx /* This puts the sign bit of the real part + in bit 4. So we can use it to index a + small array to select 0 or Inf. */ + fsincos /* cos(y) : sin(y) */ + fnstsw + testl $0x0400, %eax + jnz 5f + fldl MOX(huge_nan_null_null,%edx,1) + movl 4(%esp), %edx /* Pointer to memory for result. */ + fld %st + fstpt 12(%edx) + fstpt (%edx) + ftst + fnstsw + shll $7, %eax + andl $0x8000, %eax + orl %eax, 8(%edx) + fstp %st(0) + ftst + fnstsw + shll $7, %eax + andl $0x8000, %eax + orl %eax, 20(%edx) + fstp %st(0) + ret $4 + /* We must reduce the argument to fsincos. */ + .align ALIGNARG(4) +5: fldt MO(twopi) + fxch +6: fprem1 + fnstsw + testl $0x400, %eax + jnz 6b + fstp %st(1) + fsincos + fldl MOX(huge_nan_null_null,%edx,1) + movl 4(%esp), %edx /* Pointer to memory for result. */ + fld %st + fstpt 12(%edx) + fstpt (%edx) + ftst + fnstsw + shll $7, %eax + andl $0x8000, %eax + orl %eax, 8(%edx) + fstp %st(0) + ftst + fnstsw + shll $7, %eax + andl $0x8000, %eax + orl %eax, 20(%edx) + fstp %st(0) + ret $4 + + /* The real part is +-Inf and the imaginary part is +-0. So return + +-Inf+-0i. */ + .align ALIGNARG(4) +4: movl 4(%esp), %eax /* Pointer to memory for result. */ + fstpt 12(%eax) + shrl $5, %edx + fstp %st(0) + andl $16, %edx + fldl MOX(huge_nan_null_null,%edx,1) + fstpt (%eax) + ret $4 + + /* The real part is +-Inf, the imaginary is also is not finite. */ + .align ALIGNARG(4) +3: fstp %st(0) + fstp %st(0) /* <empty> */ + andb $0x45, %ah + andb $0x47, %dh + xorb %dh, %ah + jnz 30f + fldl MO(infinity) /* Raise invalid exception. */ + fmull MO(zero) + fstp %st(0) +30: movl %edx, %eax + shrl $5, %edx + shll $4, %eax + andl $16, %edx + andl $32, %eax + orl %eax, %edx + movl 4(%esp), %eax /* Pointer to memory for result. */ + + fldl MOX(huge_nan_null_null,%edx,1) + fldl MOX(huge_nan_null_null+8,%edx,1) + fxch + fstpt (%eax) + fstpt 12(%eax) + ret $4 + + /* The real part is NaN. */ + .align ALIGNARG(4) +20: fldl MO(infinity) /* Raise invalid exception. */ + fmull MO(zero) + fstp %st(0) +2: fstp %st(0) + fstp %st(0) + movl 4(%esp), %eax /* Pointer to memory for result. */ + fldl MO(huge_nan_null_null+8) + fld %st(0) + fstpt (%eax) + fstpt 12(%eax) + ret $4 + +END(__cexpl) +weak_alias (__cexpl, cexpl) diff --git a/sysdeps/i386/fpu/s_copysign.S b/sysdeps/i386/fpu/s_copysign.S new file mode 100644 index 0000000000..2520a94427 --- /dev/null +++ b/sysdeps/i386/fpu/s_copysign.S @@ -0,0 +1,20 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_copysign.S,v 1.4 1995/05/08 23:53:02 jtc Exp $") + +ENTRY(__copysign) + movl 16(%esp),%edx + movl 8(%esp),%eax + andl $0x80000000,%edx + andl $0x7fffffff,%eax + orl %edx,%eax + movl %eax,8(%esp) + fldl 4(%esp) + ret +END (__copysign) +weak_alias (__copysign, copysign) diff --git a/sysdeps/i386/fpu/s_copysignf.S b/sysdeps/i386/fpu/s_copysignf.S new file mode 100644 index 0000000000..57b1a6f119 --- /dev/null +++ b/sysdeps/i386/fpu/s_copysignf.S @@ -0,0 +1,20 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_copysignf.S,v 1.3 1995/05/08 23:53:25 jtc Exp $") + +ENTRY(__copysignf) + movl 8(%esp),%edx + movl 4(%esp),%eax + andl $0x80000000,%edx + andl $0x7fffffff,%eax + orl %edx,%eax + movl %eax,4(%esp) + flds 4(%esp) + ret +END (__copysignf) +weak_alias (__copysignf, copysignf) diff --git a/sysdeps/i386/fpu/s_copysignl.S b/sysdeps/i386/fpu/s_copysignl.S new file mode 100644 index 0000000000..2163e7b014 --- /dev/null +++ b/sysdeps/i386/fpu/s_copysignl.S @@ -0,0 +1,21 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__copysignl) + movl 24(%esp),%edx + movl 12(%esp),%eax + andl $0x8000,%edx + andl $0x7fff,%eax + orl %edx,%eax + movl %eax,12(%esp) + fldt 4(%esp) + ret +END (__copysignl) +weak_alias (__copysignl, copysignl) diff --git a/sysdeps/i386/fpu/s_cos.S b/sysdeps/i386/fpu/s_cos.S new file mode 100644 index 0000000000..ac8b1459d9 --- /dev/null +++ b/sysdeps/i386/fpu/s_cos.S @@ -0,0 +1,29 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_cos.S,v 1.5 1995/05/08 23:54:00 jtc Exp $") + +ENTRY(__cos) + fldl 4(%esp) + fcos + fnstsw %ax + testl $0x400,%eax + jnz 1f + ret + .align ALIGNARG(4) +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fnstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fcos + ret +END (__cos) +weak_alias (__cos, cos) diff --git a/sysdeps/i386/fpu/s_cosf.S b/sysdeps/i386/fpu/s_cosf.S new file mode 100644 index 0000000000..d8e8090639 --- /dev/null +++ b/sysdeps/i386/fpu/s_cosf.S @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_cosf.S,v 1.3 1995/05/08 23:55:16 jtc Exp $") + +/* A float's domain isn't large enough to require argument reduction. */ +ENTRY(__cosf) + flds 4(%esp) + fcos + ret +END (__cosf) +weak_alias (__cosf, cosf) diff --git a/sysdeps/i386/fpu/s_cosl.S b/sysdeps/i386/fpu/s_cosl.S new file mode 100644 index 0000000000..61c9010c99 --- /dev/null +++ b/sysdeps/i386/fpu/s_cosl.S @@ -0,0 +1,31 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__cosl) + fldt 4(%esp) + fcos + fnstsw %ax + testl $0x400,%eax + jnz 1f + ret + .align ALIGNARG(4) +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fnstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fcos + ret +END (__cosl) +weak_alias (__cosl, cosl) diff --git a/sysdeps/i386/fpu/s_exp2.S b/sysdeps/i386/fpu/s_exp2.S new file mode 100644 index 0000000000..778c0c0eb6 --- /dev/null +++ b/sysdeps/i386/fpu/s_exp2.S @@ -0,0 +1,37 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Adapted for exp2 by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + */ + +#include <machine/asm.h> + +ENTRY(__ieee754_exp2) + fldl 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 + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + fld %st + frndint /* int(x) */ + fsubr %st,%st(1) /* fract(x) */ + fxch + f2xm1 /* 2^(fract(x)) - 1 */ + fld1 + faddp /* 2^(fract(x)) */ + fscale /* e^x */ + fstp %st(1) + ret + +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +2: ret +END (__ieee754_exp2) diff --git a/sysdeps/i386/fpu/s_exp2f.S b/sysdeps/i386/fpu/s_exp2f.S new file mode 100644 index 0000000000..c2d1af1af1 --- /dev/null +++ b/sysdeps/i386/fpu/s_exp2f.S @@ -0,0 +1,37 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Adapted for exp2 by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + */ + +#include <machine/asm.h> + +ENTRY(__ieee754_exp2f) + 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 + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + fld %st + frndint /* int(x) */ + fsubr %st,%st(1) /* fract(x) */ + fxch + f2xm1 /* 2^(fract(x)) - 1 */ + fld1 + faddp /* 2^(fract(x)) */ + fscale /* e^x */ + fstp %st(1) + ret + +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +2: ret +END (__ieee754_exp2f) diff --git a/sysdeps/i386/fpu/s_exp2l.S b/sysdeps/i386/fpu/s_exp2l.S new file mode 100644 index 0000000000..fa1fdc9236 --- /dev/null +++ b/sysdeps/i386/fpu/s_exp2l.S @@ -0,0 +1,37 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Adapted for exp2 by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + */ + +#include <machine/asm.h> + +ENTRY(__ieee754_exp2l) + fldt 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 + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + fld %st + frndint /* int(x) */ + fsubr %st,%st(1) /* fract(x) */ + fxch + f2xm1 /* 2^(fract(x)) - 1 */ + fld1 + faddp /* 2^(fract(x)) */ + fscale /* e^x */ + fstp %st(1) + ret + +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +2: ret +END (__ieee754_exp2l) diff --git a/sysdeps/i386/fpu/s_expm1.S b/sysdeps/i386/fpu/s_expm1.S new file mode 100644 index 0000000000..92beaf0776 --- /dev/null +++ b/sysdeps/i386/fpu/s_expm1.S @@ -0,0 +1,88 @@ +/* ix87 specific implementation of exp(x)-1. + Copyright (C) 1996, 1997 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>. + Corrections by H.J. Lu (hjl@gnu.ai.mit.edu), 1997. + + 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(minus1,@object) +minus1: .double -1.0 + ASM_SIZE_DIRECTIVE(minus1) + 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 + movb $0x45, %ch + andb %ah, %ch + cmpb $0x40, %ch + je 3f // If +-0, jump. +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + cmpb $0x05, %ch + je 2f // If +-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 + fsubr %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) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fsubrl MO(one) // 1-2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fstp %st(1) // 1-2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fsubrp %st, %st(1) // 2^(log2(e)*x) + ret + +2: testl $0x200, %eax // Test sign. + jz 3f // If positive, jump. + fstp %st + fldl MO(minus1) // Set result to -1.0. +3: ret +END(__expm1) +weak_alias (__expm1, expm1) diff --git a/sysdeps/i386/fpu/s_expm1f.S b/sysdeps/i386/fpu/s_expm1f.S new file mode 100644 index 0000000000..45a60fe010 --- /dev/null +++ b/sysdeps/i386/fpu/s_expm1f.S @@ -0,0 +1,88 @@ +/* ix87 specific implementation of exp(x)-1. + Copyright (C) 1996, 1997 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>. + Corrections by H.J. Lu (hjl@gnu.ai.mit.edu), 1997. + + 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(minus1,@object) +minus1: .double -1.0 + ASM_SIZE_DIRECTIVE(minus1) + 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 + movb $0x45, %ch + andb %ah, %ch + cmpb $0x40, %ch + je 3f // If +-0, jump. +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + cmpb $0x05, %ch + je 2f // If +-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 + fsubr %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) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fsubrl MO(one) // 1-2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fstp %st(1) // 1-2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fsubrp %st, %st(1) // 2^(log2(e)*x) + ret + +2: testl $0x200, %eax // Test sign. + jz 3f // If positive, jump. + fstp %st + fldl MO(minus1) // Set result to -1.0. +3: ret +END(__expm1f) +weak_alias (__expm1f, expm1f) diff --git a/sysdeps/i386/fpu/s_expm1l.S b/sysdeps/i386/fpu/s_expm1l.S new file mode 100644 index 0000000000..13fa698cc7 --- /dev/null +++ b/sysdeps/i386/fpu/s_expm1l.S @@ -0,0 +1,88 @@ +/* ix87 specific implementation of exp(x)-1. + Copyright (C) 1996, 1997 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>. + Corrections by H.J. Lu (hjl@gnu.ai.mit.edu), 1997. + + 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(minus1,@object) +minus1: .double -1.0 + ASM_SIZE_DIRECTIVE(minus1) + 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 + movb $0x45, %ch + andb %ah, %ch + cmpb $0x40, %ch + je 3f // If +-0, jump. +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + cmpb $0x05, %ch + je 2f // If +-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 + fsubr %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) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fsubrl MO(one) // 1-2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fstp %st(1) // 1-2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x) + fsubrp %st, %st(1) // 2^(log2(e)*x) + ret + +2: testl $0x200, %eax // Test sign. + jz 3f // If positive, jump. + fstp %st + fldl MO(minus1) // Set result to -1.0. +3: ret +END(__expm1l) +weak_alias (__expm1l, expm1l) diff --git a/sysdeps/i386/fpu/s_fdim.S b/sysdeps/i386/fpu/s_fdim.S new file mode 100644 index 0000000000..7a1e2ffa49 --- /dev/null +++ b/sysdeps/i386/fpu/s_fdim.S @@ -0,0 +1,51 @@ +/* Compute positive difference. + Copyright (C) 1997, 1999 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__fdim) + fldl 4(%esp) // x + fldl 12(%esp) // x : y + + fucom %st(1) + fnstsw + sahf + jp 1f + + fsubrp %st, %st(1) + jc 2f + + fstp %st(0) + fldz + jmp 2f + +1: fxam + fnstsw + andb $0x45, %ah + cmpb $0x01, %ah + je 3f + + fxch +3: fstp %st(1) + +2: ret +END(__fdim) +weak_alias (__fdim, fdim) diff --git a/sysdeps/i386/fpu/s_fdimf.S b/sysdeps/i386/fpu/s_fdimf.S new file mode 100644 index 0000000000..7e58f6752b --- /dev/null +++ b/sysdeps/i386/fpu/s_fdimf.S @@ -0,0 +1,51 @@ +/* Compute positive difference. + Copyright (C) 1997, 1999 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__fdimf) + flds 4(%esp) // x + flds 8(%esp) // x : y + + fucom %st(1) + fnstsw + sahf + jp 1f + + fsubrp %st, %st(1) + jc 2f + + fstp %st(0) + fldz + jmp 2f + +1: fxam + fnstsw + andb $0x45, %ah + cmpb $0x01, %ah + je 3f + + fxch +3: fstp %st(1) + +2: ret +END(__fdimf) +weak_alias (__fdimf, fdimf) diff --git a/sysdeps/i386/fpu/s_fdiml.S b/sysdeps/i386/fpu/s_fdiml.S new file mode 100644 index 0000000000..7a8c18a5d3 --- /dev/null +++ b/sysdeps/i386/fpu/s_fdiml.S @@ -0,0 +1,51 @@ +/* Compute positive difference. + Copyright (C) 1997, 1999 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__fdiml) + fldt 4(%esp) // x + fldt 16(%esp) // x : y + + fucom %st(1) + fnstsw + sahf + jp 1f + + fsubrp %st, %st(1) + jc 2f + + fstp %st(0) + fldz + jmp 2f + +1: fxam + fnstsw + andb $0x45, %ah + cmpb $0x01, %ah + je 3f + + fxch +3: fstp %st(1) + +2: ret +END(__fdiml) +weak_alias (__fdiml, fdiml) diff --git a/sysdeps/i386/fpu/s_finite.S b/sysdeps/i386/fpu/s_finite.S new file mode 100644 index 0000000000..63c766a950 --- /dev/null +++ b/sysdeps/i386/fpu/s_finite.S @@ -0,0 +1,15 @@ +/* + * Written by Joe Keane <jgk@jgk.org>. + */ + +#include <machine/asm.h> + +ENTRY(__finite) + movl 8(%esp),%eax + movl $0xFFEFFFFF,%ecx + subl %eax,%ecx + xorl %ecx,%eax + shrl $31, %eax + ret +END (__finite) +weak_alias (__finite, finite) diff --git a/sysdeps/i386/fpu/s_finitef.S b/sysdeps/i386/fpu/s_finitef.S new file mode 100644 index 0000000000..dabb71a115 --- /dev/null +++ b/sysdeps/i386/fpu/s_finitef.S @@ -0,0 +1,15 @@ +/* + * Written by Joe Keane <jgk@jgk.org>. + */ + +#include <machine/asm.h> + +ENTRY(__finitef) + movl 4(%esp),%eax + movl $0xFF7FFFFF,%ecx + subl %eax,%ecx + xorl %ecx,%eax + shrl $31,%eax + ret +END (__finitef) +weak_alias (__finitef, finitef) diff --git a/sysdeps/i386/fpu/s_finitel.S b/sysdeps/i386/fpu/s_finitel.S new file mode 100644 index 0000000000..acc5ad4cd0 --- /dev/null +++ b/sysdeps/i386/fpu/s_finitel.S @@ -0,0 +1,14 @@ +/* + * Written by Joe Keane <jgk@jgk.org>. + */ + +#include <machine/asm.h> + +ENTRY(__finitel) + movl 12(%esp),%eax + orl $0xffff8000, %eax + incl %eax + shrl $31, %eax + ret +END (__finitel) +weak_alias (__finitel, finitel) diff --git a/sysdeps/i386/fpu/s_floor.S b/sysdeps/i386/fpu/s_floor.S new file mode 100644 index 0000000000..20a8660424 --- /dev/null +++ b/sysdeps/i386/fpu/s_floor.S @@ -0,0 +1,32 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_floor.S,v 1.4 1995/05/09 00:01:59 jtc Exp $") + +ENTRY(__floor) + fldl 4(%esp) + subl $8,%esp + + 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 $0x400,%edx /* round towards -oo */ + orl 4(%esp),%edx + andl $0xf7ff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ + + frndint /* round */ + + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp + ret +END (__floor) +weak_alias (__floor, floor) diff --git a/sysdeps/i386/fpu/s_floorf.S b/sysdeps/i386/fpu/s_floorf.S new file mode 100644 index 0000000000..eca93a2aa3 --- /dev/null +++ b/sysdeps/i386/fpu/s_floorf.S @@ -0,0 +1,32 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_floorf.S,v 1.3 1995/05/09 00:04:32 jtc Exp $") + +ENTRY(__floorf) + flds 4(%esp) + subl $8,%esp + + 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 $0x400,%edx /* round towards -oo */ + orl 4(%esp),%edx + andl $0xf7ff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ + + frndint /* round */ + + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp + ret +END (__floorf) +weak_alias (__floorf, floorf) diff --git a/sysdeps/i386/fpu/s_floorl.S b/sysdeps/i386/fpu/s_floorl.S new file mode 100644 index 0000000000..c2bf091d73 --- /dev/null +++ b/sysdeps/i386/fpu/s_floorl.S @@ -0,0 +1,33 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__floorl) + fldt 4(%esp) + subl $8,%esp + + 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 $0x400,%edx /* round towards -oo */ + orl 4(%esp),%edx + andl $0xf7ff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ + + frndint /* round */ + + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp + ret +END (__floorl) +weak_alias (__floorl, floorl) diff --git a/sysdeps/i386/fpu/s_fma.S b/sysdeps/i386/fpu/s_fma.S new file mode 100644 index 0000000000..2affafcbf5 --- /dev/null +++ b/sysdeps/i386/fpu/s_fma.S @@ -0,0 +1,31 @@ +/* Compute (X * Y) + Z as ternary operation. + Copyright (C) 1997, 1998 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__fma) + fldl 4(%esp) // x + fmull 12(%esp) // x * y + fldl 20(%esp) // z : x * y + faddp // (x * y) + z + ret +END(__fma) +weak_alias (__fma, fma) diff --git a/sysdeps/i386/fpu/s_fmaf.S b/sysdeps/i386/fpu/s_fmaf.S new file mode 100644 index 0000000000..03c6cc0f9e --- /dev/null +++ b/sysdeps/i386/fpu/s_fmaf.S @@ -0,0 +1,31 @@ +/* Compute (X * Y) + Z as ternary operation. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__fmaf) + flds 4(%esp) // x + fmuls 8(%esp) // x * y + flds 12(%esp) // z : x * y + faddp // (x * y) + z + ret +END(__fmaf) +weak_alias (__fmaf, fmaf) diff --git a/sysdeps/i386/fpu/s_fmal.S b/sysdeps/i386/fpu/s_fmal.S new file mode 100644 index 0000000000..c15fca8dab --- /dev/null +++ b/sysdeps/i386/fpu/s_fmal.S @@ -0,0 +1,32 @@ +/* Compute (X * Y) + Z as ternary operation. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__fmal) + fldt 4(%esp) // x + fldt 16(%esp) // x : y + fmulp // x * y + fldt 28(%esp) // z : x * y + faddp // (x * y) + z + ret +END(__fmal) +weak_alias (__fmal, fmal) diff --git a/sysdeps/i386/fpu/s_fmax.S b/sysdeps/i386/fpu/s_fmax.S new file mode 100644 index 0000000000..3dbaa139ad --- /dev/null +++ b/sysdeps/i386/fpu/s_fmax.S @@ -0,0 +1,44 @@ +/* Compute maximum of two numbers, regarding NaN as missing argument. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__fmax) + fldl 12(%esp) // y + fxam + fnstsw + fldl 4(%esp) // y : x + + andb $0x45, %ah + cmpb $0x01, %ah + je 1f // y == NaN + + fucom %st(1) + fnstsw + sahf + jnc 1f + + fxch %st(1) +1: fstp %st(1) + + ret +END(__fmax) +weak_alias (__fmax, fmax) diff --git a/sysdeps/i386/fpu/s_fmaxf.S b/sysdeps/i386/fpu/s_fmaxf.S new file mode 100644 index 0000000000..5dd94a43c4 --- /dev/null +++ b/sysdeps/i386/fpu/s_fmaxf.S @@ -0,0 +1,44 @@ +/* Compute maximum of two numbers, regarding NaN as missing argument. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__fmaxf) + flds 8(%esp) // y + fxam + fnstsw + flds 4(%esp) // y : x + + andb $0x45, %ah + cmpb $0x01, %ah + je 1f // y == NaN + + fucom %st(1) + fnstsw + sahf + jnc 1f + + fxch %st(1) +1: fstp %st(1) + + ret +END(__fmaxf) +weak_alias (__fmaxf, fmaxf) diff --git a/sysdeps/i386/fpu/s_fmaxl.S b/sysdeps/i386/fpu/s_fmaxl.S new file mode 100644 index 0000000000..d833bd1985 --- /dev/null +++ b/sysdeps/i386/fpu/s_fmaxl.S @@ -0,0 +1,44 @@ +/* Compute maximum of two numbers, regarding NaN as missing argument. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__fmaxl) + fldt 16(%esp) // y + fxam + fnstsw + fldt 4(%esp) // y : x + + andb $0x45, %ah + cmpb $0x01, %ah + je 1f // y == NaN + + fucom %st(1) + fnstsw + sahf + jnc 1f + + fxch %st(1) +1: fstp %st(1) + + ret +END(__fmaxl) +weak_alias (__fmaxl, fmaxl) diff --git a/sysdeps/i386/fpu/s_fmin.S b/sysdeps/i386/fpu/s_fmin.S new file mode 100644 index 0000000000..1ac46c7d21 --- /dev/null +++ b/sysdeps/i386/fpu/s_fmin.S @@ -0,0 +1,44 @@ +/* Compute minimum of two numbers, regarding NaN as missing argument. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__fmin) + fldl 4(%esp) // x + fldl 12(%esp) // x : y + + fxam + fnstsw + andb $0x45, %ah + cmpb $0x01, %ah + je 1f // y == NaN + + fucom %st(1) + fnstsw + sahf + jc 2f + +1: fxch %st(1) +2: fstp %st(1) + + ret +END(__fmin) +weak_alias (__fmin, fmin) diff --git a/sysdeps/i386/fpu/s_fminf.S b/sysdeps/i386/fpu/s_fminf.S new file mode 100644 index 0000000000..7dd346d2bc --- /dev/null +++ b/sysdeps/i386/fpu/s_fminf.S @@ -0,0 +1,44 @@ +/* Compute minimum of two numbers, regarding NaN as missing argument. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__fminf) + flds 4(%esp) // x + flds 8(%esp) // x : y + + fxam + fnstsw + andb $0x45, %ah + cmpb $0x01, %ah + je 1f // y == NaN + + fucom %st(1) + fnstsw + sahf + jc 2f + +1: fxch %st(1) +2: fstp %st(1) + + ret +END(__fminf) +weak_alias (__fminf, fminf) diff --git a/sysdeps/i386/fpu/s_fminl.S b/sysdeps/i386/fpu/s_fminl.S new file mode 100644 index 0000000000..883735d32c --- /dev/null +++ b/sysdeps/i386/fpu/s_fminl.S @@ -0,0 +1,44 @@ +/* Compute minimum of two numbers, regarding NaN as missing argument. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__fminl) + fldt 4(%esp) // x + fldt 16(%esp) // x : y + + fxam + fnstsw + andb $0x45, %ah + cmpb $0x01, %ah + je 1f // y == NaN + + fucom %st(1) + fnstsw + sahf + jc 2f + +1: fxch %st(1) +2: fstp %st(1) + + ret +END(__fminl) +weak_alias (__fminl, fminl) diff --git a/sysdeps/i386/fpu/s_frexp.S b/sysdeps/i386/fpu/s_frexp.S new file mode 100644 index 0000000000..6a05f26720 --- /dev/null +++ b/sysdeps/i386/fpu/s_frexp.S @@ -0,0 +1,82 @@ +/* ix87 specific frexp implementation for double. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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(two54,@object) +two54: .byte 0, 0, 0, 0, 0, 0, 0x50, 0x43 + ASM_SIZE_DIRECTIVE(two54) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__frexp) + movl 4(%esp), %ecx + movl 8(%esp), %eax + movl %eax, %edx + andl $0x7fffffff, %eax + orl %eax, %ecx + jz 1f + xorl %ecx, %ecx + cmpl $0x7ff00000, %eax + jae 1f + + cmpl $0x00100000, %eax + jae 2f + + fldl 4(%esp) +#ifdef PIC + call 3f +3: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-3b], %edx +#endif + fmull MO(two54) + movl $-54, %ecx + fstpl 4(%esp) + movl 8(%esp), %eax + movl %eax, %edx + andl $0x7fffffff, %eax + +2: shrl $20, %eax + andl $0x800fffff, %edx + subl $1022, %eax + orl $0x3fe00000, %edx + addl %eax, %ecx + movl %edx, 8(%esp) + + /* Store %ecx in the variable pointed to by the second argument, + get the factor from the stack and return. */ +1: movl 12(%esp), %eax + fldl 4(%esp) + movl %ecx, (%eax) + ret +END(__frexp) +weak_alias (__frexp, frexp) diff --git a/sysdeps/i386/fpu/s_frexpf.S b/sysdeps/i386/fpu/s_frexpf.S new file mode 100644 index 0000000000..1021b97aee --- /dev/null +++ b/sysdeps/i386/fpu/s_frexpf.S @@ -0,0 +1,80 @@ +/* ix87 specific frexp implementation for float. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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(two25,@object) +two25: .byte 0, 0, 0, 0x4c + ASM_SIZE_DIRECTIVE(two25) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__frexpf) + movl 4(%esp), %eax + xorl %ecx, %ecx + movl %eax, %edx + andl $0x7fffffff, %eax + jz 1f + cmpl $0x7f800000, %eax + jae 1f + + cmpl $0x00800000, %eax + jae 2f + + flds 4(%esp) +#ifdef PIC + call 3f +3: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-3b], %edx +#endif + fmuls MO(two25) + movl $-25, %ecx + fstps 4(%esp) + movl 4(%esp), %eax + movl %eax, %edx + andl $0x7fffffff, %eax + +2: shrl $23, %eax + andl $0x807fffff, %edx + subl $126, %eax + orl $0x3f000000, %edx + addl %eax, %ecx + movl %edx, 4(%esp) + + /* Store %ecx in the variable pointed to by the second argument, + get the factor from the stack and return. */ +1: movl 8(%esp), %eax + flds 4(%esp) + movl %ecx, (%eax) + ret +END(__frexpf) +weak_alias (__frexpf, frexpf) diff --git a/sysdeps/i386/fpu/s_frexpl.S b/sysdeps/i386/fpu/s_frexpl.S new file mode 100644 index 0000000000..e3019ced6d --- /dev/null +++ b/sysdeps/i386/fpu/s_frexpl.S @@ -0,0 +1,82 @@ +/* ix87 specific frexp implementation for long double. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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(two64,@object) +two64: .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43 + ASM_SIZE_DIRECTIVE(two64) + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__frexpl) + movl 4(%esp), %ecx + movl 12(%esp), %eax + orl 8(%esp), %ecx + movl %eax, %edx + andl $0x7fff, %eax + orl %eax, %ecx + jz 1f + xorl %ecx, %ecx + cmpl $0x7fff, %eax + je 1f + + cmpl $0, %eax + je 2f + + fldt 4(%esp) +#ifdef PIC + call 3f +3: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-3b], %edx +#endif + fmull MO(two64) /* It's not necessary to use a 80bit factor */ + movl $-64, %ecx + fstpt 4(%esp) + movl 12(%esp), %eax + movl %eax, %edx + andl $0x7fff, %eax + +2: andl $0x8000, %edx + subl $16382, %eax + orl $0x3ffe, %edx + addl %eax, %ecx + movl %edx, 12(%esp) + + /* Store %ecx in the variable pointed to by the second argument, + get the factor from the stack and return. */ +1: movl 16(%esp), %eax + fldt 4(%esp) + movl %ecx, (%eax) + ret +END(__frexpl) +weak_alias (__frexpl, frexpl) diff --git a/sysdeps/i386/fpu/s_ilogb.S b/sysdeps/i386/fpu/s_ilogb.S new file mode 100644 index 0000000000..36fb000c02 --- /dev/null +++ b/sysdeps/i386/fpu/s_ilogb.S @@ -0,0 +1,22 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_ilogb.S,v 1.5 1995/10/12 15:53:09 jtc Exp $") + +ENTRY(__ilogb) + fldl 4(%esp) + fxtract + pushl %eax + fstp %st + + fistpl (%esp) + fwait + popl %eax + + ret +END (__ilogb) +weak_alias (__ilogb, ilogb) diff --git a/sysdeps/i386/fpu/s_ilogbf.S b/sysdeps/i386/fpu/s_ilogbf.S new file mode 100644 index 0000000000..54f9d4647e --- /dev/null +++ b/sysdeps/i386/fpu/s_ilogbf.S @@ -0,0 +1,22 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_ilogbf.S,v 1.4 1995/10/22 20:32:43 pk Exp $") + +ENTRY(__ilogbf) + flds 4(%esp) + fxtract + pushl %eax + fstp %st + + fistpl (%esp) + fwait + popl %eax + + ret +END (__ilogbf) +weak_alias (__ilogbf, ilogbf) diff --git a/sysdeps/i386/fpu/s_ilogbl.S b/sysdeps/i386/fpu/s_ilogbl.S new file mode 100644 index 0000000000..1dad93abeb --- /dev/null +++ b/sysdeps/i386/fpu/s_ilogbl.S @@ -0,0 +1,23 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__ilogbl) + fldt 4(%esp) + fxtract + pushl %eax + fstp %st + + fistpl (%esp) + fwait + popl %eax + + ret +END (__ilogbl) +weak_alias (__ilogbl, ilogbl) diff --git a/sysdeps/i386/fpu/s_isinfl.c b/sysdeps/i386/fpu/s_isinfl.c new file mode 100644 index 0000000000..f07898fd1b --- /dev/null +++ b/sysdeps/i386/fpu/s_isinfl.c @@ -0,0 +1,36 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Change for long double by Ulrich Drepper <drepper@cygnus.com>. + * Intel i387 specific version. + * Public domain. + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: $"; +#endif + +/* + * isinfl(x) returns 1 if x is inf, -1 if x is -inf, else 0; + * no branching! + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + int __isinfl(long double x) +#else + int __isinfl(x) + long double x; +#endif +{ + int32_t se,hx,lx; + GET_LDOUBLE_WORDS(se,hx,lx,x); + /* This additional ^ 0x80000000 is necessary because in Intel's + internal representation of the implicit one is explicit. */ + lx |= (hx ^ 0x80000000) | ((se & 0x7fff) ^ 0x7fff); + lx |= -lx; + se &= 0x8000; + return ~(lx >> 31) & (1 - (se >> 14)); +} +weak_alias (__isinfl, isinfl) diff --git a/sysdeps/i386/fpu/s_isnanl.c b/sysdeps/i386/fpu/s_isnanl.c new file mode 100644 index 0000000000..6a74b956cc --- /dev/null +++ b/sysdeps/i386/fpu/s_isnanl.c @@ -0,0 +1,47 @@ +/* s_isnanl.c -- long double version for i387 of s_isnan.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 + +/* + * isnanl(x) returns 1 is x is nan, else 0; + * no branching! + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + int __isnanl(long double x) +#else + int __isnanl(x) + long double x; +#endif +{ + int32_t se,hx,lx; + GET_LDOUBLE_WORDS(se,hx,lx,x); + se = (se & 0x7fff) << 1; + /* The additional & 0x7fffffff is required because Intel's + extended format has the normally implicit 1 explicit + present. Sigh! */ + lx |= hx & 0x7fffffff; + se |= (u_int32_t)(lx|(-lx))>>31; + se = 0xfffe - se; + return (int)((u_int32_t)(se))>>16; +} +weak_alias (__isnanl, isnanl) diff --git a/sysdeps/i386/fpu/s_llrint.S b/sysdeps/i386/fpu/s_llrint.S new file mode 100644 index 0000000000..6109ec45c3 --- /dev/null +++ b/sysdeps/i386/fpu/s_llrint.S @@ -0,0 +1,34 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__llrint) + fldl 4(%esp) + subl $8, %esp + fistpll (%esp) + fwait + popl %eax + popl %edx + ret +END(__llrint) +weak_alias (__llrint, llrint) diff --git a/sysdeps/i386/fpu/s_llrintf.S b/sysdeps/i386/fpu/s_llrintf.S new file mode 100644 index 0000000000..c1e1c225c0 --- /dev/null +++ b/sysdeps/i386/fpu/s_llrintf.S @@ -0,0 +1,34 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__llrintf) + flds 4(%esp) + subl $8, %esp + fistpll (%esp) + fwait + popl %eax + popl %edx + ret +END(__llrintf) +weak_alias (__llrintf, llrintf) diff --git a/sysdeps/i386/fpu/s_llrintl.S b/sysdeps/i386/fpu/s_llrintl.S new file mode 100644 index 0000000000..d894897d5e --- /dev/null +++ b/sysdeps/i386/fpu/s_llrintl.S @@ -0,0 +1,34 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__llrintl) + fldt 4(%esp) + subl $8, %esp + fistpll (%esp) + fwait + popl %eax + popl %edx + ret +END(__llrintl) +weak_alias (__llrintl, llrintl) diff --git a/sysdeps/i386/fpu/s_log1p.S b/sysdeps/i386/fpu/s_log1p.S new file mode 100644 index 0000000000..10e8a36369 --- /dev/null +++ b/sysdeps/i386/fpu/s_log1p.S @@ -0,0 +1,62 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_log1p.S,v 1.7 1995/05/09 00:10:58 jtc Exp $") + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + /* The fyl2xp1 can only be used for values in + -1 + sqrt(2) / 2 <= x <= 1 - sqrt(2) / 2 + 0.29 is a safe value. + */ +limit: .double 0.29 +one: .double 1.0 + +/* + * Use the fyl2xp1 function when the argument is in the range -0.29 to 0.29, + * otherwise fyl2x with the needed extra computation. + */ + .text +ENTRY(__log1p) + fldln2 + + fldl 4(%esp) + +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + + fld %st + fabs +#ifdef PIC + fcompl limit@GOTOFF(%edx) +#else + fcompl limit +#endif + fnstsw + sahf + jc 2f + +#ifdef PIC + faddl one@GOTOFF(%edx) +#else + faddl one +#endif + fyl2x + ret + +2: fyl2xp1 + ret + +END (__log1p) +weak_alias (__log1p, log1p) diff --git a/sysdeps/i386/fpu/s_log1pf.S b/sysdeps/i386/fpu/s_log1pf.S new file mode 100644 index 0000000000..df9fdcbcfc --- /dev/null +++ b/sysdeps/i386/fpu/s_log1pf.S @@ -0,0 +1,62 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_log1pf.S,v 1.4 1995/05/09 00:13:05 jtc Exp $") + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + /* The fyl2xp1 can only be used for values in + -1 + sqrt(2) / 2 <= x <= 1 - sqrt(2) / 2 + 0.29 is a safe value. + */ +limit: .float 0.29 +one: .float 1.0 + +/* + * Use the fyl2xp1 function when the argument is in the range -0.29 to 0.29, + * otherwise fyl2x with the needed extra computation. + */ + .text +ENTRY(__log1pf) + fldln2 + + flds 4(%esp) + +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + + fld %st + fabs +#ifdef PIC + fcomps limit@GOTOFF(%edx) +#else + fcomps limit +#endif + fnstsw + sahf + jc 2f + +#ifdef PIC + fadds one@GOTOFF(%edx) +#else + fadds one +#endif + fyl2x + ret + +2: fyl2xp1 + ret + +END (__log1pf) +weak_alias (__log1pf, log1pf) diff --git a/sysdeps/i386/fpu/s_log1pl.S b/sysdeps/i386/fpu/s_log1pl.S new file mode 100644 index 0000000000..05a17b2831 --- /dev/null +++ b/sysdeps/i386/fpu/s_log1pl.S @@ -0,0 +1,68 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_log1p.S,v 1.7 1995/05/09 00:10:58 jtc Exp $") + +#ifdef __ELF__ + .section .rodata +#else + .text +#endif + .align ALIGNARG(4) + /* The fyl2xp1 can only be used for values in + -1 + sqrt(2) / 2 <= x <= 1 - sqrt(2) / 2 + 0.29 is a safe value. + */ +limit: .tfloat 0.29 + /* Please note: we use a double value here. Since 1.0 has + an exact representation this does not effect the accuracy + but it helps to optimize the code. */ +one: .double 1.0 + +/* + * Use the fyl2xp1 function when the argument is in the range -0.29 to 0.29, + * otherwise fyl2x with the needed extra computation. + */ + .text +ENTRY(__log1pl) + fldln2 + + fldt 4(%esp) + +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + + fld %st + fabs +#ifdef PIC + fldt limit@GOTOFF(%edx) +#else + fldt limit +#endif + fcompp + fnstsw + sahf + jnc 2f + +#ifdef PIC + faddl one@GOTOFF(%edx) +#else + faddl one +#endif + fyl2x + ret + +2: fyl2xp1 + ret + +END (__log1pl) +weak_alias (__log1pl, log1pl) diff --git a/sysdeps/i386/fpu/s_log2.S b/sysdeps/i386/fpu/s_log2.S new file mode 100644 index 0000000000..4632c96f67 --- /dev/null +++ b/sysdeps/i386/fpu/s_log2.S @@ -0,0 +1,59 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + */ + +#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) + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__log2) +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + fldl MO(one) + fldl 4(%esp) // x : 1 + fld %st // x : x : 1 + fsub %st(2), %st // x-1 : x : 1 + fld %st // x-1 : x-1 : x : 1 + fabs // |x-1| : x-1 : x : 1 + fcompl MO(limit) // x-1 : x : 1 + fnstsw // x-1 : x : 1 + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : 1 + fyl2xp1 // log(x) + ret + +2: fstp %st(0) // x : 1 + fyl2x // log(x) + ret +END (__log2) +weak_alias (__log2, log2) diff --git a/sysdeps/i386/fpu/s_log2f.S b/sysdeps/i386/fpu/s_log2f.S new file mode 100644 index 0000000000..bfdd0ef88e --- /dev/null +++ b/sysdeps/i386/fpu/s_log2f.S @@ -0,0 +1,59 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + */ + +#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) + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__log2f) +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + fldl MO(one) + flds 4(%esp) // x : 1 + fld %st // x : x : 1 + fsub %st(2), %st // x-1 : x : 1 + fld %st // x-1 : x-1 : x : 1 + fabs // |x-1| : x-1 : x : 1 + fcompl MO(limit) // x-1 : x : 1 + fnstsw // x-1 : x : 1 + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : 1 + fyl2xp1 // log(x) + ret + +2: fstp %st(0) // x : 1 + fyl2x // log(x) + ret +END (__log2f) +weak_alias (__log2f, log2f) diff --git a/sysdeps/i386/fpu/s_log2l.S b/sysdeps/i386/fpu/s_log2l.S new file mode 100644 index 0000000000..184981c4bb --- /dev/null +++ b/sysdeps/i386/fpu/s_log2l.S @@ -0,0 +1,59 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + */ + +#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) + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ + ASM_TYPE_DIRECTIVE(limit,@object) +limit: .double 0.29 + ASM_SIZE_DIRECTIVE(limit) + + +#ifdef PIC +#define MO(op) op##@GOTOFF(%edx) +#else +#define MO(op) op +#endif + + .text +ENTRY(__log2l) +#ifdef PIC + call 1f +1: popl %edx + addl $_GLOBAL_OFFSET_TABLE_+[.-1b], %edx +#endif + fldl MO(one) + fldt 4(%esp) // x : 1 + fld %st // x : x : 1 + fsub %st(2), %st // x-1 : x : 1 + fld %st // x-1 : x-1 : x : 1 + fabs // |x-1| : x-1 : x : 1 + fcompl MO(limit) // x-1 : x : 1 + fnstsw // x-1 : x : 1 + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : 1 + fyl2xp1 // log(x) + ret + +2: fstp %st(0) // x : 1 + fyl2x // log(x) + ret +END (__log2l) +weak_alias (__log2l, log2l) diff --git a/sysdeps/i386/fpu/s_logb.S b/sysdeps/i386/fpu/s_logb.S new file mode 100644 index 0000000000..f78c091c8a --- /dev/null +++ b/sysdeps/i386/fpu/s_logb.S @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_logb.S,v 1.4 1995/05/09 00:14:30 jtc Exp $") + +ENTRY(__logb) + fldl 4(%esp) + fxtract + fstp %st + ret +END (__logb) +weak_alias (__logb, logb) diff --git a/sysdeps/i386/fpu/s_logbf.S b/sysdeps/i386/fpu/s_logbf.S new file mode 100644 index 0000000000..91eb3d2925 --- /dev/null +++ b/sysdeps/i386/fpu/s_logbf.S @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_logbf.S,v 1.3 1995/05/09 00:15:12 jtc Exp $") + +ENTRY(__logbf) + flds 4(%esp) + fxtract + fstp %st + ret +END (__logbf) +weak_alias (__logbf, logbf) diff --git a/sysdeps/i386/fpu/s_logbl.S b/sysdeps/i386/fpu/s_logbl.S new file mode 100644 index 0000000000..5c9a9c1c9d --- /dev/null +++ b/sysdeps/i386/fpu/s_logbl.S @@ -0,0 +1,17 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__logbl) + fldt 4(%esp) + fxtract + fstp %st + ret +END (__logbl) +weak_alias (__logbl, logbl) diff --git a/sysdeps/i386/fpu/s_lrint.S b/sysdeps/i386/fpu/s_lrint.S new file mode 100644 index 0000000000..40956ab0d0 --- /dev/null +++ b/sysdeps/i386/fpu/s_lrint.S @@ -0,0 +1,33 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__lrint) + fldl 4(%esp) + subl $4, %esp + fistpl (%esp) + fwait + popl %eax + ret +END(__lrint) +weak_alias (__lrint, lrint) diff --git a/sysdeps/i386/fpu/s_lrintf.S b/sysdeps/i386/fpu/s_lrintf.S new file mode 100644 index 0000000000..2f49bdbabe --- /dev/null +++ b/sysdeps/i386/fpu/s_lrintf.S @@ -0,0 +1,33 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__lrintf) + flds 4(%esp) + subl $4, %esp + fistpl (%esp) + fwait + popl %eax + ret +END(__lrintf) +weak_alias (__lrintf, lrintf) diff --git a/sysdeps/i386/fpu/s_lrintl.S b/sysdeps/i386/fpu/s_lrintl.S new file mode 100644 index 0000000000..3a06c91b29 --- /dev/null +++ b/sysdeps/i386/fpu/s_lrintl.S @@ -0,0 +1,33 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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 <sysdep.h> + + .text +ENTRY(__lrintl) + fldt 4(%esp) + subl $4, %esp + fistpl (%esp) + fwait + popl %eax + ret +END(__lrintl) +weak_alias (__lrintl, lrintl) diff --git a/sysdeps/i386/fpu/s_nearbyint.S b/sysdeps/i386/fpu/s_nearbyint.S new file mode 100644 index 0000000000..65ce4f76a1 --- /dev/null +++ b/sysdeps/i386/fpu/s_nearbyint.S @@ -0,0 +1,25 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ +/* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. */ + +#include <machine/asm.h> + +ENTRY(__nearbyint) + fldl 4(%esp) + pushl %eax + pushl %ecx + fnstcw (%esp) + movl (%esp), %eax + andl $~0x20, %eax + movl %eax, 4(%esp) + fldcw 4(%esp) + frndint + fclex + fldcw (%esp) + popl %ecx + popl %eax + ret +END (__nearbyint) +weak_alias (__nearbyint, nearbyint) diff --git a/sysdeps/i386/fpu/s_nearbyintf.S b/sysdeps/i386/fpu/s_nearbyintf.S new file mode 100644 index 0000000000..090c631607 --- /dev/null +++ b/sysdeps/i386/fpu/s_nearbyintf.S @@ -0,0 +1,25 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ +/* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. */ + +#include <machine/asm.h> + +ENTRY(__nearbyintf) + flds 4(%esp) + pushl %eax + pushl %ecx + fnstcw (%esp) + movl (%esp), %eax + andl $~0x20, %eax + movl %eax, 4(%esp) + fldcw 4(%esp) + frndint + fclex + fldcw (%esp) + popl %ecx + popl %eax + ret +END (__nearbyintf) +weak_alias (__nearbyintf, nearbyintf) diff --git a/sysdeps/i386/fpu/s_nearbyintl.S b/sysdeps/i386/fpu/s_nearbyintl.S new file mode 100644 index 0000000000..2f60af8f18 --- /dev/null +++ b/sysdeps/i386/fpu/s_nearbyintl.S @@ -0,0 +1,25 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ +/* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. */ + +#include <machine/asm.h> + +ENTRY(__nearbyintl) + fldt 4(%esp) + pushl %eax + pushl %ecx + fnstcw (%esp) + movl (%esp), %eax + andl $~0x20, %eax + movl %eax, 4(%esp) + fldcw 4(%esp) + frndint + fclex + fldcw (%esp) + popl %ecx + popl %eax + ret +END (__nearbyintl) +weak_alias (__nearbyintl, nearbyintl) diff --git a/sysdeps/i386/fpu/s_nextafterl.c b/sysdeps/i386/fpu/s_nextafterl.c new file mode 100644 index 0000000000..4596c6b93c --- /dev/null +++ b/sysdeps/i386/fpu/s_nextafterl.c @@ -0,0 +1,102 @@ +/* s_nextafterl.c -- long double version of s_nextafter.c. + * Special version for i387. + * 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 + +/* IEEE functions + * nextafterl(x,y) + * return the next machine floating-point number of x in the + * direction toward y. + * Special cases: + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + long double __nextafterl(long double x, long double y) +#else + long double __nextafterl(x,y) + long double x,y; +#endif +{ + int32_t hx,hy,ix,iy; + u_int32_t lx,ly,esx,esy; + + GET_LDOUBLE_WORDS(esx,hx,lx,x); + GET_LDOUBLE_WORDS(esy,hy,ly,y); + ix = esx&0x7fff; /* |x| */ + iy = esy&0x7fff; /* |y| */ + + /* The additional &hx/&hy is required because Intel's extended format + has the normally implicit 1 explicit present. Sigh! */ + if(((ix==0x7fff)&&(((hx|lx)|-(hx|lx))&hx)>>31!=0) || /* x is nan */ + ((iy==0x7fff)&&(((hy|ly)|-(hy|ly))&hy)>>31!=0)) /* y is nan */ + return x+y; + if(x==y) return y; /* x=y, return y */ + if((ix|hx|lx)==0) { /* x == 0 */ + SET_LDOUBLE_WORDS(x,esx&0x8000,0,1);/* return +-minsubnormal */ + y = x*x; + if(y==x) return y; else return x; /* raise underflow flag */ + } + if(esx<0x8000) { /* x > 0 */ + if(ix>iy||((ix==iy) && (hx>hy||((hx==hy)&&(lx>ly))))) { + /* x > y, x -= ulp */ + if(lx==0) { + if (hx==0) esx -= 1; + hx -= 1; + } + lx -= 1; + } else { /* x < y, x += ulp */ + lx += 1; + if(lx==0) { + hx += 1; + if (hx==0) + esx += 1; + } + } + } else { /* x < 0 */ + if(esy>=0||(ix>iy||((ix==iy)&&(hx>hy||((hx==hy)&&(lx>ly)))))){ + /* x < y, x -= ulp */ + if(lx==0) { + if (hx==0) esx -= 1; + hx -= 1; + } + lx -= 1; + } else { /* x > y, x += ulp */ + lx += 1; + if(lx==0) { + hx += 1; + if (hx==0) esx += 1; + } + } + } + esy = esx&0x7fff; + if(esy==0x7fff) return x+x; /* overflow */ + if(esy==0) { /* underflow */ + y = x*x; + if(y!=x) { /* raise underflow flag */ + SET_LDOUBLE_WORDS(y,esx,hx,lx); + return y; + } + } + SET_LDOUBLE_WORDS(x,esx,hx,lx); + return x; +} +weak_alias (__nextafterl, nextafterl) diff --git a/sysdeps/i386/fpu/s_remquo.S b/sysdeps/i386/fpu/s_remquo.S new file mode 100644 index 0000000000..8ada191771 --- /dev/null +++ b/sysdeps/i386/fpu/s_remquo.S @@ -0,0 +1,36 @@ +/* + * Written by Ulrich Drepper <drepper@cygnus.com>. + * Based on e_remainder by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +ENTRY(__remquo) + fldl 12(%esp) + fldl 4(%esp) +1: fprem1 + fstsw %ax + sahf + jp 1b + fstp %st(1) + /* Compute the congruent of the quotient. */ + movl %eax, %ecx + shrl $8, %eax + shrl $12, %ecx + andl $4, %ecx + andl $3, %eax + orl %eax, %ecx + movl $0xef2960, %eax + shrl %cl, %eax + andl $3, %eax + movl 20(%esp), %ecx + movl 8(%esp), %edx + xorl 16(%esp), %edx + testl $0x80000000, %edx + jz 1f + negl %eax +1: movl %eax, (%ecx) + ret +END (__remquo) +weak_alias (__remquo, remquo) diff --git a/sysdeps/i386/fpu/s_remquof.S b/sysdeps/i386/fpu/s_remquof.S new file mode 100644 index 0000000000..f60aec9c46 --- /dev/null +++ b/sysdeps/i386/fpu/s_remquof.S @@ -0,0 +1,36 @@ +/* + * Written by Ulrich Drepper <drepper@cygnus.com>. + * Based on e_remainder by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +ENTRY(__remquof) + flds 8(%esp) + flds 4(%esp) +1: fprem1 + fstsw %ax + sahf + jp 1b + fstp %st(1) + /* Compute the congruent of the quotient. */ + movl %eax, %ecx + shrl $8, %eax + shrl $12, %ecx + andl $4, %ecx + andl $3, %eax + orl %eax, %ecx + movl $0xef2960, %eax + shrl %cl, %eax + andl $3, %eax + movl 12(%esp), %ecx + movl 4(%esp), %edx + xorl 8(%esp), %edx + testl $0x80000000, %edx + jz 1f + negl %eax +1: movl %eax, (%ecx) + ret +END (__remquof) +weak_alias (__remquof, remquof) diff --git a/sysdeps/i386/fpu/s_remquol.S b/sysdeps/i386/fpu/s_remquol.S new file mode 100644 index 0000000000..115d6e874b --- /dev/null +++ b/sysdeps/i386/fpu/s_remquol.S @@ -0,0 +1,36 @@ +/* + * Written by Ulrich Drepper <drepper@cygnus.com>. + * Based on e_remainder by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +ENTRY(__remquol) + fldt 16(%esp) + fldt 4(%esp) +1: fprem1 + fstsw %ax + sahf + jp 1b + fstp %st(1) + /* Compute the congruent of the quotient. */ + movl %eax, %ecx + shrl $8, %eax + shrl $12, %ecx + andl $4, %ecx + andl $3, %eax + orl %eax, %ecx + movl $0xef2960, %eax + shrl %cl, %eax + andl $3, %eax + movl 28(%esp), %ecx + movl 12(%esp), %edx + xorl 24(%esp), %edx + testl $0x8000, %edx + jz 1f + negl %eax +1: movl %eax, (%ecx) + ret +END (__remquol) +weak_alias (__remquol, remquol) diff --git a/sysdeps/i386/fpu/s_rint.S b/sysdeps/i386/fpu/s_rint.S new file mode 100644 index 0000000000..be36c5f0ca --- /dev/null +++ b/sysdeps/i386/fpu/s_rint.S @@ -0,0 +1,15 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_rint.S,v 1.4 1995/05/09 00:16:08 jtc Exp $") + +ENTRY(__rint) + fldl 4(%esp) + frndint + ret +END (__rint) +weak_alias (__rint, rint) diff --git a/sysdeps/i386/fpu/s_rintf.S b/sysdeps/i386/fpu/s_rintf.S new file mode 100644 index 0000000000..2b358c1cf1 --- /dev/null +++ b/sysdeps/i386/fpu/s_rintf.S @@ -0,0 +1,15 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_rintf.S,v 1.3 1995/05/09 00:17:22 jtc Exp $") + +ENTRY(__rintf) + flds 4(%esp) + frndint + ret +END (__rintf) +weak_alias (__rintf, rintf) diff --git a/sysdeps/i386/fpu/s_rintl.S b/sysdeps/i386/fpu/s_rintl.S new file mode 100644 index 0000000000..fd1ae6324e --- /dev/null +++ b/sysdeps/i386/fpu/s_rintl.S @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__rintl) + fldt 4(%esp) + frndint + ret +END (__rintl) +weak_alias (__rintl, rintl) diff --git a/sysdeps/i386/fpu/s_scalbln.c b/sysdeps/i386/fpu/s_scalbln.c new file mode 100644 index 0000000000..1009713fbc --- /dev/null +++ b/sysdeps/i386/fpu/s_scalbln.c @@ -0,0 +1,2 @@ +/* Nothing to do. This function is the same as scalbn. So we define an + alias. */ diff --git a/sysdeps/i386/fpu/s_scalblnf.c b/sysdeps/i386/fpu/s_scalblnf.c new file mode 100644 index 0000000000..5e558c3540 --- /dev/null +++ b/sysdeps/i386/fpu/s_scalblnf.c @@ -0,0 +1,2 @@ +/* Nothing to do. This function is the same as scalbnf. So we define an + alias. */ diff --git a/sysdeps/i386/fpu/s_scalblnl.c b/sysdeps/i386/fpu/s_scalblnl.c new file mode 100644 index 0000000000..cda2ec11c8 --- /dev/null +++ b/sysdeps/i386/fpu/s_scalblnl.c @@ -0,0 +1,2 @@ +/* Nothing to do. This function is the same as scalbnl. So we define an + alias. */ diff --git a/sysdeps/i386/fpu/s_scalbn.S b/sysdeps/i386/fpu/s_scalbn.S new file mode 100644 index 0000000000..ea9e25f094 --- /dev/null +++ b/sysdeps/i386/fpu/s_scalbn.S @@ -0,0 +1,19 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_scalbn.S,v 1.4 1995/05/09 00:19:06 jtc Exp $") + +ENTRY(__scalbn) + fildl 12(%esp) + fldl 4(%esp) + fscale + fstp %st(1) + ret +END (__scalbn) +weak_alias (__scalbn, scalbn) +strong_alias (__scalbn, __scalbln) +weak_alias (__scalbn, scalbln) diff --git a/sysdeps/i386/fpu/s_scalbnf.S b/sysdeps/i386/fpu/s_scalbnf.S new file mode 100644 index 0000000000..dc8cfb4296 --- /dev/null +++ b/sysdeps/i386/fpu/s_scalbnf.S @@ -0,0 +1,19 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_scalbnf.S,v 1.3 1995/05/09 00:19:59 jtc Exp $") + +ENTRY(__scalbnf) + fildl 8(%esp) + flds 4(%esp) + fscale + fstp %st(1) + ret +END (__scalbnf) +weak_alias (__scalbnf, scalbnf) +strong_alias (__scalbnf, __scalblnf) +weak_alias (__scalbnf, scalblnf) diff --git a/sysdeps/i386/fpu/s_scalbnl.S b/sysdeps/i386/fpu/s_scalbnl.S new file mode 100644 index 0000000000..295494b3d2 --- /dev/null +++ b/sysdeps/i386/fpu/s_scalbnl.S @@ -0,0 +1,20 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__scalbnl) + fildl 16(%esp) + fldt 4(%esp) + fscale + fstp %st(1) + ret +END (__scalbnl) +weak_alias (__scalbnl, scalbnl) +strong_alias (__scalbnl, __scalblnl) +weak_alias (__scalbnl, scalblnl) diff --git a/sysdeps/i386/fpu/s_significand.S b/sysdeps/i386/fpu/s_significand.S new file mode 100644 index 0000000000..4859b7ed71 --- /dev/null +++ b/sysdeps/i386/fpu/s_significand.S @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_significand.S,v 1.4 1995/05/09 00:21:47 jtc Exp $") + +ENTRY(__significand) + fldl 4(%esp) + fxtract + fstp %st(1) + ret +END (__significand) +weak_alias (__significand, significand) diff --git a/sysdeps/i386/fpu/s_significandf.S b/sysdeps/i386/fpu/s_significandf.S new file mode 100644 index 0000000000..3a2de97759 --- /dev/null +++ b/sysdeps/i386/fpu/s_significandf.S @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_significandf.S,v 1.3 1995/05/09 00:24:07 jtc Exp $") + +ENTRY(__significandf) + flds 4(%esp) + fxtract + fstp %st(1) + ret +END (__significandf) +weak_alias (__significandf, significandf) diff --git a/sysdeps/i386/fpu/s_significandl.S b/sysdeps/i386/fpu/s_significandl.S new file mode 100644 index 0000000000..e3a69cba73 --- /dev/null +++ b/sysdeps/i386/fpu/s_significandl.S @@ -0,0 +1,17 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__significandl) + fldt 4(%esp) + fxtract + fstp %st(0) + ret +END (__significandl) +weak_alias (__significandl, significandl) diff --git a/sysdeps/i386/fpu/s_sin.S b/sysdeps/i386/fpu/s_sin.S new file mode 100644 index 0000000000..eb22d7e98b --- /dev/null +++ b/sysdeps/i386/fpu/s_sin.S @@ -0,0 +1,29 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_sin.S,v 1.5 1995/05/09 00:25:54 jtc Exp $") + +ENTRY(__sin) + fldl 4(%esp) + fsin + fnstsw %ax + testl $0x400,%eax + jnz 1f + ret + .align ALIGNARG(4) +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fnstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fsin + ret +END (__sin) +weak_alias (__sin, sin) diff --git a/sysdeps/i386/fpu/s_sincos.S b/sysdeps/i386/fpu/s_sincos.S new file mode 100644 index 0000000000..fe99f42d18 --- /dev/null +++ b/sysdeps/i386/fpu/s_sincos.S @@ -0,0 +1,48 @@ +/* Compute sine and cosine of argument. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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> + +ENTRY(__sincos) + fldl 4(%esp) + fsincos + movl 12(%esp), %ecx + movl 16(%esp), %edx + fnstsw %ax + testl $0x400,%eax + jnz 1f + fstpl (%edx) + fstpl (%ecx) + ret + .align ALIGNARG(4) +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fnstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fsincos + fstpl (%edx) + fstpl (%ecx) + ret +END(__sincos) +weak_alias(__sincos, sincos) diff --git a/sysdeps/i386/fpu/s_sincosf.S b/sysdeps/i386/fpu/s_sincosf.S new file mode 100644 index 0000000000..5bb13f3c33 --- /dev/null +++ b/sysdeps/i386/fpu/s_sincosf.S @@ -0,0 +1,48 @@ +/* Compute sine and cosine of argument. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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> + +ENTRY(__sincosf) + flds 4(%esp) + fsincos + movl 8(%esp), %ecx + movl 12(%esp), %edx + fnstsw %ax + testl $0x400,%eax + jnz 1f + fstps (%edx) + fstps (%ecx) + ret + .align ALIGNARG(4) +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fnstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fsincos + fstps (%edx) + fstps (%ecx) + ret +END(__sincosf) +weak_alias(__sincosf, sincosf) diff --git a/sysdeps/i386/fpu/s_sincosl.S b/sysdeps/i386/fpu/s_sincosl.S new file mode 100644 index 0000000000..8b6694f09f --- /dev/null +++ b/sysdeps/i386/fpu/s_sincosl.S @@ -0,0 +1,48 @@ +/* Compute sine and cosine of argument. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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> + +ENTRY(__sincosl) + fldt 4(%esp) + fsincos + movl 16(%esp), %ecx + movl 20(%esp), %edx + fnstsw %ax + testl $0x400,%eax + jnz 1f + fstpt (%edx) + fstpt (%ecx) + ret + .align ALIGNARG(4) +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fnstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fsincos + fstpt (%edx) + fstpt (%ecx) + ret +END(__sincosl) +weak_alias(__sincosl, sincosl) diff --git a/sysdeps/i386/fpu/s_sinf.S b/sysdeps/i386/fpu/s_sinf.S new file mode 100644 index 0000000000..a010d60f5e --- /dev/null +++ b/sysdeps/i386/fpu/s_sinf.S @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_sinf.S,v 1.3 1995/05/09 00:27:53 jtc Exp $") + +/* A float's domain isn't large enough to require argument reduction. */ +ENTRY(__sinf) + flds 4(%esp) + fsin + ret +END (__sinf) +weak_alias (__sinf, sinf) diff --git a/sysdeps/i386/fpu/s_sinl.S b/sysdeps/i386/fpu/s_sinl.S new file mode 100644 index 0000000000..3e215de5e1 --- /dev/null +++ b/sysdeps/i386/fpu/s_sinl.S @@ -0,0 +1,31 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__sinl) + fldt 4(%esp) + fsin + fnstsw %ax + testl $0x400,%eax + jnz 1f + ret + .align ALIGNARG(4) +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fnstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fsin + ret +END (__sinl) +weak_alias (__sinl, sinl) diff --git a/sysdeps/i386/fpu/s_tan.S b/sysdeps/i386/fpu/s_tan.S new file mode 100644 index 0000000000..7b3547af4c --- /dev/null +++ b/sysdeps/i386/fpu/s_tan.S @@ -0,0 +1,30 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_tan.S,v 1.5 1995/05/09 00:30:00 jtc Exp $") + +ENTRY(__tan) + fldl 4(%esp) + fptan + fnstsw %ax + testl $0x400,%eax + jnz 1f + fstp %st(0) + ret +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fptan + fstp %st(0) + ret +END (__tan) +weak_alias (__tan, tan) diff --git a/sysdeps/i386/fpu/s_tanf.S b/sysdeps/i386/fpu/s_tanf.S new file mode 100644 index 0000000000..7a7509119b --- /dev/null +++ b/sysdeps/i386/fpu/s_tanf.S @@ -0,0 +1,17 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: s_tanf.S,v 1.3 1995/05/09 00:31:09 jtc Exp $") + +/* A float's domain isn't large enough to require argument reduction. */ +ENTRY(__tanf) + flds 4(%esp) + fptan + fstp %st(0) + ret +END (__tanf) +weak_alias (__tanf, tanf) diff --git a/sysdeps/i386/fpu/s_tanl.S b/sysdeps/i386/fpu/s_tanl.S new file mode 100644 index 0000000000..f2bdd6a605 --- /dev/null +++ b/sysdeps/i386/fpu/s_tanl.S @@ -0,0 +1,32 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <machine/asm.h> + +RCSID("$NetBSD: $") + +ENTRY(__tanl) + fldt 4(%esp) + fptan + fnstsw %ax + testl $0x400,%eax + jnz 1f + fstp %st(0) + ret +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fptan + fstp %st(0) + ret +END (__tanl) +weak_alias (__tanl, tanl) diff --git a/sysdeps/i386/fpu/s_trunc.S b/sysdeps/i386/fpu/s_trunc.S new file mode 100644 index 0000000000..3100d716a9 --- /dev/null +++ b/sysdeps/i386/fpu/s_trunc.S @@ -0,0 +1,36 @@ +/* Truncate double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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> + +ENTRY(__trunc) + fldl 4(%esp) + subl $8, %esp + fstcw 4(%esp) + movl $0xc00, %edx + orl 4(%esp), %edx + movl %edx, (%esp) + fldcw (%esp) + frndint + fldcw 4(%esp) + addl $8, %esp + ret +END(__trunc) +weak_alias (__trunc, trunc) diff --git a/sysdeps/i386/fpu/s_truncf.S b/sysdeps/i386/fpu/s_truncf.S new file mode 100644 index 0000000000..275e5f714c --- /dev/null +++ b/sysdeps/i386/fpu/s_truncf.S @@ -0,0 +1,36 @@ +/* Truncate float value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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> + +ENTRY(__truncf) + flds 4(%esp) + subl $8, %esp + fstcw 4(%esp) + movl $0xc00, %edx + orl 4(%esp), %edx + movl %edx, (%esp) + fldcw (%esp) + frndint + fldcw 4(%esp) + addl $8, %esp + ret +END(__truncf) +weak_alias (__truncf, truncf) diff --git a/sysdeps/i386/fpu/s_truncl.S b/sysdeps/i386/fpu/s_truncl.S new file mode 100644 index 0000000000..4c0bb0ce53 --- /dev/null +++ b/sysdeps/i386/fpu/s_truncl.S @@ -0,0 +1,36 @@ +/* Truncate long double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + 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> + +ENTRY(__truncl) + fldt 4(%esp) + subl $8, %esp + fstcw 4(%esp) + movl $0xc00, %edx + orl 4(%esp), %edx + movl %edx, (%esp) + fldcw (%esp) + frndint + fldcw 4(%esp) + addl $8, %esp + ret +END(__truncl) +weak_alias (__truncl, truncl) diff --git a/sysdeps/i386/fpu/t_exp.c b/sysdeps/i386/fpu/t_exp.c new file mode 100644 index 0000000000..fd37963b05 --- /dev/null +++ b/sysdeps/i386/fpu/t_exp.c @@ -0,0 +1 @@ +/* Empty. Not needed. */ |