Scalar lgamma_stirling and lbeta truncate log input to integer

[rok-cesnovar]

Description

https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/lbeta.hpp does not include <cmath> and does not declare using std::log anywhere. Which probably means that that would use stan::math::log. But, given that only declare stan::math::log(int) that means that the input gets truncated to integer. Luckily this doesnt effect vectorized versions and eigen versions.

We have the same issue in https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/lgamma_stirling.hpp where there is a missing using statement.

This came up when trying to fix #1692

cc: @martinmodrak I think you were having issue with these two functions in the past if I remember correctly, but cant find the thread/PR right now. Hopefully the fix of this wil help in any way.

Will fix in #1702

Current Version:

v3.1.0

[martinmodrak]

Thanks for looking into this, this is indeed my code (the relevant PR is #1612). Nevertheless the code is passing some quite strict tests, including with non-integer numbers, so I don't find it likely that the log is truncated. But rather safe than sorry I guess :-)

[rok-cesnovar]

Yeah, I guess if it passes some strict tests, then log is not being truncated, luckily :)

Something like this was happening on an experimental branch where we were using exp() in a custom function in the stan::math namespace and the input got truncated and days were wasted figuring this one :)

There is most likely a global using std::log statement somewhere or something like that which is saving us in this case.

[bob-carpenter]

Are they passing these tests on all compilers? What can happen is that top-level ::log can sneak in from some headers or in some compilers. Is there a way to figure out which log is getting resolved?

Also, I think we should change all of our functions specified for int and/or double and rewrite them to take a template parameter with an enable-if for is-arithmetic. Then if there's a more specific library function it'll match, but otherwise the templated version in the math lib will be the best match.

This is all in theory---I haven't tried it!

[rok-cesnovar]

Ok, got to the bottom of this. First, the most important, there is no user-facing bug in terms of precision and performance but its a bit too dangerous at least for the latter.

In lgamma_stirling we are saved by the global namespace ::log that is included by <cmath>. If instead of <cmath> we would include #include <stan/math/prim/fun/log.hpp> this would use the vectorized version of log(). If that one would not exist (or was not included pre-flatten) it would use stan::math::log(int). I went and verified all these claims on Ubuntu with g++ and clang++.

Including stan/math/prim/fun/log.hpp instead of cmath can happen pretty quickly and the question is do we lose anything if log(3.3) calls the vectorized call shown below?
Tagging @t4c1 @SteveBronder @andrjohns who would know more on this.

template <typename T, typename = require_not_eigen_vt<std::is_arithmetic, T>>
inline auto log(const T& x) {
  return apply_scalar_unary<log_fun, T>::apply(x);
}

In lbeta we are saved by the same global namespace ::log that is included by <cmath> which is included by lgamma_stirling header that is included.There is no simple way to not include cmath here so I didnt go and try that.

The same things happen for other functions that dont have their own declarations for primitive scalars.

#include <stan/math/prim/fun/cos.hpp>
#include <iostream>

int main() {
    std::cout << stan::math::cos(1.1) << std::endl;
}

This uses https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/cos.hpp#L35

Is that fine or is it any slower?

If its slower, we need to add definitions for scalars for these functions.
If its fine, we can remove all the stan::math scalar functions like the ones defined for cbrt for example https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/cbrt.hpp

[bob-carpenter]

On Feb 13, 2020, at 7:11 AM, Rok Češnovar ***@***.***> wrote: Ok, got to the bottom of this. First, the most important, there is no user-facing bug in terms of precision and performance but its a bit too dangerous at least for the latter. In lgamma_stirling we are saved by the global namespace ::log that is included by <cmath>. If instead of <cmath> we would include #include <stan/math/prim/fun/log.hpp> this would use the vectorized version of log(). If that one would not exist (or was not included pre-flatten) it would use stan::math::log(int). I went and verified all these claims on Ubuntu with g++ and clang++.

That's great. I thought some of these behaved differently w.r.t putting the C math functions in the top-level namespace.

Including stan/math/prim/fun/log.hpp instead of cmath can happen pretty quickly and the question is do we lose anything if log(3.3) calls the vectorized call shown below? Tagging @t4c1 @SteveBronder @andrjohns who would know more on this.

The code generated should be efficient and produce the correct results. apply_scalar_unary has an overload for scalars and the log_fun defines a static function so there's no constructor overhead.

template <typename T, typename = require_not_eigen_vt<std::is_arithmetic, T>> inline auto log(const T& x) { return apply_scalar_unary<log_fun, T>::apply (x); } In lbeta we are saved by the same global namespace ::log that is included by <cmath> which is included by lgamma_stirling header that is included.There is no simple way to not include cmath here so I didnt go and try that. The same things happen for other functions that dont have their own declarations for primitive scalars. #include <stan/math/prim/fun/cos.hpp> # include <iostream> int main () { std::cout << stan::math::cos(1.1 ) << std::endl; } This uses https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/cos.hpp#L35 Is that fine or is it any slower?

I don't see why it'd be slower. Compilers are really good at inlining trivial calls.

If its slower, we need to add definitions for scalars for these functions. If its fine, we can remove all the stan::math scalar functions like the ones defined for cbrt for example https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/cbrt.hpp

I don't think it's slower. Those scalar functions all got added because things wouldn't compile without them. That may have all been laid in before the vectorization, though. We've been struggling with the C math namespaces from the get go.

[rok-cesnovar]

I wasn't sure how apply_scalar_unary worked and wanted to double check. In that case we don't need any of the scalar function definitions wherever we have the vectorized definitions. Will clean that up in the updated PR.

[bob-carpenter]

Thanks so much for digging into this. All this cleanup and getting to the bottom of things that've been thorns in our side for years is making me so happy!

…

On Feb 13, 2020, at 9:06 AM, Rok Češnovar ***@***.***> wrote: I wasn't sure how apply_scalar_unary worked and wanted to double check. In that case we don't need any of the scalar function definitions wherever we have the vectorized definitions. Will clean that up in the updated PR. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[rok-cesnovar]

No problem. I strangely enjoy solving these kind of puzzles.

[martinmodrak]

Thanks for checking this and for making the codebase better.

[rok-cesnovar]

Fixed by #1711