Skip to content

Use <cmath> for lgamma - #419

Closed
MichaelChirico wants to merge 2 commits into
StochasticTree:mainfrom
MichaelChirico:std-lgamma
Closed

MichaelChirico wants to merge 2 commits into
StochasticTree:mainfrom
MichaelChirico:std-lgamma

Conversation

@MichaelChirico

Copy link
Copy Markdown
Contributor

This avoids any boost headaches and should be equivalent (up to machine precision).

Written by Gemini while pursuing a crash. If boost is involved in the crash, it's a bit uglier to debug vs. if std::lgamma is used. It's also consistent with nearby usages of std::log().

OTOH, AIUI boost::math::lgamma is slightly more accurate, so perhaps the choice was intentional. E.g. Gemini offers the following:


Where Boost is slightly more accurate:

  1. Peak ULP Error:

    • boost::math::lgamma: Maintains a peak error strictly bounded under $\approx 0.6 - 0.9$ ULP across the positive real line. Boost generates its Lanczos expansion coefficients using 128-bit quad-precision arithmetic and high-degree polynomials.
    • Standard std::lgamma (glibc / libm): Peak error typically hovers around $1.0 - 2.0$ ULPs in worst-case regions.
  2. Near Roots ($x \approx 1$ and $x \approx 2$):

    • Since $\ln(\Gamma(1)) = 0$ and $\ln(\Gamma(2)) = 0$, evaluating $\ln(\Gamma(x))$ near $x = 1$ or $x = 2$ suffers from catastrophic cancellation in standard polynomial approximations.
    • Boost.Math includes dedicated Taylor/Padé expansions around $x = 1$ and $x = 2$ (using log1p variants) to preserve 15–16 significant digits across the zero-crossings.

Does it matter for stochtree?

  • Both implementations are accurate to 15–16 decimal digits ($\approx 10^{-16}$).
  • In stochtree, lgamma(a) is called on shape parameters ($a \ge 2.0$) to compute log-marginal likelihoods. A difference of $&lt; 1$ ULP ($10^{-16}$) is negligible compared to standard floating-point addition order.

@andrewherren

Copy link
Copy Markdown
Collaborator

Thanks for this PR! While I agree that making the debug traces easier to follow is an important goal for stochtree, this change doesn't remove our need for boost entirely and it introduces potential data races in multi-threading settings. Per cppreference:

The POSIX version of lgamma is not thread-safe: each execution of the function stores the sign of the gamma function of num in the static external variable signgam.

I am going to close it, though I will merge the rest of your PRs and also likely upgrade to Eigen 5

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants