Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bugfix/issue 1250 lgamma #1255

Merged
merged 47 commits into from Sep 14, 2019
Merged
Changes from all commits
Commits
Show all changes
47 commits
Select commit Hold shift + click to select a range
5c7407d
switch to thread safe and non-blocking version of lgamma from boost
wds15 May 30, 2019
1d92c2a
change back test expectation and fix a doxygen error
wds15 May 30, 2019
4651f94
add note in doc why we avoid use of the std version of lgamma
wds15 May 30, 2019
05f3ecc
align behavior at 0 for lgamma
wds15 May 30, 2019
8fd0f77
cpplint
wds15 May 30, 2019
6b7586b
fix cpplint
wds15 May 30, 2019
9b738d2
avoid extreme test leading to NaN
wds15 May 30, 2019
1e663cd
help my memory
wds15 May 30, 2019
cce935f
reenable large number test for binomial_coefficient_log_test
wds15 May 31, 2019
b626924
lower large number test
wds15 May 31, 2019
031c206
adjust test to new lgamma precision targets
wds15 May 31, 2019
b42d27a
Merge remote-tracking branch 'origin/develop' into bugfix/issue-1250-…
wds15 Jun 20, 2019
cfe1903
switch to lgamma_r which is the reentrant c implementation
wds15 Jun 20, 2019
8d3ec7d
use math.h which makes lgamma_r hopefully available on mingw gcc
weberse2 Jun 21, 2019
a7aa430
Merge branch 'bugfix/issue-1250-lgamma' of https://github.com/stan-de…
wds15 Jun 21, 2019
f01d617
on non-POSIX systems like Windows use boost of lgamma_r
wds15 Jun 21, 2019
a2cc131
make things work on Windows
Jun 22, 2019
1ff4299
Merge branch 'bugfix/issue-1250-lgamma' of https://github.com/stan-de…
wds15 Jun 22, 2019
09e5dff
[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…
stan-buildbot Jun 22, 2019
83a37c9
make use of _REENTRANT functions
wds15 Jun 22, 2019
64a4ffa
make sure that reentrant functions are not used on migw32
wds15 Jun 22, 2019
073a0c0
Merge branch 'bugfix/issue-1250-lgamma' of https://github.com/stan-de…
wds15 Jun 22, 2019
6521f3f
enhance test-math-dependencies to detect wrong use of std::lgamma whi…
wds15 Jul 3, 2019
d5e8474
remove std::lgamma again
wds15 Jul 3, 2019
df6b2e5
Merge remote-tracking branch 'origin/develop' into bugfix/issue-1250-…
wds15 Jul 3, 2019
af13453
avoid windows config error of Jenkins
wds15 Jul 4, 2019
0082361
revert to develop Jenkinsfile due to syntax errors
wds15 Jul 8, 2019
9f320b9
Merge branch 'bugfix/issue-1250-lgamma' of https://github.com/stan-de…
weberse2 Jul 15, 2019
f0dcc6d
spelling
weberse2 Jul 15, 2019
bb3ccb3
always define _REENTRANT and simplify branching logic to only check f…
wds15 Jul 22, 2019
32d5fbe
Merge branch 'bugfix/issue-1250-lgamma' of https://github.com/stan-de…
wds15 Jul 22, 2019
f26315a
Merge remote-tracking branch 'origin/develop' into bugfix/issue-1250-…
wds15 Jul 25, 2019
1ef61d4
update README.md and add error note about _REENTRANT in case users at…
wds15 Jul 31, 2019
7acf35d
Merge remote-tracking branch 'origin/develop' into bugfix/issue-1250-…
wds15 Jul 31, 2019
714058a
cpplint
wds15 Jul 31, 2019
664c698
Merge remote-tracking branch 'origin/develop' into bugfix/issue-1250-…
wds15 Aug 12, 2019
4856838
reenable large number test
wds15 Aug 17, 2019
3429945
add test to ensure we have the right boost lgamma
wds15 Aug 17, 2019
419ad92
add some large number tests for boost math functions
wds15 Aug 17, 2019
90bb66d
backport commit https://github.com/boostorg/math/commit/fd07b121a619e…
wds15 Aug 17, 2019
1168f93
add refs to lgamma discussions
wds15 Aug 17, 2019
ff9c3d4
test only boost functions which actually use boost_policy_t
wds15 Aug 17, 2019
57f8cc9
spelling
wds15 Aug 17, 2019
5c1b85e
address review comments
wds15 Sep 8, 2019
bcf6e7b
Merge remote-tracking branch 'origin/develop' into bugfix/issue-1250-…
wds15 Sep 8, 2019
6a1af0e
Merge branch 'develop' into bugfix/issue-1250-lgamma
rok-cesnovar Sep 9, 2019
cb96cad
fixed merge
rok-cesnovar Sep 9, 2019
File filter...
Filter file types
Jump to…
Jump to file or symbol
Failed to load files and symbols.

Always

Just for now

@@ -38,7 +38,7 @@ If this is in the file `/path/to/foo/foo.cpp`, then you can compile and run this

```
> cd /path/to/foo
> clang++ -std=c++1y -I /path/to/stan-math -I /path/to/Eigen -I /path/to/boost -I /path/to/sundials foo.cpp
> clang++ -std=c++1y -I /path/to/stan-math -I /path/to/Eigen -I /path/to/boost -I /path/to/sundials -D_REENTRANT foo.cpp
> ./a.out
log normal(1 | 2, 3)=-2.07311
```
@@ -53,7 +53,7 @@ The `-I` includes provide paths pointing to the four necessary includes:
Note that the paths should *not* include the final directories `stan`, `Eigen`, or `boost` on the paths. An example of a real instantiation:

```
clang++ -std=c++1y -I ~/stan-dev/math -I ~/stan-dev/math/lib/eigen_3.3.3/ -I ~/stan-dev/math/lib/boost_1.69.0/ -I ~/stan-dev/math/lib/sundials_4.1.0/include foo.cpp
clang++ -std=c++1y -I ~/stan-dev/math -I ~/stan-dev/math/lib/eigen_3.3.3/ -I ~/stan-dev/math/lib/boost_1.69.0/ -I ~/stan-dev/math/lib/sundials_4.1.0/include -D_REENTRANT foo.cpp
```

The following directories all exist below the links given to `-I`: `~/stan-dev/math/stan` and `~/stan-dev/math/lib/eigen_3.3.3/Eigen` and `~stan-dev/math/lib/boost_1.69.0/boost` and `~stan-dev/math/lib/sundials_4.1.0/include`.
@@ -51,6 +51,15 @@ inline double lanczos13m53::lanczos_sum<double>(const double& x)
static_cast<double>(23531376880.41075968857200767445163675473L),
static_cast<double>(0u)
};

