libc/src/__support/math/exp_double_utils.h - llvm-project - Git at Google


            
            //===----------------------------------------------------------------------===//

            
            //

            
            // 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

            
            /// Shared double-precision utilities and tables for exponential functions when

            
            /// LIBC_MATH_HAS_SKIP_ACCURATE_PASS is enabled.

            
            ///

            
            //===----------------------------------------------------------------------===//


            #ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_DOUBLE_UTILS_H

            
            #define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_DOUBLE_UTILS_H


            #include "hdr/errno_macros.h"

            
            #include "src/__support/CPP/bit.h"

            
            #include "src/__support/FPUtil/FEnvImpl.h"

            
            #include "src/__support/FPUtil/FPBits.h"

            
            #include "src/__support/FPUtil/multiply_add.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/math/exp_utils.h"


            namespace LIBC_NAMESPACE_DECL {

            
            namespace math {

            
            namespace double_eval {


            #ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS


            #ifdef LIBC_MATH_HAS_SMALL_TABLES


            LIBC_INLINE_VAR constexpr int EXP_N_BITS = 3;

            
            LIBC_INLINE_VAR constexpr double EXP_SCALE = 8.0;

            
            LIBC_INLINE_VAR constexpr double EXP_SCALE_INV = 0x1.0p-3;

            
            LIBC_INLINE_VAR constexpr int EXP_MASK = 0x07;


            // Degree-5 polynomial approximation for (exp(x) - 1 - x) / x^2 on

            
            // [-2^-4 * log(2), 2^-4 * log(2)]:

            
            // Generated by Sollya with:

            
            // > display = hexadecimal!;

            
            // > prec = 500;

            
            // > P = fpminimax((exp(x) - 1 - x)/x^2, 5, [|D...|],

            
            //                 [-2^(-4) * log(2), 2^(-4) * log(2)]);

            
            // > dirtyinfnorm(exp(x) - (1 + x + x^2 * P),

            
            //                [-2^(-4) * log(2), 2^(-4) * log(2)]);

            
            //   0x1.6828...p-57 < 1.5 * 2^-57.

            
            LIBC_INLINE_VAR constexpr double EXP_P_COEFFS[6] = {

            
                0x1.000000000002ep-1,  // p_0

            
                0x1.555555555549cp-3,  // p_1

            
                0x1.5555554e95a5p-5,   // p_2

            
                0x1.11111114c1e11p-7,  // p_3

            
                0x1.6c1b70af254c4p-10, // p_4

            
                0x1.a01b62c901438p-13, // p_5

            
            };


            // Evaluate P(u) ~ (e^u - 1 - u) / u^2 using Estrin's scheme:

            
            LIBC_ALWAYS_INLINE double exp_poly(double u, double u2) {

            
              double c0 = fputil::multiply_add(u, EXP_P_COEFFS[1], EXP_P_COEFFS[0]);

            
              double c1 = fputil::multiply_add(u, EXP_P_COEFFS[3], EXP_P_COEFFS[2]);

            
              double u4 = u2 * u2;

            
              double c2 = fputil::multiply_add(u, EXP_P_COEFFS[5], EXP_P_COEFFS[4]);

            
              double c01 = fputil::multiply_add(u2, c1, c0);

            
              return fputil::multiply_add(u4, c2, c01);

            
            }


            // Lookup table for 2^(k * 2^-3) with k = 0..7.

            
            // Generated by Sollya with:

            
            // > display = hexadecimal!;

            
            // > for i from 0 to 7 do print(round(2^(i * 2^-3), D, RN), ",");

            
            LIBC_INLINE_VAR constexpr double EXP2_MID[8] = {

            
                0x1.0000000000000p0, 0x1.172b83c7d517bp0, 0x1.306fe0a31b715p0,

            
                0x1.4bfdad5362a27p0, 0x1.6a09e667f3bcdp0, 0x1.8ace5422aa0dbp0,

            
                0x1.ae89f995ad3adp0, 0x1.d5818dcfba487p0,

            
            };


            #else // !LIBC_MATH_HAS_SMALL_TABLES


            // Lookup table for 2^(k * 2^-7) with k = 0..127.

            
            // Generated by Sollya with:

            
            // > display = hexadecimal!;

            
            // > for i from 0 to 127 do print(round(2^(i * 2^-7), D, RN), ",");

            
            LIBC_INLINE_VAR constexpr double EXP2_MID[128] = {

            
                0x1.0000000000000p0, 0x1.0163da9fb3335p0, 0x1.02c9a3e778061p0,

            
                0x1.04315e86e7f85p0, 0x1.059b0d3158574p0, 0x1.0706b29ddf6dep0,

            
                0x1.0874518759bc8p0, 0x1.09e3ecac6f383p0, 0x1.0b5586cf9890fp0,

            
                0x1.0cc922b7247f7p0, 0x1.0e3ec32d3d1a2p0, 0x1.0fb66affed31bp0,

            
                0x1.11301d0125b51p0, 0x1.12abdc06c31ccp0, 0x1.1429aaea92dep0,

            
                0x1.15a98c8a58e51p0, 0x1.172b83c7d517bp0, 0x1.18af9388c8deap0,

            
                0x1.1a35beb6fcb75p0, 0x1.1bbe084045cd4p0, 0x1.1d4873168b9aap0,

            
                0x1.1ed5022fcd91dp0, 0x1.2063b88628cd6p0, 0x1.21f49917ddc96p0,

            
                0x1.2387a6e756238p0, 0x1.251ce4fb2a63fp0, 0x1.26b4565e27cddp0,

            
                0x1.284dfe1f56381p0, 0x1.29e9df51fdee1p0, 0x1.2b87fd0dad99p0,

            
                0x1.2d285a6e4030bp0, 0x1.2ecafa93e2f56p0, 0x1.306fe0a31b715p0,

            
                0x1.32170fc4cd831p0, 0x1.33c08b26416ffp0, 0x1.356c55f929ff1p0,

            
                0x1.371a7373aa9cbp0, 0x1.38cae6d05d866p0, 0x1.3a7db34e59ff7p0,

            
                0x1.3c32dc313a8e5p0, 0x1.3dea64c123422p0, 0x1.3fa4504ac801cp0,

            
                0x1.4160a21f72e2ap0, 0x1.431f5d950a897p0, 0x1.44e086061892dp0,

            
                0x1.46a41ed1d0057p0, 0x1.486a2b5c13cdp0,  0x1.4a32af0d7d3dep0,

            
                0x1.4bfdad5362a27p0, 0x1.4dcb299fddd0dp0, 0x1.4f9b2769d2ca7p0,

            
                0x1.516daa2cf6642p0, 0x1.5342b569d4f82p0, 0x1.551a4ca5d920fp0,

            
                0x1.56f4736b527dap0, 0x1.58d12d497c7fdp0, 0x1.5ab07dd485429p0,

            
                0x1.5c9268a5946b7p0, 0x1.5e76f15ad2148p0, 0x1.605e1b976dc09p0,

            
                0x1.6247eb03a5585p0, 0x1.6434634ccc32p0,  0x1.6623882552225p0,

            
                0x1.68155d44ca973p0, 0x1.6a09e667f3bcdp0, 0x1.6c012750bdabfp0,

            
                0x1.6dfb23c651a2fp0, 0x1.6ff7df9519484p0, 0x1.71f75e8ec5f74p0,

            
                0x1.73f9a48a58174p0, 0x1.75feb564267c9p0, 0x1.780694fde5d3fp0,

            
                0x1.7a11473eb0187p0, 0x1.7c1ed0130c132p0, 0x1.7e2f336cf4e62p0,

            
                0x1.80427543e1a12p0, 0x1.82589994cce13p0, 0x1.8471a4623c7adp0,

            
                0x1.868d99b4492edp0, 0x1.88ac7d98a6699p0, 0x1.8ace5422aa0dbp0,

            
                0x1.8cf3216b5448cp0, 0x1.8f1ae99157736p0, 0x1.9145b0b91ffc6p0,

            
                0x1.93737b0cdc5e5p0, 0x1.95a44cbc8520fp0, 0x1.97d829fde4e5p0,

            
                0x1.9a0f170ca07bap0, 0x1.9c49182a3f09p0,  0x1.9e86319e32323p0,

            
                0x1.a0c667b5de565p0, 0x1.a309bec4a2d33p0, 0x1.a5503b23e255dp0,

            
                0x1.a799e1330b358p0, 0x1.a9e6b5579fdbfp0, 0x1.ac36bbfd3f37ap0,

            
                0x1.ae89f995ad3adp0, 0x1.b0e07298db666p0, 0x1.b33a2b84f15fbp0,

            
                0x1.b59728de5593ap0, 0x1.b7f76f2fb5e47p0, 0x1.ba5b030a1064ap0,

            
                0x1.bcc1e904bc1d2p0, 0x1.bf2c25bd71e09p0, 0x1.c199bdd85529cp0,

            
                0x1.c40ab5fffd07ap0, 0x1.c67f12e57d14bp0, 0x1.c8f6d9406e7b5p0,

            
                0x1.cb720dcef9069p0, 0x1.cdf0b555dc3fap0, 0x1.d072d4a07897cp0,

            
                0x1.d2f87080d89f2p0, 0x1.d5818dcfba487p0, 0x1.d80e316c98398p0,

            
                0x1.da9e603db3285p0, 0x1.dd321f301b46p0,  0x1.dfc97337b9b5fp0,

            
                0x1.e264614f5a129p0, 0x1.e502ee78b3ff6p0, 0x1.e7a51fbc74c83p0,

            
                0x1.ea4afa2a490dap0, 0x1.ecf482d8e67f1p0, 0x1.efa1bee615a27p0,

            
                0x1.f252b376bba97p0, 0x1.f50765b6e454p0,  0x1.f7bfdad9cbe14p0,

            
                0x1.fa7c1819e90d8p0, 0x1.fd3c22b8f71f1p0,

            
            };


            #endif // LIBC_MATH_HAS_SMALL_TABLES


            LIBC_NOINLINE LIBC_INLINE double exp_denorm(double x, int hi, double r,

            
                                                        bool is_exp2) {

            
              using FPBits = typename fputil::FPBits<double>;

            
              if (is_exp2 && LIBC_UNLIKELY(static_cast<double>(hi) == x))

            
                return FPBits(uint64_t(1) << (hi + 1074)).get_val();


              double res = ziv_test_denorm</*SKIP_ZIV_TEST=*/true>(hi, r, 0.0, 0.0).value();

            
              if (LIBC_UNLIKELY(FPBits(res).is_normal()))

            
                return res;


              fputil::set_errno_if_required(ERANGE);

            
              fputil::raise_underflow_except_if_required<double>();

            
              return res;

            
            }


            LIBC_MATH_SMALL_TABLES_INLINE double

            
            exp2_eval(double x, double u, int hi, unsigned idx, bool is_exp2 = false) {

            
              double exp2_mid = EXP2_MID[idx];

            
              double u2 = u * u;

            
              double p = exp_poly(u, u2);

            
              // u + u^2 * p

            
              double poly_lo = fputil::multiply_add(u2, p, u);

            
              double r = fputil::multiply_add(exp2_mid, poly_lo, exp2_mid);


              if (LIBC_UNLIKELY(hi <= -1022))

            
                return exp_denorm(x, hi, r, is_exp2);


              return scale_exp2(r, hi);

            
            }


            #endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS


            } // namespace double_eval

            
            } // namespace math

            
            } // namespace LIBC_NAMESPACE_DECL


            #endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_DOUBLE_UTILS_H
	//===----------------------------------------------------------------------===//
	//
	// 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
	/// Shared double-precision utilities and tables for exponential functions when
	/// LIBC_MATH_HAS_SKIP_ACCURATE_PASS is enabled.
	///
	//===----------------------------------------------------------------------===//

	#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_DOUBLE_UTILS_H
	#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_DOUBLE_UTILS_H

	#include "hdr/errno_macros.h"
	#include "src/__support/CPP/bit.h"
	#include "src/__support/FPUtil/FEnvImpl.h"
	#include "src/__support/FPUtil/FPBits.h"
	#include "src/__support/FPUtil/multiply_add.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/math/exp_utils.h"

	namespace LIBC_NAMESPACE_DECL {
	namespace math {
	namespace double_eval {

	#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS

	#ifdef LIBC_MATH_HAS_SMALL_TABLES

	LIBC_INLINE_VAR constexpr int EXP_N_BITS = 3;
	LIBC_INLINE_VAR constexpr double EXP_SCALE = 8.0;
	LIBC_INLINE_VAR constexpr double EXP_SCALE_INV = 0x1.0p-3;
	LIBC_INLINE_VAR constexpr int EXP_MASK = 0x07;

	// Degree-5 polynomial approximation for (exp(x) - 1 - x) / x^2 on
	// [-2^-4 * log(2), 2^-4 * log(2)]:
	// Generated by Sollya with:
	// > display = hexadecimal!;
	// > prec = 500;
	// > P = fpminimax((exp(x) - 1 - x)/x^2, 5, [\|D...\|],
	// [-2^(-4) * log(2), 2^(-4) * log(2)]);
	// > dirtyinfnorm(exp(x) - (1 + x + x^2 * P),
	// [-2^(-4) * log(2), 2^(-4) * log(2)]);
	// 0x1.6828...p-57 < 1.5 * 2^-57.
	LIBC_INLINE_VAR constexpr double EXP_P_COEFFS[6] = {
	0x1.000000000002ep-1, // p_0
	0x1.555555555549cp-3, // p_1
	0x1.5555554e95a5p-5, // p_2
	0x1.11111114c1e11p-7, // p_3
	0x1.6c1b70af254c4p-10, // p_4
	0x1.a01b62c901438p-13, // p_5
	};

	// Evaluate P(u) ~ (e^u - 1 - u) / u^2 using Estrin's scheme:
	LIBC_ALWAYS_INLINE double exp_poly(double u, double u2) {
	double c0 = fputil::multiply_add(u, EXP_P_COEFFS[1], EXP_P_COEFFS[0]);
	double c1 = fputil::multiply_add(u, EXP_P_COEFFS[3], EXP_P_COEFFS[2]);
	double u4 = u2 * u2;
	double c2 = fputil::multiply_add(u, EXP_P_COEFFS[5], EXP_P_COEFFS[4]);
	double c01 = fputil::multiply_add(u2, c1, c0);
	return fputil::multiply_add(u4, c2, c01);
	}

	// Lookup table for 2^(k * 2^-3) with k = 0..7.
	// Generated by Sollya with:
	// > display = hexadecimal!;
	// > for i from 0 to 7 do print(round(2^(i * 2^-3), D, RN), ",");
	LIBC_INLINE_VAR constexpr double EXP2_MID[8] = {
	0x1.0000000000000p0, 0x1.172b83c7d517bp0, 0x1.306fe0a31b715p0,
	0x1.4bfdad5362a27p0, 0x1.6a09e667f3bcdp0, 0x1.8ace5422aa0dbp0,
	0x1.ae89f995ad3adp0, 0x1.d5818dcfba487p0,
	};

	#else // !LIBC_MATH_HAS_SMALL_TABLES

	// Lookup table for 2^(k * 2^-7) with k = 0..127.
	// Generated by Sollya with:
	// > display = hexadecimal!;
	// > for i from 0 to 127 do print(round(2^(i * 2^-7), D, RN), ",");
	LIBC_INLINE_VAR constexpr double EXP2_MID[128] = {
	0x1.0000000000000p0, 0x1.0163da9fb3335p0, 0x1.02c9a3e778061p0,
	0x1.04315e86e7f85p0, 0x1.059b0d3158574p0, 0x1.0706b29ddf6dep0,
	0x1.0874518759bc8p0, 0x1.09e3ecac6f383p0, 0x1.0b5586cf9890fp0,
	0x1.0cc922b7247f7p0, 0x1.0e3ec32d3d1a2p0, 0x1.0fb66affed31bp0,
	0x1.11301d0125b51p0, 0x1.12abdc06c31ccp0, 0x1.1429aaea92dep0,
	0x1.15a98c8a58e51p0, 0x1.172b83c7d517bp0, 0x1.18af9388c8deap0,
	0x1.1a35beb6fcb75p0, 0x1.1bbe084045cd4p0, 0x1.1d4873168b9aap0,
	0x1.1ed5022fcd91dp0, 0x1.2063b88628cd6p0, 0x1.21f49917ddc96p0,
	0x1.2387a6e756238p0, 0x1.251ce4fb2a63fp0, 0x1.26b4565e27cddp0,
	0x1.284dfe1f56381p0, 0x1.29e9df51fdee1p0, 0x1.2b87fd0dad99p0,
	0x1.2d285a6e4030bp0, 0x1.2ecafa93e2f56p0, 0x1.306fe0a31b715p0,
	0x1.32170fc4cd831p0, 0x1.33c08b26416ffp0, 0x1.356c55f929ff1p0,
	0x1.371a7373aa9cbp0, 0x1.38cae6d05d866p0, 0x1.3a7db34e59ff7p0,
	0x1.3c32dc313a8e5p0, 0x1.3dea64c123422p0, 0x1.3fa4504ac801cp0,
	0x1.4160a21f72e2ap0, 0x1.431f5d950a897p0, 0x1.44e086061892dp0,
	0x1.46a41ed1d0057p0, 0x1.486a2b5c13cdp0, 0x1.4a32af0d7d3dep0,
	0x1.4bfdad5362a27p0, 0x1.4dcb299fddd0dp0, 0x1.4f9b2769d2ca7p0,
	0x1.516daa2cf6642p0, 0x1.5342b569d4f82p0, 0x1.551a4ca5d920fp0,
	0x1.56f4736b527dap0, 0x1.58d12d497c7fdp0, 0x1.5ab07dd485429p0,
	0x1.5c9268a5946b7p0, 0x1.5e76f15ad2148p0, 0x1.605e1b976dc09p0,
	0x1.6247eb03a5585p0, 0x1.6434634ccc32p0, 0x1.6623882552225p0,
	0x1.68155d44ca973p0, 0x1.6a09e667f3bcdp0, 0x1.6c012750bdabfp0,
	0x1.6dfb23c651a2fp0, 0x1.6ff7df9519484p0, 0x1.71f75e8ec5f74p0,
	0x1.73f9a48a58174p0, 0x1.75feb564267c9p0, 0x1.780694fde5d3fp0,
	0x1.7a11473eb0187p0, 0x1.7c1ed0130c132p0, 0x1.7e2f336cf4e62p0,
	0x1.80427543e1a12p0, 0x1.82589994cce13p0, 0x1.8471a4623c7adp0,
	0x1.868d99b4492edp0, 0x1.88ac7d98a6699p0, 0x1.8ace5422aa0dbp0,
	0x1.8cf3216b5448cp0, 0x1.8f1ae99157736p0, 0x1.9145b0b91ffc6p0,
	0x1.93737b0cdc5e5p0, 0x1.95a44cbc8520fp0, 0x1.97d829fde4e5p0,
	0x1.9a0f170ca07bap0, 0x1.9c49182a3f09p0, 0x1.9e86319e32323p0,
	0x1.a0c667b5de565p0, 0x1.a309bec4a2d33p0, 0x1.a5503b23e255dp0,
	0x1.a799e1330b358p0, 0x1.a9e6b5579fdbfp0, 0x1.ac36bbfd3f37ap0,
	0x1.ae89f995ad3adp0, 0x1.b0e07298db666p0, 0x1.b33a2b84f15fbp0,
	0x1.b59728de5593ap0, 0x1.b7f76f2fb5e47p0, 0x1.ba5b030a1064ap0,
	0x1.bcc1e904bc1d2p0, 0x1.bf2c25bd71e09p0, 0x1.c199bdd85529cp0,
	0x1.c40ab5fffd07ap0, 0x1.c67f12e57d14bp0, 0x1.c8f6d9406e7b5p0,
	0x1.cb720dcef9069p0, 0x1.cdf0b555dc3fap0, 0x1.d072d4a07897cp0,
	0x1.d2f87080d89f2p0, 0x1.d5818dcfba487p0, 0x1.d80e316c98398p0,
	0x1.da9e603db3285p0, 0x1.dd321f301b46p0, 0x1.dfc97337b9b5fp0,
	0x1.e264614f5a129p0, 0x1.e502ee78b3ff6p0, 0x1.e7a51fbc74c83p0,
	0x1.ea4afa2a490dap0, 0x1.ecf482d8e67f1p0, 0x1.efa1bee615a27p0,
	0x1.f252b376bba97p0, 0x1.f50765b6e454p0, 0x1.f7bfdad9cbe14p0,
	0x1.fa7c1819e90d8p0, 0x1.fd3c22b8f71f1p0,
	};

	#endif // LIBC_MATH_HAS_SMALL_TABLES

	LIBC_NOINLINE LIBC_INLINE double exp_denorm(double x, int hi, double r,
	bool is_exp2) {
	using FPBits = typename fputil::FPBits<double>;
	if (is_exp2 && LIBC_UNLIKELY(static_cast<double>(hi) == x))
	return FPBits(uint64_t(1) << (hi + 1074)).get_val();

	double res = ziv_test_denorm</SKIP_ZIV_TEST=/true>(hi, r, 0.0, 0.0).value();
	if (LIBC_UNLIKELY(FPBits(res).is_normal()))
	return res;

	fputil::set_errno_if_required(ERANGE);
	fputil::raise_underflow_except_if_required<double>();
	return res;
	}

	LIBC_MATH_SMALL_TABLES_INLINE double
	exp2_eval(double x, double u, int hi, unsigned idx, bool is_exp2 = false) {
	double exp2_mid = EXP2_MID[idx];
	double u2 = u * u;
	double p = exp_poly(u, u2);
	// u + u^2 * p
	double poly_lo = fputil::multiply_add(u2, p, u);
	double r = fputil::multiply_add(exp2_mid, poly_lo, exp2_mid);

	if (LIBC_UNLIKELY(hi <= -1022))
	return exp_denorm(x, hi, r, is_exp2);

	return scale_exp2(r, hi);
	}

	#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS

	} // namespace double_eval
	} // namespace math
	} // namespace LIBC_NAMESPACE_DECL

	#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_DOUBLE_UTILS_H