boost/math/special_functions/owens_t.hpp

   1 // Copyright Benjamin Sobotta 2012
   2
   3 //  Use, modification and distribution are subject to the
   4 //  Boost Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_OWENS_T_HPP
   8 #define BOOST_OWENS_T_HPP
   9
  10 // Reference:
  11 // Mike Patefield, David Tandy
  12 // FAST AND ACCURATE CALCULATION OF OWEN'S T-FUNCTION
  13 // Journal of Statistical Software, 5 (5), 1-25
  14
  15 #ifdef _MSC_VER
  16 #  pragma once
  17 #endif
  18
  19 #include <boost/config/no_tr1/cmath.hpp>
  20 #include <boost/math/special_functions/erf.hpp>
  21 #include <boost/math/special_functions/expm1.hpp>
  22 #include <boost/throw_exception.hpp>
  23 #include <boost/assert.hpp>
  24 #include <boost/math/constants/constants.hpp>
  25 #include <boost/math/tools/big_constant.hpp>
  26
  27 #include <stdexcept>
  28
  29 namespace boost
  30 {
  31    namespace math
  32    {
  33       namespace detail
  34       {
  35          // owens_t_znorm1(x) = P(-oo<Z<=x)-0.5 with Z being normally distributed.
  36          template<typename RealType>
  37          inline RealType owens_t_znorm1(const RealType x)
  38          {
  39             using namespace boost::math::constants;
  40             return erf(x*one_div_root_two<RealType>())*half<RealType>();
  41          } // RealType owens_t_znorm1(const RealType x)
  42
  43          // owens_t_znorm2(x) = P(x<=Z<oo) with Z being normally distributed.
  44          template<typename RealType>
  45          inline RealType owens_t_znorm2(const RealType x)
  46          {
  47             using namespace boost::math::constants;
  48             return erfc(x*one_div_root_two<RealType>())*half<RealType>();
  49          } // RealType owens_t_znorm2(const RealType x)
  50
  51          // Auxiliary function, it computes an array key that is used to determine
  52          // the specific computation method for Owen's T and the order thereof
  53          // used in owens_t_dispatch.
  54          template<typename RealType>
  55          inline unsigned short owens_t_compute_code(const RealType h, const RealType a)
  56          {
  57             static const RealType hrange[] =
  58             {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};
  59
  60             static const RealType arange[] = {0.025, 0.09, 0.15, 0.36, 0.5, 0.9, 0.99999};
  61             /*
  62             original select array from paper:
  63             1, 1, 2,13,13,13,13,13,13,13,13,16,16,16, 9
  64             1, 2, 2, 3, 3, 5, 5,14,14,15,15,16,16,16, 9
  65             2, 2, 3, 3, 3, 5, 5,15,15,15,15,16,16,16,10
  66             2, 2, 3, 5, 5, 5, 5, 7, 7,16,16,16,16,16,10
  67             2, 3, 3, 5, 5, 6, 6, 8, 8,17,17,17,12,12,11
  68             2, 3, 5, 5, 5, 6, 6, 8, 8,17,17,17,12,12,12
  69             2, 3, 4, 4, 6, 6, 8, 8,17,17,17,17,17,12,12
  70             2, 3, 4, 4, 6, 6,18,18,18,18,17,17,17,12,12
  71             */
  72             // subtract one because the array is written in FORTRAN in mind - in C arrays start @ zero
  73             static const unsigned short select[] =
  74             {
  75                0,    0 ,   1  , 12   ,12 ,  12  , 12  , 12 ,  12  , 12  , 12  , 15  , 15 ,  15  ,  8,
  76                0  ,  1  ,  1   , 2 ,   2   , 4  ,  4  , 13 ,  13  , 14  , 14 ,  15  , 15  , 15  ,  8,
  77                1  ,  1   , 2 ,   2  ,  2  ,  4   , 4  , 14  , 14 ,  14  , 14 ,  15  , 15 ,  15  ,  9,
  78                1  ,  1   , 2 ,   4  ,  4  ,  4   , 4  ,  6  ,  6 ,  15  , 15 ,  15 ,  15 ,  15  ,  9,
  79                1  ,  2   , 2  ,  4  ,  4  ,  5   , 5  ,  7  ,  7  , 16   ,16 ,  16 ,  11 ,  11 ,  10,
  80                1  ,  2   , 4  ,  4   , 4  ,  5   , 5  ,  7  ,  7  , 16  , 16 ,  16 ,  11  , 11 ,  11,
  81                1  ,  2   , 3  ,  3  ,  5  ,  5   , 7  ,  7  , 16 ,  16  , 16 ,  16 ,  16  , 11 ,  11,
  82                1  ,  2   , 3   , 3   , 5  ,  5 ,  17  , 17  , 17 ,  17  , 16 ,  16 ,  16 ,  11 ,  11
  83             };
  84
  85             unsigned short ihint = 14, iaint = 7;
  86             for(unsigned short i = 0; i != 14; i++)
  87             {
  88                if( h <= hrange[i] )
  89                {
  90                   ihint = i;
  91                   break;
  92                }
  93             } // for(unsigned short i = 0; i != 14; i++)
  94
  95             for(unsigned short i = 0; i != 7; i++)
  96             {
  97                if( a <= arange[i] )
  98                {
  99                   iaint = i;
 100                   break;
 101                }
 102             } // for(unsigned short i = 0; i != 7; i++)
 103
 104             // interprete select array as 8x15 matrix
 105             return select[iaint*15 + ihint];
 106
 107          } // unsigned short owens_t_compute_code(const RealType h, const RealType a)
 108
 109          template<typename RealType>
 110          inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const mpl::int_<53>&)
 111          {
 112             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
 113
 114             BOOST_ASSERT(icode<18);
 115
 116             return ord[icode];
 117          } // unsigned short owens_t_get_order(const unsigned short icode, RealType, mpl::int<53> const&)
 118
 119          template<typename RealType>
 120          inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const mpl::int_<64>&)
 121         {
 122            // method ================>>>       {1, 1, 1, 1, 1,  1,  1,  1,  2,  2,  2,  3, 4,  4,  4,  4,  5, 6}
 123            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
 124
 125           BOOST_ASSERT(icode<18);
 126
 127           return ord[icode];
 128         } // unsigned short owens_t_get_order(const unsigned short icode, RealType, mpl::int<64> const&)
 129
 130          template<typename RealType, typename Policy>
 131          inline unsigned short owens_t_get_order(const unsigned short icode, RealType r, const Policy&)
 132          {
 133             typedef typename policies::precision<RealType, Policy>::type precision_type;
 134             typedef typename mpl::if_<
 135                mpl::or_<
 136                   mpl::less_equal<precision_type, mpl::int_<0> >,
 137                   mpl::greater<precision_type, mpl::int_<53> >
 138                >,
 139                mpl::int_<64>,
 140                mpl::int_<53>
 141             >::type tag_type;
 142
 143             return owens_t_get_order_imp(icode, r, tag_type());
 144          }
 145
 146          // compute the value of Owen's T function with method T1 from the reference paper
 147          template<typename RealType>
 148          inline RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m)
 149          {
 150             BOOST_MATH_STD_USING
 151             using namespace boost::math::constants;
 152
 153             const RealType hs = -h*h*half<RealType>();
 154             const RealType dhs = exp( hs );
 155             const RealType as = a*a;
 156
 157             unsigned short j=1;
 158             RealType jj = 1;
 159             RealType aj = a * one_div_two_pi<RealType>();
 160             RealType dj = expm1( hs );
 161             RealType gj = hs*dhs;
 162
 163             RealType val = atan( a ) * one_div_two_pi<RealType>();
 164
 165             while( true )
 166             {
 167                val += dj*aj/jj;
 168
 169                if( m <= j )
 170                   break;
 171
 172                j++;
 173                jj += static_cast<RealType>(2);
 174                aj *= as;
 175                dj = gj - dj;
 176                gj *= hs / static_cast<RealType>(j);
 177             } // while( true )
 178
 179             return val;
 180          } // RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m)
 181
 182          // compute the value of Owen's T function with method T2 from the reference paper
 183          template<typename RealType, class Policy>
 184          inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy&, const mpl::false_&)
 185          {
 186             BOOST_MATH_STD_USING
 187             using namespace boost::math::constants;
 188
 189             const unsigned short maxii = m+m+1;
 190             const RealType hs = h*h;
 191             const RealType as = -a*a;
 192             const RealType y = static_cast<RealType>(1) / hs;
 193
 194             unsigned short ii = 1;
 195             RealType val = 0;
 196             RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
 197             RealType z = owens_t_znorm1(ah)/h;
 198
 199             while( true )
 200             {
 201                val += z;
 202                if( maxii <= ii )
 203                {
 204                   val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
 205                   break;
 206                } // if( maxii <= ii )
 207                z = y * ( vi - static_cast<RealType>(ii) * z );
 208                vi *= as;
 209                ii += 2;
 210             } // while( true )
 211
 212             return val;
 213          } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)
 214
 215          // compute the value of Owen's T function with method T3 from the reference paper
 216          template<typename RealType>
 217          inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const mpl::int_<53>&)
 218          {
 219             BOOST_MATH_STD_USING
 220             using namespace boost::math::constants;
 221
 222       const unsigned short m = 20;
 223
 224             static const RealType c2[] =
 225             {
 226                0.99999999999999987510,
 227                -0.99999999999988796462,      0.99999999998290743652,
 228                -0.99999999896282500134,      0.99999996660459362918,
 229                -0.99999933986272476760,      0.99999125611136965852,
 230                -0.99991777624463387686,      0.99942835555870132569,
 231                -0.99697311720723000295,      0.98751448037275303682,
 232                -0.95915857980572882813,      0.89246305511006708555,
 233                -0.76893425990463999675,      0.58893528468484693250,
 234                -0.38380345160440256652,      0.20317601701045299653,
 235                -0.82813631607004984866E-01,  0.24167984735759576523E-01,
 236                -0.44676566663971825242E-02,  0.39141169402373836468E-03
 237             };
 238
 239             const RealType as = a*a;
 240             const RealType hs = h*h;
 241             const RealType y = static_cast<RealType>(1)/hs;
 242
 243             RealType ii = 1;
 244             unsigned short i = 0;
 245             RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
 246             RealType zi = owens_t_znorm1(ah)/h;
 247             RealType val = 0;
 248
 249             while( true )
 250             {
 251                BOOST_ASSERT(i < 21);
 252                val += zi*c2[i];
 253                if( m <= i ) // if( m < i+1 )
 254                {
 255                   val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
 256                   break;
 257                } // if( m < i )
 258                zi = y * (ii*zi - vi);
 259                vi *= as;
 260                ii += 2;
 261                i++;
 262             } // while( true )
 263
 264             return val;
 265          } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah)
 266
 267         // compute the value of Owen's T function with method T3 from the reference paper
 268         template<class RealType>
 269         inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const mpl::int_<64>&)
 270         {
 271           BOOST_MATH_STD_USING
 272           using namespace boost::math::constants;
 273
 274           const unsigned short m = 30;
 275
 276           static const RealType c2[] =
 277           {
 278              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999999999729978162447266851932041876728736094298092917625009873),
 279              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99999999999999999999467056379678391810626533251885323416799874878563998732905968),
 280              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999824849349313270659391127814689133077036298754586814091034842536),
 281              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999999997703859616213643405880166422891953033591551179153879839440241685),
 282              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999998394883415238173334565554173013941245103172035286759201504179038147),
 283              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999993063616095509371081203145247992197457263066869044528823599399470977),
 284              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999999797336340409464429599229870590160411238245275855903767652432017766116267),
 285              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999999999574958412069046680119051639753412378037565521359444170241346845522403274),
 286              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999933226234193375324943920160947158239076786103108097456617750134812033362048),
 287              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999188923242461073033481053037468263536806742737922476636768006622772762168467),
 288              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999992195143483674402853783549420883055129680082932629160081128947764415749728967),
 289              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999993935137206712830997921913316971472227199741857386575097250553105958772041501),
 290              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99996135597690552745362392866517133091672395614263398912807169603795088421057688716),
 291              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99979556366513946026406788969630293820987757758641211293079784585126692672425362469),
 292              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.999092789629617100153486251423850590051366661947344315423226082520411961968929483),
 293              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.996593837411918202119308620432614600338157335862888580671450938858935084316004769854),
 294              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.98910017138386127038463510314625339359073956513420458166238478926511821146316469589567),
 295              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.970078558040693314521331982203762771512160168582494513347846407314584943870399016019),
 296              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.92911438683263187495758525500033707204091967947532160289872782771388170647150321633673),
 297              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.8542058695956156057286980736842905011429254735181323743367879525470479126968822863),
 298              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.73796526033030091233118357742803709382964420335559408722681794195743240930748630755),
 299              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.58523469882837394570128599003785154144164680587615878645171632791404210655891158),
 300              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.415997776145676306165661663581868460503874205343014196580122174949645271353372263),
 301              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.2588210875241943574388730510317252236407805082485246378222935376279663808416534365),
 302              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.1375535825163892648504646951500265585055789019410617565727090346559210218472356689),
 303              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.0607952766325955730493900985022020434830339794955745989150270485056436844239206648),
 304              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0216337683299871528059836483840390514275488679530797294557060229266785853764115),
 305              BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.00593405693455186729876995814181203900550014220428843483927218267309209471516256),
 306              BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0011743414818332946510474576182739210553333860106811865963485870668929503649964142),
 307              BOOST_MATH_BIG_CONSTANT(RealType, 260, -1.489155613350368934073453260689881330166342484405529981510694514036264969925132e-4),
 308              BOOST_MATH_BIG_CONSTANT(RealType, 260, 9.072354320794357587710929507988814669454281514268844884841547607134260303118208e-6)
 309           };
 310
 311           const RealType as = a*a;
 312           const RealType hs = h*h;
 313           const RealType y = 1 / hs;
 314
 315           RealType ii = 1;
 316           unsigned short i = 0;
 317           RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();
 318           RealType zi = owens_t_znorm1(ah)/h;
 319           RealType val = 0;
 320
 321           while( true )
 322           {
 323               BOOST_ASSERT(i < 31);
 324               val += zi*c2[i];
 325               if( m <= i ) // if( m < i+1 )
 326               {
 327                 val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();
 328                 break;
 329               } // if( m < i )
 330               zi = y * (ii*zi - vi);
 331               vi *= as;
 332               ii += 2;
 333               i++;
 334           } // while( true )
 335
 336           return val;
 337         } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah)
 338
 339         template<class RealType, class Policy>
 340         inline RealType owens_t_T3(const RealType h, const RealType a, const RealType ah, const Policy&)
 341         {
 342             typedef typename policies::precision<RealType, Policy>::type precision_type;
 343             typedef typename mpl::if_<
 344                mpl::or_<
 345                   mpl::less_equal<precision_type, mpl::int_<0> >,
 346                   mpl::greater<precision_type, mpl::int_<53> >
 347                >,
 348                mpl::int_<64>,
 349                mpl::int_<53>
 350             >::type tag_type;
 351
 352             return owens_t_T3_imp(h, a, ah, tag_type());
 353         }
 354
 355          // compute the value of Owen's T function with method T4 from the reference paper
 356          template<typename RealType>
 357          inline RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m)
 358          {
 359             BOOST_MATH_STD_USING
 360             using namespace boost::math::constants;
 361
 362             const unsigned short maxii = m+m+1;
 363             const RealType hs = h*h;
 364             const RealType as = -a*a;
 365
 366             unsigned short ii = 1;
 367             RealType ai = a * exp( -hs*(static_cast<RealType>(1)-as)*half<RealType>() ) * one_div_two_pi<RealType>();
 368             RealType yi = 1;
 369             RealType val = 0;
 370
 371             while( true )
 372             {
 373                val += ai*yi;
 374                if( maxii <= ii )
 375                   break;
 376                ii += 2;
 377                yi = (static_cast<RealType>(1)-hs*yi) / static_cast<RealType>(ii);
 378                ai *= as;
 379             } // while( true )
 380
 381             return val;
 382          } // RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m)
 383
 384          // compute the value of Owen's T function with method T5 from the reference paper
 385          template<typename RealType>
 386          inline RealType owens_t_T5_imp(const RealType h, const RealType a, const mpl::int_<53>&)
 387          {
 388             BOOST_MATH_STD_USING
 389             /*
 390                NOTICE:
 391                - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre
 392                  polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre
 393                  quadrature, because T5(h,a,m) contains only x^2 terms.
 394                - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor
 395                  of 1/(2*pi) according to T5(h,a,m).
 396              */
 397
 398             const unsigned short m = 13;
 399             static const RealType pts[] = {0.35082039676451715489E-02,
 400                0.31279042338030753740E-01,  0.85266826283219451090E-01,
 401                0.16245071730812277011,      0.25851196049125434828,
 402                0.36807553840697533536,      0.48501092905604697475,
 403                0.60277514152618576821,      0.71477884217753226516,
 404                0.81475510988760098605,      0.89711029755948965867,
 405                0.95723808085944261843,      0.99178832974629703586};
 406             static const RealType wts[] = { 0.18831438115323502887E-01,
 407                0.18567086243977649478E-01,  0.18042093461223385584E-01,
 408                0.17263829606398753364E-01,  0.16243219975989856730E-01,
 409                0.14994592034116704829E-01,  0.13535474469662088392E-01,
 410                0.11886351605820165233E-01,  0.10070377242777431897E-01,
 411                0.81130545742299586629E-02,  0.60419009528470238773E-02,
 412                0.38862217010742057883E-02,  0.16793031084546090448E-02};
 413
 414             const RealType as = a*a;
 415             const RealType hs = -h*h*boost::math::constants::half<RealType>();
 416
 417             RealType val = 0;
 418             for(unsigned short i = 0; i < m; ++i)
 419             {
 420                BOOST_ASSERT(i < 13);
 421                const RealType r = static_cast<RealType>(1) + as*pts[i];
 422                val += wts[i] * exp( hs*r ) / r;
 423             } // for(unsigned short i = 0; i < m; ++i)
 424
 425             return val*a;
 426          } // RealType owens_t_T5(const RealType h, const RealType a)
 427
 428         // compute the value of Owen's T function with method T5 from the reference paper
 429         template<typename RealType>
 430         inline RealType owens_t_T5_imp(const RealType h, const RealType a, const mpl::int_<64>&)
 431         {
 432           BOOST_MATH_STD_USING
 433             /*
 434               NOTICE:
 435               - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre
 436               polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre
 437               quadrature, because T5(h,a,m) contains only x^2 terms.
 438               - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor
 439               of 1/(2*pi) according to T5(h,a,m).
 440             */
 441
 442           const unsigned short m = 19;
 443           static const RealType pts[] = {
 444                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0016634282895983227941),
 445                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.014904509242697054183),
 446                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.04103478879005817919),
 447                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.079359853513391511008),
 448                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.1288612130237615133),
 449                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.18822336642448518856),
 450                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.25586876186122962384),
 451                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.32999972011807857222),
 452                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.40864620815774761438),
 453                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.48971819306044782365),
 454                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.57106118513245543894),
 455                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.6505134942981533829),
 456                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.72596367859928091618),
 457                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.79540665919549865924),
 458                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.85699701386308739244),
 459                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.90909804422384697594),
 460                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.95032536436570154409),
 461                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.97958418733152273717),
 462                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.99610366384229088321)
 463           };
 464           static const RealType wts[] = {
 465                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012975111395684900835),
 466                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012888764187499150078),
 467                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012716644398857307844),
 468                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012459897461364705691),
 469                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012120231988292330388),
 470                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011699908404856841158),
 471                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011201723906897224448),
 472                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.010628993848522759853),
 473                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0099855296835573320047),
 474                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0092756136096132857933),
 475                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0085039700881139589055),
 476                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0076757344408814561254),
 477                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0067964187616556459109),
 478                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.005871875456524750363),
 479                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0049082589542498110071),
 480                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0039119870792519721409),
 481                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0028897090921170700834),
 482                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0018483371329504443947),
 483                BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.00079623320100438873578)
 484           };
 485
 486           const RealType as = a*a;
 487           const RealType hs = -h*h*boost::math::constants::half<RealType>();
 488
 489           RealType val = 0;
 490           for(unsigned short i = 0; i < m; ++i)
 491             {
 492               BOOST_ASSERT(i < 19);
 493               const RealType r = 1 + as*pts[i];
 494               val += wts[i] * exp( hs*r ) / r;
 495             } // for(unsigned short i = 0; i < m; ++i)
 496
 497           return val*a;
 498         } // RealType owens_t_T5(const RealType h, const RealType a)
 499
 500         template<class RealType, class Policy>
 501         inline RealType owens_t_T5(const RealType h, const RealType a, const Policy&)
 502         {
 503             typedef typename policies::precision<RealType, Policy>::type precision_type;
 504             typedef typename mpl::if_<
 505                mpl::or_<
 506                   mpl::less_equal<precision_type, mpl::int_<0> >,
 507                   mpl::greater<precision_type, mpl::int_<53> >
 508                >,
 509                mpl::int_<64>,
 510                mpl::int_<53>
 511             >::type tag_type;
 512
 513             return owens_t_T5_imp(h, a, tag_type());
 514         }
 515
 516
 517          // compute the value of Owen's T function with method T6 from the reference paper
 518          template<typename RealType>
 519          inline RealType owens_t_T6(const RealType h, const RealType a)
 520          {
 521             BOOST_MATH_STD_USING
 522             using namespace boost::math::constants;
 523
 524             const RealType normh = owens_t_znorm2( h );
 525             const RealType y = static_cast<RealType>(1) - a;
 526             const RealType r = atan2(y, static_cast<RealType>(1 + a) );
 527
 528             RealType val = normh * ( static_cast<RealType>(1) - normh ) * half<RealType>();
 529
 530             if( r != 0 )
 531                val -= r * exp( -y*h*h*half<RealType>()/r ) * one_div_two_pi<RealType>();
 532
 533             return val;
 534          } // RealType owens_t_T6(const RealType h, const RealType a, const unsigned short m)
 535
 536          template <class T, class Policy>
 537          std::pair<T, T> owens_t_T1_accelerated(T h, T a, const Policy& pol)
 538          {
 539             //
 540             // This is the same series as T1, but:
 541             // * The Taylor series for atan has been combined with that for T1,
 542             //   reducing but not eliminating cancellation error.
 543             // * The resulting alternating series is then accelerated using method 1
 544             //   from H. Cohen, F. Rodriguez Villegas, D. Zagier,
 545             //   "Convergence acceleration of alternating series", Bonn, (1991).
 546             //
 547             BOOST_MATH_STD_USING
 548             static const char* function = "boost::math::owens_t<%1%>(%1%, %1%)";
 549             T half_h_h = h * h / 2;
 550             T a_pow = a;
 551             T aa = a * a;
 552             T exp_term = exp(-h * h / 2);
 553             T one_minus_dj_sum = exp_term;
 554             T sum = a_pow * exp_term;
 555             T dj_pow = exp_term;
 556             T term = sum;
 557             T abs_err;
 558             int j = 1;
 559
 560             //
 561             // Normally with this form of series acceleration we can calculate
 562             // up front how many terms will be required - based on the assumption
 563             // that each term decreases in size by a factor of 3.  However,
 564             // that assumption does not apply here, as the underlying T1 series can
 565             // go quite strongly divergent in the early terms, before strongly
 566             // converging later.  Various "guestimates" have been tried to take account
 567             // of this, but they don't always work.... so instead set "n" to the
 568             // largest value that won't cause overflow later, and abort iteration
 569             // when the last accelerated term was small enough...
 570             //
 571             int n;
 572             try
 573             {
 574                n = itrunc(T(tools::log_max_value<T>() / 6));
 575             }
 576             catch(...)
 577             {
 578                n = (std::numeric_limits<int>::max)();
 579             }
 580             n = (std::min)(n, 1500);
 581             T d = pow(3 + sqrt(T(8)), n);
 582             d = (d + 1 / d) / 2;
 583             T b = -1;
 584             T c = -d;
 585             c = b - c;
 586             sum *= c;
 587             b = -n * n * b * 2;
 588             abs_err = ldexp(fabs(sum), -tools::digits<T>());
 589
 590             while(j < n)
 591             {
 592                a_pow *= aa;
 593                dj_pow *= half_h_h / j;
 594                one_minus_dj_sum += dj_pow;
 595                term = one_minus_dj_sum * a_pow / (2 * j + 1);
 596                c = b - c;
 597                sum += c * term;
 598                abs_err += ldexp((std::max)(T(fabs(sum)), T(fabs(c*term))), -tools::digits<T>());
 599                b = (j + n) * (j - n) * b / ((j + T(0.5)) * (j + 1));
 600                ++j;
 601                //
 602                // Include an escape route to prevent calculating too many terms:
 603                //
 604                if((j > 10) && (fabs(sum * tools::epsilon<T>()) > fabs(c * term)))
 605                   break;
 606             }
 607             abs_err += fabs(c * term);
 608             if(sum < 0)  // sum must always be positive, if it's negative something really bad has happend:
 609                policies::raise_evaluation_error(function, 0, T(0), pol);
 610             return std::pair<T, T>((sum / d) / boost::math::constants::two_pi<T>(), abs_err / sum);
 611          }
 612
 613          template<typename RealType, class Policy>
 614          inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy&, const mpl::true_&)
 615          {
 616             BOOST_MATH_STD_USING
 617             using namespace boost::math::constants;
 618
 619             const unsigned short maxii = m+m+1;
 620             const RealType hs = h*h;
 621             const RealType as = -a*a;
 622             const RealType y = static_cast<RealType>(1) / hs;
 623
 624             unsigned short ii = 1;
 625             RealType val = 0;
 626             RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
 627             RealType z = owens_t_znorm1(ah)/h;
 628             RealType last_z = fabs(z);
 629             RealType lim = policies::get_epsilon<RealType, Policy>();
 630
 631             while( true )
 632             {
 633                val += z;
 634                //
 635                // This series stops converging after a while, so put a limit
 636                // on how far we go before returning our best guess:
 637                //
 638                if((fabs(lim * val) > fabs(z)) || ((ii > maxii) && (fabs(z) > last_z)) || (z == 0))
 639                {
 640                   val *= exp( -hs*half<RealType>() ) / root_two_pi<RealType>();
 641                   break;
 642                } // if( maxii <= ii )
 643                last_z = fabs(z);
 644                z = y * ( vi - static_cast<RealType>(ii) * z );
 645                vi *= as;
 646                ii += 2;
 647             } // while( true )
 648
 649             return val;
 650          } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)
 651
 652          template<typename RealType, class Policy>
 653          inline std::pair<RealType, RealType> owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&)
 654          {
 655             //
 656             // This is the same series as T2, but with acceleration applied.
 657             // Note that we have to be *very* careful to check that nothing bad
 658             // has happened during evaluation - this series will go divergent
 659             // and/or fail to alternate at a drop of a hat! :-(
 660             //
 661             BOOST_MATH_STD_USING
 662             using namespace boost::math::constants;
 663
 664             const RealType hs = h*h;
 665             const RealType as = -a*a;
 666             const RealType y = static_cast<RealType>(1) / hs;
 667
 668             unsigned short ii = 1;
 669             RealType val = 0;
 670             RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
 671             RealType z = boost::math::detail::owens_t_znorm1(ah)/h;
 672             RealType last_z = fabs(z);
 673
 674             //
 675             // Normally with this form of series acceleration we can calculate
 676             // up front how many terms will be required - based on the assumption
 677             // that each term decreases in size by a factor of 3.  However,
 678             // that assumption does not apply here, as the underlying T1 series can
 679             // go quite strongly divergent in the early terms, before strongly
 680             // converging later.  Various "guestimates" have been tried to take account
 681             // of this, but they don't always work.... so instead set "n" to the
 682             // largest value that won't cause overflow later, and abort iteration
 683             // when the last accelerated term was small enough...
 684             //
 685             int n;
 686             try
 687             {
 688                n = itrunc(RealType(tools::log_max_value<RealType>() / 6));
 689             }
 690             catch(...)
 691             {
 692                n = (std::numeric_limits<int>::max)();
 693             }
 694             n = (std::min)(n, 1500);
 695             RealType d = pow(3 + sqrt(RealType(8)), n);
 696             d = (d + 1 / d) / 2;
 697             RealType b = -1;
 698             RealType c = -d;
 699             int s = 1;
 700
 701             for(int k = 0; k < n; ++k)
 702             {
 703                //
 704                // Check for both convergence and whether the series has gone bad:
 705                //
 706                if(
 707                   (fabs(z) > last_z)     // Series has gone divergent, abort
 708                   || (fabs(val) * tools::epsilon<RealType>() > fabs(c * s * z))  // Convergence!
 709                   || (z * s < 0)         // Series has stopped alternating - all bets are off - abort.
 710                   )
 711                {
 712                   break;
 713                }
 714                c = b - c;
 715                val += c * s * z;
 716                b = (k + n) * (k - n) * b / ((k + RealType(0.5)) * (k + 1));
 717                last_z = fabs(z);
 718                s = -s;
 719                z = y * ( vi - static_cast<RealType>(ii) * z );
 720                vi *= as;
 721                ii += 2;
 722             } // while( true )
 723             RealType err = fabs(c * z) / val;
 724             return std::pair<RealType, RealType>(val * exp( -hs*half<RealType>() ) / (d * root_two_pi<RealType>()), err);
 725          } // RealType owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&)
 726
 727          template<typename RealType, typename Policy>
 728          inline RealType T4_mp(const RealType h, const RealType a, const Policy& pol)
 729          {
 730             BOOST_MATH_STD_USING
 731
 732             const RealType hs = h*h;
 733             const RealType as = -a*a;
 734
 735             unsigned short ii = 1;
 736             RealType ai = constants::one_div_two_pi<RealType>() * a * exp( -0.5*hs*(1.0-as) );
 737             RealType yi = 1.0;
 738             RealType val = 0.0;
 739
 740             RealType lim = boost::math::policies::get_epsilon<RealType, Policy>();
 741
 742             while( true )
 743             {
 744                RealType term = ai*yi;
 745                val += term;
 746                if((yi != 0) && (fabs(val * lim) > fabs(term)))
 747                   break;
 748                ii += 2;
 749                yi = (1.0-hs*yi) / static_cast<RealType>(ii);
 750                ai *= as;
 751                if(ii > (std::min)(1500, (int)policies::get_max_series_iterations<Policy>()))
 752                   policies::raise_evaluation_error("boost::math::owens_t<%1%>", 0, val, pol);
 753             } // while( true )
 754
 755             return val;
 756          } // arg_type owens_t_T4(const arg_type h, const arg_type a, const unsigned short m)
 757
 758
 759          // This routine dispatches the call to one of six subroutines, depending on the values
 760          // of h and a.
 761          // preconditions: h >= 0, 0<=a<=1, ah=a*h
 762          //
 763          // Note there are different versions for different precisions....
 764          template<typename RealType, typename Policy>
 765          inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, mpl::int_<64> const&)
 766          {
 767             // Simple main case for 64-bit precision or less, this is as per the Patefield-Tandy paper:
 768             BOOST_MATH_STD_USING
 769             //
 770             // Handle some special cases first, these are from
 771             // page 1077 of Owen's original paper:
 772             //
 773             if(h == 0)
 774             {
 775                return atan(a) * constants::one_div_two_pi<RealType>();
 776             }
 777             if(a == 0)
 778             {
 779                return 0;
 780             }
 781             if(a == 1)
 782             {
 783                return owens_t_znorm2(RealType(-h)) * owens_t_znorm2(h) / 2;
 784             }
 785             if(a >= tools::max_value<RealType>())
 786             {
 787                return owens_t_znorm2(RealType(fabs(h)));
 788             }
 789             RealType val = 0; // avoid compiler warnings, 0 will be overwritten in any case
 790             const unsigned short icode = owens_t_compute_code(h, a);
 791             const unsigned short m = owens_t_get_order(icode, val /* just a dummy for the type */, pol);
 792             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
 793
 794             // determine the appropriate method, T1 ... T6
 795             switch( meth[icode] )
 796             {
 797             case 1: // T1
 798                val = owens_t_T1(h,a,m);
 799                break;
 800             case 2: // T2
 801                typedef typename policies::precision<RealType, Policy>::type precision_type;
 802                typedef mpl::bool_<(precision_type::value == 0) || (precision_type::value > 64)> tag_type;
 803                val = owens_t_T2(h, a, m, ah, pol, tag_type());
 804                break;
 805             case 3: // T3
 806                val = owens_t_T3(h,a,ah, pol);
 807                break;
 808             case 4: // T4
 809                val = owens_t_T4(h,a,m);
 810                break;
 811             case 5: // T5
 812                val = owens_t_T5(h,a, pol);
 813                break;
 814             case 6: // T6
 815                val = owens_t_T6(h,a);
 816                break;
 817             default:
 818                BOOST_THROW_EXCEPTION(std::logic_error("selection routine in Owen's T function failed"));
 819             }
 820             return val;
 821          }
 822
 823          template<typename RealType, typename Policy>
 824          inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const mpl::int_<65>&)
 825          {
 826             // Arbitrary precision version:
 827             BOOST_MATH_STD_USING
 828             //
 829             // Handle some special cases first, these are from
 830             // page 1077 of Owen's original paper:
 831             //
 832             if(h == 0)
 833             {
 834                return atan(a) * constants::one_div_two_pi<RealType>();
 835             }
 836             if(a == 0)
 837             {
 838                return 0;
 839             }
 840             if(a == 1)
 841             {
 842                return owens_t_znorm2(RealType(-h)) * owens_t_znorm2(h) / 2;
 843             }
 844             if(a >= tools::max_value<RealType>())
 845             {
 846                return owens_t_znorm2(RealType(fabs(h)));
 847             }
 848             // Attempt arbitrary precision code, this will throw if it goes wrong:
 849             typedef typename boost::math::policies::normalise<Policy, boost::math::policies::evaluation_error<> >::type forwarding_policy;
 850             std::pair<RealType, RealType> p1(0, tools::max_value<RealType>()), p2(0, tools::max_value<RealType>());
 851             RealType target_precision = policies::get_epsilon<RealType, Policy>() * 1000;
 852             bool have_t1(false), have_t2(false);
 853             if(ah < 3)
 854             {
 855                try
 856                {
 857                   have_t1 = true;
 858                   p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
 859                   if(p1.second < target_precision)
 860                      return p1.first;
 861                }
 862                catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
 863             }
 864             if(ah > 1)
 865             {
 866                try
 867                {
 868                   have_t2 = true;
 869                   p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
 870                   if(p2.second < target_precision)
 871                      return p2.first;
 872                }
 873                catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
 874             }
 875             //
 876             // If we haven't tried T1 yet, do it now - sometimes it succeeds and the number of iterations
 877             // is fairly low compared to T4.
 878             //
 879             if(!have_t1)
 880             {
 881                try
 882                {
 883                   have_t1 = true;
 884                   p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
 885                   if(p1.second < target_precision)
 886                      return p1.first;
 887                }
 888                catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
 889             }
 890             //
 891             // If we haven't tried T2 yet, do it now - sometimes it succeeds and the number of iterations
 892             // is fairly low compared to T4.
 893             //
 894             if(!have_t2)
 895             {
 896                try
 897                {
 898                   have_t2 = true;
 899                   p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
 900                   if(p2.second < target_precision)
 901                      return p2.first;
 902                }
 903                catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
 904             }
 905             //
 906             // OK, nothing left to do but try the most expensive option which is T4,
 907             // this is often slow to converge, but when it does converge it tends to
 908             // be accurate:
 909             try
 910             {
 911                return T4_mp(h, a, pol);
 912             }
 913             catch(const boost::math::evaluation_error&){}  // T4 may fail and throw, that's OK
 914             //
 915             // Now look back at the results from T1 and T2 and see if either gave better
 916             // results than we could get from the 64-bit precision versions.
 917             //
 918             if((std::min)(p1.second, p2.second) < 1e-20)
 919             {
 920                return p1.second < p2.second ? p1.first : p2.first;
 921             }
 922             //
 923             // We give up - no arbitrary precision versions succeeded!
 924             //
 925             return owens_t_dispatch(h, a, ah, pol, mpl::int_<64>());
 926          } // RealType owens_t_dispatch(RealType h, RealType a, RealType ah)
 927          template<typename RealType, typename Policy>
 928          inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const mpl::int_<0>&)
 929          {
 930             // We don't know what the precision is until runtime:
 931             if(tools::digits<RealType>() <= 64)
 932                return owens_t_dispatch(h, a, ah, pol, mpl::int_<64>());
 933             return owens_t_dispatch(h, a, ah, pol, mpl::int_<65>());
 934          }
 935          template<typename RealType, typename Policy>
 936          inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol)
 937          {
 938             // Figure out the precision and forward to the correct version:
 939             typedef typename policies::precision<RealType, Policy>::type precision_type;
 940             typedef typename mpl::if_c<
 941                precision_type::value == 0,
 942                mpl::int_<0>,
 943                typename mpl::if_c<
 944                   precision_type::value <= 64,
 945                   mpl::int_<64>,
 946                   mpl::int_<65>
 947                >::type
 948             >::type tag_type;
 949             return owens_t_dispatch(h, a, ah, pol, tag_type());
 950          }
 951          // compute Owen's T function, T(h,a), for arbitrary values of h and a
 952          template<typename RealType, class Policy>
 953          inline RealType owens_t(RealType h, RealType a, const Policy& pol)
 954          {
 955             BOOST_MATH_STD_USING
 956             // exploit that T(-h,a) == T(h,a)
 957             h = fabs(h);
 958
 959             // Use equation (2) in the paper to remap the arguments
 960             // such that h>=0 and 0<=a<=1 for the call of the actual
 961             // computation routine.
 962
 963             const RealType fabs_a = fabs(a);
 964             const RealType fabs_ah = fabs_a*h;
 965
 966             RealType val = 0.0; // avoid compiler warnings, 0.0 will be overwritten in any case
 967
 968             if(fabs_a <= 1)
 969             {
 970                val = owens_t_dispatch(h, fabs_a, fabs_ah, pol);
 971             } // if(fabs_a <= 1.0)
 972             else
 973             {
 974                if( h <= 0.67 )
 975                {
 976                   const RealType normh = owens_t_znorm1(h);
 977                   const RealType normah = owens_t_znorm1(fabs_ah);
 978                   val = static_cast<RealType>(1)/static_cast<RealType>(4) - normh*normah -
 979                      owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol);
 980                } // if( h <= 0.67 )
 981                else
 982                {
 983                   const RealType normh = detail::owens_t_znorm2(h);
 984                   const RealType normah = detail::owens_t_znorm2(fabs_ah);
 985                   val = constants::half<RealType>()*(normh+normah) - normh*normah -
 986                      owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol);
 987                } // else [if( h <= 0.67 )]
 988             } // else [if(fabs_a <= 1)]
 989
 990             // exploit that T(h,-a) == -T(h,a)
 991             if(a < 0)
 992             {
 993                return -val;
 994             } // if(a < 0)
 995
 996             return val;
 997          } // RealType owens_t(RealType h, RealType a)
 998
 999          template <class T, class Policy, class tag>
