1/* ix87 specific implementation of arcsinh.
2   Copyright (C) 1996 Free Software Foundation, Inc.
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
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.
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
14   Lesser General Public License for more details.
15
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.  */
20
21#include <machine/asm.h>
22
23#ifdef __ELF__
24	.section .rodata
25#else
26	.text
27#endif
28
29	.align ALIGNARG(4)
30	ASM_TYPE_DIRECTIVE(one,@object)
31one:	.double 1.0
32	ASM_SIZE_DIRECTIVE(one)
33	ASM_TYPE_DIRECTIVE(limit,@object)
34limit:	.double 0.29
35	ASM_SIZE_DIRECTIVE(limit)
36
37#ifdef PIC
38#define MO(op) op##@GOTOFF(%edx)
39#else
40#define MO(op) op
41#endif
42
43	.text
44ENTRY(__ieee754_acosh)
45	movl	8(%esp), %ecx
46	cmpl	$0x3ff00000, %ecx
47	jl	5f			// < 1 => invalid
48	fldln2				// log(2)
49	fldl	4(%esp)			// x : log(2)
50	cmpl	$0x41b00000, %ecx
51	ja	3f			// x > 2^28
52#ifdef	PIC
53	call	1f
541:	popl	%edx
55	addl	$_GLOBAL_OFFSET_TABLE_+[.-1b], %edx
56#endif
57	cmpl	$0x40000000, %ecx
58	ja	4f			// x > 2
59
60	// 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
61	fsubl	MO(one)			// x-1 : log(2)
62	fld	%st			// x-1 : x-1 : log(2)
63	fmul	%st(1)			// (x-1)^2 : x-1 : log(2)
64	fadd	%st(1)			// x-1+(x-1)^2 : x-1 : log(2)
65	fadd	%st(1)			// 2*(x-1)+(x-1)^2 : x-1 : log(2)
66	fsqrt				// sqrt(2*(x-1)+(x-1)^2) : x-1 : log(2)
67	faddp				// x-1+sqrt(2*(x-1)+(x-1)^2) : log(2)
68	fcoml	MO(limit)
69	fnstsw
70	sahf
71	ja	2f
72	fyl2xp1				// log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
73	ret
74
752:	faddl	MO(one)			// x+sqrt(2*(x-1)+(x-1)^2) : log(2)
76	fyl2x				// log(x+sqrt(2*(x-1)+(x-1)^2))
77	ret
78
79	// x > 2^28 => y = log(x) + log(2)
80	.align ALIGNARG(4)
813:	fyl2x				// log(x)
82	fldln2				// log(2) : log(x)
83	faddp				// log(x)+log(2)
84	ret
85
86	// 2^28 > x > 2 => y = log(2*x - 1/(x+sqrt(x*x-1)))
87	.align ALIGNARG(4)
884:	fld	%st			// x : x : log(2)
89	fadd	%st, %st(1)		// x : 2*x : log(2)
90	fld	%st			// x : x : 2*x : log(2)
91	fmul	%st(1)			// x^2 : x : 2*x : log(2)
92	fsubl	MO(one)			// x^2-1 : x : 2*x : log(2)
93	fsqrt				// sqrt(x^2-1) : x : 2*x : log(2)
94	faddp				// x+sqrt(x^2-1) : 2*x : log(2)
95	fdivrl	MO(one)			// 1/(x+sqrt(x^2-1)) : 2*x : log(2)
96	fsubrp				// 2*x+1/(x+sqrt(x^2)-1) : log(2)
97	fyl2x				// log(2*x+1/(x+sqrt(x^2-1)))
98	ret
99
100	// x < 1 => NaN
101	.align ALIGNARG(4)
1025:	fldz
103	fdiv	%st, %st(0)
104	ret
105END(__ieee754_acosh)
106