]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
92f5a8d4 TL |
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__ | |
f67539c2 | 40 | // nor #pragma diagnostic ignored work :( |
92f5a8d4 TL |
41 | // |
42 | #pragma GCC system_header | |
43 | #endif | |
44 | ||
7c673cae FG |
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. | |
20effc67 TL |
52 | template<typename RealType, class Policy> |
53 | inline RealType owens_t_znorm1(const RealType x, const Policy& pol) | |
7c673cae FG |
54 | { |
55 | using namespace boost::math::constants; | |
20effc67 | 56 | return boost::math::erf(x*one_div_root_two<RealType>(), pol)*half<RealType>(); |
7c673cae FG |
57 | } // RealType owens_t_znorm1(const RealType x) |
58 | ||
59 | // owens_t_znorm2(x) = P(x<=Z<oo) with Z being normally distributed. | |
20effc67 TL |
60 | template<typename RealType, class Policy> |
61 | inline RealType owens_t_znorm2(const RealType x, const Policy& pol) | |
7c673cae FG |
62 | { |
63 | using namespace boost::math::constants; | |
20effc67 | 64 | return boost::math::erfc(x*one_div_root_two<RealType>(), pol)*half<RealType>(); |
7c673cae FG |
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> | |
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 | ||
f67539c2 | 120 | // interpret select array as 8x15 matrix |
7c673cae FG |
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> | |
f67539c2 | 126 | inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const boost::integral_constant<int, 53>&) |
7c673cae FG |
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]; | |
f67539c2 | 133 | } // unsigned short owens_t_get_order(const unsigned short icode, RealType, boost::integral_constant<int, 53> const&) |
7c673cae FG |
134 | |
135 | template<typename RealType> | |
f67539c2 | 136 | inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const boost::integral_constant<int, 64>&) |
7c673cae FG |
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]; | |
f67539c2 | 144 | } // unsigned short owens_t_get_order(const unsigned short icode, RealType, boost::integral_constant<int, 64> const&) |
7c673cae FG |
145 | |
146 | template<typename RealType, typename 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; | |
f67539c2 TL |
150 | typedef boost::integral_constant<int, |
151 | precision_type::value <= 0 ? 64 : | |
152 | precision_type::value <= 53 ? 53 : 64 | |
153 | > tag_type; | |
7c673cae FG |
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> | |
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> | |
20effc67 | 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&) |
7c673cae FG |
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>(); | |
20effc67 | 209 | RealType z = owens_t_znorm1(ah, pol)/h; |
7c673cae FG |
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 | |
20effc67 TL |
228 | template<typename RealType, class Policy> |
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) | |
7c673cae FG |
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>(); | |
20effc67 | 258 | RealType zi = owens_t_znorm1(ah, pol)/h; |
7c673cae FG |
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 | |
20effc67 TL |
280 | template<class RealType, class Policy> |
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) | |
7c673cae FG |
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>(); | |
20effc67 | 330 | RealType zi = owens_t_znorm1(ah, pol)/h; |
7c673cae FG |
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> | |
20effc67 | 352 | inline RealType owens_t_T3(const RealType h, const RealType a, const RealType ah, const Policy& pol) |
7c673cae FG |
353 | { |
354 | typedef typename policies::precision<RealType, Policy>::type precision_type; | |
f67539c2 TL |
355 | typedef boost::integral_constant<int, |
356 | precision_type::value <= 0 ? 64 : | |
357 | precision_type::value <= 53 ? 53 : 64 | |
358 | > tag_type; | |
7c673cae | 359 | |
20effc67 | 360 | return owens_t_T3_imp(h, a, ah, tag_type(), pol); |
7c673cae FG |
361 | } |
362 | ||
363 | // compute the value of Owen's T function with method T4 from the reference paper | |
364 | template<typename RealType> | |
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> | |
f67539c2 | 394 | inline RealType owens_t_T5_imp(const RealType h, const RealType a, const boost::integral_constant<int, 53>&) |
7c673cae FG |
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> | |
f67539c2 | 440 | inline RealType owens_t_T5_imp(const RealType h, const RealType a, const boost::integral_constant<int, 64>&) |
7c673cae FG |
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> | |
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; | |
f67539c2 TL |
514 | typedef boost::integral_constant<int, |
515 | precision_type::value <= 0 ? 64 : | |
516 | precision_type::value <= 53 ? 53 : 64 | |
517 | > tag_type; | |
7c673cae FG |
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 | |
20effc67 TL |
524 | template<typename RealType, class Policy> |
525 | inline RealType owens_t_T6(const RealType h, const RealType a, const Policy& pol) | |
7c673cae FG |
526 | { |
527 | BOOST_MATH_STD_USING | |
528 | using namespace boost::math::constants; | |
529 | ||
20effc67 | 530 | const RealType normh = owens_t_znorm2(h, pol); |
7c673cae FG |
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> | |
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 | |
f67539c2 | 572 | // converging later. Various "guesstimates" have been tried to take account |
7c673cae FG |
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); | |
f67539c2 | 618 | if(sum < 0) // sum must always be positive, if it's negative something really bad has happened: |
7c673cae FG |
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> | |
20effc67 | 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&) |
7c673cae FG |
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>(); | |
20effc67 | 637 | RealType z = owens_t_znorm1(ah, pol)/h; |
7c673cae FG |
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> | |
20effc67 | 663 | inline std::pair<RealType, RealType> owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy& pol) |
7c673cae FG |
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>(); | |
20effc67 | 681 | RealType z = boost::math::detail::owens_t_znorm1(ah, pol)/h; |
7c673cae FG |
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 | |
f67539c2 | 690 | // converging later. Various "guesstimates" have been tried to take account |
7c673cae FG |
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> | |
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> | |
f67539c2 | 779 | inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, boost::integral_constant<int, 64> const&) |
7c673cae FG |
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 | { | |
20effc67 | 797 | return owens_t_znorm2(RealType(-h), pol) * owens_t_znorm2(h, pol) / 2; |
7c673cae FG |
798 | } |
799 | if(a >= tools::max_value<RealType>()) | |
800 | { | |
20effc67 | 801 | return owens_t_znorm2(RealType(fabs(h)), pol); |
7c673cae FG |
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; | |
f67539c2 | 816 | typedef boost::integral_constant<bool, (precision_type::value == 0) || (precision_type::value > 64)> tag_type; |
7c673cae FG |
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 | |
20effc67 | 829 | val = owens_t_T6(h,a, pol); |
7c673cae FG |
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> | |
f67539c2 | 838 | inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const boost::integral_constant<int, 65>&) |
7c673cae FG |
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 | { | |
20effc67 | 856 | return owens_t_znorm2(RealType(-h), pol) * owens_t_znorm2(h, pol) / 2; |
7c673cae FG |
857 | } |
858 | if(a >= tools::max_value<RealType>()) | |
859 | { | |
20effc67 | 860 | return owens_t_znorm2(RealType(fabs(h)), pol); |
7c673cae FG |
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 | // | |
f67539c2 | 959 | return owens_t_dispatch(h, a, ah, pol, boost::integral_constant<int, 64>()); |
7c673cae FG |
960 | } // RealType owens_t_dispatch(RealType h, RealType a, RealType ah) |
961 | template<typename RealType, typename Policy> | |
f67539c2 | 962 | inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const boost::integral_constant<int, 0>&) |
7c673cae FG |
963 | { |
964 | // We don't know what the precision is until runtime: | |
965 | if(tools::digits<RealType>() <= 64) | |
f67539c2 TL |
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>()); | |
7c673cae FG |
968 | } |
969 | template<typename RealType, typename Policy> | |
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; | |
f67539c2 TL |
974 | typedef boost::integral_constant<int, |
975 | precision_type::value <= 0 ? 0 : | |
976 | precision_type::value <= 64 ? 64 : 65 | |
977 | > tag_type; | |
978 | ||
7c673cae FG |
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> | |
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 | { | |
20effc67 TL |
1006 | const RealType normh = owens_t_znorm1(h, pol); |
1007 | const RealType normah = owens_t_znorm1(fabs_ah, pol); | |
7c673cae FG |
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 | { | |
20effc67 TL |
1013 | const RealType normh = detail::owens_t_znorm2(h, pol); |
1014 | const RealType normah = detail::owens_t_znorm2(fabs_ah, pol); | |
7c673cae FG |
1015 | val = constants::half<RealType>()*(normh+normah) - normh*normah - |
1016 | owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol); | |
1017 | } // else [if( h <= 0.67 )] | |
1018 | } // else [if(fabs_a <= 1)] | |
1019 | ||
1020 | // exploit that T(h,-a) == -T(h,a) | |
1021 | if(a < 0) | |
1022 | { | |
1023 | return -val; | |
1024 | } // if(a < 0) | |
1025 | ||
1026 | return val; | |
1027 | } // RealType owens_t(RealType h, RealType a) | |
1028 | ||
1029 | template <class T, class Policy, class tag> | |
1030 | struct owens_t_initializer | |
1031 | { | |
1032 | struct init | |
1033 | { | |
1034 | init() | |
1035 | { | |
1036 | do_init(tag()); | |
1037 | } | |
1038 | template <int N> | |
f67539c2 TL |
1039 | static void do_init(const boost::integral_constant<int, N>&){} |
1040 | static void do_init(const boost::integral_constant<int, 64>&) | |
7c673cae FG |
1041 | { |
1042 | boost::math::owens_t(static_cast<T>(7), static_cast<T>(0.96875), Policy()); | |
1043 | boost::math::owens_t(static_cast<T>(2), static_cast<T>(0.5), Policy()); | |
1044 | } | |
1045 | void force_instantiate()const{} | |
1046 | }; | |
1047 | static const init initializer; | |
1048 | static void force_instantiate() | |
1049 | { | |
1050 | initializer.force_instantiate(); | |
1051 | } | |
1052 | }; | |
1053 | ||
1054 | template <class T, class Policy, class tag> | |
1055 | const typename owens_t_initializer<T, Policy, tag>::init owens_t_initializer<T, Policy, tag>::initializer; | |
1056 | ||
1057 | } // namespace detail | |
1058 | ||
1059 | template <class T1, class T2, class Policy> | |
1060 | inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a, const Policy& pol) | |
1061 | { | |
1062 | typedef typename tools::promote_args<T1, T2>::type result_type; | |
1063 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1064 | typedef typename policies::precision<value_type, Policy>::type precision_type; | |
f67539c2 TL |
1065 | typedef boost::integral_constant<int, |
1066 | precision_type::value <= 0 ? 0 : | |
1067 | precision_type::value <= 64 ? 64 : 65 | |
1068 | > tag_type; | |
7c673cae FG |
1069 | |
1070 | detail::owens_t_initializer<result_type, Policy, tag_type>::force_instantiate(); | |
1071 | ||
1072 | return policies::checked_narrowing_cast<result_type, Policy>(detail::owens_t(static_cast<value_type>(h), static_cast<value_type>(a), pol), "boost::math::owens_t<%1%>(%1%,%1%)"); | |
1073 | } | |
1074 | ||
1075 | template <class T1, class T2> | |
1076 | inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a) | |
1077 | { | |
1078 | return owens_t(h, a, policies::policy<>()); | |
1079 | } | |
1080 | ||
1081 | ||
1082 | } // namespace math | |
1083 | } // namespace boost | |
1084 | ||
1085 | #ifdef BOOST_MSVC | |
1086 | #pragma warning(pop) | |
1087 | #endif | |
1088 | ||
1089 | #endif | |
1090 | // EOF |