static const double lim = 4.31965e+25; // By experiment, the largest x for which the SSE2 code does not go bad.

if (x > lim)
{
double z = 1 / x;
return ((((((((((((coeff[24] * z + coeff[22]) * z + coeff[20]) * z + coeff[18]) * z + coeff[16]) * z + coeff[14]) * z + coeff[12]) * z + coeff[10]) * z + coeff[8]) * z + coeff[6]) * z + coeff[4]) * z + coeff[2]) * z + coeff[0]) / ((((((((((((coeff[25] * z + coeff[23]) * z + coeff[21]) * z + coeff[19]) * z + coeff[17]) * z + coeff[15]) * z + coeff[13]) * z + coeff[11]) * z + coeff[9]) * z + coeff[7]) * z + coeff[5]) * z + coeff[3]) * z + coeff[1]);
}

__m128d vx = _mm_load1_pd(&x);
__m128d sum_even = _mm_load_pd(coeff);
__m128d sum_odd = _mm_load_pd(coeff+2);
@@ -136,6 +145,15 @@ inline double lanczos13m53::lanczos_sum_expG_scaled<double>(const double& x)
static_cast<double>(56906521.91347156388090791033559122686859L),
static_cast<double>(0u)
};

static const double lim = 4.76886e+25; // By experiment, the largest x for which the SSE2 code does not go bad.

if (x > lim)
{
double z = 1 / x;
return ((((((((((((coeff[24] * z + coeff[22]) * z + coeff[20]) * z + coeff[18]) * z + coeff[16]) * z + coeff[14]) * z + coeff[12]) * z + coeff[10]) * z + coeff[8]) * z + coeff[6]) * z + coeff[4]) * z + coeff[2]) * z + coeff[0]) / ((((((((((((coeff[25] * z + coeff[23]) * z + coeff[21]) * z + coeff[19]) * z + coeff[17]) * z + coeff[15]) * z + coeff[13]) * z + coeff[11]) * z + coeff[9]) * z + coeff[7]) * z + coeff[5]) * z + coeff[3]) * z + coeff[1]);
}

