math/special_functions/owens_t.hpp

0001 // Copyright Benjamin Sobotta 2012
0002
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_OWENS_T_HPP
0008 #define BOOST_OWENS_T_HPP
0009
0010 // Reference:
0011 // Mike Patefield, David Tandy
0012 // FAST AND ACCURATE CALCULATION OF OWEN'S T-FUNCTION
0013 // Journal of Statistical Software, 5 (5), 1-25
0014
0015 #ifdef _MSC_VER
0016 #  pragma once
0017 #endif
0018
0019 #include <boost/math/special_functions/math_fwd.hpp>
0020 #include <boost/math/special_functions/erf.hpp>
0021 #include <boost/math/special_functions/expm1.hpp>
0022 #include <boost/math/tools/throw_exception.hpp>
0023 #include <boost/math/tools/assert.hpp>
0024 #include <boost/math/constants/constants.hpp>
0025 #include <boost/math/tools/big_constant.hpp>
0026
0027 #include <stdexcept>
0028 #include <cmath>
0029
0030 #ifdef _MSC_VER
0031 #pragma warning(push)
0032 #pragma warning(disable:4127)
0033 #endif
0034
0035 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0036 //
0037 // This is the only way we can avoid
0038 // warning: non-standard suffix on floating constant [-Wpedantic]
0039 // when building with -Wall -pedantic.  Neither __extension__
0040 // nor #pragma diagnostic ignored work :(
0041 //
0042 #pragma GCC system_header
0043 #endif
0044
0045 namespace boost
0046 {
0047    namespace math
0048    {
0049       namespace detail
0050       {
0051          // owens_t_znorm1(x) = P(-oo<Z<=x)-0.5 with Z being normally distributed.
0052          template<typename RealType, class Policy>
0053          inline RealType owens_t_znorm1(const RealType x, const Policy& pol)
0054          {
0055             using namespace boost::math::constants;
0056             return boost::math::erf(x*one_div_root_two<RealType>(), pol)*half<RealType>();
0057          } // RealType owens_t_znorm1(const RealType x)
0058
0059          // owens_t_znorm2(x) = P(x<=Z<oo) with Z being normally distributed.
0060          template<typename RealType, class Policy>
0061          inline RealType owens_t_znorm2(const RealType x, const Policy& pol)
0062          {
0063             using namespace boost::math::constants;
0064             return boost::math::erfc(x*one_div_root_two<RealType>(), pol)*half<RealType>();
0065          } // RealType owens_t_znorm2(const RealType x)
0066
0067          // Auxiliary function, it computes an array key that is used to determine
0068          // the specific computation method for Owen's T and the order thereof
0069          // used in owens_t_dispatch.
0070          template<typename RealType>
0071          inline unsigned short owens_t_compute_code(const RealType h, const RealType a)
0072          {
0073             static const RealType hrange[] =
0074             { 0.02f, 0.06f, 0.09f, 0.125f, 0.26f, 0.4f,  0.6f,  1.6f,  1.7f,  2.33f,  2.4f,  3.36f, 3.4f,  4.8f };
0075
0076             static const RealType arange[] = { 0.025f, 0.09f, 0.15f, 0.36f, 0.5f, 0.9f, 0.99999f };
0077             /*
0078             original select array from paper:
0079             1, 1, 2,13,13,13,13,13,13,13,13,16,16,16, 9
0080             1, 2, 2, 3, 3, 5, 5,14,14,15,15,16,16,16, 9
0081             2, 2, 3, 3, 3, 5, 5,15,15,15,15,16,16,16,10
0082             2, 2, 3, 5, 5, 5, 5, 7, 7,16,16,16,16,16,10
0083             2, 3, 3, 5, 5, 6, 6, 8, 8,17,17,17,12,12,11
0084             2, 3, 5, 5, 5, 6, 6, 8, 8,17,17,17,12,12,12
0085             2, 3, 4, 4, 6, 6, 8, 8,17,17,17,17,17,12,12
0086             2, 3, 4, 4, 6, 6,18,18,18,18,17,17,17,12,12
0087             */
0088             // subtract one because the array is written in FORTRAN in mind - in C arrays start @ zero
0089             static const unsigned short select[] =
0090             {
0091                0,    0 ,   1  , 12   ,12 ,  12  , 12  , 12 ,  12  , 12  , 12  , 15  , 15 ,  15  ,  8,
0092                0  ,  1  ,  1   , 2 ,   2   , 4  ,  4  , 13 ,  13  , 14  , 14 ,  15  , 15  , 15  ,  8,
0093                1  ,  1   , 2 ,   2  ,  2  ,  4   , 4  , 14  , 14 ,  14  , 14 ,  15  , 15 ,  15  ,  9,
0094                1  ,  1   , 2 ,   4  ,  4  ,  4   , 4  ,  6  ,  6 ,  15  , 15 ,  15 ,  15 ,  15  ,  9,
0095                1  ,  2   , 2  ,  4  ,  4  ,  5   , 5  ,  7  ,  7  , 16   ,16 ,  16 ,  11 ,  11 ,  10,
0096                1  ,  2   , 4  ,  4   , 4  ,  5   , 5  ,  7  ,  7  , 16  , 16 ,  16 ,  11  , 11 ,  11,
0097                1  ,  2   , 3  ,  3  ,  5  ,  5   , 7  ,  7  , 16 ,  16  , 16 ,  16 ,  16  , 11 ,  11,
0098                1  ,  2   , 3   , 3   , 5  ,  5 ,  17  , 17  , 17 ,  17  , 16 ,  16 ,  16 ,  11 ,  11
0099             };
0100
0101             unsigned short ihint = 14, iaint = 7;
0102             for(unsigned short i = 0; i != 14; i++)
0103             {
0104                if( h <= hrange[i] )
0105                {
0106                   ihint = i;
0107                   break;
0108                }
0109             } // for(unsigned short i = 0; i != 14; i++)
0110
0111             for(unsigned short i = 0; i != 7; i++)
0112             {
0113                if( a <= arange[i] )
0114                {
0115                   iaint = i;
0116                   break;
0117                }
0118             } // for(unsigned short i = 0; i != 7; i++)
0119
0120             // interpret select array as 8x15 matrix
0121             return select[iaint*15 + ihint];
0122
0123          } // unsigned short owens_t_compute_code(const RealType h, const RealType a)
0124
0125          template<typename RealType>
0126          inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const std::integral_constant<int, 53>&)
0127          {
0128             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
0129
0130             BOOST_MATH_ASSERT(icode<18);
0131
0132             return ord[icode];
0133          } // unsigned short owens_t_get_order(const unsigned short icode, RealType, std::integral_constant<int, 53> const&)
0134
0135          template<typename RealType>
0136          inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const std::integral_constant<int, 64>&)
0137         {
0138            // method ================>>>       {1, 1, 1, 1, 1,  1,  1,  1,  2,  2,  2,  3, 4,  4,  4,  4,  5, 6}
0139            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
0140
0141           BOOST_MATH_ASSERT(icode<18);
0142
0143           return ord[icode];
0144         } // unsigned short owens_t_get_order(const unsigned short icode, RealType, std::integral_constant<int, 64> const&)
0145
0146          template<typename RealType, typename Policy>
0147          inline unsigned short owens_t_get_order(const unsigned short icode, RealType r, const Policy&)
0148          {
0149             typedef typename policies::precision<RealType, Policy>::type precision_type;
0150             typedef std::integral_constant<int,
0151                precision_type::value <= 0 ? 64 :
0152                precision_type::value <= 53 ? 53 : 64
0153             > tag_type;
0154
0155             return owens_t_get_order_imp(icode, r, tag_type());
0156          }
0157
0158          // compute the value of Owen's T function with method T1 from the reference paper
0159          template<typename RealType, typename Policy>
0160          inline RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m, const Policy& pol)
0161          {
0162             BOOST_MATH_STD_USING
0163             using namespace boost::math::constants;
0164
0165             const RealType hs = -h*h*half<RealType>();
0166             const RealType dhs = exp( hs );
0167             const RealType as = a*a;
0168
0169             unsigned short j=1;
0170             RealType jj = 1;
0171             RealType aj = a * one_div_two_pi<RealType>();
0172             RealType dj = boost::math::expm1( hs, pol);
0173             RealType gj = hs*dhs;
0174
0175             RealType val = atan( a ) * one_div_two_pi<RealType>();
0176
0177             while( true )
0178             {
0179                val += dj*aj/jj;
0180
0181                if( m <= j )
0182                   break;
0183
0184                j++;
0185                jj += static_cast<RealType>(2);
0186                aj *= as;
0187                dj = gj - dj;
0188                gj *= hs / static_cast<RealType>(j);
0189             } // while( true )
0190
0191             return val;
0192          } // RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m)
0193
0194          // compute the value of Owen's T function with method T2 from the reference paper
0195          template<typename RealType, class Policy>
0196          inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy& pol, const std::false_type&)
0197          {
0198             BOOST_MATH_STD_USING
0199             using namespace boost::math::constants;
0200
0201             const unsigned short maxii = m+m+1;
0202             const RealType hs = h*h;
0203             const RealType as = -a*a;
0204             const RealType y = static_cast<RealType>(1) / hs;
0205
0206             unsigned short ii = 1;
0207             RealType val = 0;
0208             RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
0209             RealType z = owens_t_znorm1(ah, pol)/h;
0210
0211             while( true )
0212             {
0213                val += z;
0214                if( maxii <= ii )
0215                {
0216                   val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
0217                   break;
0218                } // if( maxii <= ii )
0219                z = y * ( vi - static_cast<RealType>(ii) * z );
0220                vi *= as;
0221                ii += 2;
0222             } // while( true )
0223
0224             return val;
0225          } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)
0226
0227          // compute the value of Owen's T function with method T3 from the reference paper
0228          template<typename RealType, class Policy>
0229          inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const std::integral_constant<int, 53>&, const Policy& pol)
0230          {
0231             BOOST_MATH_STD_USING
0232             using namespace boost::math::constants;
0233
0234       const unsigned short m = 20;
0235
0236             static const RealType c2[] =
0237             {
0238                static_cast<RealType>(0.99999999999999987510),
0239                static_cast<RealType>(-0.99999999999988796462),      static_cast<RealType>(0.99999999998290743652),
0240                static_cast<RealType>(-0.99999999896282500134),      static_cast<RealType>(0.99999996660459362918),
0241                static_cast<RealType>(-0.99999933986272476760),      static_cast<RealType>(0.99999125611136965852),
0242                static_cast<RealType>(-0.99991777624463387686),      static_cast<RealType>(0.99942835555870132569),
0243                static_cast<RealType>(-0.99697311720723000295),      static_cast<RealType>(0.98751448037275303682),
0244                static_cast<RealType>(-0.95915857980572882813),      static_cast<RealType>(0.89246305511006708555),
0245                static_cast<RealType>(-0.76893425990463999675),      static_cast<RealType>(0.58893528468484693250),
0246                static_cast<RealType>(-0.38380345160440256652),      static_cast<RealType>(0.20317601701045299653),
0247                static_cast<RealType>(-0.82813631607004984866E-01),  static_cast<RealType>(0.24167984735759576523E-01),
0248                static_cast<RealType>(-0.44676566663971825242E-02),  static_cast<RealType>(0.39141169402373836468E-03)
0249             };
0250
0251             const RealType as = a*a;
0252             const RealType hs = h*h;
0253             const RealType y = static_cast<RealType>(1)/hs;
0254
0255             RealType ii = 1;
0256             unsigned short i = 0;
0257             RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
0258             RealType zi = owens_t_znorm1(ah, pol)/h;
0259             RealType val = 0;
0260
0261             while( true )
0262             {
0263                BOOST_MATH_ASSERT(i < 21);
0264                val += zi*c2[i];
0265                if( m <= i ) // if( m < i+1 )
0266                {
0267                   val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
0268                   break;
0269                } // if( m < i )
0270                zi = y * (ii*zi - vi);
0271                vi *= as;
0272                ii += 2;
0273                i++;
0274             } // while( true )
0275
0276             return val;
0277          } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah)
0278
0279         // compute the value of Owen's T function with method T3 from the reference paper
0280         template<class RealType, class Policy>
0281         inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const std::integral_constant<int, 64>&, const Policy& pol)
0282         {
0283           BOOST_MATH_STD_USING
0284           using namespace boost::math::constants;
0285
0286           const unsigned short m = 30;
0287
0288           static const RealType c2[] =
0289           {
0290              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999999999729978162447266851932041876728736094298092917625009873),
0291              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99999999999999999999467056379678391810626533251885323416799874878563998732905968),
0292              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999824849349313270659391127814689133077036298754586814091034842536),
0293              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999999997703859616213643405880166422891953033591551179153879839440241685),
0294              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999998394883415238173334565554173013941245103172035286759201504179038147),
0295              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999993063616095509371081203145247992197457263066869044528823599399470977),
0296              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999999797336340409464429599229870590160411238245275855903767652432017766116267),
0297              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999999999574958412069046680119051639753412378037565521359444170241346845522403274),
0298              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999933226234193375324943920160947158239076786103108097456617750134812033362048),
0299              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999188923242461073033481053037468263536806742737922476636768006622772762168467),
0300              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999992195143483674402853783549420883055129680082932629160081128947764415749728967),
0301              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999993935137206712830997921913316971472227199741857386575097250553105958772041501),
0302              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99996135597690552745362392866517133091672395614263398912807169603795088421057688716),
0303              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99979556366513946026406788969630293820987757758641211293079784585126692672425362469),
0304              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.999092789629617100153486251423850590051366661947344315423226082520411961968929483),
0305              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.996593837411918202119308620432614600338157335862888580671450938858935084316004769854),
0306              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.98910017138386127038463510314625339359073956513420458166238478926511821146316469589567),
0307              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.970078558040693314521331982203762771512160168582494513347846407314584943870399016019),
0308              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.92911438683263187495758525500033707204091967947532160289872782771388170647150321633673),
0309              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.8542058695956156057286980736842905011429254735181323743367879525470479126968822863),
0310              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.73796526033030091233118357742803709382964420335559408722681794195743240930748630755),
0311              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.58523469882837394570128599003785154144164680587615878645171632791404210655891158),
0312              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.415997776145676306165661663581868460503874205343014196580122174949645271353372263),
0313              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.2588210875241943574388730510317252236407805082485246378222935376279663808416534365),
0314              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.1375535825163892648504646951500265585055789019410617565727090346559210218472356689),
0315              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.0607952766325955730493900985022020434830339794955745989150270485056436844239206648),
0316              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0216337683299871528059836483840390514275488679530797294557060229266785853764115),
0317              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.00593405693455186729876995814181203900550014220428843483927218267309209471516256),
0318              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0011743414818332946510474576182739210553333860106811865963485870668929503649964142),
0319              BOOST_MATH_BIG_CONSTANT(RealType, 260, -1.489155613350368934073453260689881330166342484405529981510694514036264969925132e-4),
0320              BOOST_MATH_BIG_CONSTANT(RealType, 260, 9.072354320794357587710929507988814669454281514268844884841547607134260303118208e-6)
0321           };
0322
0323           const RealType as = a*a;
0324           const RealType hs = h*h;
0325           const RealType y = 1 / hs;
0326
0327           RealType ii = 1;
0328           unsigned short i = 0;
0329           RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
0330           RealType zi = owens_t_znorm1(ah, pol)/h;
0331           RealType val = 0;
0332
0333           while( true )
0334           {
0335               BOOST_MATH_ASSERT(i < 31);
0336               val += zi*c2[i];
0337               if( m <= i ) // if( m < i+1 )
0338               {
0339                 val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
0340                 break;
0341               } // if( m < i )
0342               zi = y * (ii*zi - vi);
0343               vi *= as;
0344               ii += 2;
0345               i++;
0346           } // while( true )
0347
0348           return val;
0349         } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah)
0350
0351         template<class RealType, class Policy>
0352         inline RealType owens_t_T3(const RealType h, const RealType a, const RealType ah, const Policy& pol)
0353         {
0354             typedef typename policies::precision<RealType, Policy>::type precision_type;
0355             typedef std::integral_constant<int,
0356                precision_type::value <= 0 ? 64 :
0357                precision_type::value <= 53 ? 53 : 64
0358             > tag_type;
0359
0360             return owens_t_T3_imp(h, a, ah, tag_type(), pol);
0361         }
0362
0363          // compute the value of Owen's T function with method T4 from the reference paper
0364          template<typename RealType>
0365          inline RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m)
0366          {
0367             BOOST_MATH_STD_USING
0368             using namespace boost::math::constants;
0369
0370             const unsigned short maxii = m+m+1;
0371             const RealType hs = h*h;
0372             const RealType as = -a*a;
0373
0374             unsigned short ii = 1;
0375             RealType ai = a * exp( -hs*(static_cast<RealType>(1)-as)*half<RealType>() ) * one_div_two_pi<RealType>();
0376             RealType yi = 1;
0377             RealType val = 0;
0378
0379             while( true )
0380             {
0381                val += ai*yi;
0382                if( maxii <= ii )
0383                   break;
0384                ii += 2;
0385                yi = (static_cast<RealType>(1)-hs*yi) / static_cast<RealType>(ii);
0386                ai *= as;
0387             } // while( true )
0388
0389             return val;
0390          } // RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m)
0391
0392          // compute the value of Owen's T function with method T5 from the reference paper
0393          template<typename RealType>
0394          inline RealType owens_t_T5_imp(const RealType h, const RealType a, const std::integral_constant<int, 53>&)
0395          {
0396             BOOST_MATH_STD_USING
0397             /*
0398                NOTICE:
0399                - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre
0400                  polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre
0401                  quadrature, because T5(h,a,m) contains only x^2 terms.
0402                - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor
0403                  of 1/(2*pi) according to T5(h,a,m).
0404              */
0405
0406             const unsigned short m = 13;
0407             static const RealType pts[] = {
0408                static_cast<RealType>(0.35082039676451715489E-02),
0409                static_cast<RealType>(0.31279042338030753740E-01),  static_cast<RealType>(0.85266826283219451090E-01),
0410                static_cast<RealType>(0.16245071730812277011),      static_cast<RealType>(0.25851196049125434828),
0411                static_cast<RealType>(0.36807553840697533536),      static_cast<RealType>(0.48501092905604697475),
0412                static_cast<RealType>(0.60277514152618576821),      static_cast<RealType>(0.71477884217753226516),
0413                static_cast<RealType>(0.81475510988760098605),      static_cast<RealType>(0.89711029755948965867),
0414                static_cast<RealType>(0.95723808085944261843),      static_cast<RealType>(0.99178832974629703586) };
0415             static const RealType wts[] = {
0416                static_cast<RealType>(0.18831438115323502887E-01),
0417                static_cast<RealType>(0.18567086243977649478E-01),  static_cast<RealType>(0.18042093461223385584E-01),
0418                static_cast<RealType>(0.17263829606398753364E-01),  static_cast<RealType>(0.16243219975989856730E-01),
0419                static_cast<RealType>(0.14994592034116704829E-01),  static_cast<RealType>(0.13535474469662088392E-01),
0420                static_cast<RealType>(0.11886351605820165233E-01),  static_cast<RealType>(0.10070377242777431897E-01),
0421                static_cast<RealType>(0.81130545742299586629E-02),  static_cast<RealType>(0.60419009528470238773E-02),
0422                static_cast<RealType>(0.38862217010742057883E-02),  static_cast<RealType>(0.16793031084546090448E-02) };
0423
0424             const RealType as = a*a;
0425             const RealType hs = -h*h*boost::math::constants::half<RealType>();
0426
0427             RealType val = 0;
0428             for(unsigned short i = 0; i < m; ++i)
0429             {
0430                BOOST_MATH_ASSERT(i < 13);
0431                const RealType r = static_cast<RealType>(1) + as*pts[i];
0432                val += wts[i] * exp( hs*r ) / r;
0433             } // for(unsigned short i = 0; i < m; ++i)
0434
0435             return val*a;
0436          } // RealType owens_t_T5(const RealType h, const RealType a)
0437
0438         // compute the value of Owen's T function with method T5 from the reference paper
0439         template<typename RealType>
0440         inline RealType owens_t_T5_imp(const RealType h, const RealType a, const std::integral_constant<int, 64>&)
0441         {
0442           BOOST_MATH_STD_USING
0443             /*
0444               NOTICE:
0445               - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre
0446               polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre
0447               quadrature, because T5(h,a,m) contains only x^2 terms.
0448               - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor
0449               of 1/(2*pi) according to T5(h,a,m).
0450             */
0451
0452           const unsigned short m = 19;
0453           static const RealType pts[] = {
0454                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0016634282895983227941),
0455                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.014904509242697054183),
0456                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.04103478879005817919),
0457                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.079359853513391511008),
0458                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.1288612130237615133),
0459                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.18822336642448518856),
0460                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.25586876186122962384),
0461                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.32999972011807857222),
0462                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.40864620815774761438),
0463                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.48971819306044782365),
0464                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.57106118513245543894),
0465                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.6505134942981533829),
0466                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.72596367859928091618),
0467                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.79540665919549865924),
0468                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.85699701386308739244),
0469                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.90909804422384697594),
0470                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.95032536436570154409),
0471                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.97958418733152273717),
0472                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.99610366384229088321)
0473           };
0474           static const RealType wts[] = {
0475                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012975111395684900835),
0476                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012888764187499150078),
0477                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012716644398857307844),
0478                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012459897461364705691),
0479                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012120231988292330388),
0480                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011699908404856841158),
0481                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011201723906897224448),
0482                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.010628993848522759853),
0483                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0099855296835573320047),
0484                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0092756136096132857933),
0485                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0085039700881139589055),
0486                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0076757344408814561254),
0487                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0067964187616556459109),
0488                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.005871875456524750363),
0489                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0049082589542498110071),
0490                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0039119870792519721409),
0491                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0028897090921170700834),
0492                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0018483371329504443947),
0493                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.00079623320100438873578)
0494           };
0495
0496           const RealType as = a*a;
0497           const RealType hs = -h*h*boost::math::constants::half<RealType>();
0498
0499           RealType val = 0;
0500           for(unsigned short i = 0; i < m; ++i)
0501             {
0502               BOOST_MATH_ASSERT(i < 19);
0503               const RealType r = 1 + as*pts[i];
0504               val += wts[i] * exp( hs*r ) / r;
0505             } // for(unsigned short i = 0; i < m; ++i)
0506
0507           return val*a;
0508         } // RealType owens_t_T5(const RealType h, const RealType a)
0509
0510         template<class RealType, class Policy>
0511         inline RealType owens_t_T5(const RealType h, const RealType a, const Policy&)
0512         {
0513             typedef typename policies::precision<RealType, Policy>::type precision_type;
0514             typedef std::integral_constant<int,
0515                precision_type::value <= 0 ? 64 :
0516                precision_type::value <= 53 ? 53 : 64
0517             > tag_type;
0518
0519             return owens_t_T5_imp(h, a, tag_type());
0520         }
0521
0522
0523          // compute the value of Owen's T function with method T6 from the reference paper
0524          template<typename RealType, class Policy>
0525          inline RealType owens_t_T6(const RealType h, const RealType a, const Policy& pol)
0526          {
0527             BOOST_MATH_STD_USING
0528             using namespace boost::math::constants;
0529
0530             const RealType normh = owens_t_znorm2(h, pol);
0531             const RealType y = static_cast<RealType>(1) - a;
0532             const RealType r = atan2(y, static_cast<RealType>(1 + a) );
0533
0534             RealType val = normh * ( static_cast<RealType>(1) - normh ) * half<RealType>();
0535
0536             if( r != 0 )
0537                val -= r * exp( -y*h*h*half<RealType>()/r ) * one_div_two_pi<RealType>();
0538
0539             return val;
0540          } // RealType owens_t_T6(const RealType h, const RealType a, const unsigned short m)
0541
0542          template <class T, class Policy>
0543          std::pair<T, T> owens_t_T1_accelerated(T h, T a, const Policy& pol)
0544          {
0545             //
0546             // This is the same series as T1, but:
0547             // * The Taylor series for atan has been combined with that for T1,
0548             //   reducing but not eliminating cancellation error.
0549             // * The resulting alternating series is then accelerated using method 1
0550             //   from H. Cohen, F. Rodriguez Villegas, D. Zagier,
0551             //   "Convergence acceleration of alternating series", Bonn, (1991).
0552             //
0553             BOOST_MATH_STD_USING
0554             static const char* function = "boost::math::owens_t<%1%>(%1%, %1%)";
0555             T half_h_h = h * h / 2;
0556             T a_pow = a;
0557             T aa = a * a;
0558             T exp_term = exp(-h * h / 2);
0559             T one_minus_dj_sum = exp_term;
0560             T sum = a_pow * exp_term;
0561             T dj_pow = exp_term;
0562             T term = sum;
0563             T abs_err;
0564             int j = 1;
0565
0566             //
0567             // Normally with this form of series acceleration we can calculate
0568             // up front how many terms will be required - based on the assumption
0569             // that each term decreases in size by a factor of 3.  However,
0570             // that assumption does not apply here, as the underlying T1 series can
0571             // go quite strongly divergent in the early terms, before strongly
0572             // converging later.  Various "guesstimates" have been tried to take account
0573             // of this, but they don't always work.... so instead set "n" to the
0574             // largest value that won't cause overflow later, and abort iteration
0575             // when the last accelerated term was small enough...
0576             //
0577             int n;
0578 #ifndef BOOST_NO_EXCEPTIONS
0579             try
0580             {
0581 #endif
0582                n = itrunc(T(tools::log_max_value<T>() / 6));
0583 #ifndef BOOST_NO_EXCEPTIONS
0584             }
0585             catch(...)
0586             {
0587                n = (std::numeric_limits<int>::max)();
0588             }
0589 #endif
0590             n = (std::min)(n, 1500);
0591             T d = pow(3 + sqrt(T(8)), T(n));
0592             d = (d + 1 / d) / 2;
0593             T b = -1;
0594             T c = -d;
0595             c = b - c;
0596             sum *= c;
0597             b = -n * n * b * 2;
0598             abs_err = ldexp(fabs(sum), -tools::digits<T>());
0599
0600             while(j < n)
0601             {
0602                a_pow *= aa;
0603                dj_pow *= half_h_h / j;
0604                one_minus_dj_sum += dj_pow;
0605                term = one_minus_dj_sum * a_pow / (2 * j + 1);
0606                c = b - c;
0607                sum += c * term;
0608                abs_err += ldexp((std::max)(T(fabs(sum)), T(fabs(c*term))), -tools::digits<T>());
0609                b = (j + n) * (j - n) * b / ((j + T(0.5)) * (j + 1));
0610                ++j;
0611                //
0612                // Include an escape route to prevent calculating too many terms:
0613                //
0614                if((j > 10) && (fabs(sum * tools::epsilon<T>()) > fabs(c * term)))
0615                   break;
0616             }
0617             abs_err += fabs(c * term);
0618             if(sum < 0)  // sum must always be positive, if it's negative something really bad has happened:
0619                policies::raise_evaluation_error(function, 0, T(0), pol);
0620             return std::pair<T, T>((sum / d) / boost::math::constants::two_pi<T>(), abs_err / sum);
0621          }
0622
0623          template<typename RealType, class Policy>
0624          inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy& pol, const std::true_type&)
0625          {
0626             BOOST_MATH_STD_USING
0627             using namespace boost::math::constants;
0628
0629             const unsigned short maxii = m+m+1;
0630             const RealType hs = h*h;
0631             const RealType as = -a*a;
0632             const RealType y = static_cast<RealType>(1) / hs;
0633
0634             unsigned short ii = 1;
0635             RealType val = 0;
0636             RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
0637             RealType z = owens_t_znorm1(ah, pol)/h;
0638             RealType last_z = fabs(z);
0639             RealType lim = policies::get_epsilon<RealType, Policy>();
0640
0641             while( true )
0642             {
0643                val += z;
0644                //
0645                // This series stops converging after a while, so put a limit
0646                // on how far we go before returning our best guess:
0647                //
0648                if((fabs(lim * val) > fabs(z)) || ((ii > maxii) && (fabs(z) > last_z)) || (z == 0))
0649                {
0650                   val *= exp( -hs*half<RealType>() ) / root_two_pi<RealType>();
0651                   break;
0652                } // if( maxii <= ii )
0653                last_z = fabs(z);
0654                z = y * ( vi - static_cast<RealType>(ii) * z );
0655                vi *= as;
0656                ii += 2;
0657             } // while( true )
0658
0659             return val;
0660          } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)
0661
0662          template<typename RealType, class Policy>
0663          inline std::pair<RealType, RealType> owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy& pol)
0664          {
0665             //
0666             // This is the same series as T2, but with acceleration applied.
0667             // Note that we have to be *very* careful to check that nothing bad
0668             // has happened during evaluation - this series will go divergent
0669             // and/or fail to alternate at a drop of a hat! :-(
0670             //
0671             BOOST_MATH_STD_USING
0672             using namespace boost::math::constants;
0673
0674             const RealType hs = h*h;
0675             const RealType as = -a*a;
0676             const RealType y = static_cast<RealType>(1) / hs;
0677
0678             unsigned short ii = 1;
0679             RealType val = 0;
0680             RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
0681             RealType z = boost::math::detail::owens_t_znorm1(ah, pol)/h;
0682             RealType last_z = fabs(z);
0683
0684             //
0685             // Normally with this form of series acceleration we can calculate
0686             // up front how many terms will be required - based on the assumption
0687             // that each term decreases in size by a factor of 3.  However,
0688             // that assumption does not apply here, as the underlying T1 series can
0689             // go quite strongly divergent in the early terms, before strongly
0690             // converging later.  Various "guesstimates" have been tried to take account
0691             // of this, but they don't always work.... so instead set "n" to the
0692             // largest value that won't cause overflow later, and abort iteration
0693             // when the last accelerated term was small enough...
0694             //
0695             int n;
0696 #ifndef BOOST_NO_EXCEPTIONS
0697             try
0698             {
0699 #endif
0700                n = itrunc(RealType(tools::log_max_value<RealType>() / 6));
0701 #ifndef BOOST_NO_EXCEPTIONS
0702             }
0703             catch(...)
0704             {
0705                n = (std::numeric_limits<int>::max)();
0706             }
0707 #endif
0708             n = (std::min)(n, 1500);
0709             RealType d = pow(3 + sqrt(RealType(8)), RealType(n));
0710             d = (d + 1 / d) / 2;
0711             RealType b = -1;
0712             RealType c = -d;
0713             int s = 1;
0714
0715             for(int k = 0; k < n; ++k)
0716             {
0717                //
0718                // Check for both convergence and whether the series has gone bad:
0719                //
0720                if(
0721                   (fabs(z) > last_z)     // Series has gone divergent, abort
0722                   || (fabs(val) * tools::epsilon<RealType>() > fabs(c * s * z))  // Convergence!
0723                   || (z * s < 0)         // Series has stopped alternating - all bets are off - abort.
0724                   )
0725                {
0726                   break;
0727                }
0728                c = b - c;
0729                val += c * s * z;
0730                b = (k + n) * (k - n) * b / ((k + RealType(0.5)) * (k + 1));
0731                last_z = fabs(z);
0732                s = -s;
0733                z = y * ( vi - static_cast<RealType>(ii) * z );
0734                vi *= as;
0735                ii += 2;
0736             } // while( true )
0737             RealType err = fabs(c * z) / val;
0738             return std::pair<RealType, RealType>(val * exp( -hs*half<RealType>() ) / (d * root_two_pi<RealType>()), err);
0739          } // RealType owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&)
0740
0741          template<typename RealType, typename Policy>
0742          inline RealType T4_mp(const RealType h, const RealType a, const Policy& pol)
0743          {
0744             BOOST_MATH_STD_USING
0745
0746             const RealType hs = h*h;
0747             const RealType as = -a*a;
0748
0749             unsigned short ii = 1;
0750             RealType ai = constants::one_div_two_pi<RealType>() * a * exp( -0.5*hs*(1.0-as) );
0751             RealType yi = 1.0;
0752             RealType val = 0.0;
0753
0754             RealType lim = boost::math::policies::get_epsilon<RealType, Policy>();
0755
0756             while( true )
0757             {
0758                RealType term = ai*yi;
0759                val += term;
0760                if((yi != 0) && (fabs(val * lim) > fabs(term)))
0761                   break;
0762                ii += 2;
0763                yi = (1.0-hs*yi) / static_cast<RealType>(ii);
0764                ai *= as;
0765                if(ii > (std::min)(1500, (int)policies::get_max_series_iterations<Policy>()))
0766                   policies::raise_evaluation_error("boost::math::owens_t<%1%>", 0, val, pol);
0767             } // while( true )
0768
0769             return val;
0770          } // arg_type owens_t_T4(const arg_type h, const arg_type a, const unsigned short m)
0771
0772
0773          // This routine dispatches the call to one of six subroutines, depending on the values
0774          // of h and a.
0775          // preconditions: h >= 0, 0<=a<=1, ah=a*h
0776          //
0777          // Note there are different versions for different precisions....
0778          template<typename RealType, typename Policy>
0779          inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, std::integral_constant<int, 64> const&)
0780          {
0781             // Simple main case for 64-bit precision or less, this is as per the Patefield-Tandy paper:
0782             BOOST_MATH_STD_USING
0783             //
0784             // Handle some special cases first, these are from
0785             // page 1077 of Owen's original paper:
0786             //
0787             if(h == 0)
0788             {
0789                return atan(a) * constants::one_div_two_pi<RealType>();
0790             }
0791             if(a == 0)
0792             {
0793                return 0;
0794             }
0795             if(a == 1)
0796             {
0797                return owens_t_znorm2(RealType(-h), pol) * owens_t_znorm2(h, pol) / 2;
0798             }
0799             if(a >= tools::max_value<RealType>())
0800             {
0801                return owens_t_znorm2(RealType(fabs(h)), pol);
0802             }
0803             RealType val = 0; // avoid compiler warnings, 0 will be overwritten in any case
0804             const unsigned short icode = owens_t_compute_code(h, a);
0805             const unsigned short m = owens_t_get_order(icode, val /* just a dummy for the type */, pol);
0806             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
0807
0808             // determine the appropriate method, T1 ... T6
0809             switch( meth[icode] )
0810             {
0811             case 1: // T1
0812                val = owens_t_T1(h,a,m,pol);
0813                break;
0814             case 2: // T2
0815                typedef typename policies::precision<RealType, Policy>::type precision_type;
0816                typedef std::integral_constant<bool, (precision_type::value == 0) || (precision_type::value > 64)> tag_type;
0817                val = owens_t_T2(h, a, m, ah, pol, tag_type());
0818                break;
0819             case 3: // T3
0820                val = owens_t_T3(h,a,ah, pol);
0821                break;
0822             case 4: // T4
0823                val = owens_t_T4(h,a,m);
0824                break;
0825             case 5: // T5
0826                val = owens_t_T5(h,a, pol);
0827                break;
0828             case 6: // T6
0829                val = owens_t_T6(h,a, pol);
0830                break;
0831             default:
0832                val = policies::raise_evaluation_error<RealType>("boost::math::owens_t", "selection routine in Owen's T function failed with h = %1%", h, pol);
0833             }
0834             return val;
0835          }
0836
0837          template<typename RealType, typename Policy>
0838          inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const std::integral_constant<int, 65>&)
0839          {
0840             // Arbitrary precision version:
0841             BOOST_MATH_STD_USING
0842             //
0843             // Handle some special cases first, these are from
0844             // page 1077 of Owen's original paper:
0845             //
0846             if(h == 0)
0847             {
0848                return atan(a) * constants::one_div_two_pi<RealType>();
0849             }
0850             if(a == 0)
0851             {
0852                return 0;
0853             }
0854             if(a == 1)
0855             {
0856                return owens_t_znorm2(RealType(-h), pol) * owens_t_znorm2(h, pol) / 2;
0857             }
0858             if(a >= tools::max_value<RealType>())
0859             {
0860                return owens_t_znorm2(RealType(fabs(h)), pol);
0861             }
0862             // Attempt arbitrary precision code, this will throw if it goes wrong:
0863             typedef typename boost::math::policies::normalise<Policy, boost::math::policies::evaluation_error<> >::type forwarding_policy;
0864             std::pair<RealType, RealType> p1(0, tools::max_value<RealType>()), p2(0, tools::max_value<RealType>());
0865             RealType target_precision = policies::get_epsilon<RealType, Policy>() * 1000;
0866             bool have_t1(false), have_t2(false);
0867             if(ah < 3)
0868             {
0869 #ifndef BOOST_NO_EXCEPTIONS
0870                try
0871                {
0872 #endif
0873                   have_t1 = true;
0874                   p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
0875                   if(p1.second < target_precision)
0876                      return p1.first;
0877 #ifndef BOOST_NO_EXCEPTIONS
0878                }
0879                catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
0880 #endif
0881             }
0882             if(ah > 1)
0883             {
0884 #ifndef BOOST_NO_EXCEPTIONS
0885                try
0886                {
0887 #endif
0888                   have_t2 = true;
0889                   p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
0890                   if(p2.second < target_precision)
0891                      return p2.first;
0892 #ifndef BOOST_NO_EXCEPTIONS
0893                }
0894                catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
0895 #endif
0896             }
0897             //
0898             // If we haven't tried T1 yet, do it now - sometimes it succeeds and the number of iterations
0899             // is fairly low compared to T4.
0900             //
0901             if(!have_t1)
0902             {
0903 #ifndef BOOST_NO_EXCEPTIONS
0904                try
0905                {
0906 #endif
0907                   have_t1 = true;
0908                   p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
0909                   if(p1.second < target_precision)
0910                      return p1.first;
0911 #ifndef BOOST_NO_EXCEPTIONS
0912                }
0913                catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
0914 #endif
0915             }
0916             //
0917             // If we haven't tried T2 yet, do it now - sometimes it succeeds and the number of iterations
0918             // is fairly low compared to T4.
0919             //
0920             if(!have_t2)
0921             {
0922 #ifndef BOOST_NO_EXCEPTIONS
0923                try
0924                {
0925 #endif
0926                   have_t2 = true;
0927                   p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
0928                   if(p2.second < target_precision)
0929                      return p2.first;
0930 #ifndef BOOST_NO_EXCEPTIONS
0931                }
0932                catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
0933 #endif
0934             }
0935             //
0936             // OK, nothing left to do but try the most expensive option which is T4,
0937             // this is often slow to converge, but when it does converge it tends to
0938             // be accurate:
0939 #ifndef BOOST_NO_EXCEPTIONS
0940             try
0941             {
0942 #endif
0943                return T4_mp(h, a, pol);
0944 #ifndef BOOST_NO_EXCEPTIONS
0945             }
0946             catch(const boost::math::evaluation_error&){}  // T4 may fail and throw, that's OK
0947 #endif
0948             //
0949             // Now look back at the results from T1 and T2 and see if either gave better
0950             // results than we could get from the 64-bit precision versions.
0951             //
0952             if((std::min)(p1.second, p2.second) < RealType(1e-20))
0953             {
0954                return p1.second < p2.second ? p1.first : p2.first;
0955             }
0956             //
0957             // We give up - no arbitrary precision versions succeeded!
0958             //
0959             return owens_t_dispatch(h, a, ah, pol, std::integral_constant<int, 64>());
0960          } // RealType owens_t_dispatch(RealType h, RealType a, RealType ah)
0961          template<typename RealType, typename Policy>
0962          inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const std::integral_constant<int, 0>&)
0963          {
0964             // We don't know what the precision is until runtime:
0965             if(tools::digits<RealType>() <= 64)
0966                return owens_t_dispatch(h, a, ah, pol, std::integral_constant<int, 64>());
0967             return owens_t_dispatch(h, a, ah, pol, std::integral_constant<int, 65>());
0968          }
0969          template<typename RealType, typename Policy>
0970          inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol)
0971          {
0972             // Figure out the precision and forward to the correct version:
0973             typedef typename policies::precision<RealType, Policy>::type precision_type;
0974             typedef std::integral_constant<int,
0975                precision_type::value <= 0 ? 0 :
0976                precision_type::value <= 64 ? 64 : 65
0977             > tag_type;
0978
0979             return owens_t_dispatch(h, a, ah, pol, tag_type());
0980          }
0981          // compute Owen's T function, T(h,a), for arbitrary values of h and a
0982          template<typename RealType, class Policy>
0983          inline RealType owens_t(RealType h, RealType a, const Policy& pol)
0984          {
0985             BOOST_MATH_STD_USING
0986             // exploit that T(-h,a) == T(h,a)
0987             h = fabs(h);
0988
0989             // Use equation (2) in the paper to remap the arguments
0990             // such that h>=0 and 0<=a<=1 for the call of the actual
0991             // computation routine.
0992
0993             const RealType fabs_a = fabs(a);
0994             const RealType fabs_ah = fabs_a*h;
0995
0996             RealType val = static_cast<RealType>(0.0f); // avoid compiler warnings, 0.0 will be overwritten in any case
0997
0998             if(fabs_a <= 1)
0999             {
1000                val = owens_t_dispatch(h, fabs_a, fabs_ah, pol);
1001             } // if(fabs_a <= 1.0)
1002             else
1003             {
1004                if( h <= RealType(0.67) )
1005                {
1006                   const RealType normh = owens_t_znorm1(h, pol);
1007                   const RealType normah = owens_t_znorm1(fabs_ah, pol);
1008                   val = static_cast<RealType>(1)/static_cast<RealType>(4) - normh*normah -
1009                      owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol);
1010                } // if( h <= 0.67 )
1011                else
1012                {
1013                   const RealType normh = detail::owens_t_znorm2(h, pol);
1014                   const RealType normah = detail::owens_t_znorm2(fabs_ah, pol);
1015                   val = constants::half<RealType>()*(normh+normah) - normh*normah -
1016                      owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol);
1017                } // else [if( h <= 0.67 )]
1018             } // else [if(fabs_a <= 1)]
1019
1020             // exploit that T(h,-a) == -T(h,a)
1021             if(a < 0)
1022             {
1023                return -val;
1024             } // if(a < 0)
1025
1026             return val;
1027          } // RealType owens_t(RealType h, RealType a)
1028
1029          template <class T, class Policy, class tag>
1030          struct owens_t_initializer
1031          {
1032             struct init
1033             {
1034                init()
1035                {
1036                   do_init(tag());
1037                }
1038                template <int N>
1039                static void do_init(const std::integral_constant<int, N>&){}
1040                static void do_init(const std::integral_constant<int, 64>&)
1041                {
1042                   boost::math::owens_t(static_cast<T>(7), static_cast<T>(0.96875), Policy());
1043                   boost::math::owens_t(static_cast<T>(2), static_cast<T>(0.5), Policy());
1044                }
1045                void force_instantiate()const{}
1046             };
1047             static const init initializer;
1048             static void force_instantiate()
1049             {
1050                initializer.force_instantiate();
1051             }
1052          };
1053
1054          template <class T, class Policy, class tag>
1055          const typename owens_t_initializer<T, Policy, tag>::init owens_t_initializer<T, Policy, tag>::initializer;
1056
1057       } // namespace detail
1058
1059       template <class T1, class T2, class Policy>
1060       inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a, const Policy& pol)
1061       {
1062          typedef typename tools::promote_args<T1, T2>::type result_type;
1063          typedef typename policies::evaluation<result_type, Policy>::type value_type;
1064          typedef typename policies::precision<value_type, Policy>::type precision_type;
1065          typedef std::integral_constant<int,
1066             precision_type::value <= 0 ? 0 :
1067             precision_type::value <= 64 ? 64 : 65
1068          > tag_type;
1069
1070          detail::owens_t_initializer<result_type, Policy, tag_type>::force_instantiate();
1071
1072          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%)");
1073       }
1074
1075       template <class T1, class T2>
1076       inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a)
1077       {
1078          return owens_t(h, a, policies::policy<>());
1079       }
1080
1081
1082    } // namespace math
1083 } // namespace boost
1084
1085 #ifdef _MSC_VER
1086 #pragma warning(pop)
1087 #endif
1088
1089 #endif
1090 // EOF