Paul A. Bristow wrote:
This looks too obviously lifted from NR in C (or FORTRAN) to fit the Boost license.
Yes, as I mentioned, it's from a NR FORTRAN recipe. I don't know if the conversion from FORTRAN to C++ makes it a non-derivative work. The code certainly wouldn't work in poisson_distribution.hpp as is.
But the algorithm should still be usable? (But I am not qualified to judge).
If not, then at least the fix (in the original part of my post) should be used. Right now, using [random] poisson_distribution is slow and dangerous for any mean over 750. If it is changed to use ln instead of exp, it will only be slow. That is completely unrelated to the NR code later on in my post, which is a (serious) speed optimization.
ln gamma (and poisson) is available from the Math Toolkit, but this is a random requirement and you may value speed more than accuracy?
I don't see any way to use the Math Toolkit poisson distribution code to actually generate poisson random variables. As a side note, the pdf function in Math Toolkit may have a possible problem with exp(-mean) where k is smaller than max_factorial and mean is larger than ~750. One fix would simply be to use the same code for k>max_factorial in that case as well. Joel Eidsath