169 lines · c
1//===-- Single-precision atanf float 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#ifndef LIBC_SRC___SUPPORT_MATH_ATANF_FLOAT_H10#define LIBC_SRC___SUPPORT_MATH_ATANF_FLOAT_H11 12#include "src/__support/FPUtil/FEnvImpl.h"13#include "src/__support/FPUtil/FPBits.h"14#include "src/__support/FPUtil/double_double.h"15#include "src/__support/FPUtil/multiply_add.h"16#include "src/__support/FPUtil/nearest_integer.h"17#include "src/__support/macros/config.h"18#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY19 20namespace LIBC_NAMESPACE_DECL {21 22namespace math {23 24namespace atanf_internal {25 26using fputil::FloatFloat;27// atan(i/64) with i = 0..16, generated by Sollya with:28// > for i from 0 to 16 do {29// a = round(atan(i/16), SG, RN);30// b = round(atan(i/16) - a, SG, RN);31// print("{", b, ",", a, "},");32// };33static constexpr FloatFloat ATAN_I[17] = {34 {0.0f, 0.0f},35 {-0x1.1a6042p-30f, 0x1.ff55bcp-5f},36 {-0x1.54f424p-30f, 0x1.fd5baap-4f},37 {0x1.79cb6p-28f, 0x1.7b97b4p-3f},38 {-0x1.b4dfc8p-29f, 0x1.f5b76p-3f},39 {-0x1.1f0286p-27f, 0x1.362774p-2f},40 {0x1.e4defp-30f, 0x1.6f6194p-2f},41 {0x1.e611fep-29f, 0x1.a64eecp-2f},42 {0x1.586ed4p-28f, 0x1.dac67p-2f},43 {-0x1.6499e6p-26f, 0x1.0657eap-1f},44 {0x1.7bdfd6p-26f, 0x1.1e00bap-1f},45 {-0x1.98e422p-28f, 0x1.345f02p-1f},46 {0x1.934f7p-28f, 0x1.4978fap-1f},47 {0x1.c5a6c6p-27f, 0x1.5d5898p-1f},48 {0x1.5e118cp-27f, 0x1.700a7cp-1f},49 {-0x1.1d4eb6p-26f, 0x1.819d0cp-1f},50 {-0x1.777a5cp-26f, 0x1.921fb6p-1f},51};52 53// 1 / (1 + (i/16)^2) with i = 0..16, generated by Sollya with:54// > for i from 0 to 16 do {55// a = round(1 / (1 + (i/16)^2), SG, RN);56// print(a, ",");57// };58static constexpr float ATANF_REDUCED_ARG[17] = {59 0x1.0p0f, 0x1.fe01fep-1f, 0x1.f81f82p-1f, 0x1.ee9c8p-1f,60 0x1.e1e1e2p-1f, 0x1.d272cap-1f, 0x1.c0e07p-1f, 0x1.adbe88p-1f,61 0x1.99999ap-1f, 0x1.84f00cp-1f, 0x1.702e06p-1f, 0x1.5babccp-1f,62 0x1.47ae14p-1f, 0x1.34679ap-1f, 0x1.21fb78p-1f, 0x1.107fbcp-1f,63 0x1p-1f,64};65 66// Approximating atan( u / (1 + u * k/16) )67// atan( u / (1 + u * k/16) ) / u ~ 1 - k/16 * u + (k^2/256 - 1/3) * u^2 +68// + (k/16 - (k/16)^3) * u^3 + O(u^4)69LIBC_INLINE static float atanf_eval(float u, float k_over_16) {70 // (k/16)^271 float c2 = k_over_16 * k_over_16;72 // -(k/16)^373 float c3 = fputil::multiply_add(-k_over_16, c2, k_over_16);74 float u2 = u * u;75 // (k^2/256 - 1/3) + u * (k/16 - (k/16)^3)76 float a0 = fputil::multiply_add(c3, u, c2 - 0x1.555556p-2f);77 // -k/16 + u*(k^2/256 - 1/3) + u^2 * (k/16 - (k/16)^3)78 float a1 = fputil::multiply_add(u, a0, -k_over_16);79 // u - u^2 * k/16 + u^3 * ((k^2/256 - 1/3) + u^4 * (k/16 - (k/16)^3))80 return fputil::multiply_add(u2, a1, u);81}82 83} // namespace atanf_internal84 85// There are several range reduction steps we can take for atan2(y, x) as86// follow:87 88LIBC_INLINE static float atanf(float x) {89 using namespace atanf_internal;90 using FPBits = typename fputil::FPBits<float>;91 using FPBits = typename fputil::FPBits<float>;92 93 constexpr float SIGN[2] = {1.0f, -1.0f};94 constexpr FloatFloat PI_OVER_2 = {-0x1.777a5cp-25f, 0x1.921fb6p0f};95 96 FPBits x_bits(x);97 Sign s = x_bits.sign();98 float sign = SIGN[s.is_neg()];99 uint32_t x_abs = x_bits.uintval() & 0x7fff'ffffU;100 101 // x is inf or nan, |x| <= 2^-11 or |x|= > 2^11.102 if (LIBC_UNLIKELY(x_abs <= 0x3a00'0000U || x_abs >= 0x4500'0000U)) {103 if (LIBC_UNLIKELY(x_bits.is_inf()))104 return sign * PI_OVER_2.hi;105 // atan(NaN) = NaN106 if (LIBC_UNLIKELY(x_bits.is_nan())) {107 if (x_bits.is_signaling_nan()) {108 fputil::raise_except_if_required(FE_INVALID);109 return FPBits::quiet_nan().get_val();110 }111 112 return x;113 }114 // |x| >= 2^11:115 // atan(x) = sign(x) * pi/2 - atan(1/x)116 // ~ sign(x) * pi/2 - 1/x117 if (LIBC_UNLIKELY(x_abs >= 0x4200'0000))118 return fputil::multiply_add(sign, PI_OVER_2.hi, -1.0f / x);119 // x <= 2^-11:120 // atan(x) ~ x121 return x;122 }123 124 // Range reduction steps:125 // 1) atan(x) = sign(x) * atan(|x|)126 // 2) If |x| > 1, atan(|x|) = pi/2 - atan(1/|x|)127 // 3) For 1/16 < |x| < 1 + 1/32, we find k such that: | |x| - k/16 | <= 1/32.128 // Let y = |x| - k/16, then using the angle summation formula, we have:129 // atan(|x|) = atan(k/16) + atan( (|x| - k/16) / (1 + |x| * k/16) )130 // = atan(k/16) + atan( y / (1 + (y + k/16) * k/16 )131 // = atan(k/16) + atan( y / ((1 + k^2/256) + y * k/16) )132 // 4) Let u = y / (1 + k^2/256), then we can rewritten the above as:133 // atan(|x|) = atan(k/16) + atan( u / (1 + u * k/16) )134 // ~ atan(k/16) + (u - k/16 * u^2 + (k^2/256 - 1/3) * u^3 +135 // + (k/16 - (k/16)^3) * u^4) + O(u^5)136 float x_a = cpp::bit_cast<float>(x_abs);137 // |x| > 1 + 1/32, we need to invert x, so we will perform the division in138 // float-float.139 if (x_abs > 0x3f84'0000U)140 x_a = 1.0f / x_a;141 // Perform range reduction.142 // k = nearestint(x * 16)143 float k_f = fputil::nearest_integer(x_a * 0x1.0p4f);144 unsigned idx = static_cast<unsigned>(k_f);145 float k_over_16 = k_f * 0x1.0p-4f;146 float y = x_a - k_over_16;147 // u = (x - k/16) / (1 + (k/16)^2)148 float u = y * ATANF_REDUCED_ARG[idx];149 150 // atan(x) = sign(x) * atan(|x|)151 // = sign(x) * (atan(k/16) + atan(|))152 // p ~ atan(u)153 float p = atanf_eval(u, k_over_16);154 // |x| > 1 + 1/32: q ~ (pi/2 - atan(1/|x|))155 // |x| <= 1 + 1/32: q ~ atan(|x|)156 float q = (p + ATAN_I[idx].lo) + ATAN_I[idx].hi;157 if (x_abs > 0x3f84'0000U)158 q = PI_OVER_2.hi + (PI_OVER_2.lo - q);159 // |x| > 1 + 1/32: sign(x) * (pi/2 - atan(1/|x|))160 // |x| <= 1 + 1/32: sign(x) * atan(|x|)161 return sign * q;162}163 164} // namespace math165 166} // namespace LIBC_NAMESPACE_DECL167 168#endif // LIBC_SRC___SUPPORT_MATH_ATANF_FLOAT_H169