140 lines · c
1//===-- Implementation header for atanf -------------------------*- C++ -*-===//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#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_ATANF_H10#define LLVM_LIBC_SRC___SUPPORT_MATH_ATANF_H11 12#include "inv_trigf_utils.h"13#include "src/__support/FPUtil/FPBits.h"14#include "src/__support/FPUtil/PolyEval.h"15#include "src/__support/FPUtil/except_value_utils.h"16#include "src/__support/FPUtil/multiply_add.h"17#include "src/__support/FPUtil/nearest_integer.h"18#include "src/__support/macros/config.h"19#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY20 21#if defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) && \22 defined(LIBC_MATH_HAS_INTERMEDIATE_COMP_IN_FLOAT)23 24// We use float-float implementation to reduce size.25#include "atanf_float.h"26 27#else28 29namespace LIBC_NAMESPACE_DECL {30 31namespace math {32 33LIBC_INLINE static constexpr float atanf(float x) {34 using namespace inv_trigf_utils_internal;35 using FPBits = typename fputil::FPBits<float>;36 37 constexpr double FINAL_SIGN[2] = {1.0, -1.0};38 constexpr double SIGNED_PI_OVER_2[2] = {0x1.921fb54442d18p0,39 -0x1.921fb54442d18p0};40 41 FPBits x_bits(x);42 Sign sign = x_bits.sign();43 x_bits.set_sign(Sign::POS);44 uint32_t x_abs = x_bits.uintval();45 46 // x is inf or nan, |x| < 2^-4 or |x|= > 16.47 if (LIBC_UNLIKELY(x_abs <= 0x3d80'0000U || x_abs >= 0x4180'0000U)) {48 double x_d = static_cast<double>(x);49 double const_term = 0.0;50 if (LIBC_UNLIKELY(x_abs >= 0x4180'0000)) {51 // atan(+-Inf) = +-pi/2.52 if (x_bits.is_inf()) {53 volatile double sign_pi_over_2 = SIGNED_PI_OVER_2[sign.is_neg()];54 return static_cast<float>(sign_pi_over_2);55 }56 if (x_bits.is_nan())57 return x;58 // x >= 1659 x_d = -1.0 / x_d;60 const_term = SIGNED_PI_OVER_2[sign.is_neg()];61 }62 // 0 <= x < 1/16;63 if (LIBC_UNLIKELY(x_bits.is_zero()))64 return x;65 // x <= 2^-12;66 if (LIBC_UNLIKELY(x_abs < 0x3980'0000)) {67#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)68 return fputil::multiply_add(x, -0x1.0p-25f, x);69#else70 double x_d = static_cast<double>(x);71 return static_cast<float>(fputil::multiply_add(x_d, -0x1.0p-25, x_d));72#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT73 }74 // Use Taylor polynomial:75 // atan(x) ~ x * (1 - x^2 / 3 + x^4 / 5 - x^6 / 7 + x^8 / 9 - x^10 / 11).76 constexpr double ATAN_TAYLOR[6] = {77 0x1.0000000000000p+0, -0x1.5555555555555p-2, 0x1.999999999999ap-3,78 -0x1.2492492492492p-3, 0x1.c71c71c71c71cp-4, -0x1.745d1745d1746p-4,79 };80 double x2 = x_d * x_d;81 double x4 = x2 * x2;82 double c0 = fputil::multiply_add(x2, ATAN_TAYLOR[1], ATAN_TAYLOR[0]);83 double c1 = fputil::multiply_add(x2, ATAN_TAYLOR[3], ATAN_TAYLOR[2]);84 double c2 = fputil::multiply_add(x2, ATAN_TAYLOR[5], ATAN_TAYLOR[4]);85 double p = fputil::polyeval(x4, c0, c1, c2);86 double r = fputil::multiply_add(x_d, p, const_term);87 return static_cast<float>(r);88 }89 90 // Range reduction steps:91 // 1) atan(x) = sign(x) * atan(|x|)92 // 2) If |x| > 1, atan(|x|) = pi/2 - atan(1/|x|)93 // 3) For 1/16 < x <= 1, we find k such that: |x - k/16| <= 1/32.94 // 4) Then we use polynomial approximation:95 // atan(x) ~ atan((k/16) + (x - (k/16)) * Q(x - k/16)96 // = P(x - k/16)97 double x_d = 0, const_term = 0, final_sign = 0;98 int idx = 0;99 100 if (x_abs > 0x3f80'0000U) {101 // |x| > 1, we need to invert x, so we will perform range reduction in102 // double precision.103 x_d = 1.0 / static_cast<double>(x_bits.get_val());104 double k_d = fputil::nearest_integer(x_d * 0x1.0p4);105 x_d = fputil::multiply_add(k_d, -0x1.0p-4, x_d);106 idx = static_cast<int>(k_d);107 final_sign = FINAL_SIGN[sign.is_pos()];108 // Adjust constant term of the polynomial by +- pi/2.109 const_term = fputil::multiply_add(final_sign, ATAN_COEFFS[idx][0],110 SIGNED_PI_OVER_2[sign.is_neg()]);111 } else {112 // Exceptional value:113 if (LIBC_UNLIKELY(x_abs == 0x3d8d'6b23U)) { // |x| = 0x1.1ad646p-4114 return sign.is_pos() ? fputil::round_result_slightly_down(0x1.1a6386p-4f)115 : fputil::round_result_slightly_up(-0x1.1a6386p-4f);116 }117 // Perform range reduction in single precision.118 float x_f = x_bits.get_val();119 float k_f = fputil::nearest_integer(x_f * 0x1.0p4f);120 x_f = fputil::multiply_add(k_f, -0x1.0p-4f, x_f);121 x_d = static_cast<double>(x_f);122 idx = static_cast<int>(k_f);123 final_sign = FINAL_SIGN[sign.is_neg()];124 const_term = final_sign * ATAN_COEFFS[idx][0];125 }126 127 double p = atan_eval(x_d, idx);128 double r = fputil::multiply_add(final_sign * x_d, p, const_term);129 130 return static_cast<float>(r);131}132 133} // namespace math134 135} // namespace LIBC_NAMESPACE_DECL136 137#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS138 139#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ATANF_H140