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