1
- // ===-- Half-precision asinhf16 (x) function--------------------------------===//
1
+ // ===-- Half-precision asinh (x) function -- --------------------------------===//
2
2
//
3
3
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4
4
// See https://llvm.org/LICENSE.txt for license information.
5
5
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception.
6
6
//
7
- //
8
7
// ===----------------------------------------------------------------------===//
9
8
10
9
#include " src/math/asinhf16.h"
10
+ #include " explogxf.h"
11
+ #include " hdr/fenv_macros.h"
12
+ #include " src/__support/FPUtil/cast.h"
11
13
#include " src/__support/FPUtil/except_value_utils.h"
12
- #include " src/__support/FPUtil/generic/sqrt.h"
13
14
#include " src/__support/FPUtil/multiply_add.h"
15
+ #include " src/__support/FPUtil/sqrt.h"
14
16
#include " src/__support/common.h"
15
17
#include " src/__support/macros/config.h"
16
- #include " src/__support/macros/properties/types.h"
17
- #include " src/math/generic/explogxf.h"
18
+ #include " src/__support/macros/optimization.h"
18
19
19
20
namespace LIBC_NAMESPACE_DECL {
20
21
21
- #ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
22
22
static constexpr size_t N_EXCEPTS = 8 ;
23
23
24
24
static constexpr fputil::ExceptValues<float16, N_EXCEPTS> ASINHF16_EXCEPTS{
25
25
{// (input, RZ output, RU offset, RD offset, RN offset)
26
- {0x3769 , 0x372A , 1 , 0 , 1 },
27
- {0x3B5B , 0x3A96 , 1 , 0 , 0 },
28
- {0x4B1F , 0x42B3 , 1 , 0 , 0 },
29
- {0x4C9B , 0x4336 , 1 , 0 , 1 },
30
- {0xB769 , 0xB72A , 0 , 1 , 1 },
31
- {0xBB5B , 0xBA96 , 0 , 1 , 0 },
32
- {0xCB1F , 0xC2B3 , 0 , 1 , 0 },
33
- {0xCC9B , 0xC336 , 0 , 1 , 1 }}};
34
- #endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
26
+
27
+ // x = 0x1.da4p-2, asinhf16(x) = 0x1.ca8p-2 (RZ)
28
+ {0x3769 , 0x372a , 1 , 0 , 1 },
29
+ // x = 0x1.d6cp-1, asinhf16(x) = 0x1.a58p-1 (RZ)
30
+ {0x3b5b , 0x3a96 , 1 , 0 , 0 },
31
+ // x = 0x1.c7cp+3, asinhf16(x) = 0x1.accp+1 (RZ)
32
+ {0x4b1f , 0x42b3 , 1 , 0 , 0 },
33
+ // x = 0x1.26cp+4, asinhf16(x) = 0x1.cd8p+1 (RZ)
34
+ {0x4c9b , 0x4336 , 1 , 0 , 1 },
35
+ // x = -0x1.da4p-2, asinhf16(x) = -0x1.ca8p-2 (RZ)
36
+ {0xb769 , 0xb72a , 0 , 1 , 1 },
37
+ // x = -0x1.d6cp-1, asinhf16(x) = -0x1.a58p-1 (RZ)
38
+ {0xbb5b , 0xba96 , 0 , 1 , 0 },
39
+ // x = -0x1.c7cp+3, asinhf16(x) = -0x1.accp+1 (RZ)
40
+ {0xcb1f , 0xc2b3 , 0 , 1 , 0 },
41
+ // x = -0x1.26cp+4, asinhf16(x) = -0x1.cd8p+1 (RZ)
42
+ {0xcc9b , 0xc336 , 0 , 1 , 1 }}};
35
43
36
44
LLVM_LIBC_FUNCTION (float16, asinhf16, (float16 x)) {
37
45
using FPBits = fputil::FPBits<float16>;
38
46
FPBits xbits (x);
39
47
40
- float x_d = x;
41
48
uint16_t x_u = xbits.uintval ();
42
49
uint16_t x_abs = x_u & 0x7fff ;
43
50
@@ -50,35 +57,45 @@ LLVM_LIBC_FUNCTION(float16, asinhf16, (float16 x)) {
50
57
return x;
51
58
}
52
59
53
- #ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
54
60
// Handle exceptional values
55
61
if (auto r = ASINHF16_EXCEPTS.lookup (x_u); LIBC_UNLIKELY (r.has_value ()))
56
62
return r.value ();
57
- #endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
58
63
64
+ float xf = x;
59
65
const float SIGN[2 ] = {1 .0f , -1 .0f };
60
66
float x_sign = SIGN[x_u >> 15 ];
61
67
62
68
// |x| <= 0.25
63
69
if (LIBC_UNLIKELY (x_abs <= 0x3400 )) {
64
- if (LIBC_UNLIKELY (x_abs == 0 ))
65
- return x;
66
- if (LIBC_UNLIKELY ((fputil::get_round () == FE_UPWARD) && (x_u >= 0x8401 ) &&
67
- (x_u <= 0x90E6 )))
68
- return static_cast <float16>(x_d + 0x1p-24f );
69
70
70
- float x_sq = x_d * x_d;
71
+ // when |x| < 0x1.718p-5, asinhf16(x) = x. Adjust by 1 ULP for certain
72
+ // rounding types.
73
+ if (LIBC_UNLIKELY (x_abs < 0x29c6 )) {
74
+ if (((fputil::get_round () == FE_UPWARD) ||
75
+ (fputil::get_round () == FE_TOWARDZERO)) &&
76
+ xf < 0 )
77
+ return fputil::cast<float16>(xf + 0x1p-24f );
78
+ if (((fputil::get_round () == FE_DOWNWARD) ||
79
+ (fputil::get_round () == FE_TOWARDZERO)) &&
80
+ xf > 0 )
81
+ return fputil::cast<float16>(xf - 0x1p-24f );
82
+ return fputil::cast<float16>(xf);
83
+ }
84
+
85
+ float x_sq = xf * xf;
71
86
// Generated by Sollya with:
72
- // > P = fpminimax(asinh(x)/x, [|0, 2, 4, 6, 8|], [|SG...|],[0, 2^-2]);
73
- float p = fputil::polyeval (x_sq, 1 .0f , -0x1 .555556p-3f , 0x1 .3334dep-4f ,
74
- -0x1 .6f3e2p-5f , 0x1 .51d012p-5f );
87
+ // > P = fpminimax(asinh(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 2^-2]);
88
+ // The last coefficient 0x1.bd114ep-6f has been changed to 0x1.bd114ep-5f
89
+ // for better accuracy.
90
+ float p = fputil::polyeval (x_sq, 1 .0f , -0x1 .555552p-3f , 0x1 .332f6ap-4f ,
91
+ -0x1 .6c53dep-5f , 0x1 .bd114ep -5f );
75
92
76
- return static_cast <float16>(fputil::multiply_add (x_d, p, 0 . 0f ) );
93
+ return fputil::cast <float16>(xf * p );
77
94
}
78
95
79
96
// General case: asinh(x) = ln(x + sqrt(x^2 + 1))
80
- float sqrt_term = fputil::sqrt<float >(fputil::multiply_add (x_d, x_d , 1 .0f ));
97
+ float sqrt_term = fputil::sqrt<float >(fputil::multiply_add (xf, xf , 1 .0f ));
81
98
return fputil::cast<float16>(
82
- x_sign * log_eval (fputil::multiply_add (x_d , x_sign, sqrt_term)));
99
+ x_sign * log_eval (fputil::multiply_add (xf , x_sign, sqrt_term)));
83
100
}
84
101
} // namespace LIBC_NAMESPACE_DECL
0 commit comments