diff options
Diffstat (limited to 'boost/math/special_functions/detail/bessel_j1.hpp')
| -rw-r--r-- | boost/math/special_functions/detail/bessel_j1.hpp | 17 |
1 files changed, 14 insertions, 3 deletions
diff --git a/boost/math/special_functions/detail/bessel_j1.hpp b/boost/math/special_functions/detail/bessel_j1.hpp index 09d862c240..91ecd2832d 100644 --- a/boost/math/special_functions/detail/bessel_j1.hpp +++ b/boost/math/special_functions/detail/bessel_j1.hpp @@ -166,13 +166,24 @@ T bessel_j1(T x) { T y = 8 / w; T y2 = y * y; - T z = w - 0.75f * pi<T>(); BOOST_ASSERT(sizeof(PC) == sizeof(QC)); BOOST_ASSERT(sizeof(PS) == sizeof(QS)); rc = evaluate_rational(PC, QC, y2); rs = evaluate_rational(PS, QS, y2); - factor = sqrt(2 / (w * pi<T>())); - value = factor * (rc * cos(z) - y * rs * sin(z)); + factor = 1 / (sqrt(w) * constants::root_pi<T>()); + // + // What follows is really just: + // + // T z = w - 0.75f * pi<T>(); + // value = factor * (rc * cos(z) - y * rs * sin(z)); + // + // but using the sin/cos addition rules plus constants + // for the values of sin/cos of 3PI/4 which then cancel + // out with corresponding terms in "factor". + // + T sx = sin(x); + T cx = cos(x); + value = factor * (rc * (sx - cx) + y * rs * (sx + cx)); } if (x < 0) |