[Boost-commit] svn:boost r81755 - in trunk/boost/math: distributions special_functions/detail

Subject: [Boost-commit] svn:boost r81755 - in trunk/boost/math: distributions special_functions/detail
From: john_at_[hidden]
Date: 2012-12-07 07:49:44

Author: johnmaddock
Date: 2012-12-07 07:49:42 EST (Fri, 07 Dec 2012)
New Revision: 81755
URL: http://svn.boost.org/trac/boost/changeset/81755

Log:
Fix typo in fwd.hpp.
Improve accuracy of Bessel J and Y[0,1] for large z.
Text files modified: 
   trunk/boost/math/distributions/fwd.hpp                  |     2 +-                                      
   trunk/boost/math/special_functions/detail/bessel_j0.hpp |    15 +++++++++++++--                         
   trunk/boost/math/special_functions/detail/bessel_j1.hpp |    15 +++++++++++++--                         
   trunk/boost/math/special_functions/detail/bessel_y0.hpp |    15 +++++++++++++--                         
   trunk/boost/math/special_functions/detail/bessel_y1.hpp |    14 ++++++++++++--                          
   5 files changed, 52 insertions(+), 9 deletions(-)

Modified: trunk/boost/math/distributions/fwd.hpp
==============================================================================
--- trunk/boost/math/distributions/fwd.hpp	(original)
+++ trunk/boost/math/distributions/fwd.hpp	2012-12-07 07:49:42 EST (Fri, 07 Dec 2012)
@@ -121,7 +121,7 @@
    typedef boost::math::fisher_f_distribution<Type, Policy> fisher_f;\
    typedef boost::math::gamma_distribution<Type, Policy> gamma;\
    typedef boost::math::geometric_distribution<Type, Policy> geometric;\
-   typedef boost::math::hyper_geometric_distribution<Type, Policy> hyper_geometric;\
+   typedef boost::math::hypergeometric_distribution<Type, Policy> hypergeometric;\
    typedef boost::math::inverse_chi_squared_distribution<Type, Policy> inverse_chi_squared;\
    typedef boost::math::inverse_gaussian_distribution<Type, Policy> inverse_gaussian;\
    typedef boost::math::inverse_gamma_distribution<Type, Policy> inverse_gamma;\

Modified: trunk/boost/math/special_functions/detail/bessel_j0.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/bessel_j0.hpp	(original)
+++ trunk/boost/math/special_functions/detail/bessel_j0.hpp	2012-12-07 07:49:42 EST (Fri, 07 Dec 2012)
@@ -165,13 +165,24 @@
     {
         T y = 8 / x;
         T y2 = y * y;
-        T z = x - 0.25f * 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 / (x * pi<T>()));
-        value = factor * (rc * cos(z) - y * rs * sin(z));
+        //
+        // What follows is really just:
+        //
+        // T z = x - pi/4;
+        // value = factor * (rc * cos(z) - y * rs * sin(z));
+        //
+        // But using the addition formulae for sin and cos, plus
+        // the special values for sin/cos of pi/4.
+        //
+        T sx = sin(x);
+        T cx = cos(x);
+        value = factor * (rc * (cx * constants::one_div_root_two<T>() + sx * constants::half_root_two<T>()) 
+           - y * rs * (sx * constants::one_div_root_two<T>() - cx * constants::half_root_two<T>()));
     }
 
     return value;

Modified: trunk/boost/math/special_functions/detail/bessel_j1.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/bessel_j1.hpp	(original)
+++ trunk/boost/math/special_functions/detail/bessel_j1.hpp	2012-12-07 07:49:42 EST (Fri, 07 Dec 2012)
@@ -166,13 +166,24 @@
     {
         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));
+        //
+        // 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.
+        //
+        T sx = sin(x);
+        T cx = cos(x);
+        value = factor * (rc * (cx * -constants::one_div_root_two<T>() + sx * constants::half_root_two<T>()) 
+           - y * rs * (sx * -constants::one_div_root_two<T>() - cx * constants::half_root_two<T>()));
     }
 
     if (x < 0)

Modified: trunk/boost/math/special_functions/detail/bessel_y0.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/bessel_y0.hpp	(original)
+++ trunk/boost/math/special_functions/detail/bessel_y0.hpp	2012-12-07 07:49:42 EST (Fri, 07 Dec 2012)
@@ -197,11 +197,22 @@
     {
         T y = 8 / x;
         T y2 = y * y;
-        T z = x - 0.25f * pi<T>();
         rc = evaluate_rational(PC, QC, y2);
         rs = evaluate_rational(PS, QS, y2);
         factor = sqrt(2 / (x * pi<T>()));
-        value = factor * (rc * sin(z) + y * rs * cos(z));
+        //
+        // The following code is really just:
+        //
+        // T z = x - 0.25f * pi<T>();
+        // value = factor * (rc * sin(z) + y * rs * cos(z));
+        //
+        // But using the sin/cos addition formulae and constant values for
+        // sin/cos of PI/4:
+        //
+        T sx = sin(x);
+        T cx = cos(x);
+        value = factor * (rc * (sx * constants::one_div_root_two<T>() - cx * constants::half_root_two<T>()) 
+           + y * rs * (cx * constants::one_div_root_two<T>() + sx * constants::half_root_two<T>()));
     }
 
     return value;

Modified: trunk/boost/math/special_functions/detail/bessel_y1.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/bessel_y1.hpp	(original)
+++ trunk/boost/math/special_functions/detail/bessel_y1.hpp	2012-12-07 07:49:42 EST (Fri, 07 Dec 2012)
@@ -170,11 +170,21 @@
     {
         T y = 8 / x;
         T y2 = y * y;
-        T z = x - 0.75f * pi<T>();
         rc = evaluate_rational(PC, QC, y2);
         rs = evaluate_rational(PS, QS, y2);
         factor = sqrt(2 / (x * pi<T>()));
-        value = factor * (rc * sin(z) + y * rs * cos(z));
+        //
+        // This code is really just:
+        //
+        // T z = x - 0.75f * pi<T>();
+        // value = factor * (rc * sin(z) + y * rs * cos(z));
+        //
+        // But using the sin/cos addition rules, plus constants for sin/cos of 3PI/4:
+        //
+        T sx = sin(x);
+        T cx = cos(x);
+        value = factor * (rc * (sx * -constants::one_div_root_two<T>() - cx * constants::half_root_two<T>()) 
+           + y * rs * (cx * -constants::one_div_root_two<T>() + sx * constants::half_root_two<T>()));
     }
 
     return value;