1//===-- Implementation header for lgammaf -----------------------*- 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_LGAMMAF_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMAF_H
11
12#include "src/__support/FPUtil/FEnvImpl.h"
13#include "src/__support/FPUtil/FPBits.h"
14#include "src/__support/FPUtil/NearestIntegerOperations.h"
15#include "src/__support/FPUtil/PolyEval.h"
16#include "src/__support/FPUtil/cast.h"
17#include "src/__support/FPUtil/except_value_utils.h"
18#include "src/__support/FPUtil/multiply_add.h"
19#include "src/__support/macros/config.h"
20#include "src/__support/macros/optimization.h"
21#include "src/__support/math/gamma_util.h"
22#include "src/__support/math/log.h"
23
24namespace LIBC_NAMESPACE_DECL {
25
26namespace math {
27
28namespace lgammaf_internal {
29
30// P_M2(d), d = t - 1.5, approximates lgamma(t)/((t-1)(t-2)) on [1, 2].
31// Degree 18 centered-monomial coefficients from a 250-bit Chebyshev fit,
32// max relative error 2^-51.6
33LIBC_INLINE double lgamma_m2_poly(double d) {
34 constexpr double POLY_M2[19] = {
35 0x1.eeb95b094c191p-2, -0x1.2aed059bd613dp-3, 0x1.01af62a292eb9p-4,
36 -0x1.007aa83cc458bp-5, 0x1.13342351db068p-6, -0x1.34dbda42244dep-7,
37 0x1.64f066cf96646p-8, -0x1.a5033bab00972p-9, 0x1.f814855fdfeb8p-10,
38 -0x1.31418be9f2e14p-10, 0x1.7510b1d1fa536p-11, -0x1.caccf4395927bp-12,
39 0x1.1c0d19e4fef05p-12, -0x1.67d43c4be68f7p-13, 0x1.c3d6281f9ac35p-14,
40 -0x1.dd9e09255e712p-15, 0x1.2550c8f33e895p-15, -0x1.6cd0f980284cbp-15,
41 0x1.dbd1c21b8dfa7p-16};
42 double d2 = d * d, d4 = d2 * d2, d8 = d4 * d4, d16 = d8 * d8;
43 double p01 = fputil::multiply_add(x: d, y: POLY_M2[1], z: POLY_M2[0]);
44 double p23 = fputil::multiply_add(x: d, y: POLY_M2[3], z: POLY_M2[2]);
45 double p45 = fputil::multiply_add(x: d, y: POLY_M2[5], z: POLY_M2[4]);
46 double p67 = fputil::multiply_add(x: d, y: POLY_M2[7], z: POLY_M2[6]);
47 double p89 = fputil::multiply_add(x: d, y: POLY_M2[9], z: POLY_M2[8]);
48 double p1011 = fputil::multiply_add(x: d, y: POLY_M2[11], z: POLY_M2[10]);
49 double p1213 = fputil::multiply_add(x: d, y: POLY_M2[13], z: POLY_M2[12]);
50 double p1415 = fputil::multiply_add(x: d, y: POLY_M2[15], z: POLY_M2[14]);
51 double p1617 = fputil::multiply_add(x: d, y: POLY_M2[17], z: POLY_M2[16]);
52 double q03 = fputil::multiply_add(x: d2, y: p23, z: p01);
53 double q47 = fputil::multiply_add(x: d2, y: p67, z: p45);
54 double q811 = fputil::multiply_add(x: d2, y: p1011, z: p89);
55 double q1215 = fputil::multiply_add(x: d2, y: p1415, z: p1213);
56 double q1618 = fputil::multiply_add(x: d2, y: POLY_M2[18], z: p1617);
57 double r07 = fputil::multiply_add(x: d4, y: q47, z: q03);
58 double r815 = fputil::multiply_add(x: d4, y: q1215, z: q811);
59 double s015 = fputil::multiply_add(x: d8, y: r815, z: r07);
60 return fputil::multiply_add(x: d16, y: q1618, z: s015);
61}
62
63} // namespace lgammaf_internal
64
65LIBC_INLINE float lgammaf(float x) {
66 using namespace gamma_internal;
67 using namespace lgammaf_internal;
68 using FPBits = fputil::FPBits<float>;
69
70 FPBits xbits(x);
71 uint32_t x_abs = xbits.abs().uintval();
72
73 // NaN / Inf
74 if (LIBC_UNLIKELY(x_abs >= 0x7f800000u)) {
75 if (x_abs == 0x7f800000u)
76 return FPBits::inf().get_val();
77 if (xbits.is_signaling_nan()) {
78 fputil::raise_except_if_required(FE_INVALID);
79 return FPBits::quiet_nan().get_val();
80 }
81 return x;
82 }
83
84 // +/- 0 -> +Inf pole
85 if (LIBC_UNLIKELY(x_abs == 0)) {
86 fputil::raise_except_if_required(FE_DIVBYZERO);
87 fputil::set_errno_if_required(ERANGE);
88 return FPBits::inf().get_val();
89 }
90
91 // Negative integers and lgamma(1) = lgamma(2) = 0.
92 if (LIBC_UNLIKELY(is_integer(x))) {
93 if (xbits.is_neg()) {
94 fputil::raise_except_if_required(FE_DIVBYZERO);
95 fputil::set_errno_if_required(ERANGE);
96 return FPBits::inf().get_val();
97 }
98 if (x_abs == 0x3f800000u || x_abs == 0x40000000u)
99 return FPBits::zero().get_val();
100 }
101
102 double xd = fputil::cast<double>(x);
103 double abs_xd = xd < 0.0 ? -xd : xd;
104 double lgamma_val;
105
106 // For very tiny |x| (< 2^-23), use the truncated Laurent series:
107 // lgamma(x) = -log(x) - gamma*x + O(x^2) for tiny x > 0
108 // lgamma(-y) = -log(y) + gamma*y + O(y^2) for tiny y > 0
109 // Even though gamma*|x| < 2^-25 is below 1 float ULP of |result|, it tips
110 // rounding at boundary cases. The x^2 term is at most gamma*2^-46 << 2^-50.
111 if (x_abs < 0x34000000u) { // |x| < 2^-23
112 constexpr fputil::ExceptValues<float, 4> LGAMMAF_EXCEPTS_TINY{.values: {
113 // input, toward-zero result, RU, RD, RN
114 {.input: 0x9b7679ffu, .rnd_towardzero_result: 0x4247c72cu, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
115 {.input: 0x9e88452du, .rnd_towardzero_result: 0x4236bd8bu, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
116 {.input: 0xa77a8e47u, .rnd_towardzero_result: 0x42052b94u, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
117 {.input: 0xb0e17820u, .rnd_towardzero_result: 0x41a1d37bu, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
118 }};
119 if (auto r = LGAMMAF_EXCEPTS_TINY.lookup(x_bits: xbits.uintval());
120 LIBC_UNLIKELY(r.has_value()))
121 return r.value();
122 constexpr double EULER_GAMMA = 0x1.2788cfc6fb619p-1;
123 double sign_corr = xbits.is_neg() ? EULER_GAMMA : -EULER_GAMMA;
124 return fputil::cast<float>(
125 x: fputil::multiply_add(x: sign_corr, y: abs_xd, z: -math::log(x: abs_xd)));
126 }
127
128 if (x_abs < 0x3f290000u) {
129 if (xbits.is_neg()) {
130 // Small negative: x in (-0.66015625, -2^-23). Degree-25 monomial fit
131 // of g(x) = (lgamma(x) + log(-x)) / x (smooth on [-0.66, 0]), then
132 // lgamma(x) = x*g(x) - log(-x). d = x - MID is exact.
133 constexpr double MID_SN = -0x1.5200000000000p-2;
134 constexpr double POLY_SN[26] = {
135 -0x1.cf99908bbb1d7p-1, 0x1.386ced0346de0p+0, -0x1.d0b4c2d337e0fp-1,
136 0x1.db2dcaaa022a9p-1, -0x1.1188790a072bap+0, 0x1.4e50420966212p+0,
137 -0x1.a787465391482p+0, 0x1.12d6375bce19ep+1, -0x1.6b05729f76a9cp+1,
138 0x1.e602a94d73479p+1, -0x1.48d5f84ddbdefp+2, 0x1.c0e7f7d05fe6dp+2,
139 -0x1.34d52338a0d93p+3, 0x1.ab7e04eff787dp+3, -0x1.2724d2ace1c81p+4,
140 0x1.9b17bee5ae20ap+4, -0x1.34beabd4deff5p+5, 0x1.bc510ff62dcf0p+5,
141 -0x1.7406b7594b6ddp+5, 0x1.b97b7265194a1p+5, -0x1.5e2ba44fee700p+8,
142 0x1.18310d7d43c46p+9, 0x1.59472c8ca4419p+9, -0x1.2a61147e20ef1p+10,
143 -0x1.7cfcf784946edp+11, 0x1.2436ec1277b00p+12};
144 double d = xd - MID_SN;
145 double d2 = d * d, d4 = d2 * d2, d8 = d4 * d4, d16 = d8 * d8;
146 double p01 = fputil::multiply_add(x: d, y: POLY_SN[1], z: POLY_SN[0]);
147 double p23 = fputil::multiply_add(x: d, y: POLY_SN[3], z: POLY_SN[2]);
148 double p45 = fputil::multiply_add(x: d, y: POLY_SN[5], z: POLY_SN[4]);
149 double p67 = fputil::multiply_add(x: d, y: POLY_SN[7], z: POLY_SN[6]);
150 double p89 = fputil::multiply_add(x: d, y: POLY_SN[9], z: POLY_SN[8]);
151 double p1011 = fputil::multiply_add(x: d, y: POLY_SN[11], z: POLY_SN[10]);
152 double p1213 = fputil::multiply_add(x: d, y: POLY_SN[13], z: POLY_SN[12]);
153 double p1415 = fputil::multiply_add(x: d, y: POLY_SN[15], z: POLY_SN[14]);
154 double p1617 = fputil::multiply_add(x: d, y: POLY_SN[17], z: POLY_SN[16]);
155 double p1819 = fputil::multiply_add(x: d, y: POLY_SN[19], z: POLY_SN[18]);
156 double p2021 = fputil::multiply_add(x: d, y: POLY_SN[21], z: POLY_SN[20]);
157 double p2223 = fputil::multiply_add(x: d, y: POLY_SN[23], z: POLY_SN[22]);
158 double p2425 = fputil::multiply_add(x: d, y: POLY_SN[25], z: POLY_SN[24]);
159 double q03 = fputil::multiply_add(x: d2, y: p23, z: p01);
160 double q47 = fputil::multiply_add(x: d2, y: p67, z: p45);
161 double q811 = fputil::multiply_add(x: d2, y: p1011, z: p89);
162 double q1215 = fputil::multiply_add(x: d2, y: p1415, z: p1213);
163 double q1619 = fputil::multiply_add(x: d2, y: p1819, z: p1617);
164 double q2023 = fputil::multiply_add(x: d2, y: p2223, z: p2021);
165 double r07 = fputil::multiply_add(x: d4, y: q47, z: q03);
166 double r815 = fputil::multiply_add(x: d4, y: q1215, z: q811);
167 double r1623 = fputil::multiply_add(x: d4, y: q2023, z: q1619);
168 double s015 = fputil::multiply_add(x: d8, y: r815, z: r07);
169 double s1625 = fputil::multiply_add(x: d8, y: p2425, z: r1623);
170 double poly_g = fputil::multiply_add(x: d16, y: s1625, z: s015);
171 // lgamma(x) = x*g(x) - log(-x) with single rounding via FMA.
172 lgamma_val = fputil::multiply_add(x: xd, y: poly_g, z: -math::log(x: abs_xd));
173 } else {
174 // x = 0x1.f8a754p-9f
175 if (LIBC_UNLIKELY(xbits.uintval() == 0x3b7c53aau))
176 return fputil::round_result_slightly_up(value_rn: 0x1.63acc2p+2f);
177 // Small: t = x < 0.66015625. Degree-17 monomial fit of
178 // h(t) = (lgamma(t) + log(t)) / t (smooth on [0, 0.66]).
179 // d = t - MID is exact.
180 constexpr double MID_S = 0x1.5200000000000p-2;
181 constexpr double POLY_S[18] = {
182 -0x1.5dcd7586bfd88p-2, 0x1.3f88851c787c3p-1, -0x1.cdfaca3081737p-3,
183 0x1.cfac7b321198bp-4, -0x1.0939a239f89b2p-4, 0x1.44cddc5f52a43p-5,
184 -0x1.9e4cc54acdcebp-6, 0x1.0f61691962355p-6, -0x1.6a429a6c89071p-7,
185 0x1.ea573d148429cp-8, -0x1.4f6ad0c206ae6p-8, 0x1.cedaf9fc19cb5p-9,
186 -0x1.42333c4e88f7dp-9, 0x1.c2a0e3ffc858fp-10, -0x1.32153d4cb9316p-10,
187 0x1.add42f56d3c7ap-11, -0x1.972ff2b1a91edp-11, 0x1.25f56c66e3049p-11};
188 double d = abs_xd - MID_S;
189 double d2 = d * d, d4 = d2 * d2, d8 = d4 * d4, d16 = d8 * d8;
190 double p01 = fputil::multiply_add(x: d, y: POLY_S[1], z: POLY_S[0]);
191 double p23 = fputil::multiply_add(x: d, y: POLY_S[3], z: POLY_S[2]);
192 double p45 = fputil::multiply_add(x: d, y: POLY_S[5], z: POLY_S[4]);
193 double p67 = fputil::multiply_add(x: d, y: POLY_S[7], z: POLY_S[6]);
194 double p89 = fputil::multiply_add(x: d, y: POLY_S[9], z: POLY_S[8]);
195 double p1011 = fputil::multiply_add(x: d, y: POLY_S[11], z: POLY_S[10]);
196 double p1213 = fputil::multiply_add(x: d, y: POLY_S[13], z: POLY_S[12]);
197 double p1415 = fputil::multiply_add(x: d, y: POLY_S[15], z: POLY_S[14]);
198 double p1617 = fputil::multiply_add(x: d, y: POLY_S[17], z: POLY_S[16]);
199 double q03 = fputil::multiply_add(x: d2, y: p23, z: p01);
200 double q47 = fputil::multiply_add(x: d2, y: p67, z: p45);
201 double q811 = fputil::multiply_add(x: d2, y: p1011, z: p89);
202 double q1215 = fputil::multiply_add(x: d2, y: p1415, z: p1213);
203 double r07 = fputil::multiply_add(x: d4, y: q47, z: q03);
204 double r815 = fputil::multiply_add(x: d4, y: q1215, z: q811);
205 double s015 = fputil::multiply_add(x: d8, y: r815, z: r07);
206 double poly_h = fputil::multiply_add(x: d16, y: p1617, z: s015);
207 // poly_val - log(abs_xd) with single rounding via FMA.
208 lgamma_val = fputil::multiply_add(x: abs_xd, y: poly_h, z: -math::log(x: abs_xd));
209 }
210 } else if (x_abs < 0x3f800000u) {
211 if (xbits.is_neg()) {
212 // x in (-1, -0.66015625]: Gamma(x) = Gamma(x+2)/(x(x+1)), so
213 // lgamma(x) = x(x+1)*P_M2(x+0.5) - log(-x(x+1)). x+0.5 and x+1 are
214 // exact (Sterbenz); x(x+1) = (t-1)(t-2) is the M2 prefactor.
215 double w = xd * (xd + 1.0);
216 double poly = lgamma_m2_poly(d: xd + 0.5);
217 lgamma_val = fputil::multiply_add(x: w, y: poly, z: -lg_ln(x: -w));
218 } else {
219 // M1: t in [0.66015625, 1.0). lgamma(t) = (t-1) * P_M1(d), d = t - MID.
220 // Degree-14 monomial fit of lgamma(t)/(t-1), max error 2^-51.3.
221 constexpr double MID_M1 = 0x1.a900000000000p-1;
222 constexpr double POLY_M1[15] = {
223 -0x1.75cb89aad8dc7p-1, 0x1.f958114ec790dp-1, -0x1.2b968b172a05cp-1,
224 0x1.eaf0e2c8dc49fp-2, -0x1.c6efb493bfadfp-2, 0x1.c0a621f0e593ep-2,
225 -0x1.cb2416349da49p-2, 0x1.e19328033256dp-2, -0x1.010d0c64743f5p-1,
226 0x1.162e14c80bdc9p-1, -0x1.3034dec6eb9a3p-1, 0x1.4ca601da8405ep-1,
227 -0x1.70f73638718d4p-1, 0x1.ded7569ab5a36p-1, -0x1.0fdd39040adfbp+0};
228 double d = abs_xd - MID_M1;
229 double d2 = d * d, d4 = d2 * d2, d8 = d4 * d4;
230 double p01 = fputil::multiply_add(x: d, y: POLY_M1[1], z: POLY_M1[0]);
231 double p23 = fputil::multiply_add(x: d, y: POLY_M1[3], z: POLY_M1[2]);
232 double p45 = fputil::multiply_add(x: d, y: POLY_M1[5], z: POLY_M1[4]);
233 double p67 = fputil::multiply_add(x: d, y: POLY_M1[7], z: POLY_M1[6]);
234 double p89 = fputil::multiply_add(x: d, y: POLY_M1[9], z: POLY_M1[8]);
235 double p1011 = fputil::multiply_add(x: d, y: POLY_M1[11], z: POLY_M1[10]);
236 double p1213 = fputil::multiply_add(x: d, y: POLY_M1[13], z: POLY_M1[12]);
237 double q03 = fputil::multiply_add(x: d2, y: p23, z: p01);
238 double q47 = fputil::multiply_add(x: d2, y: p67, z: p45);
239 double q811 = fputil::multiply_add(x: d2, y: p1011, z: p89);
240 double q1214 = fputil::multiply_add(x: d2, y: POLY_M1[14], z: p1213);
241 double r07 = fputil::multiply_add(x: d4, y: q47, z: q03);
242 double r814 = fputil::multiply_add(x: d4, y: q1214, z: q811);
243 double poly = fputil::multiply_add(x: d8, y: r814, z: r07);
244 lgamma_val = (abs_xd - 1.0) * poly;
245 }
246 } else if (x_abs < 0x40000000u) {
247 if (xbits.is_neg()) {
248 // x in (-2, -1): Gamma(x) = Gamma(x+3)/(x(x+1)(x+2)), so
249 // lgamma(x) = (x+1)(x+2)*P_M2(x+1.5) - log(x(x+1)(x+2)). The shifts
250 // are exact (Sterbenz); (x+1)(x+2) = (t-1)(t-2) is the M2 prefactor.
251 double u = (xd + 1.0) * (xd + 2.0);
252 double poly = lgamma_m2_poly(d: xd + 1.5);
253 lgamma_val = fputil::multiply_add(x: u, y: poly, z: -lg_ln(x: xd * u));
254 } else {
255 // M2: t in [1.0, 2.0). lgamma(t) = (t-1)*(t-2) * P_M2(t - 1.5).
256 double d = abs_xd - 0x1.8p+0;
257 lgamma_val = (abs_xd - 1.0) * (abs_xd - 2.0) * lgamma_m2_poly(d);
258 }
259 } else if (x_abs < 0x4057e000u) {
260 // M3: t in [2.0, 3.373046875). lgamma(t) = (t-2) * P_M3(d), d = t - MID.
261 // Degree-15 monomial fit of lgamma(t)/(t-2), max error 2^-49.7.
262 constexpr double MID_M3 = 0x1.57e0000000000p+1;
263 constexpr double POLY_M3[16] = {
264 0x1.3c4e36a0b4775p-1, 0x1.01f945be1325fp-2, -0x1.4203a95730c77p-5,
265 0x1.227f82db7c7a1p-7, -0x1.32f092b6ec5a3p-9, 0x1.61df821ae7829p-11,
266 -0x1.aed17eeca55e6p-13, 0x1.100a8e7dcd62cp-14, -0x1.609aa16f960a0p-16,
267 0x1.d1e293fa801bbp-18, -0x1.38b436ce1b3b9p-19, 0x1.a838eec563338p-21,
268 -0x1.1a6a9cf4aee4bp-22, 0x1.8387c2a068b06p-24, -0x1.5ee6ae8c133f7p-25,
269 0x1.f195cf3f4b24ep-27};
270 if (!xbits.is_neg()) {
271 double d = abs_xd - MID_M3;
272 double d2 = d * d, d4 = d2 * d2, d8 = d4 * d4;
273 double p01 = fputil::multiply_add(x: d, y: POLY_M3[1], z: POLY_M3[0]);
274 double p23 = fputil::multiply_add(x: d, y: POLY_M3[3], z: POLY_M3[2]);
275 double p45 = fputil::multiply_add(x: d, y: POLY_M3[5], z: POLY_M3[4]);
276 double p67 = fputil::multiply_add(x: d, y: POLY_M3[7], z: POLY_M3[6]);
277 double p89 = fputil::multiply_add(x: d, y: POLY_M3[9], z: POLY_M3[8]);
278 double p1011 = fputil::multiply_add(x: d, y: POLY_M3[11], z: POLY_M3[10]);
279 double p1213 = fputil::multiply_add(x: d, y: POLY_M3[13], z: POLY_M3[12]);
280 double p1415 = fputil::multiply_add(x: d, y: POLY_M3[15], z: POLY_M3[14]);
281 double q03 = fputil::multiply_add(x: d2, y: p23, z: p01);
282 double q47 = fputil::multiply_add(x: d2, y: p67, z: p45);
283 double q811 = fputil::multiply_add(x: d2, y: p1011, z: p89);
284 double q1215 = fputil::multiply_add(x: d2, y: p1415, z: p1213);
285 double r07 = fputil::multiply_add(x: d4, y: q47, z: q03);
286 double r815 = fputil::multiply_add(x: d4, y: q1215, z: q811);
287 double poly = fputil::multiply_add(x: d8, y: r815, z: r07);
288 lgamma_val = (abs_xd - 2.0) * poly;
289 } else {
290 // Near the regular lgamma zero at x ~= -2.7475: subtractive cancellation
291 // in the reflection formula kills precision. Use a Taylor expansion
292 // centered at the zero. Range bits in (0x402f95c2, 0x40301b93).
293 // Coefficients adopted from CORE-MATH (Sibidanov, 2023).
294 if (LIBC_UNLIKELY(x_abs > 0x402f95c2u && x_abs < 0x40301b93u)) {
295 double h = (xd + 0x1.5fb410a1bd901p+1) - 0x1.a19a96d2e6f85p-54;
296 constexpr double C[8] = {-0x1.ea12da904b18cp+0, 0x1.3267f3c265a54p+3,
297 -0x1.4185ac30cadb3p+4, 0x1.f504accc3f2e4p+5,
298 -0x1.8588444c679b4p+7, 0x1.43740491dc22p+9,
299 -0x1.12400ea23f9e6p+11, 0x1.dac829f365795p+12};
300 double h2 = h * h, h4 = h2 * h2;
301 double p01 = fputil::multiply_add(x: h, y: C[1], z: C[0]);
302 double p23 = fputil::multiply_add(x: h, y: C[3], z: C[2]);
303 double p45 = fputil::multiply_add(x: h, y: C[5], z: C[4]);
304 double p67 = fputil::multiply_add(x: h, y: C[7], z: C[6]);
305 double p03 = fputil::multiply_add(x: h2, y: p23, z: p01);
306 double p47 = fputil::multiply_add(x: h2, y: p67, z: p45);
307 lgamma_val = h * fputil::multiply_add(x: h4, y: p47, z: p03);
308 } else if (LIBC_UNLIKELY(x_abs > 0x401ceccbu && x_abs < 0x401d95cau)) {
309 // Near the regular lgamma zero at x ~= -2.3614: same issue
310 double h = (xd + 0x1.3a7fc9600f86cp+1) + 0x1.55f64f98af8dp-55;
311 constexpr double C[7] = {0x1.83fe966af535fp+0, 0x1.36eebb002f61ap+2,
312 0x1.694a60589a0b3p+0, 0x1.1718d7aedb0b5p+3,
313 0x1.733a045eca0d3p+2, 0x1.8d4297421205bp+4,
314 0x1.7feea5fb29965p+4};
315 double h2 = h * h, h4 = h2 * h2;
316 double p01 = fputil::multiply_add(x: h, y: C[1], z: C[0]);
317 double p23 = fputil::multiply_add(x: h, y: C[3], z: C[2]);
318 double p45 = fputil::multiply_add(x: h, y: C[5], z: C[4]);
319 double p46 = fputil::multiply_add(x: h2, y: C[6], z: p45);
320 double p03 = fputil::multiply_add(x: h2, y: p23, z: p01);
321 lgamma_val = h * fputil::multiply_add(x: h4, y: p46, z: p03);
322 } else if (LIBC_UNLIKELY(x_abs > 0x40492009u && x_abs < 0x404940efu)) {
323 // Near the regular lgamma zero at x ~= -3.1431: same issue
324 double h = (xd + 0x1.9260dbc9e59afp+1) + 0x1.f717cd335a7b3p-53;
325 constexpr double C[7] = {0x1.f20a65f2fac55p+2, 0x1.9d4d297715105p+4,
326 0x1.c1137124d5b21p+6, 0x1.267203d24de38p+9,
327 0x1.99a63399a0b44p+11, 0x1.2941214faaf0cp+14,
328 0x1.bb912c0c9cdd1p+16};
329 double h2 = h * h, h4 = h2 * h2;
330 double p01 = fputil::multiply_add(x: h, y: C[1], z: C[0]);
331 double p23 = fputil::multiply_add(x: h, y: C[3], z: C[2]);
332 double p45 = fputil::multiply_add(x: h, y: C[5], z: C[4]);
333 double p46 = fputil::multiply_add(x: h2, y: C[6], z: p45);
334 double p03 = fputil::multiply_add(x: h2, y: p23, z: p01);
335 lgamma_val = h * fputil::multiply_add(x: h4, y: p46, z: p03);
336 } else {
337 // x in (-3.373, -2): Gamma(x) = Gamma(x+5)/(x(x+1)(x+2)(x+3)(x+4)),
338 // t = x+5 in (1.627, 3), so lgamma(x) = (x+3)*Q_M3N(d) - log|prod|
339 // with d = x + 2.6865234375 (exact). Q_M3N is a degree-17 monomial
340 // fit of lgamma(t)/(t-2) on [1.627, 3]; the (t-2) = x+3 prefactor
341 // removes the lgamma zero at t = 2.
342 constexpr double POLY_M3N[18] = {
343 0x1.091ff92b41f07p-1, 0x1.245eed13c42b0p-2,
344 -0x1.a6c3cf8a165cfp-5, 0x1.bc8d19a3ade02p-7,
345 -0x1.1231651010470p-8, 0x1.710f57a2abe2bp-10,
346 -0x1.061ddab9ab12cp-11, 0x1.81edb077ed799p-13,
347 -0x1.235fbc13dbf04p-14, 0x1.c02e6f8fb5dbfp-16,
348 -0x1.5d75e6a94d352p-17, 0x1.1384596ce083dp-18,
349 -0x1.b8da81a716039p-20, 0x1.61ac57fe0e288p-21,
350 -0x1.077101d536ad8p-22, 0x1.a59e2356ba870p-24,
351 -0x1.0cfcf8f102a43p-24, 0x1.c141a51f827f1p-26};
352 double d = xd + 0x1.57ep+1;
353 double d2 = d * d, d4 = d2 * d2, d8 = d4 * d4, d16 = d8 * d8;
354 double p01 = fputil::multiply_add(x: d, y: POLY_M3N[1], z: POLY_M3N[0]);
355 double p23 = fputil::multiply_add(x: d, y: POLY_M3N[3], z: POLY_M3N[2]);
356 double p45 = fputil::multiply_add(x: d, y: POLY_M3N[5], z: POLY_M3N[4]);
357 double p67 = fputil::multiply_add(x: d, y: POLY_M3N[7], z: POLY_M3N[6]);
358 double p89 = fputil::multiply_add(x: d, y: POLY_M3N[9], z: POLY_M3N[8]);
359 double p1011 = fputil::multiply_add(x: d, y: POLY_M3N[11], z: POLY_M3N[10]);
360 double p1213 = fputil::multiply_add(x: d, y: POLY_M3N[13], z: POLY_M3N[12]);
361 double p1415 = fputil::multiply_add(x: d, y: POLY_M3N[15], z: POLY_M3N[14]);
362 double p1617 = fputil::multiply_add(x: d, y: POLY_M3N[17], z: POLY_M3N[16]);
363 double q03 = fputil::multiply_add(x: d2, y: p23, z: p01);
364 double q47 = fputil::multiply_add(x: d2, y: p67, z: p45);
365 double q811 = fputil::multiply_add(x: d2, y: p1011, z: p89);
366 double q1215 = fputil::multiply_add(x: d2, y: p1415, z: p1213);
367 double r07 = fputil::multiply_add(x: d4, y: q47, z: q03);
368 double r815 = fputil::multiply_add(x: d4, y: q1215, z: q811);
369 double s015 = fputil::multiply_add(x: d8, y: r815, z: r07);
370 double poly = fputil::multiply_add(x: d16, y: p1617, z: s015);
371 double pa = xd * (xd + 1.0);
372 double pb = (xd + 2.0) * (xd + 3.0);
373 double prod = (pa * pb) * (xd + 4.0);
374 double aprod = prod < 0.0 ? -prod : prod;
375 lgamma_val = fputil::multiply_add(x: xd + 3.0, y: poly, z: -lg_ln(x: aprod));
376 }
377 }
378 } else {
379 // Large: |x| >= 3.373046875. Stirling + Bernoulli correction.
380 // lgamma(x) = (x-0.5)*log(x) - x + log(2*pi)/2 + (1/x)*P(1/x^2)
381 // = (x-0.5)*(log(x)-1) + STIR_CONST + (1/x)*P(1/x^2)
382 // STIR_CONST = log(2*pi)/2 - 0.5.
383 // For huge positive x, lgamma(x) overflows float. Use a linear
384 // approximation in double that maps to the correct Inf/max_normal.
385 if (LIBC_UNLIKELY(!xbits.is_neg() && x >= 0x1.895f1cp+121f)) {
386 fputil::set_errno_if_required(ERANGE);
387 fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
388 double r = fputil::multiply_add(x: xd, y: 0x1.4d3398p+6, z: 0x1.10f35ep+103);
389 return fputil::cast<float>(x: r);
390 }
391
392 // No cancellation in (x-0.5)*(log(x)-1): log(x)-1 > 0.2 on this range.
393 // Relative error ~2^-48; anything closer to a rounding boundary is in
394 // the exceptional cases tables below.
395 double lz = lg_ln(x: abs_xd);
396 double xm = abs_xd - 0.5;
397 lgamma_val = fputil::multiply_add(x: xm, y: lz - 1.0, z: 0x1.acfe390c97d69p-2);
398
399 // For |x| >= 2^20 the 1/(12x) correction is below ~2^-47 of the result;
400 // skip it and the 1/x divide.
401 if (x_abs < 0x49800000u) {
402 double inv_x = 1.0 / abs_xd;
403 double inv_x2 = inv_x * inv_x;
404 if (x_abs > 0x44fa0000u) {
405 constexpr fputil::ExceptValues<float, 3> LGAMMAF_EXCEPTS_BERN2{.values: {
406 // input, toward-zero result, RU, RD, RN
407 {.input: 0x46541516u, .rnd_towardzero_result: 0x47e1c01bu, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
408 {.input: 0x46b16323u, .rnd_towardzero_result: 0x48483adeu, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
409 {.input: 0xc6f7e151u, .rnd_towardzero_result: 0xc89116deu, .rnd_upward_offset: 0, .rnd_downward_offset: 1, .rnd_tonearest_offset: 0},
410 }};
411 if (auto r = LGAMMAF_EXCEPTS_BERN2.lookup(x_bits: xbits.uintval());
412 LIBC_UNLIKELY(r.has_value()))
413 return r.value();
414 // |x| > 2000 -> 2-term BERN2.
415 constexpr double BERN2[2] = {0x1.5555555555555p-4,
416 -0x1.6c16bfb7c65a8p-9};
417 lgamma_val += inv_x * fputil::multiply_add(x: inv_x2, y: BERN2[1], z: BERN2[0]);
418 } else if (x_abs > 0x42920000u) {
419 // Exceptional cases of this range.
420 constexpr fputil::ExceptValues<float, 2> LGAMMAF_EXCEPTS_BERN4{.values: {
421 // input, toward-zero result, RU, RD, RN
422 {.input: 0x449acf07u, .rnd_towardzero_result: 0x45ecd680u, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
423 {.input: 0xc33139a3u, .rnd_towardzero_result: 0xc43991afu, .rnd_upward_offset: 0, .rnd_downward_offset: 1, .rnd_tonearest_offset: 0},
424 }};
425 if (auto r = LGAMMAF_EXCEPTS_BERN4.lookup(x_bits: xbits.uintval());
426 LIBC_UNLIKELY(r.has_value()))
427 return r.value();
428 // |x| > 73 -> 4-term BERN4.
429 constexpr double BERN4[4] = {
430 0x1.5555555555555p-4, -0x1.6c16c16c15f75p-9, 0x1.a01a00593b36fp-11,
431 -0x1.37e91273668efp-11};
432 double inv_x4 = inv_x2 * inv_x2;
433 double p01 = fputil::multiply_add(x: inv_x2, y: BERN4[1], z: BERN4[0]);
434 double p23 = fputil::multiply_add(x: inv_x2, y: BERN4[3], z: BERN4[2]);
435 lgamma_val += inv_x * fputil::multiply_add(x: inv_x4, y: p23, z: p01);
436 } else {
437 // Exceptional cases of this range.
438 constexpr fputil::ExceptValues<float, 2> LGAMMAF_EXCEPTS_B10{.values: {
439 // input, toward-zero result, RU, RD, RN
440 {.input: 0x42468b59u, .rnd_towardzero_result: 0x430f25a7u, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
441 {.input: 0xc134eb14u, .rnd_towardzero_result: 0xc1875615u, .rnd_upward_offset: 0, .rnd_downward_offset: 1, .rnd_tonearest_offset: 0},
442 }};
443 if (auto r = LGAMMAF_EXCEPTS_B10.lookup(x_bits: xbits.uintval());
444 LIBC_UNLIKELY(r.has_value()))
445 return r.value();
446 // |x| in (3.373, 73]: degree-10 monomial fit of
447 // h(s) = stir_resid(1/sqrt(s)) / sqrt(s), s = 1/x^2, on [0, 0.088];
448 // correction = h(s) * (1/x). Max error 2^-53.2.
449 constexpr double MID_B10 = 0x1.6880000000000p-5;
450 constexpr double POLY_B10[11] = {
451 0x1.54d6b78cee955p-4, -0x1.635a5fb0cdf9fp-9,
452 0x1.7b5253f44b255p-11, -0x1.f32907e7a7adap-12,
453 0x1.1f269b6438739p-11, -0x1.e95dc64c5042cp-11,
454 0x1.17ebce09c49a7p-9, -0x1.91599e7728747p-8,
455 0x1.5a01e1f07b127p-6, -0x1.90e956754bc6ap-4,
456 0x1.dfbed80c6f035p-2};
457 double u = inv_x2 - MID_B10;
458 double u2 = u * u, u4 = u2 * u2, u8 = u4 * u4;
459 double p01 = fputil::multiply_add(x: u, y: POLY_B10[1], z: POLY_B10[0]);
460 double p23 = fputil::multiply_add(x: u, y: POLY_B10[3], z: POLY_B10[2]);
461 double p45 = fputil::multiply_add(x: u, y: POLY_B10[5], z: POLY_B10[4]);
462 double p67 = fputil::multiply_add(x: u, y: POLY_B10[7], z: POLY_B10[6]);
463 double p89 = fputil::multiply_add(x: u, y: POLY_B10[9], z: POLY_B10[8]);
464 double q03 = fputil::multiply_add(x: u2, y: p23, z: p01);
465 double q47 = fputil::multiply_add(x: u2, y: p67, z: p45);
466 double q810 = fputil::multiply_add(x: u2, y: POLY_B10[10], z: p89);
467 double r07 = fputil::multiply_add(x: u4, y: q47, z: q03);
468 double poly = fputil::multiply_add(x: u8, y: q810, z: r07);
469 lgamma_val += inv_x * poly;
470 }
471 } else {
472 constexpr fputil::ExceptValues<float, 3> LGAMMAF_EXCEPTS_HUGE{.values: {
473 // input, toward-zero result, RU, RD, RN
474 {.input: 0x65fca09fu, .rnd_towardzero_result: 0x68cead59u, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
475 {.input: 0x716e5dd5u, .rnd_towardzero_result: 0x747e2bb9u, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
476 {.input: 0x77ac5674u, .rnd_towardzero_result: 0x7acf27b2u, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
477 }};
478 if (auto r = LGAMMAF_EXCEPTS_HUGE.lookup(x_bits: xbits.uintval());
479 LIBC_UNLIKELY(r.has_value()))
480 return r.value();
481 }
482
483 if (xbits.is_neg()) {
484 // Reflection: lgamma(x) = log(pi) - lgamma(|x|) - log(|x|) -
485 // log(|sin(pi*frac_x)|). Reusing the already-computed lz = log(|x|)
486 // keeps the multiply out of the serial sinpi -> log chain.
487 double frac_x = xd - fputil::floor(x: xd);
488 lgamma_val = (0x1.250d048e7a1bdp+0 - lgamma_val) - lz;
489 lgamma_val -= lg_ln(x: lg_sinpi(x: frac_x));
490 }
491 }
492
493 float result = fputil::cast<float>(x: lgamma_val);
494 if (LIBC_UNLIKELY(FPBits(result).is_inf())) {
495 fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
496 fputil::set_errno_if_required(ERANGE);
497 }
498 return result;
499}
500
501} // namespace math
502
503} // namespace LIBC_NAMESPACE_DECL
504
505#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMAF_H
506