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	/* Please note that we use double value for 1.0.  This number
26	   has an exact representation and so we don't get accuracy
27	   problems.  The advantage is that the code is simpler.  */
28	.type one,@object
29one:	.double 1.0
30	ASM_SIZE_DIRECTIVE(one)
31	/* It is not important that this constant is precise.  It is only
32	   a value which is known to be on the safe side for using the
33	   fyl2xp1 instruction.  */
34	.type limit,@object
35limit:	.double 0.29
36	ASM_SIZE_DIRECTIVE(limit)
37
38#ifdef PIC
39#define MO(op) op##@GOTOFF(%edx)
40#else
41#define MO(op) op
42#endif
43
44	.text
45ENTRY(__ieee754_acoshl)
46	movl	12(%esp), %ecx
47	andl	$0xffff, %ecx
48	cmpl	$0x3fff, %ecx
49	jl	5f			// < 1 => invalid
50	fldln2				// log(2)
51	fldt	4(%esp)			// x : log(2)
52	cmpl	$0x4020, %ecx
53	ja	3f			// x > 2^34
54#ifdef	PIC
55	LOAD_PIC_REG (dx)
56#endif
57	cmpl	$0x4000, %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	fabs				// acosh(1) is +0 in all rounding modes
63	fld	%st			// x-1 : x-1 : log(2)
64	fmul	%st(1)			// (x-1)^2 : x-1 : log(2)
65	fadd	%st(1)			// x-1+(x-1)^2 : x-1 : log(2)
66	fadd	%st(1)			// 2*(x-1)+(x-1)^2 : x-1 : log(2)
67	fsqrt				// sqrt(2*(x-1)+(x-1)^2) : x-1 : log(2)
68	faddp				// x-1+sqrt(2*(x-1)+(x-1)^2) : log(2)
69	fcoml	MO(limit)
70	fnstsw
71	sahf
72	ja	2f
73	fyl2xp1				// log1p(x-1+sqrt(2*(x-1)+(x-1)^2))
74	ret
75
762:	faddl	MO(one)			// x+sqrt(2*(x-1)+(x-1)^2) : log(2)
77	fyl2x				// log(x+sqrt(2*(x-1)+(x-1)^2))
78	ret
79
80	// x > 2^34 => y = log(x) + log(2)
81	.align ALIGNARG(4)
823:	fyl2x				// log(x)
83	fldln2				// log(2) : log(x)
84	faddp				// log(x)+log(2)
85	ret
86
87	// 2^34 > x > 2 => y = log(2*x - 1/(x+sqrt(x*x-1)))
88	.align ALIGNARG(4)
894:	fld	%st			// x : x : log(2)
90	fadd	%st, %st(1)		// x : 2*x : log(2)
91	fld	%st			// x : x : 2*x : log(2)
92	fmul	%st(1)			// x^2 : x : 2*x : log(2)
93	fsubl	MO(one)			// x^2-1 : x : 2*x : log(2)
94	fsqrt				// sqrt(x^2-1) : x : 2*x : log(2)
95	faddp				// x+sqrt(x^2-1) : 2*x : log(2)
96	fdivrl	MO(one)			// 1/(x+sqrt(x^2-1)) : 2*x : log(2)
97	fsubrp				// 2*x+1/(x+sqrt(x^2)-1) : log(2)
98	fyl2x				// log(2*x+1/(x+sqrt(x^2-1)))
99	ret
100
101	// x < 1 => NaN
102	.align ALIGNARG(4)
1035:	fldz
104	fdiv	%st, %st(0)
105	ret
106END(__ieee754_acoshl)
107libm_alias_finite (__ieee754_acoshl, __acoshl)
108