1//===-- Implementation header for lgammaf16 ---------------------*- C++ -*-===//
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-exception
6//
7//===----------------------------------------------------------------------===//
8
9#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMAF16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMAF16_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
16#include "src/__support/FPUtil/FEnvImpl.h"
17#include "src/__support/FPUtil/FPBits.h"
18#include "src/__support/FPUtil/NearestIntegerOperations.h"
19#include "src/__support/FPUtil/PolyEval.h"
20#include "src/__support/FPUtil/cast.h"
21#include "src/__support/FPUtil/multiply_add.h"
22#include "src/__support/macros/config.h"
23#include "src/__support/macros/optimization.h"
24#include "src/__support/macros/properties/types.h"
25#include "src/__support/math/gamma_util.h"
26
27namespace LIBC_NAMESPACE_DECL {
28
29namespace math {
30
31LIBC_INLINE float16 lgammaf16(float16 x) {
32 using namespace gamma_internal;
33 using FPBits = fputil::FPBits<float16>;
34
35 FPBits xbits(x);
36 uint16_t x_abs = xbits.abs().uintval();
37
38 // Handle NaN and Inf
39 if (LIBC_UNLIKELY(x_abs >= 0x7c00u)) {
40 if (x_abs == 0x7c00u)
41 return FPBits::inf().get_val();
42 if (xbits.is_signaling_nan()) {
43 fputil::raise_except_if_required(FE_INVALID);
44 return FPBits::quiet_nan().get_val();
45 }
46 return x;
47 }
48
49 // +-0 pole error
50 if (LIBC_UNLIKELY(x_abs == 0)) {
51 fputil::raise_except_if_required(FE_DIVBYZERO);
52 fputil::set_errno_if_required(ERANGE);
53 return FPBits::inf().get_val();
54 }
55
56 if (LIBC_UNLIKELY(is_integer(x))) {
57 if (xbits.is_neg()) {
58 // pole -> +Inf
59 fputil::raise_except_if_required(FE_DIVBYZERO);
60 fputil::set_errno_if_required(ERANGE);
61 return FPBits::inf().get_val();
62 }
63 // lgamma(1) = lgamma(2) = 0
64 if (x_abs == 0x3c00u || x_abs == 0x4000u)
65 return FPBits::zero().get_val();
66 }
67
68 double xd = fputil::cast<double>(x);
69 double abs_xd = xd < 0.0 ? -xd : xd;
70 double lgamma_val;
71
72 if (LIBC_UNLIKELY(x_abs < 0x3948u)) {
73 // Small path: |x| < 0.66015625
74 // For t = |x|: lgamma(t) = t*P(u) - log(t), where P(u) approximates
75 // h(t) = (lgamma(t) + log(t)) / t, u = (t - MID_S) / HW_S in [-1, 1]
76 // For x < 0: lgamma(-t) = log(pi) - log(sin(pi*t)) - lgamma(t) - log(t)
77 // = log(pi) - lg_ln(lg_sinpi(t)) - t*P(u)
78 // Degree-12 Chebyshev approximation for h in u; max lgamma error ~1.5e-13
79 // Coefficients generated by DCT at 13 Chebyshev nodes on
80 // [1e-5, 0.66015625].
81 // Evaluated via Clenshaw: sum_{k=0}^{12} CHEB_H[k] * T_k(u)
82 constexpr double MID_S = 0x1.52014f8b588e3p-2;
83 constexpr double HW_S = 0x1.51feb074a771dp-2;
84 constexpr double CHEB_H[13] = {
85 -0x1.6ab0351d4e02fp-2, 0x1.ac55f136e49cep-3, -0x1.9fc5280cd41bfp-7,
86 0x1.1825d37c81c1ep-10, -0x1.ae2f107c31c9ep-14, 0x1.6174bdf95e9abp-17,
87 -0x1.2e6abd22bc4ecp-20, 0x1.09b66d4f14b06p-23, -0x1.dbb938f189d8ap-27,
88 0x1.afcb748e20e4ep-30, -0x1.8c221a9d89d8ap-33, 0x1.6e708e49b3aeep-36,
89 -0x1.50accec4ec4ecp-39};
90 // Clenshaw's algorithm: p = sum_k CHEB_H[k] * T_k(u_s)
91 double u_s = (abs_xd - MID_S) / HW_S, u_s2 = 2.0 * u_s;
92 double bh1 = CHEB_H[12], bh2 = 0.0;
93 for (int k = 11; k >= 1; k--) {
94 double bh0 = fputil::multiply_add(x: u_s2, y: bh1, z: CHEB_H[k] - bh2);
95 bh2 = bh1;
96 bh1 = bh0;
97 }
98 double poly_h = fputil::multiply_add(x: u_s, y: bh1, z: CHEB_H[0] - bh2);
99 // abs_xd * poly_h = lgamma(abs_xd) + log(abs_xd)
100 double poly_val = abs_xd * poly_h;
101 if (xbits.is_neg())
102 lgamma_val = 0x1.250d048e7a1bdp+0 - lg_ln(x: lg_sinpi(x: abs_xd)) - poly_val;
103 else
104 lgamma_val = poly_val - lg_ln(x: abs_xd);
105 } else if (x_abs > 0x42BEu) { // 3.3710938
106 double log_abs_xd = lg_ln(x: abs_xd);
107 // lgamma(x) = (x-0.5)*log(x) - x + log(2*pi)/2 + Bernoulli corrections
108 // Stirling base: (x-0.5)*(log(x)-1) + (log(2*pi)-1)/2
109 lgamma_val = (abs_xd - 0.5) * (log_abs_xd - 1.0) + 0x1.acfe390c97d69p-2;
110
111 double inv_x = 1.0 / abs_xd, inv_x2 = inv_x * inv_x;
112
113 if (x_abs > 0x64AEu) {
114 // |x| > 1198 1-term Stirling-Bernoulli: 1/(12*x)
115 lgamma_val += inv_x * 0x1.5555555555555p-4;
116 } else if (x_abs > 0x549Eu) {
117 // |x| in (73.875, 1198] 2-term correction
118 // Minimax polynomial in 1/x^2 for R(x) = lgamma(x) - stirling_base(x),
119 // where stirling_base(x) = (x-0.5)*log(x) - x + log(2*pi)/2
120 // Exact Bernoulli: B_2/(1*2) = 1/12, B_4/(3*4) = -1/360
121 // Generated by Sollya with:
122 // > stir = proc(x) { return lgamma(x)-((x-0.5)*log(x)-x+log(2*pi)/2);
123 // }; > P = fpminimax(stir(x), [|1,3|], [|D...|], [73.875, 1198]);
124 constexpr double BERN2[2] = {0x1.555555547fbadp-4, -0x1.6c0fd270c465p-9};
125 lgamma_val += inv_x * fputil::polyeval(x: inv_x2, a0: BERN2[0], a: BERN2[1]);
126 } else if (x_abs > 0x4955u) {
127 // |x| in (10.664, 73.875] 4-term correction
128 // Minimax polynomial in 1/x^2 for R(x); exact Bernoulli terms are
129 // 1/12, -1/360, 1/1260, -1/1680
130 // Generated by Sollya with:
131 // > display = hexadecimal;
132 // > stir = proc(x) { return lgamma(x)-((x-0.5)*log(x)-x+log(2*pi)/2);
133 // }; > P = fpminimax(stir(x), [|1,3,5,7|], [|D...|], [10.664, 73.875]);
134 constexpr double BERN4[4] = {0x1.555555554de0bp-4, -0x1.6c16bdc45944fp-9,
135 0x1.a0077f300ecb3p-11,
136 -0x1.2e9cfff3b29c2p-11};
137 double inv_x4 = inv_x2 * inv_x2;
138 lgamma_val +=
139 inv_x * fputil::multiply_add(
140 x: inv_x4, y: fputil::multiply_add(x: inv_x2, y: BERN4[3], z: BERN4[2]),
141 z: fputil::multiply_add(x: inv_x2, y: BERN4[1], z: BERN4[0]));
142 } else {
143 // |x| in (3.373, 10.664]: 8-term correction.
144 // Degree-7 minimax polynomial in 1/x^2 for R(x), generated by Sollya
145 // with:
146 // > stir = proc(x) { return lgamma(x)-((x-0.5)*log(x)-x+log(2*pi)/2);
147 // }; > P = fpminimax(stir(x), [|1,3,5,7,9,11,13,15|], [|D...|],
148 // [3.373, 10.664]);
149 constexpr double BERN8[8] = {
150 0x1.5555555551286p-4, -0x1.6c16c0e7c4cf4p-9, 0x1.a0193267fe6f2p-11,
151 -0x1.37e87ec19cb45p-11, 0x1.b40011dfff081p-11, -0x1.c16c8946b19b6p-10,
152 0x1.e9f47ace150d8p-9, -0x1.4f5843a71a338p-8};
153 double inv_x4 = inv_x2 * inv_x2, inv_x8 = inv_x4 * inv_x4;
154 double p01 = fputil::multiply_add(x: inv_x2, y: BERN8[1], z: BERN8[0]);
155 double p23 = fputil::multiply_add(x: inv_x2, y: BERN8[3], z: BERN8[2]);
156 double p45 = fputil::multiply_add(x: inv_x2, y: BERN8[5], z: BERN8[4]);
157 double p67 = fputil::multiply_add(x: inv_x2, y: BERN8[7], z: BERN8[6]);
158 double p03 = fputil::multiply_add(x: inv_x4, y: p23, z: p01);
159 double p47 = fputil::multiply_add(x: inv_x4, y: p67, z: p45);
160 lgamma_val += inv_x * fputil::multiply_add(x: inv_x8, y: p47, z: p03);
161 }
162
163 if (xbits.is_neg()) {
164 // Reflection formula
165 // lgamma(x) = log(pi) - log|x| - log|sin(pi*frac_x)| - lgamma|x|
166 double frac_x = xd - fputil::floor(x: xd);
167 lgamma_val = 0x1.250d048e7a1bdp+0 - lgamma_val - log_abs_xd;
168 lgamma_val -= lg_ln(x: lg_sinpi(x: frac_x));
169 }
170 } else {
171 // Medium path: |x| in [0.66015625, 3.373046875]
172 // lgamma(|x|) = (|x|-1)*(|x|-2) * P(u), where P(u) approximates
173 // g(t) = lgamma(t) / ((t-1)*(t-2)), u = (t - MID_M) / HW_M in [-1, 1]
174 // Degree-20 Chebyshev approximation for g in u; max lgamma error ~1.4e-10
175 // Coefficients generated by DCT at 21 Chebyshev nodes on
176 // [0.66015625, 3.373046875]. Evaluated via Clenshaw: sum_{k=0}^{20}
177 // CHEB_G[k] * T_k(u)
178 constexpr double MID_M = 0x1.0220000000000p+1;
179 constexpr double HW_M = 0x1.5b40000000000p+0;
180 constexpr double CHEB_G[21] = {
181 0x1.d690d4de8b355p-2, -0x1.53439726f7d8bp-3, 0x1.568fc6653c026p-5,
182 -0x1.8d074f73f36efp-7, 0x1.eefebf4ca8e86p-9, -0x1.4275913d85ff5p-10,
183 0x1.b075f0d96e492p-12, -0x1.27ed35f357339p-13, 0x1.9b12cbe71a800p-15,
184 -0x1.20c949c2829b5p-16, 0x1.996eb95241555p-18, -0x1.2460d6c60f800p-19,
185 0x1.a42297aae1862p-21, -0x1.2f6c7d28230a5p-22, 0x1.b8394aac30c31p-24,
186 -0x1.409726ce1110ep-25, 0x1.d484ac42db6dbp-27, -0x1.574dae55861abp-28,
187 0x1.f782605555555p-30, -0x1.6cefb8f65aa2ep-31, 0x1.d9f268f3cf3cfp-33};
188 // Clenshaw's algorithm: poly_g = sum_k CHEB_G[k] * T_k(u_m)
189 double u_m = (abs_xd - MID_M) / HW_M, u_m2 = 2.0 * u_m;
190 double bg1 = CHEB_G[20], bg2 = 0.0;
191 for (int k = 19; k >= 1; k--) {
192 double bg0 = fputil::multiply_add(x: u_m2, y: bg1, z: CHEB_G[k] - bg2);
193 bg2 = bg1;
194 bg1 = bg0;
195 }
196 double poly_g = fputil::multiply_add(x: u_m, y: bg1, z: CHEB_G[0] - bg2);
197 lgamma_val = (abs_xd - 1.0) * (abs_xd - 2.0) * poly_g;
198
199 if (xbits.is_neg()) {
200 // Reflection
201 // (x) = log(pi) - lgamma|x| - log(|x|*|sin(pi*frac_x)|)
202 double frac_x = xd - fputil::floor(x: xd);
203 lgamma_val = 0x1.250d048e7a1bdp+0 - lgamma_val;
204 lgamma_val -= lg_ln(x: lg_sinpi(x: frac_x) * abs_xd);
205 }
206 }
207
208 float16 result = fputil::cast<float16>(x: lgamma_val);
209 if (LIBC_UNLIKELY(fputil::FPBits<float16>(result).is_inf()))
210 fputil::set_errno_if_required(ERANGE);
211 return result;
212}
213
214} // namespace math
215
216} // namespace LIBC_NAMESPACE_DECL
217
218#endif // LIBC_TYPES_HAS_FLOAT16
219
220#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMAF16_H
221