Google Git Google Git
Sign in
llvm/llvm-project/libc/refs/heads/master/./src/__support/math/powf_utils.h
blob: dc36fbffedeaa5490740e5c8d73036419ebe4ef5 [file] [edit]
//===----------------------------------------------------------------------===//
//
// 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
//
//===----------------------------------------------------------------------===//
///
/// \file
/// Common utilities and tables for single-precision powf(x, y).
///
//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POWF_UTILS_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_POWF_UTILS_H
#include "hdr/errno_macros.h"
#include "hdr/fenv_macros.h"
#include "src/__support/CPP/bit.h"
#include "src/__support/CPP/optional.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/double_double.h"
#include "src/__support/FPUtil/is_integer.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/FPUtil/rounding_mode.h"
#include "src/__support/FPUtil/triple_double.h"
#include "src/__support/common.h"
#include "src/__support/macros/attributes.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
#include "src/__support/macros/properties/cpu_features.h"
#include "src/__support/math/common_constants.h"
#if defined(LIBC_TARGET_CPU_HAS_FPU_FLOAT) || \
!defined(LIBC_MATH_HAS_SMALL_TABLES)
#include "src/__support/FPUtil/sqrt.h"
#endif // LIBC_TARGET_CPU_HAS_FPU_FLOAT || !LIBC_MATH_HAS_SMALL_TABLES
#if !defined(LIBC_MATH_HAS_SMALL_TABLES) && \
!defined(LIBC_MATH_HAS_INTERMEDIATE_COMP_IN_FLOAT)
#include "src/__support/math/exp10f.h"
#include "src/__support/math/exp2f.h"
#endif // !LIBC_MATH_HAS_SMALL_TABLES &&
// !LIBC_MATH_HAS_INTERMEDIATE_COMP_IN_FLOAT
namespace LIBC_NAMESPACE_DECL {
namespace math {
namespace powf_internal {
using fputil::DoubleDouble;
using fputil::FloatFloat;
using fputil::TripleDouble;
// Table of -log2(R[i]) represented in triple-double precision {lo, mid, hi},
// where R[i] = 2^-8 * ceil(2^8 * (1 - 2^-8) / (1 + i * 2^-7)) for i = 0..127.
//
// The least significant bit of the high part `hi` is chosen to be >= 2^-20,
// so that the total precision of (e_x + hi) is at most 8 + 20 = 28 bits
// (since |e_x| <= 150 < 2^8 fits in 8 bits).
// Because the single-precision exponent y has 24 bits of mantissa, the product
// y * (e_x + hi) spans at most 24 + 28 = 52 bits <= 53 bits, making it exact
// in double precision.
// The mid part provides 53 bits for fast evaluation, and the lo part provides
// another 53 bits for the accurate fallback, yielding >= 127 bits of precision.
//
// Generated by Sollya with:
// > display = hexadecimal;
// > prec = 200;
// > for i from 0 to 127 do {
// if (i == 0) then {
// print("{0.0, 0.0, 0.0},");
// } else if (i == 127) then {
// print("{0.0, 0.0, 1.0},");
// } else {
// r = 2^(-8) * ceil( 2^8 * (1 - 2^(-8)) / (1 + i*2^(-7)) );
// l = -log2(r);
// h = nearestint(l * 2^20) * 2^(-20);
// m = round(l - h, D, RN);
// lo = round(l - h - m, D, RN);
// print("{", lo, ", ", m, ", ", h, "},");
// };
// };
LIBC_INLINE_VAR constexpr TripleDouble LOG2_R_TD[128] = {
{0.0, 0.0, 0.0},
{0x1.84a2c615b70adp-79, -0x1.177c23362928cp-25, 0x1.72c8p-7},
{-0x1.f27b820fd03eap-76, -0x1.179e0caa9c9abp-22, 0x1.744p-6},
{-0x1.f27ef487c8f34p-77, -0x1.c6cea541f5b7p-23, 0x1.184cp-5},
{-0x1.e3f80fbc71454p-76, -0x1.66c4d4e554434p-22, 0x1.773ap-5},
{-0x1.9f8ef14d5f6eep-79, -0x1.70700a00fdd55p-24, 0x1.d6ecp-5},
{0x1.452bbce7398c1p-77, 0x1.53002a4e86631p-23, 0x1.1bb3p-4},
{-0x1.990555535afdp-81, 0x1.fcd15f101c142p-25, 0x1.4c56p-4},
{0x1.447e30ad393eep-78, 0x1.25b3eed319cedp-22, 0x1.7d6p-4},
{0x1.b7759da88a2dap-76, -0x1.4195120d8486fp-22, 0x1.960dp-4},
{0x1.cee7766ece702p-78, 0x1.45b878e27d0d9p-23, 0x1.c7b5p-4},
{-0x1.a55c745ecdc2fp-77, 0x1.770744593a4cbp-22, 0x1.f9c9p-4},
{0x1.f7ec992caa67fp-77, 0x1.c673032495d24p-22, 0x1.097ep-3},
{-0x1.433638c6ece3ep-77, -0x1.1eaa65b49696ep-22, 0x1.22dbp-3},
{0x1.58f27b6518824p-76, 0x1.b2866f2850b22p-22, 0x1.3c6f8p-3},
{-0x1.86bdcfdfd4a4cp-79, 0x1.8ee37cd2ea9d3p-25, 0x1.494f8p-3},
{-0x1.ff7044a68a7fap-80, 0x1.7e86f9c2154fbp-24, 0x1.633a8p-3},
{-0x1.aa21694561327p-81, 0x1.8e3cfc25f0ce6p-26, 0x1.7046p-3},
{-0x1.d209f2d4239c6p-87, 0x1.57f7a64ccd537p-28, 0x1.8a898p-3},
{-0x1.a55e97e60e632p-76, -0x1.a761c09fbd2aep-22, 0x1.97c2p-3},
{0x1.261179225541ep-76, 0x1.24bea9a2c66f3p-22, 0x1.b26p-3},
{-0x1.08fa30510fca9p-82, -0x1.60002ccfe43f5p-25, 0x1.bfc68p-3},
{-0x1.63ec8d56242f9p-76, 0x1.69f220e97f22cp-22, 0x1.dac2p-3},
{0x1.8bcdaf0534365p-76, -0x1.6164f64c210ep-22, 0x1.e858p-3},
{0x1.1003282896056p-78, -0x1.0c1678ae89767p-24, 0x1.01d9cp-2},
{0x1.01bcc7025fa92p-78, -0x1.f26a05c813d57p-22, 0x1.08bdp-2},
{-0x1.fe8a8648e9ebcp-80, 0x1.4d8fc561c8d44p-24, 0x1.169cp-2},
{0x1.08dfb23650c75p-79, -0x1.362ad8f7ca2dp-22, 0x1.1d984p-2},
{-0x1.f8d5a89861a5ep-79, 0x1.2b13cd6c4d042p-22, 0x1.249ccp-2},
{-0x1.a1c872983511ep-76, -0x1.1c8f11979a5dbp-22, 0x1.32cp-2},
{0x1.e8e21bff3336bp-77, 0x1.c2ab3edefe569p-23, 0x1.39de8p-2},
{0x1.fd1994fb2c4a1p-80, 0x1.7c3eca28e69cap-26, 0x1.4106p-2},
{0x1.6b94b51cf76b1p-80, -0x1.34c4e99e1c6c6p-24, 0x1.4f6fcp-2},
{-0x1.31d55da1d0f66p-76, -0x1.194a871b63619p-22, 0x1.56b24p-2},
{-0x1.378b22691e28bp-77, 0x1.e3dd5c1c885aep-23, 0x1.5dfdcp-2},
{0x1.99e302970e411p-83, -0x1.6ccf3b1129b7cp-23, 0x1.6552cp-2},
{0x1.20164a049664dp-82, -0x1.2f346e2bf924bp-23, 0x1.6cb1p-2},
{-0x1.d14aac4d864c3p-77, -0x1.fa61aaa59c1d8p-23, 0x1.7b8ap-2},
{0x1.496ab4e4b293fp-79, 0x1.90c11fd32a3abp-22, 0x1.8304cp-2},
{-0x1.d209f2d4239c6p-86, 0x1.57f7a64ccd537p-27, 0x1.8a898p-2},
{0x1.eae3326327babp-81, 0x1.249ba76fee235p-27, 0x1.9218p-2},
{0x1.fa05bddfded8cp-77, -0x1.aad2729b21ae5p-23, 0x1.99b08p-2},
{-0x1.624140d175ba2p-77, 0x1.71810a5e1818p-22, 0x1.a8ff8p-2},
{0x1.f1c5160c515c1p-81, -0x1.6172fe015e13cp-27, 0x1.b0b68p-2},
{-0x1.86a6204eec8cp-79, 0x1.5ec6c1bfbf89ap-24, 0x1.b877cp-2},
{0x1.718f761dd3915p-78, 0x1.678bf6cdedf51p-24, 0x1.c0438p-2},
{-0x1.d4ee66c3700e4p-76, 0x1.c2d45fe43895ep-22, 0x1.c819cp-2},
{-0x1.7d14533586306p-77, -0x1.9ee52ed49d71dp-22, 0x1.cffbp-2},
{0x1.5ce9fb5a7bb5bp-81, 0x1.5786af187a96bp-27, 0x1.d7e6cp-2},
{-0x1.ae6face57ad3bp-77, 0x1.3ab0dc56138c9p-23, 0x1.dfdd8p-2},
{0x1.5ac93b443d55fp-78, 0x1.fe538ab34efb5p-22, 0x1.e7df4p-2},
{0x1.f1753e0ae1e8fp-76, -0x1.e4fee07aa4b68p-22, 0x1.efec8p-2},
{0x1.cdfd4c297069bp-76, -0x1.172f32fe67287p-22, 0x1.f804cp-2},
{0x1.97a0e8f3ba742p-79, -0x1.9a83ff9ab9cc8p-22, 0x1.00144p-1},
{-0x1.800450f5b2357p-78, -0x1.68cb06cece193p-22, 0x1.042bep-1},
{-0x1.a839041241fe7p-78, 0x1.8cd71ddf82e2p-22, 0x1.08494p-1},
{0x1.ed0b8eeccca86p-78, 0x1.5e18ab2df3ae6p-22, 0x1.0c6cap-1},
{0x1.3dd41df9689b3p-79, 0x1.5dee4d9d8a273p-25, 0x1.1096p-1},
{-0x1.990555535afdp-82, 0x1.fcd15f101c142p-26, 0x1.14c56p-1},
{-0x1.1773d02c9055cp-77, -0x1.2474b0f992ba1p-23, 0x1.18faep-1},
{-0x1.4aeef330c53c1p-78, 0x1.4b5a92a606047p-24, 0x1.1d368p-1},
{0x1.8e6ff749ebacbp-77, 0x1.16186fcf54bbdp-22, 0x1.21786p-1},
{0x1.c09d761c548ebp-84, 0x1.18efabeb7d722p-27, 0x1.25c0ap-1},
{0x1.aaa73a428e1e4p-78, -0x1.e5fc7d238691dp-24, 0x1.2a0f4p-1},
{-0x1.af2f3d8b63fbap-79, 0x1.f5809faf6283cp-22, 0x1.2e644p-1},
{-0x1.af2f3d8b63fbap-79, 0x1.f5809faf6283cp-22, 0x1.2e644p-1},
{0x1.78de359f2bb88p-77, 0x1.c6e1dcd0cb449p-22, 0x1.32bfep-1},
{-0x1.415ae1a715618p-76, 0x1.76e0e8f74b4d5p-22, 0x1.37222p-1},
{-0x1.4991b5375621fp-79, -0x1.cb82c89692d99p-24, 0x1.3b8b2p-1},
{-0x1.827d37deb2236p-76, -0x1.63161c5432aebp-22, 0x1.3ffaep-1},
{0x1.9576edac01c78p-77, 0x1.458104c41b901p-22, 0x1.44716p-1},
{0x1.9576edac01c78p-77, 0x1.458104c41b901p-22, 0x1.44716p-1},
{-0x1.05a27b81e2219p-77, -0x1.cd9d0cde578d5p-22, 0x1.48efp-1},
{0x1.237616778b4bap-82, 0x1.b9884591add87p-26, 0x1.4d738p-1},
{0x1.3b7d7e5d148bbp-76, 0x1.c6042978605ffp-22, 0x1.51ff2p-1},
{-0x1.cc3f936a5977cp-79, -0x1.fc4c96b37dcf6p-22, 0x1.56922p-1},
{0x1.20164a049664dp-83, -0x1.2f346e2bf924bp-24, 0x1.5b2c4p-1},
{0x1.20164a049664dp-83, -0x1.2f346e2bf924bp-24, 0x1.5b2c4p-1},
{-0x1.a212919a92f7ap-77, 0x1.c4e4fbb68a4d1p-22, 0x1.5fcdcp-1},
{-0x1.b64b03f7230ddp-77, -0x1.9d499bd9b3226p-23, 0x1.6476ep-1},
{-0x1.1ec6379e6e3b9p-77, -0x1.f89b355ede26fp-23, 0x1.69278p-1},
{-0x1.1ec6379e6e3b9p-77, -0x1.f89b355ede26fp-23, 0x1.69278p-1},
{-0x1.4ba44c03bfbbdp-78, 0x1.53c7e319f6e92p-24, 0x1.6ddfcp-1},
{-0x1.c36fc650d030fp-77, -0x1.b291f070528c7p-22, 0x1.729fep-1},
{-0x1.69e5693a7f067p-80, 0x1.2967a451a7b48p-25, 0x1.7767cp-1},
{-0x1.69e5693a7f067p-80, 0x1.2967a451a7b48p-25, 0x1.7767cp-1},
{0x1.6598aae91499ap-76, 0x1.244fcff690fcep-22, 0x1.7c37ap-1},
{0x1.99d61ec432837p-77, 0x1.46fd97f5dc572p-23, 0x1.810fap-1},
{0x1.99d61ec432837p-77, 0x1.46fd97f5dc572p-23, 0x1.810fap-1},
{0x1.855c42078f81bp-76, -0x1.f3a7352663e5p-22, 0x1.85efep-1},
{-0x1.59408e815107p-77, 0x1.b3cda690370b5p-23, 0x1.8ad84p-1},
{-0x1.59408e815107p-77, 0x1.b3cda690370b5p-23, 0x1.8ad84p-1},
{0x1.33b318085e50ap-78, 0x1.3226b211bf1d9p-23, 0x1.8fc92p-1},
{0x1.343fe7c9cb4aep-79, 0x1.d24b136c101eep-23, 0x1.94c28p-1},
{0x1.343fe7c9cb4aep-79, 0x1.d24b136c101eep-23, 0x1.94c28p-1},
{-0x1.d19522e56fe6p-76, 0x1.7c40c7907e82ap-22, 0x1.99c48p-1},
{-0x1.23b9d8ea55c3ep-77, -0x1.e81781d97ee91p-22, 0x1.9ecf6p-1},
{-0x1.23b9d8ea55c3ep-77, -0x1.e81781d97ee91p-22, 0x1.9ecf6p-1},
{0x1.829440c24aeb6p-78, -0x1.6a77813f94e01p-22, 0x1.a3e3p-1},
{-0x1.624140d175ba2p-76, -0x1.1cfdeb43cfdp-22, 0x1.a8ffap-1},
{-0x1.624140d175ba2p-76, -0x1.1cfdeb43cfdp-22, 0x1.a8ffap-1},
{0x1.afa6f024fb045p-77, -0x1.f983f74d3138fp-23, 0x1.ae256p-1},
{-0x1.603ad3a5d326dp-78, -0x1.e278ae1a1f51fp-23, 0x1.b3546p-1},
{-0x1.603ad3a5d326dp-78, -0x1.e278ae1a1f51fp-23, 0x1.b3546p-1},
{-0x1.0c1e0e5855d6ap-77, -0x1.97552b7b5ea45p-23, 0x1.b88ccp-1},
{-0x1.0c1e0e5855d6ap-77, -0x1.97552b7b5ea45p-23, 0x1.b88ccp-1},
{0x1.c817ad56baa16p-78, -0x1.19b4f3c72c4f8p-24, 0x1.bdceap-1},
{0x1.44c47ac1bf62bp-77, 0x1.f7402d26f1a12p-23, 0x1.c31a2p-1},
{0x1.44c47ac1bf62bp-77, 0x1.f7402d26f1a12p-23, 0x1.c31a2p-1},
{-0x1.69b9465eae1e6p-78, -0x1.2056d5dd31d96p-23, 0x1.c86f8p-1},
{-0x1.69b9465eae1e6p-78, -0x1.2056d5dd31d96p-23, 0x1.c86f8p-1},
{-0x1.24a6d9d1d1904p-79, -0x1.6e46335aae723p-24, 0x1.cdcecp-1},
{-0x1.3826144575ac4p-76, -0x1.beb244c59f331p-22, 0x1.d3382p-1},
{-0x1.3826144575ac4p-76, -0x1.beb244c59f331p-22, 0x1.d3382p-1},
{0x1.dbc96b3b12b25p-81, 0x1.16c071e93fd97p-27, 0x1.d8abap-1},
{0x1.dbc96b3b12b25p-81, 0x1.16c071e93fd97p-27, 0x1.d8abap-1},
{0x1.68a8ccdbd1f33p-77, 0x1.d8175819530c2p-22, 0x1.de298p-1},
{0x1.68a8ccdbd1f33p-77, 0x1.d8175819530c2p-22, 0x1.de298p-1},
{0x1.e586711df5ea1p-79, 0x1.51bd552842c1cp-23, 0x1.e3b2p-1},
{0x1.e586711df5ea1p-79, 0x1.51bd552842c1cp-23, 0x1.e3b2p-1},
{-0x1.bc25adf042483p-79, 0x1.914e204f19d94p-22, 0x1.e9452p-1},
{-0x1.bc25adf042483p-79, 0x1.914e204f19d94p-22, 0x1.e9452p-1},
{0x1.d7d82b65c5686p-76, 0x1.c55d997da24fdp-22, 0x1.eee32p-1},
{0x1.d7d82b65c5686p-76, 0x1.c55d997da24fdp-22, 0x1.eee32p-1},
{-0x1.3f108c0857ca3p-77, -0x1.685c2d2298a6ep-22, 0x1.f48c4p-1},
{-0x1.3f108c0857ca3p-77, -0x1.685c2d2298a6ep-22, 0x1.f48c4p-1},
{-0x1.bd800bca7a221p-78, 0x1.7a4887bd74039p-22, 0x1.fa406p-1},
{0.0, 0.0, 1.0},
};
// Lookup table for 4-bit (16 entries) range reduction for float_eval:
// r(0) = 1.0f
// r(i) = 2^-5 * ceil(2^5 * (1 - 2^-5) / (1 + i * 2^-4)), i = 1..15.
// The constants are chosen so that dx = r * m_x - 1 is exact in
// single precision (with hardware FMA or 14-bit Sterbenz split), and
// -2^-5 <= dx < 2^-4.
//
// Generated by Sollya with:
// > display = hexadecimal;
// > for i from 0 to 15 do {
// if (i == 0) then {
// print("0x1.0p+0f,");
// } else if (i == 15) then {
// print("0x1.0p-1f,");
// } else {
// r = 2^(-5) * ceil( 2^5 * (1 - 2^(-5)) / (1 + i*2^(-4)) );
// print(r @ "f,");
// };
// };
LIBC_INLINE_VAR constexpr float R_16[16] = {
0x1.0p+0f, 0x1.ep-1f, 0x1.cp-1f, 0x1.bp-1f, 0x1.9p-1f, 0x1.8p-1f,
0x1.7p-1f, 0x1.6p-1f, 0x1.5p-1f, 0x1.4p-1f, 0x1.4p-1f, 0x1.3p-1f,
0x1.2p-1f, 0x1.2p-1f, 0x1.1p-1f, 0x1.0p-1f,
};
// Lookup table for 5-bit (32 entries) range reduction for double_eval fast
// path:
// r(0) = 1.0
// r(i) = 2^-6 * ceil(2^6 * (1 - 2^-6) / (1 + i * 2^-5)), i = 1..31.
LIBC_INLINE_VAR constexpr double R_32_D[32] = {
0x1.0p+0, 0x1.fp-1, 0x1.ep-1, 0x1.dp-1, 0x1.cp-1, 0x1.b8p-1, 0x1.bp-1,
0x1.ap-1, 0x1.98p-1, 0x1.9p-1, 0x1.8p-1, 0x1.78p-1, 0x1.7p-1, 0x1.68p-1,
0x1.6p-1, 0x1.58p-1, 0x1.5p-1, 0x1.5p-1, 0x1.48p-1, 0x1.4p-1, 0x1.38p-1,
0x1.38p-1, 0x1.3p-1, 0x1.28p-1, 0x1.2p-1, 0x1.2p-1, 0x1.18p-1, 0x1.18p-1,
0x1.1p-1, 0x1.1p-1, 0x1.08p-1, 0x1.0p-1,
};
// Table of -log2(R_32_D[i]) represented in double precision for double_eval.
//
// Generated by Sollya with:
// > display = hexadecimal;
// > prec = 200;
// > for i from 0 to 31 do {
// if (i == 0) then {
// print("0x0.0000000000000p+0,");
// } else if (i == 31) then {
// print("0x1.0000000000000p+0,");
// } else {
// r = 2^(-6) * ceil( 2^6 * (1 - 2^(-6)) / (1 + i*2^(-5)) );
// print(round(-log2(r), D, RN), ",");
// };
// };
LIBC_INLINE_VAR constexpr double LOG2_R_32[32] = {
0x0.0000000000000p+0, 0x1.77394c9d958d5p-5, 0x1.7d60496cfbb4cp-4,
0x1.22dadc2ab3497p-3, 0x1.8a8980abfbd32p-3, 0x1.bfc67a7fff4ccp-3,
0x1.f5fd8a9063e35p-3, 0x1.32bfee370ee68p-2, 0x1.4f6fbb2cec598p-2,
0x1.6cb0f6865c8eap-2, 0x1.a8ff971810a5ep-2, 0x1.c819dc2d45fe4p-2,
0x1.e7df5fe538ab3p-2, 0x1.042bd4b9a7c99p-1, 0x1.14c560fe68af9p-1,
0x1.25c0a0463bebp-1, 0x1.37222bb70747cp-1, 0x1.37222bb70747cp-1,
0x1.48eef19317991p-1, 0x1.5b2c3da19723bp-1, 0x1.6ddfc2a78fc63p-1,
0x1.6ddfc2a78fc63p-1, 0x1.810fa51bf65fdp-1, 0x1.94c287492c4dbp-1,
0x1.a8ff971810a5ep-1, 0x1.a8ff971810a5ep-1, 0x1.bdce9dcc96187p-1,
0x1.bdce9dcc96187p-1, 0x1.d338120a6dd9dp-1, 0x1.d338120a6dd9dp-1,
0x1.e9452c8a71028p-1, 0x1.0000000000000p+0,
};
// Table of -log2(R_16[i]) represented in FloatFloat precision {lo, hi} for
// float_eval.
// We choose the precision of the high part to be 24 - 8 = 16 bits, so that
// e_xf + LOG2_R_FF_16[i].hi
// is exact in float for |e_x| <= 150.
// For i = 15, R_16[15] = 0.5, so -log2(R_16[15]) = 1.0f is exact.
//
// Generated by Sollya with:
// > display = hexadecimal;
// > prec = 200;
// > for i from 0 to 15 do {
// if (i == 0) then {
// print("{0x0.000000p+0f, 0x0.000000p+0f},");
// } else if (i == 15) then {
// print("{0x0.000000p+0f, 0x1.000000p+0f},");
// } else {
// r = 2^(-5) * ceil( 2^5 * (1 - 2^(-5)) / (1 + i*2^(-4)) );
// l = -log2(r);
// h = round(1 + l, 17, RN) - 1;
// lo = round(l - h, SG, RN);
// print("{" @ lo @ "f, " @ h @ "f},");
// };
// };
LIBC_INLINE_VAR constexpr FloatFloat LOG2_R_FF_16[16] = {
{0x0.000000p+0f, 0x0.000000p+0f}, {0x1.25b3eep-22f, 0x1.7d6p-4f},
{0x1.80abfcp-19f, 0x1.8a88p-3f}, {-0x1.3ab7cep-18f, 0x1.f6p-3f},
{0x1.ed0cbap-19f, 0x1.6cbp-2f}, {-0x1.a39fbep-20f, 0x1.a9p-2f},
{-0x1.40358ep-19f, 0x1.e7ep-2f}, {-0x1.3e032ep-18f, 0x1.14c6p-1f},
{0x1.5db83ap-20f, 0x1.3722p-1f}, {0x1.ed0cbap-20f, 0x1.5b2cp-1f},
{0x1.ed0cbap-20f, 0x1.5b2cp-1f}, {-0x1.6b9026p-19f, 0x1.811p-1f},
{-0x1.a39fbep-19f, 0x1.a9p-1f}, {-0x1.a39fbep-19f, 0x1.a9p-1f},
{0x1.20a6dep-21f, 0x1.d338p-1f}, {0x0.000000p+0f, 0x1.000000p+0f},
};
// Upper bound for y = 150 / |log2(1 - 2^-24)|, generated by Sollya:
// > y = round(-150 / log2(1 - 2^(-24)), SG, RU);
// > y;
// 0x1.9fe368p30
// > printsingle(y);
// 0x4ecff1b4
constexpr uint32_t Y_UPPER_BOUND = 0x4ecf'f1b4;
// Lower bound for y = 2^-25 / 150, generated by Sollya:
// > y = round(2^(-25) / 150, SG, RD);
// > y;
// 0x1.b4e81ap-33
// > printsingle(y);
// 0x2f5a740d
constexpr uint32_t Y_LOWER_BOUND = 0x2f5a'740d;
LIBC_INLINE float set_overflow(Sign sign = Sign::POS) {
fputil::set_errno_if_required(ERANGE);
fputil::raise_overflow_except_if_required<float>();
using FPBits = fputil::FPBits<float>;
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
int rounding = fputil::quick_get_round();
if (rounding == FE_TOWARDZERO)
return FPBits::max_normal(sign).get_val();
if (rounding == FE_DOWNWARD)
return sign.is_neg() ? FPBits::inf(Sign::NEG).get_val()
: FPBits::max_normal(Sign::POS).get_val();
if (rounding == FE_UPWARD)
return sign.is_neg() ? FPBits::max_normal(Sign::NEG).get_val()
: FPBits::inf(Sign::POS).get_val();
#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
return sign.is_neg() ? FPBits::inf(Sign::NEG).get_val()
: FPBits::inf(Sign::POS).get_val();
}
LIBC_INLINE float set_underflow(Sign sign = Sign::POS) {
fputil::set_errno_if_required(ERANGE);
fputil::raise_underflow_except_if_required<float>();
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
using FPBits = fputil::FPBits<float>;
int rounding = fputil::quick_get_round();
if (rounding == FE_UPWARD && sign.is_pos())
return FPBits::min_subnormal(Sign::POS).get_val();
if (rounding == FE_DOWNWARD && sign.is_neg())
return FPBits::min_subnormal(Sign::NEG).get_val();
#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
return sign.is_neg() ? -0.0f : 0.0f;
}
// Fast checks for special inputs:
// x = 0, +-1, 2^k, 10, +- Inf
// y = 0, +-1, 2, 0.5
LIBC_ALWAYS_INLINE cpp::optional<float> check_special_inputs(float x, float y) {
using FPBits = fputil::FPBits<float>;
FPBits xbits(x), ybits(y);
bool x_sign = xbits.sign() == Sign::NEG;
bool y_sign = ybits.sign() == Sign::NEG;
FPBits x_abs = xbits.abs();
FPBits y_abs = ybits.abs();
// If x or y is signaling NaN
if (x_abs.is_signaling_nan() || y_abs.is_signaling_nan()) {
fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}
if (x == 1.0f || y == 0.0f)
return 1.0f;
if (x == 0.0f) {
if (y_abs.is_nan())
return y;
if (y_abs.is_inf())
return y_sign ? FPBits::inf().get_val() : 0.0f;
bool out_is_neg = x_sign && fputil::is_odd_integer(y);
if (y_sign) {
// pow(0, negative number) = Inf
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_DIVBYZERO);
return out_is_neg ? FPBits::inf(Sign::NEG).get_val()
: FPBits::inf(Sign::POS).get_val();
}
// pow(0, positive number) = 0
return out_is_neg ? -0.0f : 0.0f;
}
if (y == 1.0f)
return x;
if (y == 2.0f)
return x * x;
#ifdef LIBC_TARGET_CPU_HAS_FPU_FLOAT
if (y == 0.5f && !x_sign)
return fputil::sqrt<float>(x);
#endif // LIBC_TARGET_CPU_HAS_FPU_FLOAT
// Exact power of 2:
// pow(2^k, y) = 2^(k * y).
// This is exact if k * y is an integer.
uint32_t y_a = y_abs.uintval();
if (x_abs.is_normal() && xbits.get_mantissa() == 0 && y_a > Y_LOWER_BOUND &&
y_a < Y_UPPER_BOUND) {
if (!x_sign || fputil::is_integer(y)) {
Sign out_sign =
(x_sign && fputil::is_odd_integer(y)) ? Sign::NEG : Sign::POS;
int e_x = xbits.get_exponent();
double ey = static_cast<double>(e_x) * static_cast<double>(y);
double ey_int = fputil::nearest_integer(ey);
if (ey_int == ey) {
if (ey > 127.0)
return set_overflow(out_sign);
if (ey < -149.0)
return set_underflow(out_sign);
int k = static_cast<int>(ey);
if (k >= -126)
return FPBits::create_value(out_sign, static_cast<uint32_t>(k + 127),
0)
.get_val();
return FPBits::create_value(out_sign, 0, 1U << (k + 149)).get_val();
}
}
}
#if !defined(LIBC_MATH_HAS_SMALL_TABLES) && \
!defined(LIBC_MATH_HAS_INTERMEDIATE_COMP_IN_FLOAT)
if (x == 2.0f)
return math::exp2f(y);
if (x == 10.0f)
return math::exp10f(y);
// For x = 2^(+- 2^n):
// x^y = 2^(+- 2^n * y).
if (!x_sign && x_abs.is_normal() && xbits.get_mantissa() == 0 &&
y_a > Y_LOWER_BOUND && y_a < Y_UPPER_BOUND) {
int e_x = xbits.get_exponent();
uint32_t abs_e = static_cast<uint32_t>(e_x > 0 ? e_x : -e_x);
if (cpp::has_single_bit(abs_e)) {
float hi = static_cast<float>(e_x) * y;
if (hi >= 128.0f)
return set_overflow();
if (hi < -150.0f)
return set_underflow();
if (hi >= -126.0f)
return math::exp2f(hi);
}
}
#endif // !LIBC_MATH_HAS_SMALL_TABLES &&
// !LIBC_MATH_HAS_INTERMEDIATE_COMP_IN_FLOAT
return cpp::nullopt;
}
// Filters out extreme input ranges, infinities, NaNs, normalizes denormal
// inputs, and handles negative bases. Returns a value if the result is
// determined, or nullopt if regular evaluation should proceed (in which case x,
// y, ex, sign may be updated).
template <typename SignType = uint64_t>
LIBC_ALWAYS_INLINE cpp::optional<float>
check_exceptional_cases(float &x, float &y, int &ex, SignType &sign) {
using FloatBits = fputil::FPBits<float>;
FloatBits xbits(x), ybits(y);
bool x_sign = xbits.sign() == Sign::NEG;
bool y_sign = ybits.sign() == Sign::NEG;
FloatBits x_abs = xbits.abs();
FloatBits y_abs = ybits.abs();
uint32_t x_a = x_abs.uintval();
uint32_t y_a = y_abs.uintval();
// If x or y is signaling NaN
if (x_abs.is_signaling_nan() || y_abs.is_signaling_nan()) {
fputil::raise_except_if_required(FE_INVALID);
return FloatBits::quiet_nan().get_val();
}
if (x == 1.0f || y == 0.0f)
return 1.0f;
// Extreme |y| >= Y_UPPER_BOUND: |y * log2(x)| >= 150 for all x != 1.
if (y_a >= Y_UPPER_BOUND) {
if (x_abs.is_nan())
return x;
if (ybits.is_nan())
return y;
if (x_a == FloatBits::one().uintval())
return 1.0f;
bool is_overflow = (x_a < FloatBits::one().uintval()) == y_sign;
if (y_abs.is_inf() || x_a == FloatBits::inf().uintval())
return is_overflow ? FloatBits::inf().get_val() : 0.0f;
return is_overflow ? set_overflow() : set_underflow();
}
// y is finite and non-zero.
if (x_a > FloatBits::inf().uintval()) {
// x is NaN.
return x;
}
// Handle negative base x < 0 (including -0.0f and -Inf):
// - If y is an integer: (-|x|)^y = (-1)^y * |x|^y.
// - If y is not an integer and x is finite non-zero: (-|x|)^y is undefined
// in real numbers; raises FE_INVALID, sets errno to EDOM, and returns
// quiet NaN.
if (x_sign) {
if (fputil::is_integer(y)) {
if (fputil::is_odd_integer(y)) {
if constexpr (sizeof(SignType) == 8)
sign = 0x8000'0000'0000'0000ULL;
else
sign = 0x8000'0000U;
}
} else if (x != 0.0f && x_a != FloatBits::inf().uintval()) {
// pow( negative, non-integer ) = NaN
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_INVALID);
return FloatBits::quiet_nan().get_val();
}
x = -x;
}
if (x_a == FloatBits::inf().uintval()) {
// pow(+-Inf, y):
// - y < 0: returns +-0.0 (negative if x = -Inf and y is odd integer).
// - y > 0: returns +-Inf (negative if x = -Inf and y is odd integer).
if (y_sign)
return (sign != 0) ? -0.0f : 0.0f;
return (sign != 0) ? FloatBits::inf(Sign::NEG).get_val()
: FloatBits::inf(Sign::POS).get_val();
}
if (x == 0.0f) {
if (y_sign) {
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_DIVBYZERO);
return (sign != 0) ? FloatBits::inf(Sign::NEG).get_val()
: FloatBits::inf(Sign::POS).get_val();
}
return (sign != 0) ? -0.0f : 0.0f;
}
if (y_a <= Y_LOWER_BOUND) {
#if defined(LIBC_MATH_HAS_NO_EXCEPT) && \
defined(LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY)
return 1.0f;
#else
volatile float one = 1.0f;
volatile float eps = ((x_a < FloatBits::one().uintval()) == !y_sign)
? -0x1.0p-50f
: 0x1.0p-50f;
return one + eps;
#endif
}
// Normalize denormal inputs.
if (x_a < FloatBits::min_normal().uintval()) {
int shift = cpp::countl_zero(x_a) - 8;
#ifdef LIBC_TARGET_CPU_HAS_FPU_DOUBLE
if (!y_sign &&
static_cast<double>(y) * static_cast<double>(125 + shift) >= 150.0)
return set_underflow(sign == 0 ? Sign::POS : Sign::NEG);
#else
uint32_t log2_bound = static_cast<uint32_t>(125 + shift);
if (!y_sign &&
(y_a >= 0x4000'0000U || (y_a >= FloatBits::one().uintval() &&
(y_a - 0x3f00'0000U) * log2_bound >=
(150U << FloatBits::FRACTION_LEN))))
return set_underflow(sign == 0 ? Sign::POS : Sign::NEG);
#endif
ex -= shift;
x = cpp::bit_cast<float>(x_a << shift);
}
return cpp::nullopt;
}
#ifndef LIBC_MATH_HAS_SMALL_TABLES
// Check if x^y is an exact rounding boundary case (a 25-bit dyadic float,
// which is either an exact 24-bit float or an exact halfway midpoint).
//
// Reference:
// Lauter, C. and Lefevre, V., "Rounding Boundary Cases for Values of the
// Exponential and Power Functions," IEEE Trans. Comput. 58(8):1063-1074.
//
// 1. x = 2^e: x^y = 2^(e * y) is in F_25 iff e * y is an integer.
// 2. x = m * 2^e, y = n * 2^f with m, n odd integers:
// x^y in F_25 implies:
// - 0 <= y <= 15, n <= 15.
// - f >= -3.
// - e * y is an integer.
// - When f < 0 (i.e. y = n / 2^(-f)): m^(2^f) is an integer, verified
// by taking repeated square roots -f times.
// - The resulting significand (m^(2^f))^n <= 2^25.
//
// Returns true and sets exact_m and exact_exp if x^y is an exact boundary.
LIBC_INLINE bool is_exact_rounding_boundary(float x, float y, uint32_t &exact_m,
int &exact_exp) {
using FPBits = fputil::FPBits<float>;
FPBits xbits(x);
int x_e = 0;
uint32_t x_mant = xbits.get_mantissa();
if (LIBC_UNLIKELY(xbits.get_biased_exponent() == 0)) {
if (x_mant != 0) {
int shift = cpp::countl_zero(x_mant) - 8;
x_mant = (x_mant << shift) & FPBits::FRACTION_MASK;
x_e = -126 - shift;
}
} else {
x_e = xbits.get_exponent();
}
// Case 1: x = 2^e.
if (x_mant == 0) {
double e = static_cast<double>(x_e);
double ey = e * static_cast<double>(y);
if (fputil::nearest_integer(ey) == ey) {
exact_m = 1;
if (ey > 150.0)
exact_exp = 200;
else if (ey < -200.0)
exact_exp = -200;
else
exact_exp = static_cast<int>(ey);
return true;
}
return false;
}
// Case 2: x is not a power of 2.
if (y < 0.0f || y > 15.0f)
return false;
// Decompose y = n * 2^f with n odd.
FPBits ybits(y);
uint32_t y_mant = ybits.get_mantissa() | (1U << FPBits::FRACTION_LEN);
int y_exp = ybits.get_exponent() - static_cast<int>(FPBits::FRACTION_LEN);
int tz_y = cpp::countr_zero(y_mant);
uint32_t n = y_mant >> tz_y;
int f = y_exp + tz_y;
if (n > 15 || f < -3)
return false;
// Decompose x = m * 2^e with m odd.
uint32_t full_x_mant = x_mant | (1U << FPBits::FRACTION_LEN);
int tz_x = cpp::countr_zero(full_x_mant);
uint32_t m = full_x_mant >> tz_x;
int e = x_e - static_cast<int>(FPBits::FRACTION_LEN) + tz_x;
if (f < 0) {
// Non-integer exponent: y = n / 2^(-f) with -f in {1, 2, 3}.
double ey = static_cast<double>(e) * static_cast<double>(y);
if (fputil::nearest_integer(ey) != ey)
return false;
// Check if m^(2^f) is an integer by taking repeated square roots -f times.
int count = -f;
uint32_t cur = m;
for (int i = 0; i < count; ++i) {
uint32_t s =
static_cast<uint32_t>(fputil::sqrt<double>(static_cast<double>(cur)));
if (s * s != cur)
return false;
cur = s;
}
// Compute res = cur^n and check if res <= 2^25.
uint32_t res = 1;
for (uint32_t i = 0; i < n; ++i) {
if (res > (1U << 25) / cur)
return false;
res *= cur;
}
exact_m = res;
exact_exp = static_cast<int>(ey);
return true;
} else {
// Integer exponent: f >= 0, so y is an integer.
int y_int = static_cast<int>(y);
uint32_t res = 1;
for (int i = 0; i < y_int; ++i) {
if (res > (1U << 25) / m)
return false;
res *= m;
}
exact_m = res;
exact_exp = e * y_int;
return true;
}
}
#endif // !LIBC_MATH_HAS_SMALL_TABLES
} // namespace powf_internal
} // namespace math
} // namespace LIBC_NAMESPACE_DECL
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POWF_UTILS_H