https://svn.boost.org/trac10/ticket/9104 <p> The way ellint_e_imp reduces its input angle into the range |phi| &lt;= pi/2 has a bug when executing in double precision mode on intel chips (as opposed to the default double extended-precision mode) that results in errors as great as 100%. </p> <p> The relevant lines are at boost/math/special_functions/ellint_2.hpp: </p> <pre class="wiki">T rphi = boost::math::tools::fmod_workaround(phi, T(constants::pi&lt;T&gt;() / 2)); T m = floor((2 * phi) / constants::pi&lt;T&gt;()); </pre><p> Compiled with optimizations on both g++ and clang++, fmod_workaround() essentially becomes a "fprem" instruction, and (2 * phi) / pi() becomes an "fdivr" instruction. "fprem" ignores the precision mode and always produces exact results. "fdivr," however, gives less accurate results in plain double precision mode. This leads to the following error when, for example, phi = M_PI (slightly less than boost's more accurate pi()): </p> <p> In plain double precision, rphi is computed as slightly less than pi()/2 (correct since M_PI &lt; pi()). However, m = 2 (incorrect, rounded result). The end result is m * pi()/2 + rphi is close to 3pi/2 instead of M_PI, causing the computed integral to be 50% too large. In extended precision, the "fdivr" is accurate enough to yield m = 1 correctly. </p> <p> There are two ways I can see of avoiding this: </p> <pre class="wiki">T rphi = boost::math::tools::fmod_workaround(phi, T(constants::pi&lt;T&gt;() / 2)); T m = (phi - rphi) / T(constants::pi&lt;T&gt;() / 2); </pre><p> and the less verbose: </p> <pre class="wiki">T m = floor(phi / T(constants::pi&lt;T&gt;() / 2)); T rphi = phi - m * T(constants::pi&lt;T&gt;() / 2); </pre><p> Both would ensure rphi and m are consistent with the original phi. Note: in the second formula for m, T(constants::pi&lt;T&gt;() / 2) is folded into a constant that is exactly pi()/2, whereas the floor() argument in the original version generated extra code to double phi. </p> en-us /htdocs/site/boost.png https://svn.boost.org/trac10/ticket/9104 Trac 1.4.3 John Maddock Sun, 15 Sep 2013 16:56:44 GMT https://svn.boost.org/trac10/ticket/9104#comment:1 https://svn.boost.org/trac10/ticket/9104#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">fixed</span> </li> </ul> <p> (In <a class="changeset" href="https://svn.boost.org/trac10/changeset/85678" title="Suppress warning in fraction.hpp. Fix internal consistency of argument ...">[85678]</a>) Suppress warning in fraction.hpp. Fix internal consistency of argument reduction in elliptic integrals when the argument is very close to a multiple of PI/2. Fixes <a class="closed ticket" href="https://svn.boost.org/trac10/ticket/9107" title="#9107: Bugs: unused parameter warning in fraction with clang (closed: fixed)">#9107</a>. Fixes <a class="closed ticket" href="https://svn.boost.org/trac10/ticket/9104" title="#9104: Bugs: boost::math::ellint_2 bug in x86_64 double precision (closed: fixed)">#9104</a>. </p> Ticket John Maddock Sat, 28 Sep 2013 16:19:56 GMT