Chris@16: // Copyright Benjamin Sobotta 2012
Chris@16: 
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_OWENS_T_HPP
Chris@16: #define BOOST_OWENS_T_HPP
Chris@16: 
Chris@16: // Reference:
Chris@16: // Mike Patefield, David Tandy
Chris@16: // FAST AND ACCURATE CALCULATION OF OWEN'S T-FUNCTION
Chris@16: // Journal of Statistical Software, 5 (5), 1-25
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #  pragma once
Chris@16: #endif
Chris@16: 
Chris@101: #include <boost/math/special_functions/math_fwd.hpp>
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: #include <boost/math/special_functions/erf.hpp>
Chris@16: #include <boost/math/special_functions/expm1.hpp>
Chris@16: #include <boost/throw_exception.hpp>
Chris@16: #include <boost/assert.hpp>
Chris@16: #include <boost/math/constants/constants.hpp>
Chris@16: #include <boost/math/tools/big_constant.hpp>
Chris@16: 
Chris@16: #include <stdexcept>
Chris@16: 
Chris@101: #ifdef BOOST_MSVC
Chris@101: #pragma warning(push)
Chris@101: #pragma warning(disable:4127)
Chris@101: #endif
Chris@101: 
Chris@16: namespace boost
Chris@16: {
Chris@16:    namespace math
Chris@16:    {
Chris@16:       namespace detail
Chris@16:       {
Chris@16:          // owens_t_znorm1(x) = P(-oo<Z<=x)-0.5 with Z being normally distributed.
Chris@16:          template<typename RealType>
Chris@16:          inline RealType owens_t_znorm1(const RealType x)
Chris@16:          {
Chris@16:             using namespace boost::math::constants;
Chris@16:             return erf(x*one_div_root_two<RealType>())*half<RealType>();
Chris@16:          } // RealType owens_t_znorm1(const RealType x)
Chris@16: 
Chris@16:          // owens_t_znorm2(x) = P(x<=Z<oo) with Z being normally distributed.
Chris@16:          template<typename RealType>
Chris@16:          inline RealType owens_t_znorm2(const RealType x)
Chris@16:          {
Chris@16:             using namespace boost::math::constants;
Chris@16:             return erfc(x*one_div_root_two<RealType>())*half<RealType>();
Chris@16:          } // RealType owens_t_znorm2(const RealType x)
Chris@16: 
Chris@16:          // Auxiliary function, it computes an array key that is used to determine
Chris@16:          // the specific computation method for Owen's T and the order thereof
Chris@16:          // used in owens_t_dispatch.
Chris@16:          template<typename RealType>
Chris@16:          inline unsigned short owens_t_compute_code(const RealType h, const RealType a)
Chris@16:          {
Chris@16:             static const RealType hrange[] =
Chris@16:             {0.02, 0.06, 0.09, 0.125, 0.26, 0.4,  0.6,  1.6,  1.7,  2.33,  2.4,  3.36, 3.4,  4.8};
Chris@16: 
Chris@16:             static const RealType arange[] = {0.025, 0.09, 0.15, 0.36, 0.5, 0.9, 0.99999};
Chris@16:             /*
Chris@16:             original select array from paper:
Chris@16:             1, 1, 2,13,13,13,13,13,13,13,13,16,16,16, 9
Chris@16:             1, 2, 2, 3, 3, 5, 5,14,14,15,15,16,16,16, 9
Chris@16:             2, 2, 3, 3, 3, 5, 5,15,15,15,15,16,16,16,10
Chris@16:             2, 2, 3, 5, 5, 5, 5, 7, 7,16,16,16,16,16,10
Chris@16:             2, 3, 3, 5, 5, 6, 6, 8, 8,17,17,17,12,12,11
Chris@16:             2, 3, 5, 5, 5, 6, 6, 8, 8,17,17,17,12,12,12
Chris@16:             2, 3, 4, 4, 6, 6, 8, 8,17,17,17,17,17,12,12
Chris@16:             2, 3, 4, 4, 6, 6,18,18,18,18,17,17,17,12,12
Chris@16:             */                  
Chris@16:             // subtract one because the array is written in FORTRAN in mind - in C arrays start @ zero
Chris@16:             static const unsigned short select[] =
Chris@16:             {
Chris@16:                0,    0 ,   1  , 12   ,12 ,  12  , 12  , 12 ,  12  , 12  , 12  , 15  , 15 ,  15  ,  8,
Chris@16:                0  ,  1  ,  1   , 2 ,   2   , 4  ,  4  , 13 ,  13  , 14  , 14 ,  15  , 15  , 15  ,  8,
Chris@16:                1  ,  1   , 2 ,   2  ,  2  ,  4   , 4  , 14  , 14 ,  14  , 14 ,  15  , 15 ,  15  ,  9,
Chris@16:                1  ,  1   , 2 ,   4  ,  4  ,  4   , 4  ,  6  ,  6 ,  15  , 15 ,  15 ,  15 ,  15  ,  9,
Chris@16:                1  ,  2   , 2  ,  4  ,  4  ,  5   , 5  ,  7  ,  7  , 16   ,16 ,  16 ,  11 ,  11 ,  10,
Chris@16:                1  ,  2   , 4  ,  4   , 4  ,  5   , 5  ,  7  ,  7  , 16  , 16 ,  16 ,  11  , 11 ,  11,
Chris@16:                1  ,  2   , 3  ,  3  ,  5  ,  5   , 7  ,  7  , 16 ,  16  , 16 ,  16 ,  16  , 11 ,  11,
Chris@16:                1  ,  2   , 3   , 3   , 5  ,  5 ,  17  , 17  , 17 ,  17  , 16 ,  16 ,  16 ,  11 ,  11
Chris@16:             };
Chris@16: 
Chris@16:             unsigned short ihint = 14, iaint = 7;
Chris@16:             for(unsigned short i = 0; i != 14; i++)
Chris@16:             {
Chris@16:                if( h <= hrange[i] )
Chris@16:                {
Chris@16:                   ihint = i;
Chris@16:                   break;
Chris@16:                }
Chris@16:             } // for(unsigned short i = 0; i != 14; i++)
Chris@16: 
Chris@16:             for(unsigned short i = 0; i != 7; i++)
Chris@16:             {
Chris@16:                if( a <= arange[i] )
Chris@16:                {
Chris@16:                   iaint = i;
Chris@16:                   break;
Chris@16:                }
Chris@16:             } // for(unsigned short i = 0; i != 7; i++)
Chris@16: 
Chris@16:             // interprete select array as 8x15 matrix
Chris@16:             return select[iaint*15 + ihint];
Chris@16: 
Chris@16:          } // unsigned short owens_t_compute_code(const RealType h, const RealType a)
Chris@16: 
Chris@16:          template<typename RealType>
Chris@16:          inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const mpl::int_<53>&)
Chris@16:          {
Chris@16:             static const unsigned short ord[] = {2, 3, 4, 5, 7, 10, 12, 18, 10, 20, 30, 0, 4, 7, 8, 20, 0, 0}; // 18 entries
Chris@16: 
Chris@16:             BOOST_ASSERT(icode<18);
Chris@16: 
Chris@16:             return ord[icode];
Chris@16:          } // unsigned short owens_t_get_order(const unsigned short icode, RealType, mpl::int<53> const&)
Chris@16: 
Chris@16:          template<typename RealType>
Chris@16:          inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const mpl::int_<64>&)
Chris@16:         {
Chris@16:            // method ================>>>       {1, 1, 1, 1, 1,  1,  1,  1,  2,  2,  2,  3, 4,  4,  4,  4,  5, 6}
Chris@16:            static const unsigned short ord[] = {3, 4, 5, 6, 8, 11, 13, 19, 10, 20, 30,  0, 7, 10, 11, 23,  0, 0}; // 18 entries
Chris@16: 
Chris@16:           BOOST_ASSERT(icode<18);
Chris@16: 
Chris@16:           return ord[icode];
Chris@16:         } // unsigned short owens_t_get_order(const unsigned short icode, RealType, mpl::int<64> const&)
Chris@16: 
Chris@16:          template<typename RealType, typename Policy>
Chris@16:          inline unsigned short owens_t_get_order(const unsigned short icode, RealType r, const Policy&)
Chris@16:          {
Chris@16:             typedef typename policies::precision<RealType, Policy>::type precision_type;
Chris@16:             typedef typename mpl::if_<
Chris@16:                mpl::or_<
Chris@16:                   mpl::less_equal<precision_type, mpl::int_<0> >,
Chris@16:                   mpl::greater<precision_type, mpl::int_<53> >
Chris@16:                >,
Chris@16:                mpl::int_<64>,
Chris@16:                mpl::int_<53>
Chris@16:             >::type tag_type;
Chris@16: 
Chris@16:             return owens_t_get_order_imp(icode, r, tag_type());
Chris@16:          }
Chris@16: 
Chris@16:          // compute the value of Owen's T function with method T1 from the reference paper
Chris@101:          template<typename RealType, typename Policy>
Chris@101:          inline RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m, const Policy& pol)
Chris@16:          {
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             using namespace boost::math::constants;
Chris@16: 
Chris@16:             const RealType hs = -h*h*half<RealType>();
Chris@16:             const RealType dhs = exp( hs );
Chris@16:             const RealType as = a*a;
Chris@16: 
Chris@16:             unsigned short j=1;
Chris@16:             RealType jj = 1;
Chris@16:             RealType aj = a * one_div_two_pi<RealType>();
Chris@101:             RealType dj = boost::math::expm1( hs, pol);
Chris@16:             RealType gj = hs*dhs;
Chris@16: 
Chris@16:             RealType val = atan( a ) * one_div_two_pi<RealType>();
Chris@16: 
Chris@16:             while( true )
Chris@16:             {
Chris@16:                val += dj*aj/jj;
Chris@16: 
Chris@16:                if( m <= j )
Chris@16:                   break;
Chris@16: 
Chris@16:                j++;
Chris@16:                jj += static_cast<RealType>(2);
Chris@16:                aj *= as;
Chris@16:                dj = gj - dj;
Chris@16:                gj *= hs / static_cast<RealType>(j);
Chris@16:             } // while( true )
Chris@16: 
Chris@16:             return val;
Chris@16:          } // RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m)
Chris@16: 
Chris@16:          // compute the value of Owen's T function with method T2 from the reference paper
Chris@16:          template<typename RealType, class Policy>
Chris@16:          inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy&, const mpl::false_&)
Chris@16:          {
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             using namespace boost::math::constants;
Chris@16: 
Chris@16:             const unsigned short maxii = m+m+1;
Chris@16:             const RealType hs = h*h;
Chris@16:             const RealType as = -a*a;
Chris@16:             const RealType y = static_cast<RealType>(1) / hs;
Chris@16: 
Chris@16:             unsigned short ii = 1;
Chris@16:             RealType val = 0;
Chris@16:             RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
Chris@16:             RealType z = owens_t_znorm1(ah)/h;
Chris@16: 
Chris@16:             while( true )
Chris@16:             {
Chris@16:                val += z;
Chris@16:                if( maxii <= ii )
Chris@16:                {
Chris@16:                   val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
Chris@16:                   break;
Chris@16:                } // if( maxii <= ii )
Chris@16:                z = y * ( vi - static_cast<RealType>(ii) * z );
Chris@16:                vi *= as;
Chris@16:                ii += 2;
Chris@16:             } // while( true )
Chris@16: 
Chris@16:             return val;
Chris@16:          } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)
Chris@16: 
Chris@16:          // compute the value of Owen's T function with method T3 from the reference paper
Chris@16:          template<typename RealType>
Chris@16:          inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const mpl::int_<53>&)
Chris@16:          {
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             using namespace boost::math::constants;
Chris@16: 
Chris@16:       const unsigned short m = 20;
Chris@16: 
Chris@16:             static const RealType c2[] =
Chris@16:             {
Chris@16:                0.99999999999999987510,
Chris@16:                -0.99999999999988796462,      0.99999999998290743652,
Chris@16:                -0.99999999896282500134,      0.99999996660459362918,
Chris@16:                -0.99999933986272476760,      0.99999125611136965852,
Chris@16:                -0.99991777624463387686,      0.99942835555870132569,
Chris@16:                -0.99697311720723000295,      0.98751448037275303682,
Chris@16:                -0.95915857980572882813,      0.89246305511006708555,
Chris@16:                -0.76893425990463999675,      0.58893528468484693250,
Chris@16:                -0.38380345160440256652,      0.20317601701045299653,
Chris@16:                -0.82813631607004984866E-01,  0.24167984735759576523E-01,
Chris@16:                -0.44676566663971825242E-02,  0.39141169402373836468E-03
Chris@16:             };
Chris@16: 
Chris@16:             const RealType as = a*a;
Chris@16:             const RealType hs = h*h;
Chris@16:             const RealType y = static_cast<RealType>(1)/hs;
Chris@16: 
Chris@16:             RealType ii = 1;
Chris@16:             unsigned short i = 0;
Chris@16:             RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
Chris@16:             RealType zi = owens_t_znorm1(ah)/h;
Chris@16:             RealType val = 0;
Chris@16: 
Chris@16:             while( true )
Chris@16:             {
Chris@16:                BOOST_ASSERT(i < 21);
Chris@16:                val += zi*c2[i];
Chris@16:                if( m <= i ) // if( m < i+1 )
Chris@16:                {
Chris@16:                   val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
Chris@16:                   break;
Chris@16:                } // if( m < i )
Chris@16:                zi = y * (ii*zi - vi);
Chris@16:                vi *= as;
Chris@16:                ii += 2;
Chris@16:                i++;
Chris@16:             } // while( true )
Chris@16: 
Chris@16:             return val;
Chris@16:          } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah)
Chris@16: 
Chris@16:         // compute the value of Owen's T function with method T3 from the reference paper
Chris@16:         template<class RealType>
Chris@16:         inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const mpl::int_<64>&)
Chris@16:         {
Chris@16:           BOOST_MATH_STD_USING
Chris@16:           using namespace boost::math::constants;
Chris@16:           
Chris@16:           const unsigned short m = 30;
Chris@16: 
Chris@16:           static const RealType c2[] =
Chris@16:           {
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999999999729978162447266851932041876728736094298092917625009873),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99999999999999999999467056379678391810626533251885323416799874878563998732905968),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999824849349313270659391127814689133077036298754586814091034842536),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999999997703859616213643405880166422891953033591551179153879839440241685),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999998394883415238173334565554173013941245103172035286759201504179038147),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999993063616095509371081203145247992197457263066869044528823599399470977),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999999797336340409464429599229870590160411238245275855903767652432017766116267),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999999999574958412069046680119051639753412378037565521359444170241346845522403274),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999933226234193375324943920160947158239076786103108097456617750134812033362048),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999188923242461073033481053037468263536806742737922476636768006622772762168467),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999992195143483674402853783549420883055129680082932629160081128947764415749728967),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999993935137206712830997921913316971472227199741857386575097250553105958772041501),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99996135597690552745362392866517133091672395614263398912807169603795088421057688716),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99979556366513946026406788969630293820987757758641211293079784585126692672425362469),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.999092789629617100153486251423850590051366661947344315423226082520411961968929483),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.996593837411918202119308620432614600338157335862888580671450938858935084316004769854),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.98910017138386127038463510314625339359073956513420458166238478926511821146316469589567),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.970078558040693314521331982203762771512160168582494513347846407314584943870399016019),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.92911438683263187495758525500033707204091967947532160289872782771388170647150321633673),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.8542058695956156057286980736842905011429254735181323743367879525470479126968822863),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.73796526033030091233118357742803709382964420335559408722681794195743240930748630755),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.58523469882837394570128599003785154144164680587615878645171632791404210655891158),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.415997776145676306165661663581868460503874205343014196580122174949645271353372263),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.2588210875241943574388730510317252236407805082485246378222935376279663808416534365),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.1375535825163892648504646951500265585055789019410617565727090346559210218472356689),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.0607952766325955730493900985022020434830339794955745989150270485056436844239206648),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0216337683299871528059836483840390514275488679530797294557060229266785853764115),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.00593405693455186729876995814181203900550014220428843483927218267309209471516256),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0011743414818332946510474576182739210553333860106811865963485870668929503649964142),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, -1.489155613350368934073453260689881330166342484405529981510694514036264969925132e-4),
Chris@16:              BOOST_MATH_BIG_CONSTANT(RealType, 260, 9.072354320794357587710929507988814669454281514268844884841547607134260303118208e-6)
Chris@16:           };
Chris@16: 
Chris@16:           const RealType as = a*a;
Chris@16:           const RealType hs = h*h;
Chris@16:           const RealType y = 1 / hs;
Chris@16: 
Chris@16:           RealType ii = 1;
Chris@16:           unsigned short i = 0;
Chris@16:           RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
Chris@16:           RealType zi = owens_t_znorm1(ah)/h;
Chris@16:           RealType val = 0;
Chris@16: 
Chris@16:           while( true )
Chris@16:           {
Chris@16:               BOOST_ASSERT(i < 31);
Chris@16:               val += zi*c2[i];
Chris@16:               if( m <= i ) // if( m < i+1 )
Chris@16:               {
Chris@16:                 val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
Chris@16:                 break;
Chris@16:               } // if( m < i )
Chris@16:               zi = y * (ii*zi - vi);
Chris@16:               vi *= as;
Chris@16:               ii += 2;
Chris@16:               i++;
Chris@16:           } // while( true )
Chris@16: 
Chris@16:           return val;
Chris@16:         } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah)
Chris@16: 
Chris@16:         template<class RealType, class Policy>
Chris@16:         inline RealType owens_t_T3(const RealType h, const RealType a, const RealType ah, const Policy&)
Chris@16:         {
Chris@16:             typedef typename policies::precision<RealType, Policy>::type precision_type;
Chris@16:             typedef typename mpl::if_<
Chris@16:                mpl::or_<
Chris@16:                   mpl::less_equal<precision_type, mpl::int_<0> >,
Chris@16:                   mpl::greater<precision_type, mpl::int_<53> >
Chris@16:                >,
Chris@16:                mpl::int_<64>,
Chris@16:                mpl::int_<53>
Chris@16:             >::type tag_type;
Chris@16: 
Chris@16:             return owens_t_T3_imp(h, a, ah, tag_type());
Chris@16:         }
Chris@16: 
Chris@16:          // compute the value of Owen's T function with method T4 from the reference paper
Chris@16:          template<typename RealType>
Chris@16:          inline RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m)
Chris@16:          {
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             using namespace boost::math::constants;
Chris@16: 
Chris@16:             const unsigned short maxii = m+m+1;
Chris@16:             const RealType hs = h*h;
Chris@16:             const RealType as = -a*a;
Chris@16: 
Chris@16:             unsigned short ii = 1;
Chris@16:             RealType ai = a * exp( -hs*(static_cast<RealType>(1)-as)*half<RealType>() ) * one_div_two_pi<RealType>();
Chris@16:             RealType yi = 1;
Chris@16:             RealType val = 0;
Chris@16: 
Chris@16:             while( true )
Chris@16:             {
Chris@16:                val += ai*yi;
Chris@16:                if( maxii <= ii )
Chris@16:                   break;
Chris@16:                ii += 2;
Chris@16:                yi = (static_cast<RealType>(1)-hs*yi) / static_cast<RealType>(ii);
Chris@16:                ai *= as;
Chris@16:             } // while( true )
Chris@16: 
Chris@16:             return val;
Chris@16:          } // RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m)
Chris@16: 
Chris@16:          // compute the value of Owen's T function with method T5 from the reference paper
Chris@16:          template<typename RealType>
Chris@16:          inline RealType owens_t_T5_imp(const RealType h, const RealType a, const mpl::int_<53>&)
Chris@16:          {
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             /*
Chris@16:                NOTICE:
Chris@16:                - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre
Chris@16:                  polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre
Chris@16:                  quadrature, because T5(h,a,m) contains only x^2 terms.
Chris@16:                - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor
Chris@16:                  of 1/(2*pi) according to T5(h,a,m).
Chris@16:              */
Chris@16: 
Chris@16:             const unsigned short m = 13;
Chris@16:             static const RealType pts[] = {0.35082039676451715489E-02,
Chris@16:                0.31279042338030753740E-01,  0.85266826283219451090E-01,
Chris@16:                0.16245071730812277011,      0.25851196049125434828,
Chris@16:                0.36807553840697533536,      0.48501092905604697475,
Chris@16:                0.60277514152618576821,      0.71477884217753226516,
Chris@16:                0.81475510988760098605,      0.89711029755948965867,
Chris@16:                0.95723808085944261843,      0.99178832974629703586};
Chris@16:             static const RealType wts[] = { 0.18831438115323502887E-01,
Chris@16:                0.18567086243977649478E-01,  0.18042093461223385584E-01,
Chris@16:                0.17263829606398753364E-01,  0.16243219975989856730E-01,
Chris@16:                0.14994592034116704829E-01,  0.13535474469662088392E-01,
Chris@16:                0.11886351605820165233E-01,  0.10070377242777431897E-01,
Chris@16:                0.81130545742299586629E-02,  0.60419009528470238773E-02,
Chris@16:                0.38862217010742057883E-02,  0.16793031084546090448E-02};
Chris@16: 
Chris@16:             const RealType as = a*a;
Chris@16:             const RealType hs = -h*h*boost::math::constants::half<RealType>();
Chris@16: 
Chris@16:             RealType val = 0;
Chris@16:             for(unsigned short i = 0; i < m; ++i)
Chris@16:             {
Chris@16:                BOOST_ASSERT(i < 13);
Chris@16:                const RealType r = static_cast<RealType>(1) + as*pts[i];
Chris@16:                val += wts[i] * exp( hs*r ) / r;
Chris@16:             } // for(unsigned short i = 0; i < m; ++i)
Chris@16: 
Chris@16:             return val*a;
Chris@16:          } // RealType owens_t_T5(const RealType h, const RealType a)
Chris@16: 
Chris@16:         // compute the value of Owen's T function with method T5 from the reference paper
Chris@16:         template<typename RealType>
Chris@16:         inline RealType owens_t_T5_imp(const RealType h, const RealType a, const mpl::int_<64>&)
Chris@16:         {
Chris@16:           BOOST_MATH_STD_USING
Chris@16:             /*
Chris@16:               NOTICE:
Chris@16:               - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre
Chris@16:               polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre
Chris@16:               quadrature, because T5(h,a,m) contains only x^2 terms.
Chris@16:               - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor
Chris@16:               of 1/(2*pi) according to T5(h,a,m).
Chris@16:             */
Chris@16: 
Chris@16:           const unsigned short m = 19;
Chris@16:           static const RealType pts[] = {
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0016634282895983227941),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.014904509242697054183),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.04103478879005817919),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.079359853513391511008),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.1288612130237615133),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.18822336642448518856),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.25586876186122962384),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.32999972011807857222),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.40864620815774761438),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.48971819306044782365),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.57106118513245543894),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.6505134942981533829),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.72596367859928091618),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.79540665919549865924),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.85699701386308739244),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.90909804422384697594),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.95032536436570154409),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.97958418733152273717),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.99610366384229088321)
Chris@16:           };
Chris@16:           static const RealType wts[] = {
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012975111395684900835),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012888764187499150078),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012716644398857307844),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012459897461364705691),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012120231988292330388),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011699908404856841158),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011201723906897224448),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.010628993848522759853),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0099855296835573320047),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0092756136096132857933),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0085039700881139589055),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0076757344408814561254),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0067964187616556459109),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.005871875456524750363),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0049082589542498110071),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0039119870792519721409),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0028897090921170700834),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0018483371329504443947),
Chris@16:                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.00079623320100438873578)
Chris@16:           };
Chris@16: 
Chris@16:           const RealType as = a*a;
Chris@16:           const RealType hs = -h*h*boost::math::constants::half<RealType>();
Chris@16: 
Chris@16:           RealType val = 0;
Chris@16:           for(unsigned short i = 0; i < m; ++i)
Chris@16:             {
Chris@16:               BOOST_ASSERT(i < 19);
Chris@16:               const RealType r = 1 + as*pts[i];
Chris@16:               val += wts[i] * exp( hs*r ) / r;
Chris@16:             } // for(unsigned short i = 0; i < m; ++i)
Chris@16: 
Chris@16:           return val*a;
Chris@16:         } // RealType owens_t_T5(const RealType h, const RealType a)
Chris@16: 
Chris@16:         template<class RealType, class Policy>
Chris@16:         inline RealType owens_t_T5(const RealType h, const RealType a, const Policy&)
Chris@16:         {
Chris@16:             typedef typename policies::precision<RealType, Policy>::type precision_type;
Chris@16:             typedef typename mpl::if_<
Chris@16:                mpl::or_<
Chris@16:                   mpl::less_equal<precision_type, mpl::int_<0> >,
Chris@16:                   mpl::greater<precision_type, mpl::int_<53> >
Chris@16:                >,
Chris@16:                mpl::int_<64>,
Chris@16:                mpl::int_<53>
Chris@16:             >::type tag_type;
Chris@16: 
Chris@16:             return owens_t_T5_imp(h, a, tag_type());
Chris@16:         }
Chris@16: 
Chris@16: 
Chris@16:          // compute the value of Owen's T function with method T6 from the reference paper
Chris@16:          template<typename RealType>
Chris@16:          inline RealType owens_t_T6(const RealType h, const RealType a)
Chris@16:          {
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             using namespace boost::math::constants;
Chris@16: 
Chris@16:             const RealType normh = owens_t_znorm2( h );
Chris@16:             const RealType y = static_cast<RealType>(1) - a;
Chris@16:             const RealType r = atan2(y, static_cast<RealType>(1 + a) );
Chris@16: 
Chris@16:             RealType val = normh * ( static_cast<RealType>(1) - normh ) * half<RealType>();
Chris@16: 
Chris@16:             if( r != 0 )
Chris@16:                val -= r * exp( -y*h*h*half<RealType>()/r ) * one_div_two_pi<RealType>();
Chris@16: 
Chris@16:             return val;
Chris@16:          } // RealType owens_t_T6(const RealType h, const RealType a, const unsigned short m)
Chris@16: 
Chris@16:          template <class T, class Policy>
Chris@16:          std::pair<T, T> owens_t_T1_accelerated(T h, T a, const Policy& pol)
Chris@16:          {
Chris@16:             //
Chris@16:             // This is the same series as T1, but:
Chris@16:             // * The Taylor series for atan has been combined with that for T1, 
Chris@16:             //   reducing but not eliminating cancellation error.
Chris@16:             // * The resulting alternating series is then accelerated using method 1
Chris@16:             //   from H. Cohen, F. Rodriguez Villegas, D. Zagier, 
Chris@16:             //   "Convergence acceleration of alternating series", Bonn, (1991).
Chris@16:             //
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             static const char* function = "boost::math::owens_t<%1%>(%1%, %1%)";
Chris@16:             T half_h_h = h * h / 2;
Chris@16:             T a_pow = a;
Chris@16:             T aa = a * a;
Chris@16:             T exp_term = exp(-h * h / 2);
Chris@16:             T one_minus_dj_sum = exp_term; 
Chris@16:             T sum = a_pow * exp_term;
Chris@16:             T dj_pow = exp_term;
Chris@16:             T term = sum;
Chris@16:             T abs_err;
Chris@16:             int j = 1;
Chris@16: 
Chris@16:             //
Chris@16:             // Normally with this form of series acceleration we can calculate
Chris@16:             // up front how many terms will be required - based on the assumption
Chris@16:             // that each term decreases in size by a factor of 3.  However,
Chris@16:             // that assumption does not apply here, as the underlying T1 series can 
Chris@16:             // go quite strongly divergent in the early terms, before strongly
Chris@16:             // converging later.  Various "guestimates" have been tried to take account
Chris@16:             // of this, but they don't always work.... so instead set "n" to the 
Chris@16:             // largest value that won't cause overflow later, and abort iteration
Chris@16:             // when the last accelerated term was small enough...
Chris@16:             //
Chris@16:             int n;
Chris@16:             try
Chris@16:             {
Chris@16:                n = itrunc(T(tools::log_max_value<T>() / 6));
Chris@16:             }
Chris@16:             catch(...)
Chris@16:             {
Chris@16:                n = (std::numeric_limits<int>::max)();
Chris@16:             }
Chris@16:             n = (std::min)(n, 1500);
Chris@16:             T d = pow(3 + sqrt(T(8)), n);
Chris@16:             d = (d + 1 / d) / 2;
Chris@16:             T b = -1;
Chris@16:             T c = -d;
Chris@16:             c = b - c;
Chris@16:             sum *= c;
Chris@16:             b = -n * n * b * 2;
Chris@16:             abs_err = ldexp(fabs(sum), -tools::digits<T>());
Chris@16: 
Chris@16:             while(j < n)
Chris@16:             {
Chris@16:                a_pow *= aa;
Chris@16:                dj_pow *= half_h_h / j;
Chris@16:                one_minus_dj_sum += dj_pow;
Chris@16:                term = one_minus_dj_sum * a_pow / (2 * j + 1);
Chris@16:                c = b - c;
Chris@16:                sum += c * term;
Chris@16:                abs_err += ldexp((std::max)(T(fabs(sum)), T(fabs(c*term))), -tools::digits<T>());
Chris@16:                b = (j + n) * (j - n) * b / ((j + T(0.5)) * (j + 1));
Chris@16:                ++j;
Chris@16:                //
Chris@16:                // Include an escape route to prevent calculating too many terms:
Chris@16:                //
Chris@16:                if((j > 10) && (fabs(sum * tools::epsilon<T>()) > fabs(c * term)))
Chris@16:                   break;
Chris@16:             }
Chris@16:             abs_err += fabs(c * term);
Chris@16:             if(sum < 0)  // sum must always be positive, if it's negative something really bad has happend:
Chris@16:                policies::raise_evaluation_error(function, 0, T(0), pol);
Chris@16:             return std::pair<T, T>((sum / d) / boost::math::constants::two_pi<T>(), abs_err / sum);
Chris@16:          }
Chris@16: 
Chris@16:          template<typename RealType, class Policy>
Chris@16:          inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy&, const mpl::true_&)
Chris@16:          {
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             using namespace boost::math::constants;
Chris@16: 
Chris@16:             const unsigned short maxii = m+m+1;
Chris@16:             const RealType hs = h*h;
Chris@16:             const RealType as = -a*a;
Chris@16:             const RealType y = static_cast<RealType>(1) / hs;
Chris@16: 
Chris@16:             unsigned short ii = 1;
Chris@16:             RealType val = 0;
Chris@16:             RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
Chris@16:             RealType z = owens_t_znorm1(ah)/h;
Chris@16:             RealType last_z = fabs(z);
Chris@16:             RealType lim = policies::get_epsilon<RealType, Policy>();
Chris@16: 
Chris@16:             while( true )
Chris@16:             {
Chris@16:                val += z;
Chris@16:                //
Chris@16:                // This series stops converging after a while, so put a limit
Chris@16:                // on how far we go before returning our best guess:
Chris@16:                //
Chris@16:                if((fabs(lim * val) > fabs(z)) || ((ii > maxii) && (fabs(z) > last_z)) || (z == 0))
Chris@16:                {
Chris@16:                   val *= exp( -hs*half<RealType>() ) / root_two_pi<RealType>();
Chris@16:                   break;
Chris@16:                } // if( maxii <= ii )
Chris@16:                last_z = fabs(z);
Chris@16:                z = y * ( vi - static_cast<RealType>(ii) * z );
Chris@16:                vi *= as;
Chris@16:                ii += 2;
Chris@16:             } // while( true )
Chris@16: 
Chris@16:             return val;
Chris@16:          } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)
Chris@16: 
Chris@16:          template<typename RealType, class Policy>
Chris@16:          inline std::pair<RealType, RealType> owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&)
Chris@16:          {
Chris@16:             //
Chris@16:             // This is the same series as T2, but with acceleration applied.
Chris@16:             // Note that we have to be *very* careful to check that nothing bad
Chris@16:             // has happened during evaluation - this series will go divergent
Chris@16:             // and/or fail to alternate at a drop of a hat! :-(
Chris@16:             //
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             using namespace boost::math::constants;
Chris@16: 
Chris@16:             const RealType hs = h*h;
Chris@16:             const RealType as = -a*a;
Chris@16:             const RealType y = static_cast<RealType>(1) / hs;
Chris@16: 
Chris@16:             unsigned short ii = 1;
Chris@16:             RealType val = 0;
Chris@16:             RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
Chris@16:             RealType z = boost::math::detail::owens_t_znorm1(ah)/h;
Chris@16:             RealType last_z = fabs(z);
Chris@16: 
Chris@16:             //
Chris@16:             // Normally with this form of series acceleration we can calculate
Chris@16:             // up front how many terms will be required - based on the assumption
Chris@16:             // that each term decreases in size by a factor of 3.  However,
Chris@16:             // that assumption does not apply here, as the underlying T1 series can 
Chris@16:             // go quite strongly divergent in the early terms, before strongly
Chris@16:             // converging later.  Various "guestimates" have been tried to take account
Chris@16:             // of this, but they don't always work.... so instead set "n" to the 
Chris@16:             // largest value that won't cause overflow later, and abort iteration
Chris@16:             // when the last accelerated term was small enough...
Chris@16:             //
Chris@16:             int n;
Chris@16:             try
Chris@16:             {
Chris@16:                n = itrunc(RealType(tools::log_max_value<RealType>() / 6));
Chris@16:             }
Chris@16:             catch(...)
Chris@16:             {
Chris@16:                n = (std::numeric_limits<int>::max)();
Chris@16:             }
Chris@16:             n = (std::min)(n, 1500);
Chris@16:             RealType d = pow(3 + sqrt(RealType(8)), n);
Chris@16:             d = (d + 1 / d) / 2;
Chris@16:             RealType b = -1;
Chris@16:             RealType c = -d;
Chris@16:             int s = 1;
Chris@16: 
Chris@16:             for(int k = 0; k < n; ++k)
Chris@16:             {
Chris@16:                //
Chris@16:                // Check for both convergence and whether the series has gone bad:
Chris@16:                //
Chris@16:                if(
Chris@16:                   (fabs(z) > last_z)     // Series has gone divergent, abort
Chris@16:                   || (fabs(val) * tools::epsilon<RealType>() > fabs(c * s * z))  // Convergence!
Chris@16:                   || (z * s < 0)         // Series has stopped alternating - all bets are off - abort.
Chris@16:                   )
Chris@16:                {
Chris@16:                   break;
Chris@16:                }
Chris@16:                c = b - c;
Chris@16:                val += c * s * z;
Chris@16:                b = (k + n) * (k - n) * b / ((k + RealType(0.5)) * (k + 1));
Chris@16:                last_z = fabs(z);
Chris@16:                s = -s;
Chris@16:                z = y * ( vi - static_cast<RealType>(ii) * z );
Chris@16:                vi *= as;
Chris@16:                ii += 2;
Chris@16:             } // while( true )
Chris@16:             RealType err = fabs(c * z) / val;
Chris@16:             return std::pair<RealType, RealType>(val * exp( -hs*half<RealType>() ) / (d * root_two_pi<RealType>()), err);
Chris@16:          } // RealType owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&)
Chris@16: 
Chris@16:          template<typename RealType, typename Policy>
Chris@16:          inline RealType T4_mp(const RealType h, const RealType a, const Policy& pol)
Chris@16:          {
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             
Chris@16:             const RealType hs = h*h;
Chris@16:             const RealType as = -a*a;
Chris@16: 
Chris@16:             unsigned short ii = 1;
Chris@16:             RealType ai = constants::one_div_two_pi<RealType>() * a * exp( -0.5*hs*(1.0-as) );
Chris@16:             RealType yi = 1.0;
Chris@16:             RealType val = 0.0;
Chris@16: 
Chris@16:             RealType lim = boost::math::policies::get_epsilon<RealType, Policy>();
Chris@16: 
Chris@16:             while( true )
Chris@16:             {
Chris@16:                RealType term = ai*yi;
Chris@16:                val += term;
Chris@16:                if((yi != 0) && (fabs(val * lim) > fabs(term)))
Chris@16:                   break;
Chris@16:                ii += 2;
Chris@16:                yi = (1.0-hs*yi) / static_cast<RealType>(ii);
Chris@16:                ai *= as;
Chris@16:                if(ii > (std::min)(1500, (int)policies::get_max_series_iterations<Policy>()))
Chris@16:                   policies::raise_evaluation_error("boost::math::owens_t<%1%>", 0, val, pol);
Chris@16:             } // while( true )
Chris@16: 
Chris@16:             return val;
Chris@16:          } // arg_type owens_t_T4(const arg_type h, const arg_type a, const unsigned short m)
Chris@16: 
Chris@16: 
Chris@16:          // This routine dispatches the call to one of six subroutines, depending on the values
Chris@16:          // of h and a.
Chris@16:          // preconditions: h >= 0, 0<=a<=1, ah=a*h
Chris@16:          //
Chris@16:          // Note there are different versions for different precisions....
Chris@16:          template<typename RealType, typename Policy>
Chris@16:          inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, mpl::int_<64> const&)
Chris@16:          {
Chris@16:             // Simple main case for 64-bit precision or less, this is as per the Patefield-Tandy paper:
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             //
Chris@16:             // Handle some special cases first, these are from
Chris@16:             // page 1077 of Owen's original paper:
Chris@16:             //
Chris@16:             if(h == 0)
Chris@16:             {
Chris@16:                return atan(a) * constants::one_div_two_pi<RealType>();
Chris@16:             }
Chris@16:             if(a == 0)
Chris@16:             {
Chris@16:                return 0;
Chris@16:             }
Chris@16:             if(a == 1)
Chris@16:             {
Chris@16:                return owens_t_znorm2(RealType(-h)) * owens_t_znorm2(h) / 2;
Chris@16:             }
Chris@16:             if(a >= tools::max_value<RealType>())
Chris@16:             {
Chris@16:                return owens_t_znorm2(RealType(fabs(h)));
Chris@16:             }
Chris@16:             RealType val = 0; // avoid compiler warnings, 0 will be overwritten in any case
Chris@16:             const unsigned short icode = owens_t_compute_code(h, a);
Chris@16:             const unsigned short m = owens_t_get_order(icode, val /* just a dummy for the type */, pol);
Chris@16:             static const unsigned short meth[] = {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6}; // 18 entries
Chris@16: 
Chris@16:             // determine the appropriate method, T1 ... T6
Chris@16:             switch( meth[icode] )
Chris@16:             {
Chris@16:             case 1: // T1
Chris@101:                val = owens_t_T1(h,a,m,pol);
Chris@16:                break;
Chris@16:             case 2: // T2
Chris@16:                typedef typename policies::precision<RealType, Policy>::type precision_type;
Chris@16:                typedef mpl::bool_<(precision_type::value == 0) || (precision_type::value > 64)> tag_type;
Chris@16:                val = owens_t_T2(h, a, m, ah, pol, tag_type());
Chris@16:                break;
Chris@16:             case 3: // T3
Chris@16:                val = owens_t_T3(h,a,ah, pol);
Chris@16:                break;
Chris@16:             case 4: // T4
Chris@16:                val = owens_t_T4(h,a,m);
Chris@16:                break;
Chris@16:             case 5: // T5
Chris@16:                val = owens_t_T5(h,a, pol);
Chris@16:                break;
Chris@16:             case 6: // T6
Chris@16:                val = owens_t_T6(h,a);
Chris@16:                break;
Chris@16:             default:
Chris@16:                BOOST_THROW_EXCEPTION(std::logic_error("selection routine in Owen's T function failed"));
Chris@16:             }
Chris@16:             return val;
Chris@16:          }
Chris@16: 
Chris@16:          template<typename RealType, typename Policy>
Chris@16:          inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const mpl::int_<65>&)
Chris@16:          {
Chris@16:             // Arbitrary precision version:
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             //
Chris@16:             // Handle some special cases first, these are from
Chris@16:             // page 1077 of Owen's original paper:
Chris@16:             //
Chris@16:             if(h == 0)
Chris@16:             {
Chris@16:                return atan(a) * constants::one_div_two_pi<RealType>();
Chris@16:             }
Chris@16:             if(a == 0)
Chris@16:             {
Chris@16:                return 0;
Chris@16:             }
Chris@16:             if(a == 1)
Chris@16:             {
Chris@16:                return owens_t_znorm2(RealType(-h)) * owens_t_znorm2(h) / 2;
Chris@16:             }
Chris@16:             if(a >= tools::max_value<RealType>())
Chris@16:             {
Chris@16:                return owens_t_znorm2(RealType(fabs(h)));
Chris@16:             }
Chris@16:             // Attempt arbitrary precision code, this will throw if it goes wrong:
Chris@16:             typedef typename boost::math::policies::normalise<Policy, boost::math::policies::evaluation_error<> >::type forwarding_policy;
Chris@16:             std::pair<RealType, RealType> p1(0, tools::max_value<RealType>()), p2(0, tools::max_value<RealType>());
Chris@16:             RealType target_precision = policies::get_epsilon<RealType, Policy>() * 1000;
Chris@16:             bool have_t1(false), have_t2(false);
Chris@16:             if(ah < 3)
Chris@16:             {
Chris@16:                try
Chris@16:                {
Chris@16:                   have_t1 = true;
Chris@16:                   p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
Chris@16:                   if(p1.second < target_precision)
Chris@16:                      return p1.first;
Chris@16:                }
Chris@16:                catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
Chris@16:             }
Chris@16:             if(ah > 1)
Chris@16:             {
Chris@16:                try
Chris@16:                {
Chris@16:                   have_t2 = true;
Chris@16:                   p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
Chris@16:                   if(p2.second < target_precision)
Chris@16:                      return p2.first;
Chris@16:                }
Chris@16:                catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
Chris@16:             }
Chris@16:             //
Chris@16:             // If we haven't tried T1 yet, do it now - sometimes it succeeds and the number of iterations
Chris@16:             // is fairly low compared to T4.
Chris@16:             //
Chris@16:             if(!have_t1)
Chris@16:             {
Chris@16:                try
Chris@16:                {
Chris@16:                   have_t1 = true;
Chris@16:                   p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
Chris@16:                   if(p1.second < target_precision)
Chris@16:                      return p1.first;
Chris@16:                }
Chris@16:                catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
Chris@16:             }
Chris@16:             //
Chris@16:             // If we haven't tried T2 yet, do it now - sometimes it succeeds and the number of iterations
Chris@16:             // is fairly low compared to T4.
Chris@16:             //
Chris@16:             if(!have_t2)
Chris@16:             {
Chris@16:                try
Chris@16:                {
Chris@16:                   have_t2 = true;
Chris@16:                   p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
Chris@16:                   if(p2.second < target_precision)
Chris@16:                      return p2.first;
Chris@16:                }
Chris@16:                catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
Chris@16:             }
Chris@16:             //
Chris@16:             // OK, nothing left to do but try the most expensive option which is T4,
Chris@16:             // this is often slow to converge, but when it does converge it tends to
Chris@16:             // be accurate:
Chris@16:             try
Chris@16:             {
Chris@16:                return T4_mp(h, a, pol);
Chris@16:             }
Chris@16:             catch(const boost::math::evaluation_error&){}  // T4 may fail and throw, that's OK
Chris@16:             //
Chris@16:             // Now look back at the results from T1 and T2 and see if either gave better
Chris@16:             // results than we could get from the 64-bit precision versions.
Chris@16:             //
Chris@16:             if((std::min)(p1.second, p2.second) < 1e-20)
Chris@16:             {
Chris@16:                return p1.second < p2.second ? p1.first : p2.first;
Chris@16:             }
Chris@16:             //
Chris@16:             // We give up - no arbitrary precision versions succeeded!
Chris@16:             //
Chris@16:             return owens_t_dispatch(h, a, ah, pol, mpl::int_<64>());
Chris@16:          } // RealType owens_t_dispatch(RealType h, RealType a, RealType ah)
Chris@16:          template<typename RealType, typename Policy>
Chris@16:          inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const mpl::int_<0>&)
Chris@16:          {
Chris@16:             // We don't know what the precision is until runtime:
Chris@16:             if(tools::digits<RealType>() <= 64)
Chris@16:                return owens_t_dispatch(h, a, ah, pol, mpl::int_<64>());
Chris@16:             return owens_t_dispatch(h, a, ah, pol, mpl::int_<65>());
Chris@16:          }
Chris@16:          template<typename RealType, typename Policy>
Chris@16:          inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol)
Chris@16:          {
Chris@16:             // Figure out the precision and forward to the correct version:
Chris@16:             typedef typename policies::precision<RealType, Policy>::type precision_type;
Chris@16:             typedef typename mpl::if_c<
Chris@16:                precision_type::value == 0,
Chris@16:                mpl::int_<0>,
Chris@16:                typename mpl::if_c<
Chris@16:                   precision_type::value <= 64,
Chris@16:                   mpl::int_<64>,
Chris@16:                   mpl::int_<65>
Chris@16:                >::type
Chris@16:             >::type tag_type;
Chris@16:             return owens_t_dispatch(h, a, ah, pol, tag_type());
Chris@16:          }
Chris@16:          // compute Owen's T function, T(h,a), for arbitrary values of h and a
Chris@16:          template<typename RealType, class Policy>
Chris@16:          inline RealType owens_t(RealType h, RealType a, const Policy& pol)
Chris@16:          {
Chris@16:             BOOST_MATH_STD_USING
Chris@16:             // exploit that T(-h,a) == T(h,a)
Chris@16:             h = fabs(h);
Chris@16: 
Chris@16:             // Use equation (2) in the paper to remap the arguments
Chris@16:             // such that h>=0 and 0<=a<=1 for the call of the actual
Chris@16:             // computation routine.
Chris@16: 
Chris@16:             const RealType fabs_a = fabs(a);
Chris@16:             const RealType fabs_ah = fabs_a*h;
Chris@16: 
Chris@16:             RealType val = 0.0; // avoid compiler warnings, 0.0 will be overwritten in any case
Chris@16: 
Chris@16:             if(fabs_a <= 1)
Chris@16:             {
Chris@16:                val = owens_t_dispatch(h, fabs_a, fabs_ah, pol);
Chris@16:             } // if(fabs_a <= 1.0)
Chris@16:             else 
Chris@16:             {
Chris@16:                if( h <= 0.67 )
Chris@16:                {
Chris@16:                   const RealType normh = owens_t_znorm1(h);
Chris@16:                   const RealType normah = owens_t_znorm1(fabs_ah);
Chris@16:                   val = static_cast<RealType>(1)/static_cast<RealType>(4) - normh*normah -
Chris@16:                      owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol);
Chris@16:                } // if( h <= 0.67 )
Chris@16:                else
Chris@16:                {
Chris@16:                   const RealType normh = detail::owens_t_znorm2(h);
Chris@16:                   const RealType normah = detail::owens_t_znorm2(fabs_ah);
Chris@16:                   val = constants::half<RealType>()*(normh+normah) - normh*normah -
Chris@16:                      owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol);
Chris@16:                } // else [if( h <= 0.67 )]
Chris@16:             } // else [if(fabs_a <= 1)]
Chris@16: 
Chris@16:             // exploit that T(h,-a) == -T(h,a)
Chris@16:             if(a < 0)
Chris@16:             {
Chris@16:                return -val;
Chris@16:             } // if(a < 0)
Chris@16: 
Chris@16:             return val;
Chris@16:          } // RealType owens_t(RealType h, RealType a)
Chris@16: 
Chris@16:          template <class T, class Policy, class tag>
Chris@16:          struct owens_t_initializer
Chris@16:          {
Chris@16:             struct init
Chris@16:             {
Chris@16:                init()
Chris@16:                {
Chris@16:                   do_init(tag());
Chris@16:                }
Chris@16:                template <int N>
Chris@16:                static void do_init(const mpl::int_<N>&){}
Chris@16:                static void do_init(const mpl::int_<64>&)
Chris@16:                {
Chris@16:                   boost::math::owens_t(static_cast<T>(7), static_cast<T>(0.96875), Policy());
Chris@16:                   boost::math::owens_t(static_cast<T>(2), static_cast<T>(0.5), Policy());
Chris@16:                }
Chris@16:                void force_instantiate()const{}
Chris@16:             };
Chris@16:             static const init initializer;
Chris@16:             static void force_instantiate()
Chris@16:             {
Chris@16:                initializer.force_instantiate();
Chris@16:             }
Chris@16:          };
Chris@16: 
Chris@16:          template <class T, class Policy, class tag>
Chris@16:          const typename owens_t_initializer<T, Policy, tag>::init owens_t_initializer<T, Policy, tag>::initializer;
Chris@16: 
Chris@16:       } // namespace detail
Chris@16: 
Chris@16:       template <class T1, class T2, class Policy>
Chris@16:       inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a, const Policy& pol)
Chris@16:       {
Chris@16:          typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:          typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:          typedef typename policies::precision<value_type, Policy>::type precision_type;
Chris@16:          typedef typename mpl::if_c<
Chris@16:                precision_type::value == 0,
Chris@16:                mpl::int_<0>,
Chris@16:                typename mpl::if_c<
Chris@16:                   precision_type::value <= 64,
Chris@16:                   mpl::int_<64>,
Chris@16:                   mpl::int_<65>
Chris@16:                >::type
Chris@16:             >::type tag_type;
Chris@16: 
Chris@16:          detail::owens_t_initializer<result_type, Policy, tag_type>::force_instantiate();
Chris@16:             
Chris@16:          return policies::checked_narrowing_cast<result_type, Policy>(detail::owens_t(static_cast<value_type>(h), static_cast<value_type>(a), pol), "boost::math::owens_t<%1%>(%1%,%1%)");
Chris@16:       }
Chris@16: 
Chris@16:       template <class T1, class T2>
Chris@16:       inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a)
Chris@16:       {
Chris@16:          return owens_t(h, a, policies::policy<>());
Chris@16:       }
Chris@16: 
Chris@16: 
Chris@16:    } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@101: #ifdef BOOST_MSVC
Chris@101: #pragma warning(pop)
Chris@101: #endif
Chris@101: 
Chris@16: #endif
Chris@16: // EOF