Allow real argument to order of modified_bessel_x_kind

[bbbales2]

Newest Description

The derivatives are here https://dlmf.nist.gov/10.15 (functions themselves are here: https://dlmf.nist.gov/10.2)

The Stan functions are these: https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/bessel_first_kind.hpp and https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/bessel_second_kind.hpp

All the stuff below was the process of getting to this issue, so I'm leaving it, though it's probably confusing and should just be ignored.

New Description (Outdated)

Edit: I'm making a new description for this issue. The old description is below (so the comments make sense).

If we can make our var work with Boost Real (requirements here: https://www.boost.org/doc/libs/1_74_0/libs/math/doc/html/math_toolkit/real_concepts.html), then we can maybe autodiff the Boost modified_bessel_x_kind functions and allow real arguments to modified_bessel_first_kind and modified_bessel_second_kind.

Not 100% sure this would work, but it's worth experimenting with.

Old Description (Outdated)

The bessel functions in the math library only support integer orders: https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/bessel_first_kind.hpp https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/bessel_second_kind.hpp

I was doing doc for the log_modified_bessel_first_kind and it has a non-integer first argument: https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/log_modified_bessel_first_kind.hpp

Are we limited to integer arguments in the first couple things for autodiff reasons? Or should we switch those to doubles? Boost apparently supports non-integer arguments (though it uses a different algorithm for them).

Also there are std implementations now: https://en.cppreference.com/w/cpp/numeric/special_functions/cyl_bessel_k https://en.cppreference.com/w/cpp/numeric/special_functions/cyl_neumann

Is there any reason to not try to switch to non-integer order arguments and use to the std implementation? I feel like I remember a performance problem in threading caused by switching a fancy math function like this once, so I wanted to ask before I opened a pull.

@bgoodri, @wds15

Current Version:

v3.3.0

[bgoodri]

Autodiffing the Bessel functions in Boost with respect to the order does not work because our var types are not considered Real types by Boost. See

https://discourse.mc-stan.org/t/trigamma-from-boost/10546/15?u=bgoodri

The versions in the C++ standard were added with C++17, which we do not utilize yet.

The derivatives with respect to order are "known" in "closed" form utilizing a lot of hypergeometric functions that are now in Boost. See the CDF workbook linked from

https://blog.wolfram.com/2016/05/16/new-derivatives-of-the-bessel-functions-have-been-discovered-with-the-help-of-the-wolfram-language/

[bbbales2]

Oh I thought it was just the floats that were C++17. I see now it's everything there that is c++17.

@SteveBronder do you know about the Real types in Boost thing? Apparently last time it was tried our vars don't quite fit the bill. I think this is the docs page: https://www.boost.org/doc/libs/1_74_0/libs/math/doc/html/math_toolkit/real_concepts.html

Do you think those are requirements var can meet? Nothing was obviously impossible to me (though there are functions that don't have derivatives, and there's probably some hidden trickiness)

This is nice and low priority but if it looks like we can't do that then at least we can have that written down somewhere (and if it looks like we can do that, we can at least make a todo).

[SteveBronder]

I've never seen the real type in boost thing. But it looks like we satisfy the basic stuff, is it possible to support itrunc or iround?

[bbbales2]

Well we can certainly implement it in Math. Don't need to expose this at the language. It'll be Bad Business autodiffing through such a thing, but same is true for floor and ceil.

[bbbales2]

I tried to do the Boost reals thing and got segfaults after I implemented the few functions it takes to compile.

Then I Googled for the derivatives to do the autodiff manually and found this: https://dlmf.nist.gov/10.2, https://dlmf.nist.gov/10.15. It's probably easier to just code up the autodiff manually than do the Boost Real thing. Updating this again.

[bgoodri]

The derivatives with respect to the order in the DLMF are not that numerically well-behaved. If you were going to go that route, the ones in https://blog.wolfram.com/2016/05/16/new-derivatives-of-the-bessel-functions-have-been-discovered-with-the-help-of-the-wolfram-language/ are possibly better. But the inability of Boost to accept var types as Real numbers (when it accepts almost everything else) is a problem that goes beyond just the Bessel functions.

[bbbales2]

But the inability of Boost to accept var types as Real numbers

I got nervous when I saw the segfault.

It reminds me of how Chris got segfaults in the initial complex number stuff. Like, uninitialized vars just lead to segfaults cuz they don't have varis.

I didn't confirmed, but that's what I fear, and that would be a Big Problem (TM).