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 <libm-alias-ldouble.h>
20#include <machine/asm.h>
21
22 .section .rodata
23
24 .align ALIGNARG(4)
25 .type huge,@object
26huge: .quad 0x89b634e7456ffa1d /* 1e+4930 */
27 .short 0x7ff8
28 ASM_SIZE_DIRECTIVE(huge)
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 .type 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 .type 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(__asinhl)
51 movl 12(%esp), %ecx
52 movl $0x7fff, %eax
53 andl %ecx, %eax
54 andl $0x8000, %ecx
55 movl %eax, %edx
56 orl $0xffff8000, %edx
57 incl %edx
58 jz 7f // x in ħInf or NaN
59 xorl %ecx, 12(%esp)
60 fldt 4(%esp) // |x|
61 cmpl $0x3fde, %eax
62 jb 2f // |x| < 2^-34
63 fldln2 // log(2) : |x|
64 cmpl $0x4020, %eax
65 fxch // |x| : log(2)
66 ja 3f // |x| > 2^34
67#ifdef PIC
68 LOAD_PIC_REG (dx)
69#endif
70 cmpl $0x4000, %eax
71 ja 5f // |x| > 2
72
73 // 2^-34 <= |x| <= 2 => y = sign(x)*log1p(|x|+|x|^2/(1+sqrt(1+|x|^2)))
74 fld %st // |x| : |x| : log(2)
75 fmul %st(1) // |x|^2 : |x| : log(2)
76 fld %st // |x|^2 : |x|^2 : |x| : log(2)
77 faddl MO(one) // 1+|x|^2 : |x|^2 : |x| : log(2)
78 fsqrt // sqrt(1+|x|^2) : |x|^2 : |x| : log(2)
79 faddl MO(one) // 1+sqrt(1+|x|^2) : |x|^2 : |x| : log(2)
80 fdivrp // |x|^2/(1+sqrt(1+|x|^2)) : |x| : log(2)
81 faddp // |x|+|x|^2/(1+sqrt(1+|x|^2)) : log(2)
82 fcoml MO(limit)
83 fnstsw
84 sahf
85 ja 6f
86 fyl2xp1
87 jecxz 4f
88 fchs
894: ret
90
917: fldt 4(%esp)
92 fadd %st
93 ret
94
956: faddl MO(one)
96 fyl2x
97 jecxz 4f
98 fchs
994: ret
100
101 // |x| < 2^-34 => y = x (inexact iff |x| != 0.0)
102 .align ALIGNARG(4)
1032:
104#ifdef PIC
105 LOAD_PIC_REG (dx)
106#endif
107 jecxz 4f
108 fchs // x
1094: fld %st // x : x
110 fldt MO(huge) // huge : x : x
111 faddp // huge+x : x
112 fstp %st(0) // x
113 cmpl $0x0001, %eax
114 jae 8f
115 fld %st(0)
116 fmul %st(0)
117 fstp %st(0)
1188: ret
119
120 // |x| > 2^34 => y = sign(x) * (log(|x|) + log(2))
121 .align ALIGNARG(4)
1223: fyl2x // log(|x|)
123 fldln2 // log(2) : log(|x|)
124 faddp // log(|x|)+log(2)
125 jecxz 4f
126 fchs
1274: ret
128
129 // |x| > 2 => y = sign(x) * log(2*|x| + 1/(|x|+sqrt(x*x+1)))
130 .align ALIGNARG(4)
1315: fld %st // |x| : |x| : log(2)
132 fadd %st, %st(1) // |x| : 2*|x| : log(2)
133 fld %st // |x| : |x| : 2*|x| : log(2)
134 fmul %st(1) // |x|^2 : |x| : 2*|x| : log(2)
135 faddl MO(one) // 1+|x|^2 : |x| : 2*|x| : log(2)
136 fsqrt // sqrt(1+|x|^2) : |x| : 2*|x| : log(2)
137 faddp // |x|+sqrt(1+|x|^2) : 2*|x| : log(2)
138 fdivrl MO(one) // 1/(|x|+sqrt(1+|x|^2)) : 2*|x| : log(2)
139 faddp // 2*|x|+1/(|x|+sqrt(1+|x|^2)) : log(2)
140 fyl2x // log(2*|x|+1/(|x|+sqrt(1+|x|^2)))
141 jecxz 4f
142 fchs
1434: ret
144END(__asinhl)
145libm_alias_ldouble (__asinh, asinh)
146

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