Skip to content

Commit 6899538

Browse files
ckormanyosckormanyos
authored
Merge pull request #551 from cppalliance/finish_gammas
fix #381 and fix #431 and fix #539 via finish gammas
2 parents 49fb1cc + 5904b37 commit 6899538

File tree

9 files changed

+453
-348
lines changed

9 files changed

+453
-348
lines changed

‎include/boost/decimal/charconv.hpp‎

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -318,7 +318,9 @@ BOOST_DECIMAL_CONSTEXPR auto to_chars_scientific_impl(char* first, char* last, c
318318

319319
// Insert the exponent character
320320
*first++ = 'e';
321-
const int abs_exp {std::abs(exp)};
321+
322+
const int abs_exp { (exp < 0) ? -exp : exp };
323+
322324
if (exp < 0)
323325
{
324326
*first++ = '-';
@@ -383,7 +385,9 @@ BOOST_DECIMAL_CONSTEXPR auto to_chars_fixed_impl(char* first, char* last, const
383385

384386
if (integer_digits < 0)
385387
{
386-
num_leading_zeros = std::abs(integer_digits);
388+
const int abs_integer_digits { (integer_digits < 0) ? -integer_digits : integer_digits };
389+
390+
num_leading_zeros = abs_integer_digits;
387391
integer_digits = 0;
388392
append_leading_zeros = true;
389393
}
@@ -683,12 +687,14 @@ BOOST_DECIMAL_CONSTEXPR auto to_chars_hex_impl(char* first, char* last, const Ta
683687
*first++ = '+';
684688
}
685689

686-
if (std::abs(exp) < 10)
690+
const int abs_exp { (exp < 0) ? -exp : exp };
691+
692+
if (abs_exp < 10)
687693
{
688694
*first++ = '0';
689695
}
690696

691-
return to_chars_integer_impl<std::uint32_t, std::uint32_t>(first, last, static_cast<std::uint32_t>(std::abs(exp)), 10);
697+
return to_chars_integer_impl<std::uint32_t, std::uint32_t>(first, last, static_cast<std::uint32_t>(abs_exp), 10);
692698
}
693699

694700
template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE TargetDecimalType>

‎include/boost/decimal/detail/cmath/impl/lgamma_impl.hpp‎

Lines changed: 109 additions & 231 deletions
Large diffs are not rendered by default.

‎include/boost/decimal/detail/cmath/impl/tgamma_impl.hpp‎

Lines changed: 151 additions & 57 deletions
Large diffs are not rendered by default.

‎include/boost/decimal/detail/cmath/lgamma.hpp‎

Lines changed: 19 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -7,16 +7,12 @@
77
#define BOOST_DECIMAL_DETAIL_CMATH_LGAMMA_HPP
88

99
#include <boost/decimal/fwd.hpp> // NOLINT(llvm-include-order)
10-
#include <boost/decimal/detail/type_traits.hpp>
11-
#include <boost/decimal/detail/config.hpp>
1210
#include <boost/decimal/detail/cmath/impl/lgamma_impl.hpp>
13-
#include <boost/decimal/detail/cmath/impl/taylor_series_result.hpp>
14-
#include <boost/decimal/detail/cmath/abs.hpp>
11+
#include <boost/decimal/detail/cmath/impl/tgamma_impl.hpp>
1512
#include <boost/decimal/detail/cmath/log.hpp>
1613
#include <boost/decimal/detail/cmath/sin.hpp>
17-
#include <boost/decimal/detail/cmath/sqrt.hpp>
18-
#include <boost/decimal/detail/cmath/tgamma.hpp>
19-
#include <boost/decimal/detail/concepts.hpp>
14+
#include <boost/decimal/detail/config.hpp>
15+
#include <boost/decimal/detail/type_traits.hpp>
2016
#include <boost/decimal/numbers.hpp>
2117

2218
#ifndef BOOST_DECIMAL_BUILD_MODULE
@@ -78,32 +74,34 @@ constexpr auto lgamma_impl(T x) noexcept
7874
}
7975
else
8076
{
81-
constexpr T half { 5, -1 };
82-
constexpr T asymp_cutoff { 8, 0 };
77+
constexpr int asymp_cutoff
78+
{
79+
std::numeric_limits<T>::digits10 < 10 ? T { 2, 1 }
80+
: std::numeric_limits<T>::digits10 < 20 ? T { 5, 1 }
81+
: T { 1, 2 }
82+
};
8383

84-
if (x < half)
84+
if (x < T { 1, -1 })
8585
{
8686
// Perform the Taylor series expansion.
87-
//result = detail::taylor_series_result(x, coefficient_table);
88-
result = detail::lgamma_taylor_series_expansion(x);
89-
90-
result = x * fma(result, x, -numbers::egamma_v<T>);
9187

92-
result -= log(x);
88+
result = (x * fma(detail::lgamma_taylor_series_expansion(x), x, -numbers::egamma_v<T>))
89+
- log(x);
9390
}
94-
else if (x < asymp_cutoff)
91+
else if (x < T { asymp_cutoff })
9592
{
9693
result = log(tgamma(x));
9794
}
9895
else
9996
{
97+
// Perform the Laurent series expansion. Do note, however, that
98+
// the coefficients of the Laurent asymptotic expansion are exactly
99+
// the same as those used for the tgamma() function.
100100

101-
// Perform the Laurent series expansion.
102-
const auto one_over_x = one / x;
103-
104-
result = detail::lgamma_laurent_series_expansion(one_over_x);
101+
constexpr T half { 5, -1 };
105102

106-
result = (x * (log(x) - one)) + log(result / sqrt(x));
103+
result = (((x - half) * (log(x)) - x))
104+
+ log(detail::tgamma_series_expansion_asymp(one / x));
107105
}
108106
}
109107
}