__m128d vx = _mm_load1_pd(&x);
__m128d sum_even = _mm_load_pd(coeff);
__m128d sum_odd = _mm_load_pd(coeff+2);
@@ -277,7 +277,11 @@ T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
T zgh = static_cast<T>(z + Lanczos::g() - boost::math::constants::half<T>());
result = log(zgh) - 1;
result *= z - 0.5f;
result += log(Lanczos::lanczos_sum_expG_scaled(z));
//
// Only add on the lanczos sum part if we're going to need it:
//
if(result * tools::epsilon<T>() < 20)
result += log(Lanczos::lanczos_sum_expG_scaled(z));
}

if(sign)
@@ -133,6 +133,8 @@ ifeq ($(OS),Linux)
endif
endif

## makes reentrant version lgamma_r available from cmath
CXXFLAGS_OS += -D_REENTRANT

################################################################################
# Setup OpenCL
@@ -128,5 +128,9 @@ test-math-dependencies:
@$(call math_warnings, 'C++14', 'stan/math',\
'<boost/utility/enable_if.hpp>',\
'Replace \<boost/utility/enable_if.hpp\> with \<type_traits\>.')
## Check that we do not use non-reentrant safe functions from std
This conversation was marked as resolved by syclik

This comment has been minimized.

Copy link
@syclik

syclik Jul 15, 2019

Member

THANK YOU!

@$(call math_warnings, 'thread-reentrant-safe', 'stan/math',\
'std::lgamma',\
'Replace std::lgamma with reentrant safe version boost::math::lgamma or lgamma_r on supporting platforms.')
## Print all warnings
@$(call print_math_warnings)
@@ -19,7 +19,8 @@ namespace math {
using boost_policy_t = boost::math::policies::policy<
boost::math::policies::overflow_error<
boost::math::policies::errno_on_error>,
boost::math::policies::pole_error<boost::math::policies::errno_on_error>>;
boost::math::policies::pole_error<boost::math::policies::errno_on_error>,
boost::math::policies::promote_double<false> >;

} // namespace math
} // namespace stan
@@ -1,8 +1,34 @@
#ifndef STAN_MATH_PRIM_SCAL_FUN_LGAMMA_HPP
#define STAN_MATH_PRIM_SCAL_FUN_LGAMMA_HPP

#include <stan/math/prim/meta.hpp>
/*
* The lgamma implementation in stan-math is based on either the
* reentrant safe lgamma_r implementation from C or the
* boost::math::lgamma implementation. The reentrant safe lgamma_r
* implementation is preferred as it is faster when compared to the
* boost version. However, the reentrant safe lgamma_r C function is
* not available with MinGW compilers used on Windows. Therefore, the
* boost version is used on Windows with MinGW compilers as fall
* back. For details on the speed evaluations, please refer to
* https://github.com/stan-dev/math/pull/1255 .
*/
#if !__MINGW32__
// _REENTRANT must be defined during compilation to ensure that cmath
// exports the reentrant safe lgamma_r version.
#if !_REENTRANT
#error \
"stan-math requires _REENTRANT being defined during compilation" \
"to make lgamma_r available."
#endif
#include <cmath>
#else
// MinGW compilers on Windows do not provide the reentrant lgamma_r
// such that we fall back to boost whenever we are on MinGW.
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/scal/fun/boost_policy.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <limits>
#endif

