1/* ix87 specific implementation of arcsinh.
2   Copyright (C) 1996, 1997 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	/* 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)
34one:	.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)
40limit:	.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
50ENTRY(__ieee754_acoshl)
51	movl	12(%esp), %ecx
52	andl	$0xffff, %ecx
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
60	call	1f
611:	popl	%edx
62	addl	$_GLOBAL_OFFSET_TABLE_+[.-1b], %edx
63#endif
64	cmpl	$0x4000, %ecx
65	ja	4f			// x > 2
66
67	// 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
68	fsubl	MO(one)			// x-1 : log(2)
69	fld	%st			// x-1 : x-1 : log(2)
70	fmul	%st(1)			// (x-1)^2 : x-1 : log(2)
71	fadd	%st(1)			// x-1+(x-1)^2 : x-1 : log(2)
72	fadd	%st(1)			// 2*(x-1)+(x-1)^2 : x-1 : log(2)
73	fsqrt				// sqrt(2*(x-1)+(x-1)^2) : x-1 : log(2)
74	faddp				// x-1+sqrt(2*(x-1)+(x-1)^2) : log(2)
75	fcoml	MO(limit)
76	fnstsw
77	sahf
78	ja	2f
79	fyl2xp1				// log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
80	ret
81
822:	faddl	MO(one)			// x+sqrt(2*(x-1)+(x-1)^2) : log(2)
83	fyl2x				// log(x+sqrt(2*(x-1)+(x-1)^2))
84	ret
85
86	// x > 2^34 => y = log(x) + log(2)
87	.align ALIGNARG(4)
883:	fyl2x				// log(x)
89	fldln2				// log(2) : log(x)
90	faddp				// log(x)+log(2)
91	ret
92
93	// 2^34 > x > 2 => y = log(2*x - 1/(x+sqrt(x*x-1)))
94	.align ALIGNARG(4)
954:	fld	%st			// x : x : log(2)
96	fadd	%st, %st(1)		// x : 2*x : log(2)
97	fld	%st			// x : x : 2*x : log(2)
98	fmul	%st(1)			// x^2 : x : 2*x : log(2)
99	fsubl	MO(one)			// x^2-1 : x : 2*x : log(2)
100	fsqrt				// sqrt(x^2-1) : x : 2*x : log(2)
101	faddp				// x+sqrt(x^2-1) : 2*x : log(2)
102	fdivrl	MO(one)			// 1/(x+sqrt(x^2-1)) : 2*x : log(2)
103	fsubrp				// 2*x+1/(x+sqrt(x^2)-1) : log(2)
104	fyl2x				// log(2*x+1/(x+sqrt(x^2-1)))
105	ret
106
107	// x < 1 => NaN
108	.align ALIGNARG(4)
1095:	fldz
110	fdiv	%st, %st(0)
111	ret
112END(__ieee754_acoshl)
113