1/* ix87 specific implementation of arcsinh.
2 Copyright (C) 1996-2024 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

source code of glibc/sysdeps/i386/fpu/e_acoshl.S