@@ -93,7 +93,38 @@ struct lgamma_taylor_series_imp
9393 + decimal64 { UINT64_C (384'615'390'346'751'857 ), -19 }, // x^26
9494 }};
9595
96- // TODO: add d128_coeffs
96+ static constexpr std::array<decimal128, 25 > d128_coeffs =
97+ {{
98+ // Use a Taylor series expansion of the logarithm of the gamma function.
99+ // N[Series[Log[Gamma[x]], {x, 0, 26}], 32]
100+ // log(1/x)
101+ // -EulerGamma // * x
102+ decimal128 { UINT64_C (822'467'033'424'113'218 ), -18 }, // x^2
103+ - decimal128 { UINT64_C (400'685'634'386'531'428 ), -18 }, // x^3
104+ + decimal128 { UINT64_C (270'580'808'427'784'548 ), -18 }, // x^4
105+ - decimal128 { UINT64_C (207'385'551'028'673'985 ), -18 }, // x^5
106+ + decimal128 { UINT64_C (169'557'176'997'408'190 ), -18 }, // x^6
107+ - decimal128 { UINT64_C (144'049'896'768'846'118 ), -18 }, // x^7
108+ + decimal128 { UINT64_C (125'509'669'524'743'042 ), -18 }, // x^8
109+ - decimal128 { UINT64_C (111'334'265'869'564'690 ), -18 }, // x^9
110+ + decimal128 { UINT64_C (100'099'457'512'781'809 ), -18 }, // x^10
111+ - decimal128 { UINT64_C (909'540'171'458'290'422 ), -19 }, // x^11
112+ + decimal128 { UINT64_C (833'538'405'461'090'040 ), -19 }, // x^12
113+ - decimal128 { UINT64_C (769'325'164'113'521'915 ), -19 }, // x^13
114+ + decimal128 { UINT64_C (714'329'462'953'613'361 ), -19 }, // x^14
115+ - decimal128 { UINT64_C (666'687'058'824'204'680 ), -19 }, // x^15
116+ + decimal128 { UINT64_C (625'009'551'412'130'407 ), -19 }, // x^16
117+ - decimal128 { UINT64_C (588'239'786'586'845'823 ), -19 }, // x^17
118+ + decimal128 { UINT64_C (555'557'676'274'036'111 ), -19 }, // x^18
119+ - decimal128 { UINT64_C (526'316'793'796'166'607 ), -19 }, // x^19
120+ + decimal128 { UINT64_C (500'000'476'981'016'936 ), -19 }, // x^20
121+ - decimal128 { UINT64_C (476'190'703'301'422'280 ), -19 }, // x^21
122+ + decimal128 { UINT64_C (454'545'562'932'046'694 ), -19 }, // x^22
123+ - decimal128 { UINT64_C (434'782'660'530'402'594 ), -19 }, // x^23
124+ + decimal128 { UINT64_C (416'666'691'503'412'105 ), -19 }, // x^24
125+ - decimal128 { UINT64_C (400'000'011'921'401'406 ), -19 }, // x^25
126+ + decimal128 { UINT64_C (384'615'390'346'751'857 ), -19 }, // x^26
127+ }};
97128};
98129
99130template <bool b>
@@ -166,6 +197,40 @@ struct lgamma_laurent_series_imp
166197 - decimal64 { UINT64_C (328'565'400'725'242'178 ), + 2 - 18 }, // * x^-25
167198 + decimal64 { UINT64_C (549'592'170'046'323'711 ), + 4 - 18 }, // * x^-26
168199 }};
200+
201+ static constexpr std::array<decimal128, 26 > d128_coeffs = {{
202+ // Use a Laurent infinite-series expansion of the logarithm of the
203+ // gamma function divided by the square root.
204+ // N[Series[Log[Gamma[x]], {x, Infinity, 26}], 32]
205+ // log( e^-x * x^x * sqrt(x) * [Series below] )
206+ // See also Wolfram Alpha(R): https://www.wolframalpha.com/input?i=Series%5BLog%5BGamma%5Bx%5D%5D%2C+%7Bx%2C+Infinity%2C+4%7D%5D
207+ decimal128 { UINT64_C (250'662'827'463'100'050 ), + 1 - 18 }, // * x^-1
208+ + decimal128 { UINT64_C (208'885'689'552'583'375 ), - 0 - 18 }, // * x^-2
209+ + decimal128 { UINT64_C (870'357'039'802'430'730 ), - 2 - 18 }, // * x^-3
210+ - decimal128 { UINT64_C (672'109'047'402'988'175 ), - 2 - 18 }, // * x^-4
211+ - decimal128 { UINT64_C (575'201'238'110'171'235 ), - 3 - 18 }, // * x^-5
212+ + decimal128 { UINT64_C (196'529'488'158'320'306 ), - 2 - 18 }, // * x^-6
213+ + decimal128 { UINT64_C (174'782'521'204'559'121 ), - 3 - 18 }, // * x^-7
214+ - decimal128 { UINT64_C (148'434'113'515'827'614 ), - 2 - 18 }, // * x^-8
215+ - decimal128 { UINT64_C (129'637'573'211'255'432 ), - 3 - 18 }, // * x^-9
216+ + decimal128 { UINT64_C (210'431'122'975'320'637 ), - 2 - 18 }, // * x^-10
217+ + decimal128 { UINT64_C (180'599'945'655'550'436 ), - 3 - 18 }, // * x^-11
218+ - decimal128 { UINT64_C (479'878'567'054'634'606 ), - 2 - 18 }, // * x^-12
219+ - decimal128 { UINT64_C (407'367'859'381'525'183 ), - 3 - 18 }, // * x^-13
220+ + decimal128 { UINT64_C (160'508'503'319'445'960 ), - 1 - 18 }, // * x^-14
221+ + decimal128 { UINT64_C (135'399'228'015'909'411 ), - 2 - 18 }, // * x^-15
222+ - decimal128 { UINT64_C (740'154'212'684'273'819 ), - 1 - 18 }, // * x^-16
223+ - decimal128 { UINT64_C (622'080'867'880'877'866 ), - 2 - 18 }, // * x^-17
224+ + decimal128 { UINT64_C (450'040'333'856'250'984 ), - 0 - 18 }, // * x^-18
225+ + decimal128 { UINT64_C (377'400'786'521'707'440 ), - 1 - 18 }, // * x^-19
226+ - decimal128 { UINT64_C (348'872'797'304'123'310 ), + 1 - 18 }, // * x^-20
227+ - decimal128 { UINT64_C (292'138'192'227'179'792 ), - 0 - 18 }, // * x^-21
228+ + decimal128 { UINT64_C (335'837'691'649'553'084 ), + 2 - 18 }, // * x^-22
229+ + decimal128 { UINT64_C (280'944'016'052'174'395 ), + 1 - 18 }, // * x^-23
230+ - decimal128 { UINT64_C (393'042'854'585'987'930 ), + 3 - 18 }, // * x^-24
231+ - decimal128 { UINT64_C (328'565'400'725'242'178 ), + 2 - 18 }, // * x^-25
232+ + decimal128 { UINT64_C (549'592'170'046'323'711 ), + 4 - 18 }, // * x^-26
233+ }};
169234};
170235
171236#if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900)
@@ -176,12 +241,18 @@ constexpr std::array<decimal32, 25> lgamma_taylor_series_imp<b>::d32_coeffs;
176241template <bool b>
177242constexpr std::array<decimal64, 25 > lgamma_taylor_series_imp<b>::d64_coeffs;
178243
244+ template <bool b>
245+ constexpr std::array<decimal128, 25 > lgamma_taylor_series_imp<b>::d128_coeffs;
246+
179247template <bool b>
180248constexpr std::array<decimal32, 26 > lgamma_laurent_series_imp<b>::d32_coeffs;
181249
182250template <bool b>
183251constexpr std::array<decimal64, 26 > lgamma_laurent_series_imp<b>::d64_coeffs;
184252
253+ template <bool b>
254+ constexpr std::array<decimal128, 26 > lgamma_laurent_series_imp<b>::d128_coeffs;
255+
185256#endif
186257
187258} // namespace lgamma_detail
@@ -190,14 +261,7 @@ using lgamma_taylor_series_table = lgamma_detail::lgamma_taylor_series_imp<true>
190261using lgamma_laurent_series_table = lgamma_detail::lgamma_laurent_series_imp<true >;
191262
192263template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T>
193- constexpr auto lgamma_taylor_series_expansion (T x) noexcept
194- {
195- // LCOV_EXCL_START
196- static_assert (!std::is_same<T, decimal128>::value, " Decimal128 has not been implemented" );
197- static_cast <void >(x);
198- return T{1 ,0 ,true };
199- // LCOV_EXCL_STOP
200- }
264+ constexpr auto lgamma_taylor_series_expansion (T x) noexcept ;
201265
202266template <>
203267constexpr auto lgamma_taylor_series_expansion<decimal32>(decimal32 x) noexcept
@@ -211,16 +275,15 @@ constexpr auto lgamma_taylor_series_expansion<decimal64>(decimal64 x) noexcept
211275 return taylor_series_result (x, lgamma_taylor_series_table::d64_coeffs);
212276}
213277
214- template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T >
215- constexpr auto lgamma_laurent_series_expansion (T x) noexcept
278+ template <>
279+ constexpr auto lgamma_taylor_series_expansion<decimal128>(decimal128 x) noexcept
216280{
217- // LCOV_EXCL_START
218- static_assert (!std::is_same<T, decimal128>::value, " Decimal128 has not been implemented" );
219- static_cast <void >(x);
220- return T{1 ,0 ,true };
221- // LCOV_EXCL_STOP
281+ return taylor_series_result (x, lgamma_taylor_series_table::d128_coeffs);
222282}
223283
284+ template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T>
285+ constexpr auto lgamma_laurent_series_expansion (T x) noexcept ;
286+
224287template <>
225288constexpr auto lgamma_laurent_series_expansion<decimal32>(decimal32 x) noexcept
226289{
@@ -233,6 +296,12 @@ constexpr auto lgamma_laurent_series_expansion<decimal64>(decimal64 x) noexcept
233296 return taylor_series_result (x, lgamma_laurent_series_table::d64_coeffs);
234297}
235298
299+ template <>
300+ constexpr auto lgamma_laurent_series_expansion<decimal128>(decimal128 x) noexcept
301+ {
302+ return taylor_series_result (x, lgamma_laurent_series_table::d128_coeffs);
303+ }
304+
236305} // namespace detail
237306} // namespace decimal
238307} // namespace boost
0 commit comments