‎include/boost/decimal/detail/cmath/pow.hpp‎

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -245,12 +245,10 @@ constexpr auto pow(T x, T a) noexcept
245245
}
246246
else if (fpc_a == FP_INFINITE)
247247
{
248-
using std::abs;
249-
250248
result =
251249
(
252-
(abs(x) < one) ? (signbit(a) ? std::numeric_limits<T>::infinity() : zero)
253-
: (abs(x) > one) ? (signbit(a) ? zero : std::numeric_limits<T>::infinity())
250+
(fabs(x) < one) ? (signbit(a) ? std::numeric_limits<T>::infinity() : zero)
251+
: (fabs(x) > one) ? (signbit(a) ? zero : std::numeric_limits<T>::infinity())
254252
: one
255253
);
256254
}

‎include/boost/decimal/detail/cmath/rint.hpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,10 @@ constexpr auto rint_impl(T1& sig, T2 exp, bool sign)
3434
using RoundType = std::conditional_t<std::is_same<T1, std::uint32_t>::value, decimal32,
3535
std::conditional_t<std::is_same<T1, std::uint64_t>::value, decimal64, decimal128>>;
3636

37-
sig /= detail::pow10(static_cast<T1>(std::abs(exp) - 1));
37+
const T2 abs_exp { (exp < T2(0)) ? -exp : exp };
38+
39+
sig /= detail::pow10(static_cast<T1>(abs_exp - 1));
40+
3841
detail::fenv_round<RoundType>(sig, sign);
3942
}
4043

‎include/boost/decimal/detail/cmath/tgamma.hpp‎

