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