[thirdparty/glibc.git] / sysdeps / i386 / fpu / e_acoshl.S

/* ix87 specific implementation of arcsinh.
   Copyright (C) 1996, 1997, 2005, 2011 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 Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 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
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with the GNU C Library; 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.cst8,"aM",@progbits,8
#else
	.text
#endif

	.p2align 3
	/* 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
	LOAD_PIC_REG (dx)
#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)
strong_alias (__ieee754_acoshl, __acoshl_finite)
Commit	Line	Data
d38cd08c	1	/* ix87 specific implementation of arcsinh.
0ac5ae23	2	Copyright (C) 1996, 1997, 2005, 2011 Free Software Foundation, Inc.
d38cd08c UD	3	This file is part of the GNU C Library.
	4	Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
	5
	6	The GNU C Library is free software; you can redistribute it and/or
41bdb6e2 AJ	7	modify it under the terms of the GNU Lesser General Public
	8	License as published by the Free Software Foundation; either
	9	version 2.1 of the License, or (at your option) any later version.
d38cd08c UD	10
	11	The GNU C Library is distributed in the hope that it will be useful,
	12	but WITHOUT ANY WARRANTY; without even the implied warranty of
	13	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
41bdb6e2	14	Lesser General Public License for more details.
d38cd08c	15
41bdb6e2 AJ	16	You should have received a copy of the GNU Lesser General Public
	17	License along with the GNU C Library; if not, write to the Free
	18	Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
	19	02111-1307 USA. */
d38cd08c UD	20
	21	#include <machine/asm.h>
	22
	23	#ifdef __ELF__
0ac5ae23	24	.section .rodata.cst8,"aM",@progbits,8
d38cd08c UD	25	#else
	26	.text
	27	#endif
	28
0ac5ae23	29	.p2align 3
d38cd08c UD	30	/* Please note that we use double value for 1.0. This number
	31	has an exact representation and so we don't get accuracy
	32	problems. The advantage is that the code is simpler. */
	33	ASM_TYPE_DIRECTIVE(one,@object)
	34	one: .double 1.0
	35	ASM_SIZE_DIRECTIVE(one)
	36	/* It is not important that this constant is precise. It is only
	37	a value which is known to be on the safe side for using the
	38	fyl2xp1 instruction. */
	39	ASM_TYPE_DIRECTIVE(limit,@object)
	40	limit: .double 0.29
	41	ASM_SIZE_DIRECTIVE(limit)
	42
	43	#ifdef PIC
	44	#define MO(op) op##@GOTOFF(%edx)
	45	#else
	46	#define MO(op) op
	47	#endif
	48
	49	.text
	50	ENTRY(__ieee754_acoshl)
	51	movl 12(%esp), %ecx
ceb2d9aa	52	andl $0xffff, %ecx
d38cd08c UD	53	cmpl $0x3fff, %ecx
	54	jl 5f // < 1 => invalid
	55	fldln2 // log(2)
	56	fldt 4(%esp) // x : log(2)
	57	cmpl $0x4020, %ecx
	58	ja 3f // x > 2^34
	59	#ifdef PIC
fee732e5	60	LOAD_PIC_REG (dx)
d38cd08c UD	61	#endif
	62	cmpl $0x4000, %ecx
	63	ja 4f // x > 2
	64
	65	// 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
	66	fsubl MO(one) // x-1 : log(2)
	67	fld %st // x-1 : x-1 : log(2)
	68	fmul %st(1) // (x-1)^2 : x-1 : log(2)
	69	fadd %st(1) // x-1+(x-1)^2 : x-1 : log(2)
	70	fadd %st(1) // 2*(x-1)+(x-1)^2 : x-1 : log(2)
	71	fsqrt // sqrt(2*(x-1)+(x-1)^2) : x-1 : log(2)
	72	faddp // x-1+sqrt(2*(x-1)+(x-1)^2) : log(2)
	73	fcoml MO(limit)
	74	fnstsw
	75	sahf
	76	ja 2f
	77	fyl2xp1 // log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
	78	ret
	79
	80	2: faddl MO(one) // x+sqrt(2*(x-1)+(x-1)^2) : log(2)
	81	fyl2x // log(x+sqrt(2*(x-1)+(x-1)^2))
	82	ret
	83
	84	// x > 2^34 => y = log(x) + log(2)
	85	.align ALIGNARG(4)
	86	3: fyl2x // log(x)
	87	fldln2 // log(2) : log(x)
	88	faddp // log(x)+log(2)
	89	ret
	90
	91	// 2^34 > x > 2 => y = log(2x - 1/(x+sqrt(xx-1)))
	92	.align ALIGNARG(4)
	93	4: fld %st // x : x : log(2)
	94	fadd %st, %st(1) // x : 2*x : log(2)
	95	fld %st // x : x : 2*x : log(2)
	96	fmul %st(1) // x^2 : x : 2*x : log(2)
	97	fsubl MO(one) // x^2-1 : x : 2*x : log(2)
	98	fsqrt // sqrt(x^2-1) : x : 2*x : log(2)
	99	faddp // x+sqrt(x^2-1) : 2*x : log(2)
	100	fdivrl MO(one) // 1/(x+sqrt(x^2-1)) : 2*x : log(2)
	101	fsubrp // 2*x+1/(x+sqrt(x^2)-1) : log(2)
	102	fyl2x // log(2*x+1/(x+sqrt(x^2-1)))
	103	ret
	104
	105	// x < 1 => NaN
	106	.align ALIGNARG(4)
	107	5: fldz
	108	fdiv %st, %st(0)
	109	ret
	110	END(__ieee754_acoshl)
0ac5ae23	111	strong_alias (__ieee754_acoshl, __acoshl_finite)