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