namespace stan {
namespace math {
@@ -33,7 +59,16 @@ namespace math {
* @return natural logarithm of the gamma function applied to
* argument
*/
inline double lgamma(double x) { return std::lgamma(x); }
inline double lgamma(double x) {
#if !__MINGW32__
int sign = 1;
return ::lgamma_r(x, &sign);
#else
if (unlikely(x == 0.0))
return std::numeric_limits<double>::infinity();
return boost::math::lgamma(x, boost_policy_t());
#endif
}

/**
* Return the natural logarithm of the gamma function applied
@@ -43,7 +78,16 @@ inline double lgamma(double x) { return std::lgamma(x); }
* @return natural logarithm of the gamma function applied to
* argument
*/
inline double lgamma(int x) { return std::lgamma(x); }
inline double lgamma(int x) {
#if !__MINGW32__
int sign = 1;
return ::lgamma_r(x, &sign);
#else
if (unlikely(x == 0.0))
return std::numeric_limits<double>::infinity();
return boost::math::lgamma(x, boost_policy_t());
#endif
}

} // namespace math
} // namespace stan
@@ -2,6 +2,7 @@
#include <boost/math/special_functions/fpclassify.hpp>
#include <gtest/gtest.h>
#include <limits>
#include <cmath>

template <typename T_N, typename T_n>
void test_binom_coefficient(const T_N& N, const T_n& n) {
@@ -28,6 +29,7 @@ TEST(MathFunctions, binomial_coefficient_log) {

test_binom_coefficient(1e9, 1e5);
test_binom_coefficient(1e50, 1e45);
test_binom_coefficient(1e20, 1e15);
}

TEST(MathFunctions, binomial_coefficient_log_nan) {
@@ -15,4 +15,6 @@ TEST(MathFunctions, digamma_nan) {
EXPECT_PRED1(boost::math::isnan<double>, stan::math::digamma(nan));

EXPECT_PRED1(boost::math::isnan<double>, stan::math::digamma(-1));

EXPECT_PRED1(boost::math::isnormal<double>, stan::math::digamma(1.0E50));
}
@@ -12,4 +12,6 @@ TEST(MathFunctions, falling_factorial) {
EXPECT_FLOAT_EQ(120, falling_factorial(5.0, 4));
EXPECT_THROW(falling_factorial(1, -4), std::domain_error);
EXPECT_THROW(falling_factorial(nan, 1), std::domain_error);
// see comments in test/unit/math/prim/scal/fun/lgamma_test.cpp
EXPECT_PRED1(boost::math::isnormal<double>, falling_factorial(1.0E30, 5));
}
@@ -53,7 +53,7 @@ TEST(MathFunctions, inc_beta_dda) {
EXPECT_FLOAT_EQ(0.0, inc_beta_dda(large_a, large_b, small_z, digamma(large_a),
digamma(large_a + large_b)))
<< "reasonable values for a, b, x";
EXPECT_FLOAT_EQ(-4.0856207e-14,
EXPECT_FLOAT_EQ(-4.2743586e-14,
inc_beta_dda(large_a, large_b, mid_z, digamma(large_a),
digamma(large_a + large_b)))
<< "reasonable values for a, b, x";
@@ -12,4 +12,17 @@ TEST(MathFunctions, lgamma_nan) {
double nan = std::numeric_limits<double>::quiet_NaN();
EXPECT_PRED1(boost::math::isnan<double>, stan::math::lgamma(nan));
EXPECT_PRED1(boost::math::isinf<double>, stan::math::lgamma(0));
// up to boost 1.70.0 the boost::math::lgamma return NaN for large
// arguments (large is greater than sqrt(largest double of 1E308)
// when used with the stan::math::boost_policy_t which avoids
// propagation of input arguments from double to long double
// internally to boost::math::lgamma. See
// https://github.com/boostorg/math/issues/242 for this. The
// stan::math::lgamma implementation is based on the
// boost::math::lgamma only when MinGW compilers are used. Thus, to
// ensure that boost::math::lgamma contains the needed bugfixes we
// test here specifically the boost::math::lgamma by testing for a
// finite return for a large argument.
EXPECT_PRED1(boost::math::isnormal<double>,
boost::math::lgamma(1.0E50, stan::math::boost_policy_t()));
}
@@ -12,4 +12,6 @@ TEST(MathFunctions, rising_factorial) {
EXPECT_FLOAT_EQ(360, rising_factorial(3.0, 4));
EXPECT_THROW(rising_factorial(1, -4), std::domain_error);
EXPECT_THROW(rising_factorial(nan, 1), std::domain_error);
// see comments in test/unit/math/prim/scal/fun/lgamma_test.cpp
EXPECT_PRED1(boost::math::isnormal<double>, rising_factorial(1.0E30, 5));
}
ProTip! Use n and p to navigate between commits in a pull request.
You can’t perform that action at this time.