Lines changed: 21 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -79,35 +79,38 @@ constexpr auto tgamma_impl(T x) noexcept
7979
}
8080
else
8181
{
82-
const auto x_is_gt_one = (x > one);
82+
constexpr int asymp_cutoff
83+
{
84+
std::numeric_limits<T>::digits10 < 10 ? T { 2, 1 }
85+
: std::numeric_limits<T>::digits10 < 20 ? T { 5, 1 }
86+
: T { 1, 2 }
87+
};
8388

84-
auto r = one;
89+
if (x < T { asymp_cutoff })
90+
{
91+
T r { 1 };
8592

86-
auto z = x;
93+
T z { x };
8794

88-
if (x_is_gt_one)
89-
{
90-
// Use scaling for arguments greater than one.
91-
// TODO(ckormanyos) The work of upscaling can become excessive for large argument.
92-
// TODO(ckormanyos) For large argument (above a cutoff) use asymptotic expansion.
95+
// Use small-argument Taylor series expansion
96+
// (with scaling for arguments greater than one).
9397

9498
for(auto k = 1; k <= nx; ++k)
9599
{
96-
r = r * (z - k);
100+
r *= (z - k);
97101
}
98102

99-
z = z - nx;
100-
}
103+
z -= nx;
101104

102-
result = detail::tgamma_series_expansion(z);
103-
result = one / (z * fma(result, z, one));
104-
105-
if (x_is_gt_one)
105+
result = r / (z * fma(detail::tgamma_series_expansion(z), z, one));
106+
}
107+
else
106108
{
107-
// Downscale: From above, when scaling for arguments greater than one.
108-
// TODO(ckormanyos) See related notes above regarding scaling for large argument.
109+
// Use large-argument asymptotic expansion.
110+
111+
const T prefix { exp(-x) * pow(x, x - T { 5, -1 }) };
109112

110-
result *= r;
113+
result = prefix * detail::tgamma_series_expansion_asymp(one / x);
111114
}
112115
}
113116
}

‎test/test_lgamma.cpp‎

Lines changed: 110 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -53,17 +53,31 @@ namespace local
5353

5454
auto result_is_ok = bool { };
5555

56+
NumericType delta { };
57+
5658
if(b == static_cast<NumericType>(0))
5759
{
58-
result_is_ok = (fabs(a - b) < tol);
60+
delta = fabs(a - b); // LCOV_EXCL_LINE
61+
62+
result_is_ok = (delta < tol); // LCOV_EXCL_LINE
5963
}
6064
else
6165
{
62-
const auto delta = fabs(1 - (a / b));
66+
delta = fabs(1 - (a / b));
6367

6468
result_is_ok = (delta < tol);
6569
}
6670

71+
// LCOV_EXCL_START
72+
if (!result_is_ok)
73+
{
74+
std::cerr << std::setprecision(std::numeric_limits<NumericType>::digits10) << "a: " << a
75+
<< "\nb: " << b
76+
<< "\ndelta: " << delta
77+
<< "\ntol: " << tol << std::endl;
78+
}
79+
// LCOV_EXCL_STOP
80+
6781
return result_is_ok;
6882
}
6983

@@ -343,6 +357,76 @@ namespace local
343357

344358
return result_is_ok;
345359
}
360+
361+
auto test_lgamma_128(const int tol_factor) -> bool
362+
{
363+
using decimal_type = boost::decimal::decimal128;
364+
365+
using str_ctrl_array_type = std::array<const char*, 21U>;
366+
367+
const str_ctrl_array_type ctrl_strings =
368+
{{
369+
// Table[N[Log[Gamma[(100 n + 10 n + 1)/100]], 36], {n, 0, 20, 1}]
370+
"4.59947987804202172251394541100874809",
371+
"-0.0540386340818523935917550731681660590",
372+
"0.102418994503958632699253052937769400",
373+
"0.997464457272922372053206167365619618",
374+
"2.32975308729902926366841147898554568",
375+
"3.97393485962892204454162289923259731",
376+
"5.86078226284320941736299492331704683",
377+
"7.94629710737608673894522918391574878",
378+
"10.2000180598708079077541397082157801",
379+
"12.5995970196581001223988397360238569",
380+
"15.1279348557753769796480417309555140",
381+
"17.7715247207270252174494824518277843",
382+
"20.5194267289921636545853277538118495",
383+
"23.3625991972628192905542283017434866",
384+
"26.2934437604612886905683626700964367",
385+
"29.3054851520909703460851836031022652",
386+
"32.3931392288932864648032352458271232",
387+
"35.5515407902467593660096806067295711",
388+
"38.7764130861225208432040187016879672",
389+
"42.0639671128620105436453477946728445",
390+
"45.4108226536777051814945280596645578",
391+
}};
392+
393+
std::array<decimal_type, std::tuple_size<str_ctrl_array_type>::value> tg_values { };
394+
std::array<decimal_type, std::tuple_size<str_ctrl_array_type>::value> ctrl_values { };
395+
396+
int nx { 0 };
397+
398+
bool result_is_ok { true };
399+
400+
const decimal_type my_tol { std::numeric_limits<decimal_type>::epsilon() * static_cast<decimal_type>(tol_factor) };
401+
402+
for(auto i = static_cast<std::size_t>(UINT8_C(0)); i < std::tuple_size<str_ctrl_array_type>::value; ++i)
403+
{
404+
const decimal_type x_arg =
405+
decimal_type
406+
{
407+
decimal_type { 1, 2 } * nx
408+
+ decimal_type { 1, 1 } * nx
409+
+ 1
410+
}
411+
/ decimal_type { 1, 2 };
412+
413+
++nx;
414+
415+
tg_values[i] = lgamma(x_arg);
416+
417+
static_cast<void>
418+
(
419+
from_chars(ctrl_strings[i], ctrl_strings[i] + std::strlen(ctrl_strings[i]), ctrl_values[i])
420+
);
421+
422+
const auto result_lgamma_is_ok = is_close_fraction(tg_values[i], ctrl_values[i], my_tol);
423+
424+
result_is_ok = (result_lgamma_is_ok && result_is_ok);
425+
}
426+
427+
return result_is_ok;
428+
}
429+
346430
} // namespace local
347431