1000          struct owens_t_initializer
1001          {
1002             struct init
1003             {
1004                init()
1005                {
1006                   do_init(tag());
1007                }
1008                template <int N>
1009                static void do_init(const mpl::int_<N>&){}
1010                static void do_init(const mpl::int_<64>&)
1011                {
1012                   boost::math::owens_t(static_cast<T>(7), static_cast<T>(0.96875), Policy());
1013                   boost::math::owens_t(static_cast<T>(2), static_cast<T>(0.5), Policy());
1014                }
1015                void force_instantiate()const{}
1016             };
1017             static const init initializer;
1018             static void force_instantiate()
1019             {
1020                initializer.force_instantiate();
1021             }
1022          };
1023
1024          template <class T, class Policy, class tag>
1025          const typename owens_t_initializer<T, Policy, tag>::init owens_t_initializer<T, Policy, tag>::initializer;
1026
1027       } // namespace detail
1028
1029       template <class T1, class T2, class Policy>
1030       inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a, const Policy& pol)
1031       {
1032          typedef typename tools::promote_args<T1, T2>::type result_type;
1033          typedef typename policies::evaluation<result_type, Policy>::type value_type;
1034          typedef typename policies::precision<value_type, Policy>::type precision_type;
1035          typedef typename mpl::if_c<
1036                precision_type::value == 0,
1037                mpl::int_<0>,
1038                typename mpl::if_c<
1039                   precision_type::value <= 64,
1040                   mpl::int_<64>,
1041                   mpl::int_<65>
1042                >::type
1043             >::type tag_type;
1044
1045          detail::owens_t_initializer<result_type, Policy, tag_type>::force_instantiate();
1046
1047          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%)");
1048       }
1049
1050       template <class T1, class T2>
1051       inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a)
1052       {
1053          return owens_t(h, a, policies::policy<>());
1054       }
1055
1056
1057    } // namespace math
1058 } // namespace boost
1059
1060 #endif
1061 // EOF