I would like to use an implementation of BesselK that accepts the order (v) as real (parameter).

I found this thread Vignette for interfacing with external C++ code via rstan, that includes exactly the BesselK example.

I was just wandering if from 2016 there was an alternative method to include it in rstan (from R C++ code for example)

(The goal is to enable coding of distributions such as SICHEL https://www.rdocumentation.org/packages/gamlss.dist/versions/5.1-1/topics/SICHEL)

That thread also put forward concerns about it’s numerical stability for differentiation, any reccomendations?

Thanks