115 lines · cpp
1//===-- Half-precision tanpif 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-exception6//7//===----------------------------------------------------------------------===//8 9#include "src/math/tanpif16.h"10#include "hdr/errno_macros.h"11#include "hdr/fenv_macros.h"12#include "src/__support/FPUtil/FEnvImpl.h"13#include "src/__support/FPUtil/FPBits.h"14#include "src/__support/FPUtil/cast.h"15#include "src/__support/FPUtil/except_value_utils.h"16#include "src/__support/FPUtil/multiply_add.h"17#include "src/__support/macros/optimization.h"18#include "src/__support/math/sincosf16_utils.h"19 20namespace LIBC_NAMESPACE_DECL {21 22#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS23constexpr size_t N_EXCEPTS = 21;24 25constexpr fputil::ExceptValues<float16, N_EXCEPTS> TANPIF16_EXCEPTS{{26 // (input, RZ output, RU offset, RD offset, RN offset)27 {0x07f2, 0x0e3d, 1, 0, 0}, {0x086a, 0x0eee, 1, 0, 1},28 {0x08db, 0x0fa0, 1, 0, 0}, {0x094c, 0x1029, 1, 0, 0},29 {0x0b10, 0x118c, 1, 0, 0}, {0x1ce0, 0x23a8, 1, 0, 1},30 {0x1235, 0x18e0, 1, 0, 0}, {0x2579, 0x2c4e, 1, 0, 0},31 {0x28b2, 0x2f68, 1, 0, 1}, {0x2a43, 0x30f4, 1, 0, 1},32 {0x31b7, 0x3907, 1, 0, 0}, {0x329d, 0x3a12, 1, 0, 1},33 {0x34f1, 0x3dd7, 1, 0, 0}, {0x3658, 0x41ee, 1, 0, 0},34 {0x38d4, 0xc1ee, 0, 1, 0}, {0x3d96, 0x41ee, 1, 0, 0},35 {0x3e6a, 0xc1ee, 0, 1, 0}, {0x40cb, 0x41ee, 1, 0, 0},36 {0x4135, 0xc1ee, 0, 1, 0}, {0x42cb, 0x41ee, 1, 0, 0},37 {0x4335, 0xc1ee, 0, 1, 0},38}};39#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS40 41LLVM_LIBC_FUNCTION(float16, tanpif16, (float16 x)) {42 using namespace sincosf16_internal;43 using FPBits = typename fputil::FPBits<float16>;44 FPBits xbits(x);45 46 uint16_t x_u = xbits.uintval();47 uint16_t x_abs = x_u & 0x7fff;48 49 // Handle exceptional values50 if (LIBC_UNLIKELY(x_abs <= 0x4335)) {51 if (LIBC_UNLIKELY(x_abs == 0U))52 return x;53 54#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS55 bool x_sign = x_u >> 15;56 57 if (auto r = TANPIF16_EXCEPTS.lookup_odd(x_abs, x_sign);58 LIBC_UNLIKELY(r.has_value()))59 return r.value();60#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS61 }62 63 // Numbers greater or equal to 2^10 are integers, or infinity, or NaN64 if (LIBC_UNLIKELY(x_abs >= 0x6400)) {65 // Check for NaN or infinity values66 if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {67 if (xbits.is_signaling_nan()) {68 fputil::raise_except_if_required(FE_INVALID);69 return FPBits::quiet_nan().get_val();70 }71 // is inf72 if (x_abs == 0x7c00) {73 fputil::set_errno_if_required(EDOM);74 fputil::raise_except_if_required(FE_INVALID);75 }76 77 return x + FPBits::quiet_nan().get_val();78 }79 80 return FPBits::zero(xbits.sign()).get_val();81 }82 // Range reduction:83 // For |x| > 1/32, we perform range reduction as follows:84 // Find k and y such that:85 // x = (k + y) * 1/3286 // k is an integer87 // |y| < 0.588 //89 // This is done by performing:90 // k = round(x * 32)91 // y = x * 32 - k92 //93 // Once k and y are computed, we then deduce the answer by the formula:94 // tan(x) = sin(x) / cos(x)95 // = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k)96 float xf = x;97 float sin_k, cos_k, sin_y, cosm1_y;98 sincospif16_eval(xf, sin_k, cos_k, sin_y, cosm1_y);99 100 if (LIBC_UNLIKELY(sin_y == 0 && cos_k == 0)) {101 fputil::set_errno_if_required(EDOM);102 fputil::raise_except_if_required(FE_DIVBYZERO);103 104 int16_t x_mp5_u = static_cast<int16_t>(x - 0.5);105 return ((x_mp5_u & 0x1) ? -1 : 1) * FPBits::inf().get_val();106 }107 108 using fputil::multiply_add;109 return fputil::cast<float16>(110 multiply_add(sin_y, cos_k, multiply_add(cosm1_y, sin_k, sin_k)) /111 multiply_add(sin_y, -sin_k, multiply_add(cosm1_y, cos_k, cos_k)));112}113 114} // namespace LIBC_NAMESPACE_DECL115