powf.cpp | Explore in Territory

//===-- Single-precision x^y function -------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//

#include "src/math/powf.h"
#include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
#include "src/__support/CPP/bit.h"
#include "src/__support/CPP/optional.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/double_double.h"
#include "src/__support/FPUtil/except_value_utils.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/FPUtil/rounding_mode.h"
#include "src/__support/FPUtil/sqrt.h" // Speedup for powf(x, 1/2) = sqrtf(x)
#include "src/__support/common.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY

#include "exp10f_impl.h" // Speedup for powf(10, y) = exp10f(y)
#include "exp2f_impl.h"  // Speedup for powf(2, y) = exp2f(y)

namespace LIBC_NAMESPACE_DECL {

DoubleDouble;
TripleDouble;

namespace {

#ifdef LIBC_TARGET_CPU_HAS_FMA
constexpr uint64_t ERR = 64;
#else
constexpr uint64_t ERR = …;
#endif // LIBC_TARGET_CPU_HAS_FMA

// We choose the precision of the high part to be 53 - 24 - 8, so that when
//   y * (e_x + LOG2_R_DD[i].hi) is exact.
// Generated by Sollya with:
// > for i from 0 to 127 do {
//     r = 2^-8 * ceil(2^8 * (1 - 2^-8) / (1 + i * 2^-7) );
//     a = -log2(r);
//     b = round(1 + a, 53 - 24 - 8, RN) - 1;
//     c = round(a - b, D, RN);
//     d = round(a - b - c, D, RN);
//     print("{", d, ",", c, ", ", b, "},");
//    };
static constexpr TripleDouble LOG2_R_TD[128] = …;

// Look up table for the second range reduction step:
// Generated by Sollya with:
// > for i from -64 to 128 do {
//     r = 2^-16 * nearestint(2^16 / (1 + i * 2^-14) );
//     a = -log2(r);
//     b = round(a, D, RN);
//     c = round(a - b, D, RN);
//     print("{", c, ", ", b, "},");
//    };
static constexpr DoubleDouble LOG2_R2_DD[] = …;

LIBC_INLINE bool is_odd_integer(float x) { … }

LIBC_INLINE bool is_integer(float x) { … }

LIBC_INLINE bool larger_exponent(double a, double b) { … }

// Calculate 2^(y * log2(x)) in double-double precision.
// At this point we can reuse the following values:
//   idx_x: index for extra precision of log2 for the middle part of log2(x).
//   dx: the reduced argument for log2(x)
//   y6: 2^6 * y.
//   lo6_hi: the high part of 2^6 * (y - (hi + mid))
//   exp2_hi_mid: high part of 2^(hi + mid)
double powf_double_double(int idx_x, double dx, double y6, double lo6_hi,
                          const DoubleDouble &exp2_hi_mid) { … }

} // namespace

LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) { … }

} // namespace LIBC_NAMESPACE_DECL
llvm/libc/src/math/generic/powf.cpp
