@@ -45,23 +45,11 @@ BOOST_MATH_GPU_ENABLED RealType cdf_imp(const cauchy_distribution<RealType, Poli
4545 //
4646 // This calculates the cdf of the Cauchy distribution and/or its complement.
4747 //
48- // The usual formula for the Cauchy cdf is:
48+ // This implementation uses the formula
4949 //
50- // cdf = 0.5 + atan( x)/pi
50+ // cdf = atan2(1, - x)/pi
5151 //
52- // But that suffers from cancellation error as x -> -INF.
53- //
54- // Recall that for x < 0:
55- //
56- // atan(x) = -pi/2 - atan(1/x)
57- //
58- // Substituting into the above we get:
59- //
60- // CDF = -atan(1/x)/pi ; x < 0
61- //
62- // So the procedure is to calculate the cdf for -fabs(x)
63- // using the above formula, and then subtract from 1 when required
64- // to get the result.
52+ // where x is the standardized (i.e. shifted and scaled) domain variable.
6553 //
6654 BOOST_MATH_STD_USING // for ADL of std functions
6755 constexpr auto function = " boost::math::cdf(cauchy<%1%>&, %1%)" ;
@@ -99,13 +87,8 @@ BOOST_MATH_GPU_ENABLED RealType cdf_imp(const cauchy_distribution<RealType, Poli
9987 { // Catches x == NaN
10088 return result;
10189 }
102- RealType mx = -fabs ((x - location) / scale); // scale is > 0
103- if (mx > -tools::epsilon<RealType>() / 8 )
104- { // special case first: x extremely close to location.
105- return static_cast <RealType>(0 .5f );
106- }
107- result = -atan (1 / mx) / constants::pi<RealType>();
108- return (((x > location) != complement) ? 1 - result : result);
90+ RealType x_std = static_cast <RealType>((complement) ? 1 : -1 )*(x - location) / scale;
91+ return atan2 (static_cast <RealType>(1 ), x_std) / constants::pi<RealType>();
10992} // cdf
11093
11194template <class RealType , class Policy >
0 commit comments