1 /* Function tanhf vectorized with SSE4.
2 Copyright (C) 2021-2024 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
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.
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.
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/. */
20 * ALGORITHM DESCRIPTION:
22 * NOTE: Since the hyperbolic tangent function is odd
23 * (tanh(x) = -tanh(-x)), below algorithm deals with the absolute
24 * value of the argument |x|: tanh(x) = sign(x) * tanh(|x|)
26 * We use a table lookup method to compute tanh(|x|).
27 * The basic idea is to split the input range into a number of subintervals
28 * and to approximate tanh(.) with a polynomial on each of them.
30 * IEEE SPECIAL CONDITIONS:
31 * x = [+, -]0, r = [+, -]0
39 * We handle special values in a callout function, aside from main path
40 * computations. "Special" for this algorithm are:
41 * INF, NAN, |x| > HUGE_THRESHOLD
44 * Main path computations are organized as follows:
45 * Actually we split the interval [0, SATURATION_THRESHOLD)
46 * into a number of subintervals. On each subinterval we approximate tanh(.)
47 * with a minimax polynomial of pre-defined degree. Polynomial coefficients
48 * are computed beforehand and stored in table. We also use
52 * here B depends on subinterval and is used to make argument
54 * We also add large fake interval [SATURATION_THRESHOLD, HUGE_THRESHOLD],
55 * where 1.0 + 0.0*y + 0.0*y^2 ... coefficients are stored - just to
56 * preserve main path computation logic but return 1.0 for all arguments.
58 * Hence reconstruction looks as follows:
59 * we extract proper polynomial and range reduction coefficients
60 * (Pj and B), corresponding to subinterval, to which |x| belongs,
63 * r := sign(x) * (P0 + P1 * y + ... + Pn * y^n)
65 * NOTE: we use multiprecision technique to multiply and sum the first
66 * K terms of the polynomial. So Pj, j = 0..K are stored in
67 * table each as a pair of target precision numbers (Pj and PLj) to
68 * achieve wider than target precision.
76 /* tanhf data tables for avx2 and sse4 implementations defined here.
78 #define ONLY_DECL_OFFSET
79 #include "svml_s_tanhf_rodata.S"
81 .section .text.sse4, "ax", @progbits
82 ENTRY(_ZGVbN4v_tanhf_sse4)
83 /* Save copy of input in xmm12. */
86 /* Here huge arguments, INF and NaNs are filtered out to callout. */
87 movdqu TANHF_DATA(_iExpMantMask)(%rip), %xmm3
91 /* Selection of arguments between [0, 0x04280000] into xmm3. */
93 /* Save xmm3 for special values check at end. */
95 psubd TANHF_DATA(_iMinIdxOfsMask)(%rip), %xmm3
97 pminsd TANHF_DATA(_iMaxIdxMask)(%rip), %xmm3
104 pshufd $0x0e, %xmm3, %xmm3
109 movaps TANHF_DATA(_sAbsMask)(%rip), %xmm1
112 leaq TANHF_DATA(_lookupTable)(%rip), %rax
113 movups (%rdx, %rax), %xmm2
114 movups (%rcx, %rax), %xmm6
117 * small table specific variables *
122 unpckhpd %xmm6, %xmm2
124 cvtps2pd %xmm0, %xmm6
126 cvtps2pd %xmm0, %xmm0
128 movups 16(%rdx, %rax), %xmm5
129 movups 16(%rsi, %rax), %xmm13
132 movaps %xmm13, %xmm11
134 movups 16(%rcx, %rax), %xmm7
135 movups 16(%rdi, %rax), %xmm3
137 unpckhpd %xmm7, %xmm5
138 unpckhpd %xmm3, %xmm13
143 movlhps %xmm7, %xmm10
144 movlhps %xmm3, %xmm11
154 movups (%rsi, %rax), %xmm2
155 movups (%rdi, %rax), %xmm7
159 unpckhpd %xmm7, %xmm2
170 cvtpd2ps %xmm0, %xmm2
171 cvtpd2ps %xmm6, %xmm0
177 /* xmm8 contains mask of special values. */
178 pcmpgtd TANHF_DATA(_iExpMask)(%rip), %xmm8
183 /* Go to special inputs processing branch */
184 jne L(SPECIAL_VALUES_BRANCH)
185 # LOE rbx rbp r12 r13 r14 r15 xmm0
186 /* No stack restoration on the fastpath. */
189 /* Cold case. edx has 1s where there was a special value that
190 needs to be handled by a tanhf call. Optimize for code size
191 more so than speed here. */
192 L(SPECIAL_VALUES_BRANCH):
193 # LOE rbx rdx rbp r12 r13 r14 r15 xmm0 xmm12
194 /* Stack coming in 16-byte aligned. Set 8-byte misaligned so on
195 call entry will be 16-byte aligned. */
197 cfi_def_cfa_offset(64)
198 movups %xmm0, 24(%rsp)
199 movups %xmm12, 40(%rsp)
201 /* Use rbx/rbp for callee save registers as they get short
202 encoding for many instructions (as compared with r12/r13). */
207 /* edx has 1s where there was a special value that needs to be handled
210 L(SPECIAL_VALUES_LOOP):
211 # LOE rbx rbp r12 r13 r14 r15
212 /* use rbp as index for special value that is saved across calls to
213 tanhf. We technically don't need a callee save register here as offset
214 to rsp is always [0, 12] so we can restore rsp by realigning to 64.
215 Essentially the tradeoff is 1 extra save/restore vs 2 extra instructions
220 /* Scalar math function call to process special input. */
221 movss 40(%rsp, %rbp, 4), %xmm0
223 /* No good way to avoid the store-forwarding fault this will cause on
224 return. `lfence` avoids the SF fault but at greater cost as it
225 serialized stack/callee save restoration. */
226 movss %xmm0, 24(%rsp, %rbp, 4)
230 jnz L(SPECIAL_VALUES_LOOP)
231 # LOE r12 r13 r14 r15
232 /* All results have been written to 24(%rsp). */
233 movups 24(%rsp), %xmm0
239 cfi_def_cfa_offset(8)
241 END(_ZGVbN4v_tanhf_sse4)