1 //===-- Half-precision tan(x) function ------------------------------------===// 2 // 3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. 4 // See https://llvm.org/LICENSE.txt for license information. 5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception. 6 // 7 //===----------------------------------------------------------------------===// 8 9 #include "src/math/tanf16.h" 10 #include "hdr/errno_macros.h" 11 #include "hdr/fenv_macros.h" 12 #include "sincosf16_utils.h" 13 #include "src/__support/FPUtil/FEnvImpl.h" 14 #include "src/__support/FPUtil/FPBits.h" 15 #include "src/__support/FPUtil/cast.h" 16 #include "src/__support/FPUtil/except_value_utils.h" 17 #include "src/__support/FPUtil/multiply_add.h" 18 #include "src/__support/macros/optimization.h" 19 20 namespace LIBC_NAMESPACE_DECL { 21 22 constexpr size_t N_EXCEPTS = 9; 23 24 constexpr fputil::ExceptValues<float16, N_EXCEPTS> TANF16_EXCEPTS{{ 25 // (input, RZ output, RU offset, RD offset, RN offset) 26 {0x2894, 0x2894, 1, 0, 1}, 27 {0x3091, 0x3099, 1, 0, 0}, 28 {0x3098, 0x30a0, 1, 0, 0}, 29 {0x55ed, 0x3911, 1, 0, 0}, 30 {0x607b, 0xc638, 0, 1, 1}, 31 {0x674e, 0x3b7d, 1, 0, 0}, 32 {0x6807, 0x4014, 1, 0, 1}, 33 {0x6f4d, 0xbe19, 0, 1, 1}, 34 {0x7330, 0xcb62, 0, 1, 0}, 35 }}; 36 37 LLVM_LIBC_FUNCTION(float16, tanf16, (float16 x)) { 38 using FPBits = fputil::FPBits<float16>; 39 FPBits xbits(x); 40 41 uint16_t x_u = xbits.uintval(); 42 uint16_t x_abs = x_u & 0x7fff; 43 bool x_sign = x_u >> 15; 44 float xf = x; 45 46 // Handle exceptional values 47 if (auto r = TANF16_EXCEPTS.lookup_odd(x_abs, x_sign); 48 LIBC_UNLIKELY(r.has_value())) 49 return r.value(); 50 51 // |x| <= 0x1.d1p-5 52 if (LIBC_UNLIKELY(x_abs <= 0x2b44)) { 53 // |x| <= 0x1.398p-11 54 if (LIBC_UNLIKELY(x_abs <= 0x10e6)) { 55 // tan(+/-0) = +/-0 56 if (LIBC_UNLIKELY(x_abs == 0)) 57 return x; 58 59 int rounding = fputil::quick_get_round(); 60 61 // Exhaustive tests show that, when: 62 // x > 0, and rounding upward or 63 // x < 0, and rounding downward then, 64 // tan(x) = x * 2^-11 + x 65 if ((xbits.is_pos() && rounding == FE_UPWARD) || 66 (xbits.is_neg() && rounding == FE_DOWNWARD)) 67 return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf)); 68 return x; 69 } 70 71 float xsq = xf * xf; 72 73 // Degree-6 minimax odd polynomial of tan(x) generated by Sollya with: 74 // > P = fpminimax(tan(x)/x, [|0, 2, 4, 6|], [|1, SG...|], [0, pi/32]); 75 float result = fputil::polyeval(xsq, 0x1p0f, 0x1.555556p-2f, 0x1.110ee4p-3f, 76 0x1.be80f6p-5f); 77 78 return fputil::cast<float16>(xf * result); 79 } 80 81 // tan(+/-inf) = NaN, and tan(NaN) = NaN 82 if (LIBC_UNLIKELY(x_abs >= 0x7c00)) { 83 // x = +/-inf 84 if (x_abs == 0x7c00) { 85 fputil::set_errno_if_required(EDOM); 86 fputil::raise_except_if_required(FE_INVALID); 87 } 88 89 return x + FPBits::quiet_nan().get_val(); 90 } 91 92 // Range reduction: 93 // For |x| > pi/32, we perform range reduction as follows: 94 // Find k and y such that: 95 // x = (k + y) * pi/32; 96 // k is an integer, |y| < 0.5 97 // 98 // This is done by performing: 99 // k = round(x * 32/pi) 100 // y = x * 32/pi - k 101 // 102 // Once k and y are computed, we then deduce the answer by the formula: 103 // tan(x) = sin(x) / cos(x) 104 // = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k) 105 float sin_k, cos_k, sin_y, cosm1_y; 106 sincosf16_eval(xf, sin_k, cos_k, sin_y, cosm1_y); 107 108 // Note that, cosm1_y = cos_y - 1: 109 using fputil::multiply_add; 110 return fputil::cast<float16>( 111 multiply_add(sin_y, cos_k, multiply_add(cosm1_y, sin_k, sin_k)) / 112 multiply_add(sin_y, -sin_k, multiply_add(cosm1_y, cos_k, cos_k))); 113 } 114 115 } // namespace LIBC_NAMESPACE_DECL 116