Hi @codeplea,

I made use of this function in the MRtrix3 software to deal with compilation on systems where our primary mathematics library dependency, Eigen, is missing the corresponding function. Its use is in converting t-values / F-values to Z-scores prior to statistical enhancement within a generalised statistical inference framework that operates on different forms of neuroimaging data.

This code was recently flagged by Thread Sanitiser. Our code generates a non-parametric null distribution by evaluating random permutations of the data across multiple threads. The issue is that std::lgamma() is not thread-safe. It evaluates the logarithm of the absolute value of the gamma function. In order to not lose information where the gamma function is negative, it stores the sign of the gamma function in a static external variable. There may be in some environments an alternative implementation lgamma_r(), which takes as a second input a pointer to an integer where that sign is instead written, which can therefore be used in a thread safe manner.

Can I suggest:

• Explicitly test that both inputs a and b to the regularised incomplete beta function are positive (ie. I presume you're not handling analytic continuation)?
• Add a note in the documentation about relevance of threading race conditions (ie. if I'm understanding correctly, these should be irrelevant since the inputs to std::lgamma() are always positive and therefore the sign of the gamma function is always positive)?

Cheers
Rob

Why make such a choice to use a return 1.0 / 0.0;

As the title, I need it in c++, but I not very good to write it...