348432
auto main() -> int
@@ -353,7 +437,7 @@ auto main() -> int
353437
using decimal_type = boost::decimal::decimal64;
354438
using float_type = double;
355439

356-
const auto result_special_issue385_is_ok = local::test_special_issue385<decimal_type, float_type>(4096);
440+
const auto result_special_issue385_is_ok = local::test_special_issue385<decimal_type, float_type>(4096);
357441

358442
BOOST_TEST(result_special_issue385_is_ok);
359443

@@ -397,15 +481,26 @@ auto main() -> int
397481
using decimal_type = boost::decimal::decimal64;
398482
using float_type = double;
399483

400-
const auto result_tgamma_is_ok = local::test_lgamma<decimal_type, float_type>(4096, 2.1L, 123.4L);
484+
const auto result_tgamma_is_ok = local::test_lgamma<decimal_type, float_type>(8192, 0.1L, 0.9L);
401485

402486
BOOST_TEST(result_tgamma_is_ok);
403487

404488
result_is_ok = (result_tgamma_is_ok && result_is_ok);
405489
}
406490

407491
{
408-
const auto result_neg32_is_ok = local::test_lgamma_neg32(1024);
492+
using decimal_type = boost::decimal::decimal64;
493+
using float_type = double;
494+
495+
const auto result_tgamma_is_ok = local::test_lgamma<decimal_type, float_type>(8192, 1.1L, 123.4L);
496+
497+
BOOST_TEST(result_tgamma_is_ok);
498+
499+
result_is_ok = (result_tgamma_is_ok && result_is_ok);
500+
}
501+
502+
{
503+
const auto result_neg32_is_ok = local::test_lgamma_neg32(2048);
409504

410505
BOOST_TEST(result_neg32_is_ok);
411506

@@ -423,6 +518,16 @@ auto main() -> int
423518
result_is_ok = (result_edge_is_ok && result_is_ok);
424519
}
425520

521+
{
522+
// TODO(ckormanyos) Can we reduce the tolerance on lgamma()-128?
523+
// Can the approximation be done a bit better (Lanczos from Math?).
524+
const auto result_lgamma128_is_ok = local::test_lgamma_128(16'000'000);
525+
526+
BOOST_TEST(result_lgamma128_is_ok);
527+
528+
result_is_ok = (result_lgamma128_is_ok && result_is_ok);
529+
}
530+
426531
result_is_ok = ((boost::report_errors() == 0) && result_is_ok);
427532

428533
return (result_is_ok ? 0 : -1);

0 commit comments

Comments
 (0)