[Boost-commit] svn:boost r78062 - in sandbox/math: boost/math/special_functions libs/math/test

Subject: [Boost-commit] svn:boost r78062 - in sandbox/math: boost/math/special_functions libs/math/test
From: john_at_[hidden]
Date: 2012-04-18 06:55:59

Author: johnmaddock
Date: 2012-04-18 06:55:57 EDT (Wed, 18 Apr 2012)
New Revision: 78062
URL: http://svn.boost.org/trac/boost/changeset/78062

Log:
Add tentative arbitrary precision version of T2.
Add special cases from Owen's original paper as test cases.
Text files modified: 
   sandbox/math/boost/math/special_functions/owens_t.hpp |    47 +++++++++++++++++++++++++++++++++++++-- 
   sandbox/math/libs/math/test/test_owens_t.cpp          |    21 +++++++++++++++++                       
   2 files changed, 64 insertions(+), 4 deletions(-)

Modified: sandbox/math/boost/math/special_functions/owens_t.hpp
==============================================================================
--- sandbox/math/boost/math/special_functions/owens_t.hpp	(original)
+++ sandbox/math/boost/math/special_functions/owens_t.hpp	2012-04-18 06:55:57 EDT (Wed, 18 Apr 2012)
@@ -180,8 +180,8 @@
          } // RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m)
 
          // compute the value of Owen's T function with method T2 from the reference paper
-         template<typename RealType>
-         inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)
+         template<typename RealType, class Policy>
+         inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy&, const mpl::false_&)
          {
             BOOST_MATH_STD_USING
             using namespace boost::math::constants;
@@ -593,6 +593,45 @@
             return (sum / d) / boost::math::constants::two_pi<T>();
          }
 
+         template<typename RealType, class Policy>
+         inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy& pol, const mpl::true_&)
+         {
+            BOOST_MATH_STD_USING
+            using namespace boost::math::constants;
+
+            const unsigned short maxii = m+m+1;
+            const RealType hs = h*h;
+            const RealType as = -a*a;
+            const RealType y = static_cast<RealType>(1) / hs;
+
+            unsigned short ii = 1;
+            RealType val = 0;
+            RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
+            RealType z = owens_t_znorm1(ah)/h;
+            RealType last_z = fabs(z);
+            RealType lim = policies::get_epsilon<RealType, Policy>();
+
+            while( true )
+            {
+               val += z;
+               //
+               // This series stops converging after a while, so put a limit
+               // on how far we go before returning our best guess:
+               //
+               if((fabs(lim * val) > fabs(z)) || ((ii > maxii) && (fabs(z) > last_z)) || (z == 0))
+               {
+                  val *= exp( -hs*half<RealType>() ) / root_two_pi<RealType>();
+                  break;
+               } // if( maxii <= ii )
+               last_z = fabs(z);
+               z = y * ( vi - static_cast<RealType>(ii) * z );
+               vi *= as;
+               ii += 2;
+            } // while( true )
+
+            return val;
+         } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)
+
          // This routine dispatches the call to one of six subroutines, depending on the values
          // of h and a.
          // preconditions: h >= 0, 0<=a<=1, ah=a*h
@@ -630,7 +669,9 @@
                val = owens_t_T1(h,a,m);
                break;
             case 2: // T2
-               val = owens_t_T2(h,a,m,ah);
+               typedef typename policies::precision<RealType, Policy>::type precision_type;
+               typedef mpl::bool_<(precision_type::value == 0) || (precision_type::value > 64)> tag_type;
+               val = owens_t_T2(h, a, m, ah, pol, tag_type());
                break;
             case 3: // T3
                val = owens_t_T3(h,a,ah, pol);

Modified: sandbox/math/libs/math/test/test_owens_t.cpp
==============================================================================
--- sandbox/math/libs/math/test/test_owens_t.cpp	(original)
+++ sandbox/math/libs/math/test/test_owens_t.cpp	2012-04-18 06:55:57 EDT (Wed, 18 Apr 2012)
@@ -32,6 +32,7 @@
 
 #include <boost/math/special_functions/owens_t.hpp> // for owens_t function.
 using boost::math::owens_t;
+#include <boost/math/distributions/normal.hpp>
 
 #include <boost/test/test_exec_monitor.hpp> // for test_main.
 #include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE and  BOOST_CHECK_CLOSE_fraction.
@@ -132,7 +133,8 @@
   cout << "Tolerance = " << tolerance << "." << endl;
 
   using  ::boost::math::owens_t;
-  using namespace std; // ADL of std names.
+  using ::boost::math::normal_distribution;
+  BOOST_MATH_STD_USING // ADL of std names.
 
   // Checks of six sub-methods T1 to T6.
   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.0625L), static_cast<RealType>(0.25L)), static_cast<RealType>(3.89119302347013668966224771378e-2L), tolerance);  // T1
@@ -160,6 +162,23 @@
   BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(0.5L), static_cast<RealType>(2L)), -owens_t(static_cast<RealType>(0.5L), static_cast<RealType>(-2L)));
   BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(0.5L), static_cast<RealType>(2L)), -owens_t(static_cast<RealType>(-0.5L), static_cast<RealType>(-2L)));
 
+  // Special relations from Owen's original paper:
+  BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(0.5), static_cast<RealType>(0)), static_cast<RealType>(0));
+  BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(10), static_cast<RealType>(0)), static_cast<RealType>(0));
+  BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(10000), static_cast<RealType>(0)), static_cast<RealType>(0));
+
+  BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0), static_cast<RealType>(2L)), atan(static_cast<RealType>(2L)) / (boost::math::constants::pi<RealType>() * 2), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0), static_cast<RealType>(0.5L)), atan(static_cast<RealType>(0.5L)) / (boost::math::constants::pi<RealType>() * 2), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0), static_cast<RealType>(2000L)), atan(static_cast<RealType>(2000L)) / (boost::math::constants::pi<RealType>() * 2), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(5), static_cast<RealType>(1)), cdf(normal_distribution<RealType>(), 5) * cdf(complement(normal_distribution<RealType>(), 5)) / 2, tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.125), static_cast<RealType>(1)), cdf(normal_distribution<RealType>(), 0.125) * cdf(complement(normal_distribution<RealType>(), 0.125)) / 2, tolerance);
+  if(std::numeric_limits<RealType>::has_infinity)
+  {
+    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.125), std::numeric_limits<RealType>::infinity()), cdf(complement(normal_distribution<RealType>(), 0.125)) / 2, tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(5), std::numeric_limits<RealType>::infinity()), cdf(complement(normal_distribution<RealType>(), 5)) / 2, tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(-0.125), std::numeric_limits<RealType>::infinity()), cdf(normal_distribution<RealType>(), -0.125) / 2, tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(-5), std::numeric_limits<RealType>::infinity()), cdf(normal_distribution<RealType>(), -5) / 2, tolerance);
+  }
 } // template <class RealType>void test_spots(RealType)
 
 template <class RealType> // Any floating-point type RealType.