1/* ix87 specific implementation of arcsinh.
2   Copyright (C) 1996-2022 Free Software Foundation, Inc.
3   This file is part of the GNU C Library.
4
5   The GNU C Library is free software; you can redistribute it and/or
6   modify it under the terms of the GNU Lesser General Public
7   License as published by the Free Software Foundation; either
8   version 2.1 of the License, or (at your option) any later version.
9
10   The GNU C Library is distributed in the hope that it will be useful,
11   but WITHOUT ANY WARRANTY; without even the implied warranty of
12   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13   Lesser General Public License for more details.
14
15   You should have received a copy of the GNU Lesser General Public
16   License along with the GNU C Library; if not, see
17   <https://www.gnu.org/licenses/>.  */
18
19#include <machine/asm.h>
20#include <libm-alias-finite.h>
21
22	.section .rodata.cst8,"aM",@progbits,8
23
24	.p2align 3
25	.type one,@object
26one:	.double 1.0
27	ASM_SIZE_DIRECTIVE(one)
28	.type limit,@object
29limit:	.double 0.29
30	ASM_SIZE_DIRECTIVE(limit)
31
32#ifdef PIC
33#define MO(op) op##@GOTOFF(%edx)
34#else
35#define MO(op) op
36#endif
37
38	.text
39ENTRY(__ieee754_acoshf)
40	movl	4(%esp), %ecx
41	cmpl	$0x3f800000, %ecx
42	jl	5f			// < 1 => invalid
43	fldln2				// log(2)
44	flds	4(%esp)			// x : log(2)
45	cmpl	$0x47000000, %ecx
46	ja	3f			// x > 2^14
47#ifdef	PIC
48	LOAD_PIC_REG (dx)
49#endif
50	cmpl	$0x40000000, %ecx
51	ja	4f			// x > 2
52
53	// 1 <= x <= 2 => y = log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
54	fsubl	MO(one)			// x-1 : log(2)
55	fabs				// acosh(1) is +0 in all rounding modes
56	fld	%st			// x-1 : x-1 : log(2)
57	fmul	%st(1)			// (x-1)^2 : x-1 : log(2)
58	fadd	%st(1)			// x-1+(x-1)^2 : x-1 : log(2)
59	fadd	%st(1)			// 2*(x-1)+(x-1)^2 : x-1 : log(2)
60	fsqrt				// sqrt(2*(x-1)+(x-1)^2) : x-1 : log(2)
61	faddp				// x-1+sqrt(2*(x-1)+(x-1)^2) : log(2)
62	fcoml	MO(limit)
63	fnstsw
64	sahf
65	ja	2f
66	fyl2xp1				// log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
67	ret
68
692:	faddl	MO(one)			// x+sqrt(2*(x-1)+(x-1)^2) : log(2)
70	fyl2x				// log(x+sqrt(2*(x-1)+(x-1)^2))
71	ret
72
73	// x > 2^14 => y = log(x) + log(2)
74	.align ALIGNARG(4)
753:	fyl2x				// log(x)
76	fldln2				// log(2) : log(x)
77	faddp				// log(x)+log(2)
78	ret
79
80	// 2^28 > x > 2 => y = log(2*x - 1/(x+sqrt(x*x-1)))
81	.align ALIGNARG(4)
824:	fld	%st			// x : x : log(2)
83	fadd	%st, %st(1)		// x : 2*x : log(2)
84	fld	%st			// x : x : 2*x : log(2)
85	fmul	%st(1)			// x^2 : x : 2*x : log(2)
86	fsubl	MO(one)			// x^2-1 : x : 2*x : log(2)
87	fsqrt				// sqrt(x^2-1) : x : 2*x : log(2)
88	faddp				// x+sqrt(x^2-1) : 2*x : log(2)
89	fdivrl	MO(one)			// 1/(x+sqrt(x^2-1)) : 2*x : log(2)
90	fsubrp				// 2*x+1/(x+sqrt(x^2)-1) : log(2)
91	fyl2x				// log(2*x+1/(x+sqrt(x^2-1)))
92	ret
93
94	// x < 1 (or -NaN) => NaN
95	.align ALIGNARG(4)
965:	flds	4(%esp)
97	fsub	%st
98	fdiv	%st, %st(0)
99	ret
100END(__ieee754_acoshf)
101libm_alias_finite (__ieee754_acoshf, __acoshf)
102