https://svn.boost.org/trac10/ticket/5433 <p> In the beta.hpp header for the beta function (i.e., *not* the beta distribution) around line 131, you have: </p> <blockquote> <p> if((c == a) &amp;&amp; (b &lt; tools::epsilon&lt;T&gt;())) </p> <blockquote> <p> return boost::math::tgamma(b, pol); </p> </blockquote> <p> else if((c == b) &amp;&amp; (a &lt; tools::epsilon&lt;T&gt;())) </p> <blockquote> <p> return boost::math::tgamma(a, pol); </p> </blockquote> </blockquote> <p> Don't you actually want if b&lt;epsilon, tgamma(a)? beta(a,b) reduces to (approximately) gamma(a)b^{-a for a&gt;&gt;b, which seems like the case for which you are trying to account?}</p> <p> Cheers, </p> <p> Eric </p> en-us /htdocs/site/boost.png https://svn.boost.org/trac10/ticket/5433 Trac 1.4.3 John Maddock Thu, 07 Apr 2011 10:08:18 GMT https://svn.boost.org/trac10/ticket/5433#comment:1 https://svn.boost.org/trac10/ticket/5433#comment:1 <ul> <li><strong>status</strong> <span class="trac-field-old">new</span> → <span class="trac-field-new">closed</span> </li> <li><strong>resolution</strong> → <span class="trac-field-new">invalid</span> </li> </ul> <p> I believe the code to be correct: </p> <p> if a &gt;&gt; b such that a + b == a (at the precision we're working to), then: </p> <p> gamma(a) / gamma(a+b) == 1 </p> <p> so: </p> <p> beta(a, b) = gamma(a) * gamma(b) / gamma(a+b) ~ gamma(b) ~ 1/b </p> Ticket Eric Butter <egbutter@…> Thu, 07 Apr 2011 15:37:02 GMT