1068 lines · plain
1// Copyright Benjamin Sobotta 20122 3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_OWENS_T_HPP8#define BOOST_OWENS_T_HPP9 10// Reference:11// Mike Patefield, David Tandy12// FAST AND ACCURATE CALCULATION OF OWEN'S T-FUNCTION13// Journal of Statistical Software, 5 (5), 1-2514 15#ifdef _MSC_VER16# pragma once17#endif18 19#include <boost/math/special_functions/math_fwd.hpp>20#include <boost/math/special_functions/erf.hpp>21#include <boost/math/special_functions/expm1.hpp>22#include <boost/math/tools/throw_exception.hpp>23#include <boost/math/tools/assert.hpp>24#include <boost/math/constants/constants.hpp>25#include <boost/math/tools/big_constant.hpp>26 27#include <stdexcept>28#include <cmath>29 30#ifdef _MSC_VER31#pragma warning(push)32#pragma warning(disable:4127)33#endif34 35#if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)36//37// This is the only way we can avoid38// 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_header43#endif44 45namespace boost46{47 namespace math48 {49 namespace detail50 {51 // owens_t_znorm1(x) = P(-oo<Z<=x)-0.5 with Z being normally distributed.52 template<typename RealType, class Policy>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>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 determine68 // the specific computation method for Owen's T and the order thereof69 // used in owens_t_dispatch.70 template<typename RealType>71 inline unsigned short owens_t_compute_code(const RealType h, const RealType a)72 {73 // LCOV_EXCL_START74 static const RealType hrange[] =75 { 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 };76 77 static const RealType arange[] = { 0.025f, 0.09f, 0.15f, 0.36f, 0.5f, 0.9f, 0.99999f };78 /*79 original select array from paper:80 1, 1, 2,13,13,13,13,13,13,13,13,16,16,16, 981 1, 2, 2, 3, 3, 5, 5,14,14,15,15,16,16,16, 982 2, 2, 3, 3, 3, 5, 5,15,15,15,15,16,16,16,1083 2, 2, 3, 5, 5, 5, 5, 7, 7,16,16,16,16,16,1084 2, 3, 3, 5, 5, 6, 6, 8, 8,17,17,17,12,12,1185 2, 3, 5, 5, 5, 6, 6, 8, 8,17,17,17,12,12,1286 2, 3, 4, 4, 6, 6, 8, 8,17,17,17,17,17,12,1287 2, 3, 4, 4, 6, 6,18,18,18,18,17,17,17,12,1288 */ 89 // subtract one because the array is written in FORTRAN in mind - in C arrays start @ zero90 static const unsigned short select[] =91 {92 0, 0 , 1 , 12 ,12 , 12 , 12 , 12 , 12 , 12 , 12 , 15 , 15 , 15 , 8,93 0 , 1 , 1 , 2 , 2 , 4 , 4 , 13 , 13 , 14 , 14 , 15 , 15 , 15 , 8,94 1 , 1 , 2 , 2 , 2 , 4 , 4 , 14 , 14 , 14 , 14 , 15 , 15 , 15 , 9,95 1 , 1 , 2 , 4 , 4 , 4 , 4 , 6 , 6 , 15 , 15 , 15 , 15 , 15 , 9,96 1 , 2 , 2 , 4 , 4 , 5 , 5 , 7 , 7 , 16 ,16 , 16 , 11 , 11 , 10,97 1 , 2 , 4 , 4 , 4 , 5 , 5 , 7 , 7 , 16 , 16 , 16 , 11 , 11 , 11,98 1 , 2 , 3 , 3 , 5 , 5 , 7 , 7 , 16 , 16 , 16 , 16 , 16 , 11 , 11,99 1 , 2 , 3 , 3 , 5 , 5 , 17 , 17 , 17 , 17 , 16 , 16 , 16 , 11 , 11100 };101 // LCOV_EXCL_STOP102 103 unsigned short ihint = 14, iaint = 7;104 for(unsigned short i = 0; i != 14; i++)105 {106 if( h <= hrange[i] )107 {108 ihint = i;109 break;110 }111 } // for(unsigned short i = 0; i != 14; i++)112 113 for(unsigned short i = 0; i != 7; i++)114 {115 if( a <= arange[i] )116 {117 iaint = i;118 break;119 }120 } // for(unsigned short i = 0; i != 7; i++)121 122 // interpret select array as 8x15 matrix123 BOOST_MATH_ASSERT(iaint * 15 + ihint < (int)(sizeof(select) / sizeof(select[0])));124 return select[iaint*15 + ihint];125 126 } // unsigned short owens_t_compute_code(const RealType h, const RealType a)127 128 template<typename RealType>129 inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const std::integral_constant<int, 53>&)130 {131 // LCOV_EXCL_START132 static const unsigned short ord[] = {2, 3, 4, 5, 7, 10, 12, 18, 10, 20, 30, 0, 4, 7, 8, 20, 0, 0}; // 18 entries133 // LCOV_EXCL_STOP134 135 BOOST_MATH_ASSERT(icode<18);136 137 return ord[icode];138 } // unsigned short owens_t_get_order(const unsigned short icode, RealType, std::integral_constant<int, 53> const&)139 140 template<typename RealType>141 inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const std::integral_constant<int, 64>&)142 {143 // method ================>>> {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6}144 // LCOV_EXCL_START145 static const unsigned short ord[] = {3, 4, 5, 6, 8, 11, 13, 19, 10, 20, 30, 0, 7, 10, 11, 23, 0, 0}; // 18 entries146 // LCOV_EXCL_STOP147 148 BOOST_MATH_ASSERT(icode<18);149 150 return ord[icode];151 } // unsigned short owens_t_get_order(const unsigned short icode, RealType, std::integral_constant<int, 64> const&)152 153 template<typename RealType, typename Policy>154 inline unsigned short owens_t_get_order(const unsigned short icode, RealType r, const Policy&)155 {156 typedef typename policies::precision<RealType, Policy>::type precision_type;157 typedef std::integral_constant<int,158 precision_type::value <= 0 ? 64 :159 precision_type::value <= 53 ? 53 : 64160 > tag_type;161 162 return owens_t_get_order_imp(icode, r, tag_type());163 }164 165 // compute the value of Owen's T function with method T1 from the reference paper166 template<typename RealType, typename Policy>167 inline RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m, const Policy& pol)168 {169 BOOST_MATH_STD_USING170 using namespace boost::math::constants;171 172 const RealType hs = -h*h*half<RealType>();173 const RealType dhs = exp( hs );174 const RealType as = a*a;175 176 unsigned short j=1;177 RealType jj = 1;178 RealType aj = a * one_div_two_pi<RealType>();179 RealType dj = boost::math::expm1( hs, pol);180 RealType gj = hs*dhs;181 182 RealType val = atan( a ) * one_div_two_pi<RealType>();183 184 while( true )185 {186 val += dj*aj/jj;187 188 if( m <= j )189 break;190 191 j++;192 jj += static_cast<RealType>(2);193 aj *= as;194 dj = gj - dj;195 gj *= hs / static_cast<RealType>(j);196 } // while( true )197 198 return val;199 } // RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m)200 201 // compute the value of Owen's T function with method T2 from the reference paper202 template<typename RealType, class Policy>203 inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy& pol, const std::false_type&)204 {205 BOOST_MATH_STD_USING206 using namespace boost::math::constants;207 208 const unsigned short maxii = m+m+1;209 const RealType hs = h*h;210 const RealType as = -a*a;211 const RealType y = static_cast<RealType>(1) / hs;212 213 unsigned short ii = 1;214 RealType val = 0;215 RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();216 RealType z = owens_t_znorm1(ah, pol)/h;217 218 while( true )219 {220 val += z;221 if( maxii <= ii )222 {223 val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();224 break;225 } // if( maxii <= ii )226 z = y * ( vi - static_cast<RealType>(ii) * z );227 vi *= as;228 ii += 2;229 } // while( true )230 231 return val;232 } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)233 234 // compute the value of Owen's T function with method T3 from the reference paper235 template<typename RealType, class Policy>236 inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const std::integral_constant<int, 53>&, const Policy& pol)237 {238 BOOST_MATH_STD_USING239 using namespace boost::math::constants;240 241 const unsigned short m = 20;242 243 // LCOV_EXCL_START244 static const RealType c2[] =245 {246 static_cast<RealType>(0.99999999999999987510),247 static_cast<RealType>(-0.99999999999988796462), static_cast<RealType>(0.99999999998290743652),248 static_cast<RealType>(-0.99999999896282500134), static_cast<RealType>(0.99999996660459362918),249 static_cast<RealType>(-0.99999933986272476760), static_cast<RealType>(0.99999125611136965852),250 static_cast<RealType>(-0.99991777624463387686), static_cast<RealType>(0.99942835555870132569),251 static_cast<RealType>(-0.99697311720723000295), static_cast<RealType>(0.98751448037275303682),252 static_cast<RealType>(-0.95915857980572882813), static_cast<RealType>(0.89246305511006708555),253 static_cast<RealType>(-0.76893425990463999675), static_cast<RealType>(0.58893528468484693250),254 static_cast<RealType>(-0.38380345160440256652), static_cast<RealType>(0.20317601701045299653),255 static_cast<RealType>(-0.82813631607004984866E-01), static_cast<RealType>(0.24167984735759576523E-01),256 static_cast<RealType>(-0.44676566663971825242E-02), static_cast<RealType>(0.39141169402373836468E-03)257 };258 // LCOV_EXCL_STOP259 260 const RealType as = a*a;261 const RealType hs = h*h;262 const RealType y = static_cast<RealType>(1)/hs;263 264 RealType ii = 1;265 unsigned short i = 0;266 RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();267 RealType zi = owens_t_znorm1(ah, pol)/h;268 RealType val = 0;269 270 while( true )271 {272 BOOST_MATH_ASSERT(i < 21);273 val += zi*c2[i];274 if( m <= i ) // if( m < i+1 )275 {276 val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();277 break;278 } // if( m < i )279 zi = y * (ii*zi - vi);280 vi *= as;281 ii += 2;282 i++;283 } // while( true )284 285 return val;286 } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah)287 288 // compute the value of Owen's T function with method T3 from the reference paper289 template<class RealType, class Policy>290 inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const std::integral_constant<int, 64>&, const Policy& pol)291 {292 BOOST_MATH_STD_USING293 using namespace boost::math::constants;294 295 const unsigned short m = 30;296 297 // LCOV_EXCL_START298 static const RealType c2[] =299 {300 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999999999729978162447266851932041876728736094298092917625009873),301 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99999999999999999999467056379678391810626533251885323416799874878563998732905968),302 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999824849349313270659391127814689133077036298754586814091034842536),303 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999999997703859616213643405880166422891953033591551179153879839440241685),304 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999998394883415238173334565554173013941245103172035286759201504179038147),305 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999993063616095509371081203145247992197457263066869044528823599399470977),306 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999999797336340409464429599229870590160411238245275855903767652432017766116267),307 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999999999574958412069046680119051639753412378037565521359444170241346845522403274),308 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999933226234193375324943920160947158239076786103108097456617750134812033362048),309 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999188923242461073033481053037468263536806742737922476636768006622772762168467),310 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999992195143483674402853783549420883055129680082932629160081128947764415749728967),311 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999993935137206712830997921913316971472227199741857386575097250553105958772041501),312 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99996135597690552745362392866517133091672395614263398912807169603795088421057688716),313 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99979556366513946026406788969630293820987757758641211293079784585126692672425362469),314 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.999092789629617100153486251423850590051366661947344315423226082520411961968929483),315 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.996593837411918202119308620432614600338157335862888580671450938858935084316004769854),316 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.98910017138386127038463510314625339359073956513420458166238478926511821146316469589567),317 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.970078558040693314521331982203762771512160168582494513347846407314584943870399016019),318 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.92911438683263187495758525500033707204091967947532160289872782771388170647150321633673),319 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.8542058695956156057286980736842905011429254735181323743367879525470479126968822863),320 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.73796526033030091233118357742803709382964420335559408722681794195743240930748630755),321 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.58523469882837394570128599003785154144164680587615878645171632791404210655891158),322 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.415997776145676306165661663581868460503874205343014196580122174949645271353372263),323 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.2588210875241943574388730510317252236407805082485246378222935376279663808416534365),324 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.1375535825163892648504646951500265585055789019410617565727090346559210218472356689),325 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.0607952766325955730493900985022020434830339794955745989150270485056436844239206648),326 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0216337683299871528059836483840390514275488679530797294557060229266785853764115),327 BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.00593405693455186729876995814181203900550014220428843483927218267309209471516256),328 BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0011743414818332946510474576182739210553333860106811865963485870668929503649964142),329 BOOST_MATH_BIG_CONSTANT(RealType, 260, -1.489155613350368934073453260689881330166342484405529981510694514036264969925132e-4),330 BOOST_MATH_BIG_CONSTANT(RealType, 260, 9.072354320794357587710929507988814669454281514268844884841547607134260303118208e-6)331 };332 // LCOV_EXCL_STOP333 334 const RealType as = a*a;335 const RealType hs = h*h;336 const RealType y = 1 / hs;337 338 RealType ii = 1;339 unsigned short i = 0;340 RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>();341 RealType zi = owens_t_znorm1(ah, pol)/h;342 RealType val = 0;343 344 while( true )345 {346 BOOST_MATH_ASSERT(i < 31);347 val += zi*c2[i];348 if( m <= i ) // if( m < i+1 )349 {350 val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>();351 break;352 } // if( m < i )353 zi = y * (ii*zi - vi);354 vi *= as;355 ii += 2;356 i++;357 } // while( true )358 359 return val;360 } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah)361 362 template<class RealType, class Policy>363 inline RealType owens_t_T3(const RealType h, const RealType a, const RealType ah, const Policy& pol)364 {365 typedef typename policies::precision<RealType, Policy>::type precision_type;366 typedef std::integral_constant<int,367 precision_type::value <= 0 ? 64 :368 precision_type::value <= 53 ? 53 : 64369 > tag_type;370 371 return owens_t_T3_imp(h, a, ah, tag_type(), pol);372 }373 374 // compute the value of Owen's T function with method T4 from the reference paper375 template<typename RealType>376 inline RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m)377 {378 BOOST_MATH_STD_USING379 using namespace boost::math::constants;380 381 const unsigned short maxii = m+m+1;382 const RealType hs = h*h;383 const RealType as = -a*a;384 385 unsigned short ii = 1;386 RealType ai = a * exp( -hs*(static_cast<RealType>(1)-as)*half<RealType>() ) * one_div_two_pi<RealType>();387 RealType yi = 1;388 RealType val = 0;389 390 while( true )391 {392 val += ai*yi;393 if( maxii <= ii )394 break;395 ii += 2;396 yi = (static_cast<RealType>(1)-hs*yi) / static_cast<RealType>(ii);397 ai *= as;398 } // while( true )399 400 return val;401 } // RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m)402 403 // compute the value of Owen's T function with method T5 from the reference paper404 template<typename RealType>405 inline RealType owens_t_T5_imp(const RealType h, const RealType a, const std::integral_constant<int, 53>&)406 {407 BOOST_MATH_STD_USING408 /*409 NOTICE:410 - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre411 polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre412 quadrature, because T5(h,a,m) contains only x^2 terms.413 - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor414 of 1/(2*pi) according to T5(h,a,m).415 */416 417 const unsigned short m = 13;418 // LCOV_EXCL_START419 static const RealType pts[] = {420 static_cast<RealType>(0.35082039676451715489E-02),421 static_cast<RealType>(0.31279042338030753740E-01), static_cast<RealType>(0.85266826283219451090E-01),422 static_cast<RealType>(0.16245071730812277011), static_cast<RealType>(0.25851196049125434828),423 static_cast<RealType>(0.36807553840697533536), static_cast<RealType>(0.48501092905604697475),424 static_cast<RealType>(0.60277514152618576821), static_cast<RealType>(0.71477884217753226516),425 static_cast<RealType>(0.81475510988760098605), static_cast<RealType>(0.89711029755948965867),426 static_cast<RealType>(0.95723808085944261843), static_cast<RealType>(0.99178832974629703586) };427 static const RealType wts[] = { 428 static_cast<RealType>(0.18831438115323502887E-01),429 static_cast<RealType>(0.18567086243977649478E-01), static_cast<RealType>(0.18042093461223385584E-01),430 static_cast<RealType>(0.17263829606398753364E-01), static_cast<RealType>(0.16243219975989856730E-01),431 static_cast<RealType>(0.14994592034116704829E-01), static_cast<RealType>(0.13535474469662088392E-01),432 static_cast<RealType>(0.11886351605820165233E-01), static_cast<RealType>(0.10070377242777431897E-01),433 static_cast<RealType>(0.81130545742299586629E-02), static_cast<RealType>(0.60419009528470238773E-02),434 static_cast<RealType>(0.38862217010742057883E-02), static_cast<RealType>(0.16793031084546090448E-02) };435 436 const RealType as = a*a;437 const RealType hs = -h*h*boost::math::constants::half<RealType>();438 // LCOV_EXCL_STOP439 440 RealType val = 0;441 for(unsigned short i = 0; i < m; ++i)442 {443 BOOST_MATH_ASSERT(i < 13);444 const RealType r = static_cast<RealType>(1) + as*pts[i];445 val += wts[i] * exp( hs*r ) / r;446 } // for(unsigned short i = 0; i < m; ++i)447 448 return val*a;449 } // RealType owens_t_T5(const RealType h, const RealType a)450 451 // compute the value of Owen's T function with method T5 from the reference paper452 template<typename RealType>453 inline RealType owens_t_T5_imp(const RealType h, const RealType a, const std::integral_constant<int, 64>&)454 {455 BOOST_MATH_STD_USING456 /*457 NOTICE:458 - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre459 polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre460 quadrature, because T5(h,a,m) contains only x^2 terms.461 - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor462 of 1/(2*pi) according to T5(h,a,m).463 */464 465 const unsigned short m = 19;466 // LCOV_EXCL_START467 static const RealType pts[] = {468 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0016634282895983227941),469 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.014904509242697054183),470 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.04103478879005817919),471 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.079359853513391511008),472 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.1288612130237615133),473 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.18822336642448518856),474 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.25586876186122962384),475 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.32999972011807857222),476 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.40864620815774761438),477 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.48971819306044782365),478 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.57106118513245543894),479 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.6505134942981533829),480 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.72596367859928091618),481 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.79540665919549865924),482 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.85699701386308739244),483 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.90909804422384697594),484 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.95032536436570154409),485 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.97958418733152273717),486 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.99610366384229088321)487 };488 static const RealType wts[] = {489 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012975111395684900835),490 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012888764187499150078),491 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012716644398857307844),492 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012459897461364705691),493 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012120231988292330388),494 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011699908404856841158),495 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011201723906897224448),496 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.010628993848522759853),497 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0099855296835573320047),498 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0092756136096132857933),499 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0085039700881139589055),500 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0076757344408814561254),501 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0067964187616556459109),502 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.005871875456524750363),503 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0049082589542498110071),504 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0039119870792519721409),505 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0028897090921170700834),506 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0018483371329504443947),507 BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.00079623320100438873578)508 };509 // LCOV_EXCL_STOP510 511 const RealType as = a*a;512 const RealType hs = -h*h*boost::math::constants::half<RealType>();513 514 RealType val = 0;515 for(unsigned short i = 0; i < m; ++i)516 {517 BOOST_MATH_ASSERT(i < 19);518 const RealType r = 1 + as*pts[i];519 val += wts[i] * exp( hs*r ) / r;520 } // for(unsigned short i = 0; i < m; ++i)521 522 return val*a;523 } // RealType owens_t_T5(const RealType h, const RealType a)524 525 template<class RealType, class Policy>526 inline RealType owens_t_T5(const RealType h, const RealType a, const Policy&)527 {528 typedef typename policies::precision<RealType, Policy>::type precision_type;529 typedef std::integral_constant<int,530 precision_type::value <= 0 ? 64 :531 precision_type::value <= 53 ? 53 : 64532 > tag_type;533 534 return owens_t_T5_imp(h, a, tag_type());535 }536 537 538 // compute the value of Owen's T function with method T6 from the reference paper539 template<typename RealType, class Policy>540 inline RealType owens_t_T6(const RealType h, const RealType a, const Policy& pol)541 {542 BOOST_MATH_STD_USING543 using namespace boost::math::constants;544 545 const RealType normh = owens_t_znorm2(h, pol);546 const RealType y = static_cast<RealType>(1) - a;547 const RealType r = atan2(y, static_cast<RealType>(1 + a) );548 549 RealType val = normh * ( static_cast<RealType>(1) - normh ) * half<RealType>();550 551 if( r != 0 )552 val -= r * exp( -y*h*h*half<RealType>()/r ) * one_div_two_pi<RealType>();553 554 return val;555 } // RealType owens_t_T6(const RealType h, const RealType a, const unsigned short m)556 557 template <class T, class Policy>558 std::pair<T, T> owens_t_T1_accelerated(T h, T a, const Policy& pol)559 {560 //561 // This is the same series as T1, but:562 // * The Taylor series for atan has been combined with that for T1, 563 // reducing but not eliminating cancellation error.564 // * The resulting alternating series is then accelerated using method 1565 // from H. Cohen, F. Rodriguez Villegas, D. Zagier, 566 // "Convergence acceleration of alternating series", Bonn, (1991).567 //568 BOOST_MATH_STD_USING569 static const char* function = "boost::math::owens_t<%1%>(%1%, %1%)";570 T half_h_h = h * h / 2;571 T a_pow = a;572 T aa = a * a;573 T exp_term = exp(-h * h / 2);574 T one_minus_dj_sum = exp_term; 575 T sum = a_pow * exp_term;576 T dj_pow = exp_term;577 T term = sum;578 T abs_err;579 int j = 1;580 581 //582 // Normally with this form of series acceleration we can calculate583 // up front how many terms will be required - based on the assumption584 // that each term decreases in size by a factor of 3. However,585 // that assumption does not apply here, as the underlying T1 series can 586 // go quite strongly divergent in the early terms, before strongly587 // converging later. Various "guesstimates" have been tried to take account588 // of this, but they don't always work.... so instead set "n" to the 589 // largest value that won't cause overflow later, and abort iteration590 // when the last accelerated term was small enough...591 //592 int n;593#ifndef BOOST_MATH_NO_EXCEPTIONS594 try595 {596#endif597 n = itrunc(T(tools::log_max_value<T>() / 6));598#ifndef BOOST_MATH_NO_EXCEPTIONS599 }600 catch(...)601 {602 n = (std::numeric_limits<int>::max)();603 }604#endif605 n = (std::min)(n, 1500);606 T d = pow(3 + sqrt(T(8)), T(n));607 d = (d + 1 / d) / 2;608 T b = -1;609 T c = -d;610 c = b - c;611 sum *= c;612 b = -n * n * b * 2;613 abs_err = ldexp(fabs(sum), -tools::digits<T>());614 615 while(j < n)616 {617 a_pow *= aa;618 dj_pow *= half_h_h / j;619 one_minus_dj_sum += dj_pow;620 term = one_minus_dj_sum * a_pow / (2 * j + 1);621 c = b - c;622 sum += c * term;623 abs_err += ldexp((std::max)(T(fabs(sum)), T(fabs(c*term))), -tools::digits<T>());624 b = (j + n) * (j - n) * b / ((j + T(0.5)) * (j + 1));625 ++j;626 //627 // Include an escape route to prevent calculating too many terms:628 //629 if((j > 10) && (fabs(sum * tools::epsilon<T>()) > fabs(c * term)))630 break;631 }632 abs_err += fabs(c * term);633 if(sum < 0) // sum must always be positive, if it's negative something really bad has happened:634 policies::raise_evaluation_error(function, 0, T(0), pol);635 return std::pair<T, T>((sum / d) / boost::math::constants::two_pi<T>(), abs_err / sum);636 }637 638 template<typename RealType, class Policy>639 inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy& pol, const std::true_type&)640 {641 BOOST_MATH_STD_USING642 using namespace boost::math::constants;643 644 const unsigned short maxii = m+m+1;645 const RealType hs = h*h;646 const RealType as = -a*a;647 const RealType y = static_cast<RealType>(1) / hs;648 649 unsigned short ii = 1;650 RealType val = 0;651 RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();652 RealType z = owens_t_znorm1(ah, pol)/h;653 RealType last_z = fabs(z);654 RealType lim = policies::get_epsilon<RealType, Policy>();655 656 while( true )657 {658 val += z;659 //660 // This series stops converging after a while, so put a limit661 // on how far we go before returning our best guess:662 //663 if((fabs(lim * val) > fabs(z)) || ((ii > maxii) && (fabs(z) > last_z)) || (z == 0))664 {665 val *= exp( -hs*half<RealType>() ) / root_two_pi<RealType>();666 break;667 } // if( maxii <= ii )668 last_z = fabs(z);669 z = y * ( vi - static_cast<RealType>(ii) * z );670 vi *= as;671 ii += 2;672 } // while( true )673 674 return val;675 } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)676 677 template<typename RealType, class Policy>678 inline std::pair<RealType, RealType> owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy& pol)679 {680 //681 // This is the same series as T2, but with acceleration applied.682 // Note that we have to be *very* careful to check that nothing bad683 // has happened during evaluation - this series will go divergent684 // and/or fail to alternate at a drop of a hat! :-(685 //686 BOOST_MATH_STD_USING687 using namespace boost::math::constants;688 689 const RealType hs = h*h;690 const RealType as = -a*a;691 const RealType y = static_cast<RealType>(1) / hs;692 693 unsigned short ii = 1;694 RealType val = 0;695 RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();696 RealType z = boost::math::detail::owens_t_znorm1(ah, pol)/h;697 RealType last_z = fabs(z);698 699 //700 // Normally with this form of series acceleration we can calculate701 // up front how many terms will be required - based on the assumption702 // that each term decreases in size by a factor of 3. However,703 // that assumption does not apply here, as the underlying T1 series can 704 // go quite strongly divergent in the early terms, before strongly705 // converging later. Various "guesstimates" have been tried to take account706 // of this, but they don't always work.... so instead set "n" to the 707 // largest value that won't cause overflow later, and abort iteration708 // when the last accelerated term was small enough...709 //710 int n;711#ifndef BOOST_MATH_NO_EXCEPTIONS712 try713 {714#endif715 n = itrunc(RealType(tools::log_max_value<RealType>() / 6));716#ifndef BOOST_MATH_NO_EXCEPTIONS717 }718 catch(...)719 {720 n = (std::numeric_limits<int>::max)();721 }722#endif723 n = (std::min)(n, 1500);724 RealType d = pow(3 + sqrt(RealType(8)), RealType(n));725 d = (d + 1 / d) / 2;726 RealType b = -1;727 RealType c = -d;728 int s = 1;729 730 for(int k = 0; k < n; ++k)731 {732 //733 // Check for both convergence and whether the series has gone bad:734 //735 if(736 (fabs(z) > last_z) // Series has gone divergent, abort737 || (fabs(val) * tools::epsilon<RealType>() > fabs(c * s * z)) // Convergence!738 || (z * s < 0) // Series has stopped alternating - all bets are off - abort.739 )740 {741 break;742 }743 c = b - c;744 val += c * s * z;745 b = (k + n) * (k - n) * b / ((k + RealType(0.5)) * (k + 1));746 last_z = fabs(z);747 s = -s;748 z = y * ( vi - static_cast<RealType>(ii) * z );749 vi *= as;750 ii += 2;751 } // while( true )752 RealType err = fabs(c * z) / val;753 return std::pair<RealType, RealType>(val * exp( -hs*half<RealType>() ) / (d * root_two_pi<RealType>()), err);754 } // RealType owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&)755 756 template<typename RealType, typename Policy>757 inline RealType T4_mp(const RealType h, const RealType a, const Policy& pol)758 {759 BOOST_MATH_STD_USING760 761 const RealType hs = h*h;762 const RealType as = -a*a;763 764 unsigned short ii = 1;765 RealType ai = constants::one_div_two_pi<RealType>() * a * exp( -0.5*hs*(1.0-as) );766 RealType yi = 1.0;767 RealType val = 0.0;768 769 RealType lim = boost::math::policies::get_epsilon<RealType, Policy>();770 771 while( true )772 {773 RealType term = ai*yi;774 val += term;775 if((yi != 0) && (fabs(val * lim) > fabs(term)))776 break;777 ii += 2;778 yi = (1.0-hs*yi) / static_cast<RealType>(ii);779 ai *= as;780 if(ii > (std::min)(1500, (int)policies::get_max_series_iterations<Policy>()))781 policies::raise_evaluation_error("boost::math::owens_t<%1%>", 0, val, pol);782 } // while( true )783 784 return val;785 } // arg_type owens_t_T4(const arg_type h, const arg_type a, const unsigned short m)786 787 788 // This routine dispatches the call to one of six subroutines, depending on the values789 // of h and a.790 // preconditions: h >= 0, 0<=a<=1, ah=a*h791 //792 // Note there are different versions for different precisions....793 template<typename RealType, typename Policy>794 inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, std::integral_constant<int, 64> const&)795 {796 // Simple main case for 64-bit precision or less, this is as per the Patefield-Tandy paper:797 BOOST_MATH_STD_USING798 //799 // Handle some special cases first, these are from800 // page 1077 of Owen's original paper:801 //802 if(h == 0)803 {804 return atan(a) * constants::one_div_two_pi<RealType>();805 }806 if(a == 0)807 {808 return 0;809 }810 if(a == 1)811 {812 return owens_t_znorm2(RealType(-h), pol) * owens_t_znorm2(h, pol) / 2;813 }814 // Rationale: when a>1 we call this routine with 1/a:815 BOOST_MATH_ASSERT(a <= 1);816 RealType val = 0; // avoid compiler warnings, 0 will be overwritten in any case817 const unsigned short icode = owens_t_compute_code(h, a);818 const unsigned short m = owens_t_get_order(icode, val /* just a dummy for the type */, pol);819 static const unsigned short meth[] = {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6}; // 18 entries820 BOOST_MATH_ASSERT(icode < sizeof(meth) / sizeof(meth[0]));821 822 // determine the appropriate method, T1 ... T6823 switch( meth[icode] )824 {825 case 1: // T1826 val = owens_t_T1(h,a,m,pol);827 break;828 case 2: // T2829 typedef typename policies::precision<RealType, Policy>::type precision_type;830 typedef std::integral_constant<bool, (precision_type::value == 0) || (precision_type::value > 64)> tag_type;831 val = owens_t_T2(h, a, m, ah, pol, tag_type());832 break;833 case 3: // T3834 val = owens_t_T3(h,a,ah, pol);835 break;836 case 4: // T4837 val = owens_t_T4(h,a,m);838 break;839 case 5: // T5840 val = owens_t_T5(h,a, pol);841 break;842 case 6: // T6843 val = owens_t_T6(h,a, pol);844 break;845 }846 return val;847 }848 849 template<typename RealType, typename Policy>850 inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const std::integral_constant<int, 65>&)851 {852 // Arbitrary precision version:853 BOOST_MATH_STD_USING854 //855 // Handle some special cases first, these are from856 // page 1077 of Owen's original paper:857 //858 if(h == 0)859 {860 return atan(a) * constants::one_div_two_pi<RealType>();861 }862 if(a == 0)863 {864 return 0;865 }866 if(a == 1)867 {868 return owens_t_znorm2(RealType(-h), pol) * owens_t_znorm2(h, pol) / 2;869 }870 if(a >= tools::max_value<RealType>())871 {872 return owens_t_znorm2(RealType(fabs(h)), pol);873 }874 // Attempt arbitrary precision code, this will throw if it goes wrong:875 typedef typename boost::math::policies::normalise<Policy, boost::math::policies::evaluation_error<> >::type forwarding_policy;876 std::pair<RealType, RealType> p1(0, tools::max_value<RealType>()), p2(0, tools::max_value<RealType>());877 RealType target_precision = policies::get_epsilon<RealType, Policy>() * 1000;878 bool have_t1(false), have_t2(false);879 if(ah < 3)880 {881#ifndef BOOST_MATH_NO_EXCEPTIONS882 try883 {884#endif885 have_t1 = true;886 p1 = owens_t_T1_accelerated(h, a, forwarding_policy());887 if(p1.second < target_precision)888 return p1.first;889#ifndef BOOST_MATH_NO_EXCEPTIONS890 }891 catch(const boost::math::evaluation_error&){} // T1 may fail and throw, that's OK892#endif893 }894 if(ah > 1)895 {896#ifndef BOOST_MATH_NO_EXCEPTIONS897 try898 {899#endif900 have_t2 = true;901 p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());902 if(p2.second < target_precision)903 return p2.first;904#ifndef BOOST_MATH_NO_EXCEPTIONS905 }906 catch(const boost::math::evaluation_error&){} // T2 may fail and throw, that's OK907#endif908 }909 //910 // If we haven't tried T1 yet, do it now - sometimes it succeeds and the number of iterations911 // is fairly low compared to T4.912 //913 if(!have_t1)914 {915#ifndef BOOST_MATH_NO_EXCEPTIONS916 try917 {918#endif919 have_t1 = true;920 p1 = owens_t_T1_accelerated(h, a, forwarding_policy());921 if(p1.second < target_precision)922 return p1.first;923#ifndef BOOST_MATH_NO_EXCEPTIONS924 }925 catch(const boost::math::evaluation_error&){} // T1 may fail and throw, that's OK926#endif927 }928 //929 // If we haven't tried T2 yet, do it now - sometimes it succeeds and the number of iterations930 // is fairly low compared to T4.931 //932 if(!have_t2)933 {934#ifndef BOOST_MATH_NO_EXCEPTIONS935 try936 {937#endif938 have_t2 = true;939 p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());940 if(p2.second < target_precision)941 return p2.first;942#ifndef BOOST_MATH_NO_EXCEPTIONS943 }944 catch(const boost::math::evaluation_error&){} // T2 may fail and throw, that's OK945#endif946 }947 //948 // OK, nothing left to do but try the most expensive option which is T4,949 // this is often slow to converge, but when it does converge it tends to950 // be accurate:951#ifndef BOOST_MATH_NO_EXCEPTIONS952 try953 {954#endif955 return T4_mp(h, a, pol);956#ifndef BOOST_MATH_NO_EXCEPTIONS957 }958 catch(const boost::math::evaluation_error&){} // T4 may fail and throw, that's OK959#endif960 //961 // Now look back at the results from T1 and T2 and see if either gave better962 // results than we could get from the 64-bit precision versions.963 //964 if((std::min)(p1.second, p2.second) < RealType(1e-20))965 {966 return p1.second < p2.second ? p1.first : p2.first;967 }968 //969 // We give up - no arbitrary precision versions succeeded!970 //971 return owens_t_dispatch(h, a, ah, pol, std::integral_constant<int, 64>());972 } // RealType owens_t_dispatch(RealType h, RealType a, RealType ah)973 template<typename RealType, typename Policy>974 inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const std::integral_constant<int, 0>&)975 {976 // We don't know what the precision is until runtime:977 if(tools::digits<RealType>() <= 64)978 return owens_t_dispatch(h, a, ah, pol, std::integral_constant<int, 64>());979 return owens_t_dispatch(h, a, ah, pol, std::integral_constant<int, 65>());980 }981 template<typename RealType, typename Policy>982 inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol)983 {984 // Figure out the precision and forward to the correct version:985 typedef typename policies::precision<RealType, Policy>::type precision_type;986 typedef std::integral_constant<int,987 precision_type::value <= 0 ? 0 :988 precision_type::value <= 64 ? 64 : 65989 > tag_type;990 991 return owens_t_dispatch(h, a, ah, pol, tag_type());992 }993 // compute Owen's T function, T(h,a), for arbitrary values of h and a994 template<typename RealType, class Policy>995 inline RealType owens_t(RealType h, RealType a, const Policy& pol)996 {997 BOOST_MATH_STD_USING998 // exploit that T(-h,a) == T(h,a)999 h = fabs(h);1000 1001 // Use equation (2) in the paper to remap the arguments1002 // such that h>=0 and 0<=a<=1 for the call of the actual1003 // computation routine.1004 1005 const RealType fabs_a = fabs(a);1006 const RealType fabs_ah = fabs_a*h;1007 1008 RealType val = static_cast<RealType>(0.0f); // avoid compiler warnings, 0.0 will be overwritten in any case1009 1010 if(fabs_a <= 1)1011 {1012 val = owens_t_dispatch(h, fabs_a, fabs_ah, pol);1013 } // if(fabs_a <= 1.0)1014 else 1015 {1016 if( h <= RealType(0.67) )1017 {1018 const RealType normh = owens_t_znorm1(h, pol);1019 const RealType normah = owens_t_znorm1(fabs_ah, pol);1020 val = static_cast<RealType>(1)/static_cast<RealType>(4) - normh*normah -1021 owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol);1022 } // if( h <= 0.67 )1023 else1024 {1025 const RealType normh = detail::owens_t_znorm2(h, pol);1026 const RealType normah = detail::owens_t_znorm2(fabs_ah, pol);1027 val = constants::half<RealType>()*(normh+normah) - normh*normah -1028 owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol);1029 } // else [if( h <= 0.67 )]1030 } // else [if(fabs_a <= 1)]1031 1032 // exploit that T(h,-a) == -T(h,a)1033 if(a < 0)1034 {1035 return -val;1036 } // if(a < 0)1037 1038 return val;1039 } // RealType owens_t(RealType h, RealType a)1040 1041 } // namespace detail1042 1043 template <class T1, class T2, class Policy>1044 inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a, const Policy& pol)1045 {1046 typedef typename tools::promote_args<T1, T2>::type result_type;1047 typedef typename policies::evaluation<result_type, Policy>::type value_type;1048 1049 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%)");1050 }1051 1052 template <class T1, class T2>1053 inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a)1054 {1055 return owens_t(h, a, policies::policy<>());1056 }1057 1058 1059 } // namespace math1060} // namespace boost1061 1062#ifdef _MSC_VER1063#pragma warning(pop)1064#endif1065 1066#endif1067// EOF1068