| //===-- Implementation header for lgammaf16 ---------------------*- C++ -*-===// |
| // |
| // 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 |
| // |
| //===----------------------------------------------------------------------===// |
| |
| #ifndef LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMAF16_H |
| #define LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMAF16_H |
| |
| #include "include/llvm-libc-macros/float16-macros.h" |
| |
| #ifdef LIBC_TYPES_HAS_FLOAT16 |
| |
| #include "src/__support/CPP/bit.h" |
| #include "src/__support/FPUtil/FEnvImpl.h" |
| #include "src/__support/FPUtil/FPBits.h" |
| #include "src/__support/FPUtil/NearestIntegerOperations.h" |
| #include "src/__support/FPUtil/PolyEval.h" |
| #include "src/__support/FPUtil/cast.h" |
| #include "src/__support/FPUtil/multiply_add.h" |
| #include "src/__support/macros/config.h" |
| #include "src/__support/macros/optimization.h" |
| #include "src/__support/macros/properties/types.h" |
| |
| namespace LIBC_NAMESPACE_DECL { |
| |
| namespace math { |
| |
| namespace lgammaf16_internal { |
| |
| LIBC_INLINE constexpr bool is_integer(float16 x) { |
| using FPBits = fputil::FPBits<float16>; |
| FPBits xbits(x); |
| uint16_t x_u = xbits.uintval(); |
| unsigned x_e = static_cast<unsigned>(xbits.get_biased_exponent()); |
| unsigned lsb = static_cast<unsigned>( |
| cpp::countr_zero(static_cast<uint16_t>(x_u | FPBits::EXP_MASK))); |
| constexpr unsigned UNIT_EXPONENT = |
| static_cast<unsigned>(FPBits::EXP_BIAS + FPBits::FRACTION_LEN); |
| return x_e + lsb >= UNIT_EXPONENT; |
| } |
| |
| // sin(pi*x) for x in (0, 1), written as (0.25 - u^2) * P(u^2) where u = x-0.5. |
| // Coefficients for P are a degree-7 polynomial in u^2 |
| // approximating sin(pi*(u+0.5)) / (pi*(0.25-u^2)) |
| // with max error ~2^{-55}, generated by Sollya with: |
| // > P = fpminimax(sin(pi*(x+0.5))/(pi*(0.25-x^2)), |
| // [|0,2,4,6,8,10,12,14|], [|D...|], [-0.5, 0.5]); |
| LIBC_INLINE constexpr double lg_sinpi(double x) { |
| constexpr double COEFFS[8] = {0x000000000000001p+2, -0x1.de9e64df22ea4p+1, |
| 0x1.472be122401f8p+0, -0x1.d4fcd82df91bp-3, |
| 0x1.9f05c97e0aab2p-6, -0x1.f3091c427b611p-10, |
| 0x1.b22c9bfdca547p-14, -0x1.15484325ef569p-18}; |
| double u = x - 0.5; |
| double u2 = u * u, u4 = u2 * u2, u8 = u4 * u4; |
| double p01 = fputil::multiply_add(u2, COEFFS[1], COEFFS[0]); |
| double p23 = fputil::multiply_add(u2, COEFFS[3], COEFFS[2]); |
| double p45 = fputil::multiply_add(u2, COEFFS[5], COEFFS[4]); |
| double p67 = fputil::multiply_add(u2, COEFFS[7], COEFFS[6]); |
| double p03 = fputil::multiply_add(u4, p23, p01); |
| double p47 = fputil::multiply_add(u4, p67, p45); |
| return (0.25 - u2) * fputil::multiply_add(u8, p47, p03); |
| } |
| |
| // Natural logarithm of x (x > 0), using 16-entry table + degree-7 polynomial. |
| // Range reduction: x = 2^e * m where m in [1, 2), decomposed as |
| // m = (1+i/16)*(1+z) so log(x) = e*log(2) + log(1+i/16) + log(1+z) |
| // = e*log(2) + IL[i] + z*P(z). |
| // P approximates log(1+z)/z on z in [-1/16, 1/16] with max |
| // error ~2^{-54}, generated by Sollya with: |
| // > P = fpminimax(log(1+x)/x, [|0,1,2,3,4,5,6,7|], [|D...|], [-1/16, 1/16]); |
| LIBC_INLINE double lg_ln(double x) { |
| using FPBits = fputil::FPBits<double>; |
| uint64_t u = FPBits(x).uintval(); |
| int e = static_cast<int>(FPBits(x).get_biased_exponent()) - 0x3ff; |
| |
| // Coefficients for log(1 + z)/z on z in [-1/16, 1/16] |
| constexpr double COEFFS[8] = {0x1.fffffffffff24p-1, -0x1.ffffffffd1d67p-2, |
| 0x1.55555537802dep-2, -0x1.ffffeca81b866p-3, |
| 0x1.999611761d772p-3, -0x1.54f3e581b61bfp-3, |
| 0x1.1e642b4cb5143p-3, -0x1.9115a5af1e1edp-4}; |
| // IL[i] = log(1 + i/16) for i = 0..15 |
| constexpr double IL[16] = { |
| 0x1.59caeec280116p-57, 0x1.f0a30c01162aap-5, 0x1.e27076e2af2ebp-4, |
| 0x1.5ff3070a793d6p-3, 0x1.c8ff7c79a9a2p-3, 0x1.1675cababa60fp-2, |
| 0x1.4618bc21c5ec2p-2, 0x1.739d7f6bbd007p-2, 0x1.9f323ecbf984dp-2, |
| 0x1.c8ff7c79a9a21p-2, 0x1.f128f5faf06ecp-2, 0x1.0be72e4252a83p-1, |
| 0x1.1e85f5e7040d1p-1, 0x1.307d7334f10bep-1, 0x1.41d8fe84672afp-1, |
| 0x1.52a2d265bc5abp-1}; |
| // IX[i] = 1 / (1 + i/16) for i = 0..15 |
| constexpr double IX[16] = { |
| 0x000000000000001p+0, 0x1.e1e1e1e1e1e1ep-1, 0x1.c71c71c71c71cp-1, |
| 0x1.af286bca1af28p-1, 0x1.999999999999ap-1, 0x1.8618618618618p-1, |
| 0x1.745d1745d1746p-1, 0x1.642c8590b2164p-1, 0x1.5555555555555p-1, |
| 0x1.47ae147ae147bp-1, 0x1.3b13b13b13b14p-1, 0x1.2f684bda12f68p-1, |
| 0x1.2492492492492p-1, 0x1.1a7b9611a7b96p-1, 0x1.1111111111111p-1, |
| 0x1.0842108421084p-1}; |
| |
| int i = static_cast<int>((u >> 48) & 0xf); |
| // Reduce to mantissa in [1, 2) |
| uint64_t mant_u = (u & (~uint64_t(0) >> 12)) | (uint64_t(0x3ff) << 52); |
| double mant = FPBits(mant_u).get_val(); |
| double z = IX[i] * mant - 1.0, z2 = z * z, z4 = z2 * z2; |
| double q01 = fputil::multiply_add(z, COEFFS[1], COEFFS[0]); |
| double q23 = fputil::multiply_add(z, COEFFS[3], COEFFS[2]); |
| double q45 = fputil::multiply_add(z, COEFFS[5], COEFFS[4]); |
| double q67 = fputil::multiply_add(z, COEFFS[7], COEFFS[6]); |
| double q03 = fputil::multiply_add(z2, q23, q01); |
| double q47 = fputil::multiply_add(z2, q67, q45); |
| return e * 0x1.62e42fefa39efp-1 + IL[i] + |
| z * fputil::multiply_add(z4, q47, q03); |
| } |
| |
| } // namespace lgammaf16_internal |
| |
| LIBC_INLINE float16 lgammaf16(float16 x) { |
| using namespace lgammaf16_internal; |
| using FPBits = fputil::FPBits<float16>; |
| |
| FPBits xbits(x); |
| uint16_t x_abs = xbits.abs().uintval(); |
| |
| // Handle NaN and Inf |
| if (LIBC_UNLIKELY(x_abs >= 0x7c00u)) { |
| if (x_abs == 0x7c00u) |
| return FPBits::inf().get_val(); |
| if (xbits.is_signaling_nan()) { |
| fputil::raise_except_if_required(FE_INVALID); |
| return FPBits::quiet_nan().get_val(); |
| } |
| return x; |
| } |
| |
| // +-0 pole error |
| if (LIBC_UNLIKELY(x_abs == 0)) { |
| fputil::raise_except_if_required(FE_DIVBYZERO); |
| fputil::set_errno_if_required(ERANGE); |
| return FPBits::inf().get_val(); |
| } |
| |
| if (LIBC_UNLIKELY(is_integer(x))) { |
| if (xbits.is_neg()) { |
| // pole -> +Inf |
| fputil::raise_except_if_required(FE_DIVBYZERO); |
| fputil::set_errno_if_required(ERANGE); |
| return FPBits::inf().get_val(); |
| } |
| // lgamma(1) = lgamma(2) = 0 |
| if (x_abs == 0x3c00u || x_abs == 0x4000u) |
| return FPBits::zero().get_val(); |
| } |
| |
| double xd = fputil::cast<double>(x); |
| double abs_xd = xd < 0.0 ? -xd : xd; |
| double lgamma_val; |
| |
| // |x| = 0.66015625 |
| if (LIBC_UNLIKELY(x_abs < 0x3948u)) { |
| // Small path: |x| < 0.66015625 |
| // For t = |x|: lgamma(t) = t*P(u) - log(t), where P(u) approximates |
| // h(t) = (lgamma(t) + log(t)) / t, u = (t - MID_S) / HW_S in [-1, 1] |
| // For x < 0: lgamma(-t) = log(pi) - log(sin(pi*t)) - lgamma(t) - log(t) |
| // = log(pi) - lg_ln(lg_sinpi(t)) - t*P(u) |
| // Degree-12 Chebyshev approximation for h in u; max lgamma error ~1.5e-13 |
| // Coefficients generated by DCT at 13 Chebyshev nodes on |
| // [1e-5, 0.66015625]. |
| // Evaluated via Clenshaw: sum_{k=0}^{12} CHEB_H[k] * T_k(u) |
| constexpr double MID_S = 0x1.52014f8b588e3p-2; |
| constexpr double HW_S = 0x1.51feb074a771dp-2; |
| constexpr double CHEB_H[13] = { |
| -0x1.6ab0351d4e02fp-2, 0x1.ac55f136e49cep-3, -0x1.9fc5280cd41bfp-7, |
| 0x1.1825d37c81c1ep-10, -0x1.ae2f107c31c9ep-14, 0x1.6174bdf95e9abp-17, |
| -0x1.2e6abd22bc4ecp-20, 0x1.09b66d4f14b06p-23, -0x1.dbb938f189d8ap-27, |
| 0x1.afcb748e20e4ep-30, -0x1.8c221a9d89d8ap-33, 0x1.6e708e49b3aeep-36, |
| -0x1.50accec4ec4ecp-39}; |
| // Clenshaw's algorithm: p = sum_k CHEB_H[k] * T_k(u_s) |
| double u_s = (abs_xd - MID_S) / HW_S, u_s2 = 2.0 * u_s; |
| double bh1 = CHEB_H[12], bh2 = 0.0; |
| for (int k = 11; k >= 1; k--) { |
| double bh0 = fputil::multiply_add(u_s2, bh1, CHEB_H[k] - bh2); |
| bh2 = bh1; |
| bh1 = bh0; |
| } |
| double poly_h = fputil::multiply_add(u_s, bh1, CHEB_H[0] - bh2); |
| // abs_xd * poly_h = lgamma(abs_xd) + log(abs_xd) |
| double poly_val = abs_xd * poly_h; |
| if (xbits.is_neg()) |
| lgamma_val = 0x1.250d048e7a1bdp+0 - lg_ln(lg_sinpi(abs_xd)) - poly_val; |
| else |
| lgamma_val = poly_val - lg_ln(abs_xd); |
| } else if (x_abs > 0x42BEu) { // 3.3710938 |
| double log_abs_xd = lg_ln(abs_xd); |
| // lgamma(x) = (x-0.5)*log(x) - x + log(2*pi)/2 + Bernoulli corrections |
| // Stirling base: (x-0.5)*(log(x)-1) + (log(2*pi)-1)/2 |
| lgamma_val = (abs_xd - 0.5) * (log_abs_xd - 1.0) + 0x1.acfe390c97d69p-2; |
| |
| double inv_x = 1.0 / abs_xd, inv_x2 = inv_x * inv_x; |
| |
| if (x_abs > 0x64AEu) { |
| // |x| > 1198 1-term Stirling-Bernoulli: 1/(12*x) |
| lgamma_val += inv_x * 0x1.5555555555555p-4; |
| } else if (x_abs > 0x549Eu) { |
| // |x| in (73.875, 1198] 2-term correction |
| // Minimax polynomial in 1/x^2 for R(x) = lgamma(x) - stirling_base(x), |
| // where stirling_base(x) = (x-0.5)*log(x) - x + log(2*pi)/2 |
| // Exact Bernoulli: B_2/(1*2) = 1/12, B_4/(3*4) = -1/360 |
| // Generated by Sollya with: |
| // > stir = proc(x) { return lgamma(x)-((x-0.5)*log(x)-x+log(2*pi)/2); |
| // }; > P = fpminimax(stir(x), [|1,3|], [|D...|], [73.875, 1198]); |
| constexpr double BERN2[2] = {0x1.555555547fbadp-4, -0x1.6c0fd270c465p-9}; |
| lgamma_val += inv_x * fputil::polyeval(inv_x2, BERN2[0], BERN2[1]); |
| } else if (x_abs > 0x4955u) { |
| // |x| in (10.664, 73.875] 4-term correction |
| // Minimax polynomial in 1/x^2 for R(x); exact Bernoulli terms are |
| // 1/12, -1/360, 1/1260, -1/1680 |
| // Generated by Sollya with: |
| // > display = hexadecimal; |
| // > stir = proc(x) { return lgamma(x)-((x-0.5)*log(x)-x+log(2*pi)/2); |
| // }; > P = fpminimax(stir(x), [|1,3,5,7|], [|D...|], [10.664, 73.875]); |
| constexpr double BERN4[4] = {0x1.555555554de0bp-4, -0x1.6c16bdc45944fp-9, |
| 0x1.a0077f300ecb3p-11, |
| -0x1.2e9cfff3b29c2p-11}; |
| double inv_x4 = inv_x2 * inv_x2; |
| lgamma_val += |
| inv_x * fputil::multiply_add( |
| inv_x4, fputil::multiply_add(inv_x2, BERN4[3], BERN4[2]), |
| fputil::multiply_add(inv_x2, BERN4[1], BERN4[0])); |
| } else { |
| // |x| in (3.373, 10.664]: 8-term correction. |
| // Degree-7 minimax polynomial in 1/x^2 for R(x), generated by Sollya |
| // with: |
| // > stir = proc(x) { return lgamma(x)-((x-0.5)*log(x)-x+log(2*pi)/2); |
| // }; > P = fpminimax(stir(x), [|1,3,5,7,9,11,13,15|], [|D...|], |
| // [3.373, 10.664]); |
| constexpr double BERN8[8] = { |
| 0x1.5555555551286p-4, -0x1.6c16c0e7c4cf4p-9, 0x1.a0193267fe6f2p-11, |
| -0x1.37e87ec19cb45p-11, 0x1.b40011dfff081p-11, -0x1.c16c8946b19b6p-10, |
| 0x1.e9f47ace150d8p-9, -0x1.4f5843a71a338p-8}; |
| double inv_x4 = inv_x2 * inv_x2, inv_x8 = inv_x4 * inv_x4; |
| double p01 = fputil::multiply_add(inv_x2, BERN8[1], BERN8[0]); |
| double p23 = fputil::multiply_add(inv_x2, BERN8[3], BERN8[2]); |
| double p45 = fputil::multiply_add(inv_x2, BERN8[5], BERN8[4]); |
| double p67 = fputil::multiply_add(inv_x2, BERN8[7], BERN8[6]); |
| double p03 = fputil::multiply_add(inv_x4, p23, p01); |
| double p47 = fputil::multiply_add(inv_x4, p67, p45); |
| lgamma_val += inv_x * fputil::multiply_add(inv_x8, p47, p03); |
| } |
| |
| if (xbits.is_neg()) { |
| // Reflection formula |
| // lgamma(x) = log(pi) - log|x| - log|sin(pi*frac_x)| - lgamma|x| |
| double frac_x = xd - fputil::floor(xd); |
| lgamma_val = 0x1.250d048e7a1bdp+0 - lgamma_val - log_abs_xd; |
| lgamma_val -= lg_ln(lg_sinpi(frac_x)); |
| } |
| } else { |
| // Medium path: |x| in [0.66015625, 3.373046875] |
| // lgamma(|x|) = (|x|-1)*(|x|-2) * P(u), where P(u) approximates |
| // g(t) = lgamma(t) / ((t-1)*(t-2)), u = (t - MID_M) / HW_M in [-1, 1] |
| // Degree-20 Chebyshev approximation for g in u; max lgamma error ~1.4e-10 |
| // Coefficients generated by DCT at 21 Chebyshev nodes on |
| // [0.66015625, 3.373046875]. Evaluated via Clenshaw: sum_{k=0}^{20} |
| // CHEB_G[k] * T_k(u) |
| constexpr double MID_M = 0x1.0220000000000p+1; |
| constexpr double HW_M = 0x1.5b40000000000p+0; |
| constexpr double CHEB_G[21] = { |
| 0x1.d690d4de8b355p-2, -0x1.53439726f7d8bp-3, 0x1.568fc6653c026p-5, |
| -0x1.8d074f73f36efp-7, 0x1.eefebf4ca8e86p-9, -0x1.4275913d85ff5p-10, |
| 0x1.b075f0d96e492p-12, -0x1.27ed35f357339p-13, 0x1.9b12cbe71a800p-15, |
| -0x1.20c949c2829b5p-16, 0x1.996eb95241555p-18, -0x1.2460d6c60f800p-19, |
| 0x1.a42297aae1862p-21, -0x1.2f6c7d28230a5p-22, 0x1.b8394aac30c31p-24, |
| -0x1.409726ce1110ep-25, 0x1.d484ac42db6dbp-27, -0x1.574dae55861abp-28, |
| 0x1.f782605555555p-30, -0x1.6cefb8f65aa2ep-31, 0x1.d9f268f3cf3cfp-33}; |
| // Clenshaw's algorithm: poly_g = sum_k CHEB_G[k] * T_k(u_m) |
| double u_m = (abs_xd - MID_M) / HW_M, u_m2 = 2.0 * u_m; |
| double bg1 = CHEB_G[20], bg2 = 0.0; |
| for (int k = 19; k >= 1; k--) { |
| double bg0 = fputil::multiply_add(u_m2, bg1, CHEB_G[k] - bg2); |
| bg2 = bg1; |
| bg1 = bg0; |
| } |
| double poly_g = fputil::multiply_add(u_m, bg1, CHEB_G[0] - bg2); |
| lgamma_val = (abs_xd - 1.0) * (abs_xd - 2.0) * poly_g; |
| |
| if (xbits.is_neg()) { |
| // Reflection |
| // (x) = log(pi) - lgamma|x| - log(|x|*|sin(pi*frac_x)|) |
| double frac_x = xd - fputil::floor(xd); |
| lgamma_val = 0x1.250d048e7a1bdp+0 - lgamma_val; |
| lgamma_val -= lg_ln(lg_sinpi(frac_x) * abs_xd); |
| } |
| } |
| |
| float16 result = fputil::cast<float16>(lgamma_val); |
| if (LIBC_UNLIKELY(fputil::FPBits<float16>(result).is_inf())) |
| fputil::set_errno_if_required(ERANGE); |
| return result; |
| } |
| |
| } // namespace math |
| |
| } // namespace LIBC_NAMESPACE_DECL |
| |
| #endif // LIBC_TYPES_HAS_FLOAT16 |
| |
| #endif // LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMAF16_H |
