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 <libm-alias-ldouble.h>
20#include <machine/asm.h>
21
22 .section .rodata
23
24 .align ALIGNARG(4)
25 .type huge,@object
26huge: .tfloat 1e+4930
27 ASM_SIZE_DIRECTIVE(huge)
28 .align ALIGNARG(4)
29 /* Please note that we use double value for 1.0. This number
30 has an exact representation and so we don't get accuracy
31 problems. The advantage is that the code is simpler. */
32 .type one,@object
33one: .double 1.0
34 ASM_SIZE_DIRECTIVE(one)
35 /* It is not important that this constant is precise. It is only
36 a value which is known to be on the safe side for using the
37 fyl2xp1 instruction. */
38 .type limit,@object
39limit: .double 0.29
40 ASM_SIZE_DIRECTIVE(limit)
41
42#ifdef PIC
43#define MO(op) op##@GOTOFF(%edx)
44#else
45#define MO(op) op
46#endif
47
48 .text
49ENTRY(__asinhl)
50 movl 12(%esp), %ecx
51 movl $0x7fff, %eax
52 andl %ecx, %eax
53 andl $0x8000, %ecx
54 movl %eax, %edx
55 orl $0xffff8000, %edx
56 incl %edx
57 jz 7f // x in ħInf or NaN
58 xorl %ecx, 12(%esp)
59 fldt 4(%esp) // |x|
60 cmpl $0x3fde, %eax
61 jb 2f // |x| < 2^-34
62 fldln2 // log(2) : |x|
63 cmpl $0x4020, %eax
64 fxch // |x| : log(2)
65 ja 3f // |x| > 2^34
66#ifdef PIC
67 LOAD_PIC_REG (dx)
68#endif
69 cmpl $0x4000, %eax
70 ja 5f // |x| > 2
71
72 // 2^-34 <= |x| <= 2 => y = sign(x)*log1p(|x|+|x|^2/(1+sqrt(1+|x|^2)))
73 fld %st // |x| : |x| : log(2)
74 fmul %st(1) // |x|^2 : |x| : log(2)
75 fld %st // |x|^2 : |x|^2 : |x| : log(2)
76 faddl MO(one) // 1+|x|^2 : |x|^2 : |x| : log(2)
77 fsqrt // sqrt(1+|x|^2) : |x|^2 : |x| : log(2)
78 faddl MO(one) // 1+sqrt(1+|x|^2) : |x|^2 : |x| : log(2)
79 fdivrp // |x|^2/(1+sqrt(1+|x|^2)) : |x| : log(2)
80 faddp // |x|+|x|^2/(1+sqrt(1+|x|^2)) : log(2)
81 fcoml MO(limit)
82 fnstsw
83 sahf
84 ja 6f
85 fyl2xp1
86 jecxz 4f
87 fchs
884: ret
89
907: fldt 4(%esp)
91 fadd %st
92 ret
93
946: faddl MO(one)
95 fyl2x
96 jecxz 4f
97 fchs
984: ret
99
100 // |x| < 2^-34 => y = x (inexact iff |x| != 0.0)
101 .align ALIGNARG(4)
1022:
103#ifdef PIC
104 LOAD_PIC_REG (dx)
105#endif
106 jecxz 4f
107 fchs // x
1084: fld %st // x : x
109 fldt MO(huge) // huge : x : x
110 faddp // huge+x : x
111 fstp %st(0) // x
112 cmpl $0x0001, %eax
113 jae 8f
114 fld %st(0)
115 fmul %st(0)
116 fstp %st(0)
1178: ret
118
119 // |x| > 2^34 => y = sign(x) * (log(|x|) + log(2))
120 .align ALIGNARG(4)
1213: fyl2x // log(|x|)
122 fldln2 // log(2) : log(|x|)
123 faddp // log(|x|)+log(2)
124 jecxz 4f
125 fchs
1264: ret
127
128 // |x| > 2 => y = sign(x) * log(2*|x| + 1/(|x|+sqrt(x*x+1)))
129 .align ALIGNARG(4)
1305: fld %st // |x| : |x| : log(2)
131 fadd %st, %st(1) // |x| : 2*|x| : log(2)
132 fld %st // |x| : |x| : 2*|x| : log(2)
133 fmul %st(1) // |x|^2 : |x| : 2*|x| : log(2)
134 faddl MO(one) // 1+|x|^2 : |x| : 2*|x| : log(2)
135 fsqrt // sqrt(1+|x|^2) : |x| : 2*|x| : log(2)
136 faddp // |x|+sqrt(1+|x|^2) : 2*|x| : log(2)
137 fdivrl MO(one) // 1/(|x|+sqrt(1+|x|^2)) : 2*|x| : log(2)
138 faddp // 2*|x|+1/(|x|+sqrt(1+|x|^2)) : log(2)
139 fyl2x // log(2*|x|+1/(|x|+sqrt(1+|x|^2)))
140 jecxz 4f
141 fchs
1424: ret
143END(__asinhl)
144libm_alias_ldouble (__asinh, asinh)
145

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