]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // (C) Copyright John Maddock 2006. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #ifndef BOOST_MATH_SPECIAL_ERF_HPP | |
7 | #define BOOST_MATH_SPECIAL_ERF_HPP | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #endif | |
12 | ||
13 | #include <boost/math/special_functions/math_fwd.hpp> | |
14 | #include <boost/math/tools/config.hpp> | |
15 | #include <boost/math/special_functions/gamma.hpp> | |
16 | #include <boost/math/tools/roots.hpp> | |
17 | #include <boost/math/policies/error_handling.hpp> | |
18 | #include <boost/math/tools/big_constant.hpp> | |
19 | ||
92f5a8d4 TL |
20 | #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128) |
21 | // | |
22 | // This is the only way we can avoid | |
23 | // warning: non-standard suffix on floating constant [-Wpedantic] | |
24 | // when building with -Wall -pedantic. Neither __extension__ | |
f67539c2 | 25 | // nor #pragma diagnostic ignored work :( |
92f5a8d4 TL |
26 | // |
27 | #pragma GCC system_header | |
28 | #endif | |
29 | ||
7c673cae FG |
30 | namespace boost{ namespace math{ |
31 | ||
32 | namespace detail | |
33 | { | |
34 | ||
35 | // | |
36 | // Asymptotic series for large z: | |
37 | // | |
38 | template <class T> | |
39 | struct erf_asympt_series_t | |
40 | { | |
41 | erf_asympt_series_t(T z) : xx(2 * -z * z), tk(1) | |
42 | { | |
43 | BOOST_MATH_STD_USING | |
44 | result = -exp(-z * z) / sqrt(boost::math::constants::pi<T>()); | |
45 | result /= z; | |
46 | } | |
47 | ||
48 | typedef T result_type; | |
49 | ||
50 | T operator()() | |
51 | { | |
52 | BOOST_MATH_STD_USING | |
53 | T r = result; | |
54 | result *= tk / xx; | |
55 | tk += 2; | |
56 | if( fabs(r) < fabs(result)) | |
57 | result = 0; | |
58 | return r; | |
59 | } | |
60 | private: | |
61 | T result; | |
62 | T xx; | |
63 | int tk; | |
64 | }; | |
65 | // | |
66 | // How large z has to be in order to ensure that the series converges: | |
67 | // | |
68 | template <class T> | |
69 | inline float erf_asymptotic_limit_N(const T&) | |
70 | { | |
71 | return (std::numeric_limits<float>::max)(); | |
72 | } | |
1e59de90 | 73 | inline float erf_asymptotic_limit_N(const std::integral_constant<int, 24>&) |
7c673cae FG |
74 | { |
75 | return 2.8F; | |
76 | } | |
1e59de90 | 77 | inline float erf_asymptotic_limit_N(const std::integral_constant<int, 53>&) |
7c673cae FG |
78 | { |
79 | return 4.3F; | |
80 | } | |
1e59de90 | 81 | inline float erf_asymptotic_limit_N(const std::integral_constant<int, 64>&) |
7c673cae FG |
82 | { |
83 | return 4.8F; | |
84 | } | |
1e59de90 | 85 | inline float erf_asymptotic_limit_N(const std::integral_constant<int, 106>&) |
7c673cae FG |
86 | { |
87 | return 6.5F; | |
88 | } | |
1e59de90 | 89 | inline float erf_asymptotic_limit_N(const std::integral_constant<int, 113>&) |
7c673cae FG |
90 | { |
91 | return 6.8F; | |
92 | } | |
93 | ||
94 | template <class T, class Policy> | |
95 | inline T erf_asymptotic_limit() | |
96 | { | |
97 | typedef typename policies::precision<T, Policy>::type precision_type; | |
1e59de90 | 98 | typedef std::integral_constant<int, |
f67539c2 TL |
99 | precision_type::value <= 0 ? 0 : |
100 | precision_type::value <= 24 ? 24 : | |
101 | precision_type::value <= 53 ? 53 : | |
102 | precision_type::value <= 64 ? 64 : | |
103 | precision_type::value <= 113 ? 113 : 0 | |
104 | > tag_type; | |
7c673cae FG |
105 | return erf_asymptotic_limit_N(tag_type()); |
106 | } | |
107 | ||
20effc67 TL |
108 | template <class T> |
109 | struct erf_series_near_zero | |
110 | { | |
111 | typedef T result_type; | |
112 | T term; | |
113 | T zz; | |
114 | int k; | |
115 | erf_series_near_zero(const T& z) : term(z), zz(-z * z), k(0) {} | |
116 | ||
117 | T operator()() | |
118 | { | |
119 | T result = term / (2 * k + 1); | |
120 | term *= zz / ++k; | |
121 | return result; | |
122 | } | |
123 | }; | |
124 | ||
125 | template <class T, class Policy> | |
126 | T erf_series_near_zero_sum(const T& x, const Policy& pol) | |
127 | { | |
128 | // | |
129 | // We need Kahan summation here, otherwise the errors grow fairly quickly. | |
130 | // This method is *much* faster than the alternatives even so. | |
131 | // | |
132 | erf_series_near_zero<T> sum(x); | |
1e59de90 | 133 | std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); |
20effc67 TL |
134 | T result = constants::two_div_root_pi<T>() * tools::kahan_sum_series(sum, tools::digits<T>(), max_iter); |
135 | policies::check_series_iterations<T>("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol); | |
136 | return result; | |
137 | } | |
138 | ||
7c673cae FG |
139 | template <class T, class Policy, class Tag> |
140 | T erf_imp(T z, bool invert, const Policy& pol, const Tag& t) | |
141 | { | |
142 | BOOST_MATH_STD_USING | |
143 | ||
144 | BOOST_MATH_INSTRUMENT_CODE("Generic erf_imp called"); | |
145 | ||
146 | if(z < 0) | |
147 | { | |
148 | if(!invert) | |
149 | return -erf_imp(T(-z), invert, pol, t); | |
150 | else | |
151 | return 1 + erf_imp(T(-z), false, pol, t); | |
152 | } | |
153 | ||
154 | T result; | |
155 | ||
156 | if(!invert && (z > detail::erf_asymptotic_limit<T, Policy>())) | |
157 | { | |
158 | detail::erf_asympt_series_t<T> s(z); | |
1e59de90 | 159 | std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); |
7c673cae FG |
160 | result = boost::math::tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, 1); |
161 | policies::check_series_iterations<T>("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol); | |
162 | } | |
163 | else | |
164 | { | |
165 | T x = z * z; | |
20effc67 | 166 | if(z < 1.3f) |
7c673cae FG |
167 | { |
168 | // Compute P: | |
20effc67 TL |
169 | // This is actually good for z p to 2 or so, but the cutoff given seems |
170 | // to be the best compromise. Performance wise, this is way quicker than anything else... | |
171 | result = erf_series_near_zero_sum(z, pol); | |
7c673cae | 172 | } |
b32b8144 FG |
173 | else if(x > 1 / tools::epsilon<T>()) |
174 | { | |
175 | // http://functions.wolfram.com/06.27.06.0006.02 | |
176 | invert = !invert; | |
177 | result = exp(-x) / (constants::root_pi<T>() * z); | |
178 | } | |
7c673cae FG |
179 | else |
180 | { | |
181 | // Compute Q: | |
182 | invert = !invert; | |
183 | result = z * exp(-x); | |
b32b8144 | 184 | result /= boost::math::constants::root_pi<T>(); |
7c673cae FG |
185 | result *= upper_gamma_fraction(T(0.5f), x, policies::get_epsilon<T, Policy>()); |
186 | } | |
187 | } | |
188 | if(invert) | |
189 | result = 1 - result; | |
190 | return result; | |
191 | } | |
192 | ||
193 | template <class T, class Policy> | |
1e59de90 | 194 | T erf_imp(T z, bool invert, const Policy& pol, const std::integral_constant<int, 53>& t) |
7c673cae FG |
195 | { |
196 | BOOST_MATH_STD_USING | |
197 | ||
198 | BOOST_MATH_INSTRUMENT_CODE("53-bit precision erf_imp called"); | |
199 | ||
92f5a8d4 TL |
200 | if ((boost::math::isnan)(z)) |
201 | return policies::raise_denorm_error("boost::math::erf<%1%>(%1%)", "Expected a finite argument but got %1%", z, pol); | |
202 | ||
7c673cae FG |
203 | if(z < 0) |
204 | { | |
205 | if(!invert) | |
206 | return -erf_imp(T(-z), invert, pol, t); | |
207 | else if(z < -0.5) | |
208 | return 2 - erf_imp(T(-z), invert, pol, t); | |
209 | else | |
210 | return 1 + erf_imp(T(-z), false, pol, t); | |
211 | } | |
212 | ||
213 | T result; | |
214 | ||
215 | // | |
216 | // Big bunch of selection statements now to pick | |
217 | // which implementation to use, | |
218 | // try to put most likely options first: | |
219 | // | |
220 | if(z < 0.5) | |
221 | { | |
222 | // | |
223 | // We're going to calculate erf: | |
224 | // | |
225 | if(z < 1e-10) | |
226 | { | |
227 | if(z == 0) | |
228 | { | |
229 | result = T(0); | |
230 | } | |
231 | else | |
232 | { | |
233 | static const T c = BOOST_MATH_BIG_CONSTANT(T, 53, 0.003379167095512573896158903121545171688); | |
234 | result = static_cast<T>(z * 1.125f + z * c); | |
235 | } | |
236 | } | |
237 | else | |
238 | { | |
239 | // Maximum Deviation Found: 1.561e-17 | |
240 | // Expected Error Term: 1.561e-17 | |
241 | // Maximum Relative Change in Control Points: 1.155e-04 | |
242 | // Max Error found at double precision = 2.961182e-17 | |
243 | ||
244 | static const T Y = 1.044948577880859375f; | |
245 | static const T P[] = { | |
246 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.0834305892146531832907), | |
247 | BOOST_MATH_BIG_CONSTANT(T, 53, -0.338165134459360935041), | |
248 | BOOST_MATH_BIG_CONSTANT(T, 53, -0.0509990735146777432841), | |
249 | BOOST_MATH_BIG_CONSTANT(T, 53, -0.00772758345802133288487), | |
250 | BOOST_MATH_BIG_CONSTANT(T, 53, -0.000322780120964605683831), | |
251 | }; | |
252 | static const T Q[] = { | |
253 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.0), | |
254 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.455004033050794024546), | |
255 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.0875222600142252549554), | |
256 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.00858571925074406212772), | |
257 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.000370900071787748000569), | |
258 | }; | |
259 | T zz = z * z; | |
260 | result = z * (Y + tools::evaluate_polynomial(P, zz) / tools::evaluate_polynomial(Q, zz)); | |
261 | } | |
262 | } | |
1e59de90 | 263 | else if(invert ? (z < 28) : (z < 5.93f)) |
7c673cae FG |
264 | { |
265 | // | |
266 | // We'll be calculating erfc: | |
267 | // | |
268 | invert = !invert; | |
269 | if(z < 1.5f) | |
270 | { | |
271 | // Maximum Deviation Found: 3.702e-17 | |
272 | // Expected Error Term: 3.702e-17 | |
273 | // Maximum Relative Change in Control Points: 2.845e-04 | |
274 | // Max Error found at double precision = 4.841816e-17 | |
275 | static const T Y = 0.405935764312744140625f; | |
276 | static const T P[] = { | |
277 | BOOST_MATH_BIG_CONSTANT(T, 53, -0.098090592216281240205), | |
278 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.178114665841120341155), | |
279 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.191003695796775433986), | |
280 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.0888900368967884466578), | |
281 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.0195049001251218801359), | |
282 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.00180424538297014223957), | |
283 | }; | |
284 | static const T Q[] = { | |
285 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.0), | |
286 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.84759070983002217845), | |
287 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.42628004845511324508), | |
288 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.578052804889902404909), | |
289 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.12385097467900864233), | |
290 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.0113385233577001411017), | |
291 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.337511472483094676155e-5), | |
292 | }; | |
293 | BOOST_MATH_INSTRUMENT_VARIABLE(Y); | |
294 | BOOST_MATH_INSTRUMENT_VARIABLE(P[0]); | |
295 | BOOST_MATH_INSTRUMENT_VARIABLE(Q[0]); | |
296 | BOOST_MATH_INSTRUMENT_VARIABLE(z); | |
297 | result = Y + tools::evaluate_polynomial(P, T(z - 0.5)) / tools::evaluate_polynomial(Q, T(z - 0.5)); | |
298 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
299 | result *= exp(-z * z) / z; | |
300 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
301 | } | |
302 | else if(z < 2.5f) | |
303 | { | |
304 | // Max Error found at double precision = 6.599585e-18 | |
305 | // Maximum Deviation Found: 3.909e-18 | |
306 | // Expected Error Term: 3.909e-18 | |
307 | // Maximum Relative Change in Control Points: 9.886e-05 | |
308 | static const T Y = 0.50672817230224609375f; | |
309 | static const T P[] = { | |
310 | BOOST_MATH_BIG_CONSTANT(T, 53, -0.0243500476207698441272), | |
311 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.0386540375035707201728), | |
312 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.04394818964209516296), | |
313 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.0175679436311802092299), | |
314 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.00323962406290842133584), | |
315 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.000235839115596880717416), | |
316 | }; | |
317 | static const T Q[] = { | |
318 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.0), | |
319 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.53991494948552447182), | |
320 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.982403709157920235114), | |
321 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.325732924782444448493), | |
322 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.0563921837420478160373), | |
323 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.00410369723978904575884), | |
324 | }; | |
325 | result = Y + tools::evaluate_polynomial(P, T(z - 1.5)) / tools::evaluate_polynomial(Q, T(z - 1.5)); | |
92f5a8d4 TL |
326 | T hi, lo; |
327 | int expon; | |
328 | hi = floor(ldexp(frexp(z, &expon), 26)); | |
329 | hi = ldexp(hi, expon - 26); | |
330 | lo = z - hi; | |
331 | T sq = z * z; | |
332 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
333 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
334 | } |
335 | else if(z < 4.5f) | |
336 | { | |
337 | // Maximum Deviation Found: 1.512e-17 | |
338 | // Expected Error Term: 1.512e-17 | |
339 | // Maximum Relative Change in Control Points: 2.222e-04 | |
340 | // Max Error found at double precision = 2.062515e-17 | |
341 | static const T Y = 0.5405750274658203125f; | |
342 | static const T P[] = { | |
343 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.00295276716530971662634), | |
344 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.0137384425896355332126), | |
345 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.00840807615555585383007), | |
346 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.00212825620914618649141), | |
347 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.000250269961544794627958), | |
348 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.113212406648847561139e-4), | |
349 | }; | |
350 | static const T Q[] = { | |
351 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.0), | |
352 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.04217814166938418171), | |
353 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.442597659481563127003), | |
354 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.0958492726301061423444), | |
355 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.0105982906484876531489), | |
356 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.000479411269521714493907), | |
357 | }; | |
358 | result = Y + tools::evaluate_polynomial(P, T(z - 3.5)) / tools::evaluate_polynomial(Q, T(z - 3.5)); | |
92f5a8d4 TL |
359 | T hi, lo; |
360 | int expon; | |
361 | hi = floor(ldexp(frexp(z, &expon), 26)); | |
362 | hi = ldexp(hi, expon - 26); | |
363 | lo = z - hi; | |
364 | T sq = z * z; | |
365 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
366 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
367 | } |
368 | else | |
369 | { | |
370 | // Max Error found at double precision = 2.997958e-17 | |
371 | // Maximum Deviation Found: 2.860e-17 | |
372 | // Expected Error Term: 2.859e-17 | |
373 | // Maximum Relative Change in Control Points: 1.357e-05 | |
374 | static const T Y = 0.5579090118408203125f; | |
375 | static const T P[] = { | |
376 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.00628057170626964891937), | |
377 | BOOST_MATH_BIG_CONSTANT(T, 53, 0.0175389834052493308818), | |
378 | BOOST_MATH_BIG_CONSTANT(T, 53, -0.212652252872804219852), | |
379 | BOOST_MATH_BIG_CONSTANT(T, 53, -0.687717681153649930619), | |
380 | BOOST_MATH_BIG_CONSTANT(T, 53, -2.5518551727311523996), | |
381 | BOOST_MATH_BIG_CONSTANT(T, 53, -3.22729451764143718517), | |
382 | BOOST_MATH_BIG_CONSTANT(T, 53, -2.8175401114513378771), | |
383 | }; | |
384 | static const T Q[] = { | |
385 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.0), | |
386 | BOOST_MATH_BIG_CONSTANT(T, 53, 2.79257750980575282228), | |
387 | BOOST_MATH_BIG_CONSTANT(T, 53, 11.0567237927800161565), | |
388 | BOOST_MATH_BIG_CONSTANT(T, 53, 15.930646027911794143), | |
389 | BOOST_MATH_BIG_CONSTANT(T, 53, 22.9367376522880577224), | |
390 | BOOST_MATH_BIG_CONSTANT(T, 53, 13.5064170191802889145), | |
391 | BOOST_MATH_BIG_CONSTANT(T, 53, 5.48409182238641741584), | |
392 | }; | |
393 | result = Y + tools::evaluate_polynomial(P, T(1 / z)) / tools::evaluate_polynomial(Q, T(1 / z)); | |
92f5a8d4 TL |
394 | T hi, lo; |
395 | int expon; | |
396 | hi = floor(ldexp(frexp(z, &expon), 26)); | |
397 | hi = ldexp(hi, expon - 26); | |
398 | lo = z - hi; | |
399 | T sq = z * z; | |
400 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
401 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
402 | } |
403 | } | |
404 | else | |
405 | { | |
406 | // | |
407 | // Any value of z larger than 28 will underflow to zero: | |
408 | // | |
409 | result = 0; | |
410 | invert = !invert; | |
411 | } | |
412 | ||
413 | if(invert) | |
414 | { | |
415 | result = 1 - result; | |
416 | } | |
417 | ||
418 | return result; | |
1e59de90 | 419 | } // template <class T, class Lanczos>T erf_imp(T z, bool invert, const Lanczos& l, const std::integral_constant<int, 53>& t) |
7c673cae FG |
420 | |
421 | ||
422 | template <class T, class Policy> | |
1e59de90 | 423 | T erf_imp(T z, bool invert, const Policy& pol, const std::integral_constant<int, 64>& t) |
7c673cae FG |
424 | { |
425 | BOOST_MATH_STD_USING | |
426 | ||
427 | BOOST_MATH_INSTRUMENT_CODE("64-bit precision erf_imp called"); | |
428 | ||
429 | if(z < 0) | |
430 | { | |
431 | if(!invert) | |
432 | return -erf_imp(T(-z), invert, pol, t); | |
433 | else if(z < -0.5) | |
434 | return 2 - erf_imp(T(-z), invert, pol, t); | |
435 | else | |
436 | return 1 + erf_imp(T(-z), false, pol, t); | |
437 | } | |
438 | ||
439 | T result; | |
440 | ||
441 | // | |
442 | // Big bunch of selection statements now to pick which | |
443 | // implementation to use, try to put most likely options | |
444 | // first: | |
445 | // | |
446 | if(z < 0.5) | |
447 | { | |
448 | // | |
449 | // We're going to calculate erf: | |
450 | // | |
451 | if(z == 0) | |
452 | { | |
453 | result = 0; | |
454 | } | |
455 | else if(z < 1e-10) | |
456 | { | |
457 | static const T c = BOOST_MATH_BIG_CONSTANT(T, 64, 0.003379167095512573896158903121545171688); | |
458 | result = z * 1.125 + z * c; | |
459 | } | |
460 | else | |
461 | { | |
462 | // Max Error found at long double precision = 1.623299e-20 | |
463 | // Maximum Deviation Found: 4.326e-22 | |
464 | // Expected Error Term: -4.326e-22 | |
465 | // Maximum Relative Change in Control Points: 1.474e-04 | |
466 | static const T Y = 1.044948577880859375f; | |
467 | static const T P[] = { | |
468 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0834305892146531988966), | |
469 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.338097283075565413695), | |
470 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.0509602734406067204596), | |
471 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.00904906346158537794396), | |
472 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.000489468651464798669181), | |
473 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.200305626366151877759e-4), | |
474 | }; | |
475 | static const T Q[] = { | |
476 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), | |
477 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.455817300515875172439), | |
478 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0916537354356241792007), | |
479 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0102722652675910031202), | |
480 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.000650511752687851548735), | |
481 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.189532519105655496778e-4), | |
482 | }; | |
483 | result = z * (Y + tools::evaluate_polynomial(P, T(z * z)) / tools::evaluate_polynomial(Q, T(z * z))); | |
484 | } | |
485 | } | |
1e59de90 | 486 | else if(invert ? (z < 110) : (z < 6.6f)) |
7c673cae FG |
487 | { |
488 | // | |
489 | // We'll be calculating erfc: | |
490 | // | |
491 | invert = !invert; | |
492 | if(z < 1.5) | |
493 | { | |
494 | // Max Error found at long double precision = 3.239590e-20 | |
495 | // Maximum Deviation Found: 2.241e-20 | |
496 | // Expected Error Term: -2.241e-20 | |
497 | // Maximum Relative Change in Control Points: 5.110e-03 | |
498 | static const T Y = 0.405935764312744140625f; | |
499 | static const T P[] = { | |
500 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.0980905922162812031672), | |
501 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.159989089922969141329), | |
502 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.222359821619935712378), | |
503 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.127303921703577362312), | |
504 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0384057530342762400273), | |
505 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.00628431160851156719325), | |
506 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.000441266654514391746428), | |
507 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.266689068336295642561e-7), | |
508 | }; | |
509 | static const T Q[] = { | |
510 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), | |
511 | BOOST_MATH_BIG_CONSTANT(T, 64, 2.03237474985469469291), | |
512 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.78355454954969405222), | |
513 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.867940326293760578231), | |
514 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.248025606990021698392), | |
515 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0396649631833002269861), | |
516 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.00279220237309449026796), | |
517 | }; | |
518 | result = Y + tools::evaluate_polynomial(P, T(z - 0.5f)) / tools::evaluate_polynomial(Q, T(z - 0.5f)); | |
92f5a8d4 TL |
519 | T hi, lo; |
520 | int expon; | |
521 | hi = floor(ldexp(frexp(z, &expon), 32)); | |
522 | hi = ldexp(hi, expon - 32); | |
523 | lo = z - hi; | |
524 | T sq = z * z; | |
525 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
526 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
527 | } |
528 | else if(z < 2.5) | |
529 | { | |
530 | // Max Error found at long double precision = 3.686211e-21 | |
531 | // Maximum Deviation Found: 1.495e-21 | |
532 | // Expected Error Term: -1.494e-21 | |
533 | // Maximum Relative Change in Control Points: 1.793e-04 | |
534 | static const T Y = 0.50672817230224609375f; | |
535 | static const T P[] = { | |
536 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.024350047620769840217), | |
537 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0343522687935671451309), | |
538 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0505420824305544949541), | |
539 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0257479325917757388209), | |
540 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.00669349844190354356118), | |
541 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.00090807914416099524444), | |
542 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.515917266698050027934e-4), | |
543 | }; | |
544 | static const T Q[] = { | |
545 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), | |
546 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.71657861671930336344), | |
547 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.26409634824280366218), | |
548 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.512371437838969015941), | |
549 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.120902623051120950935), | |
550 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0158027197831887485261), | |
551 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.000897871370778031611439), | |
552 | }; | |
553 | result = Y + tools::evaluate_polynomial(P, T(z - 1.5f)) / tools::evaluate_polynomial(Q, T(z - 1.5f)); | |
92f5a8d4 TL |
554 | T hi, lo; |
555 | int expon; | |
556 | hi = floor(ldexp(frexp(z, &expon), 32)); | |
557 | hi = ldexp(hi, expon - 32); | |
558 | lo = z - hi; | |
559 | T sq = z * z; | |
560 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
561 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
562 | } |
563 | else if(z < 4.5) | |
564 | { | |
565 | // Maximum Deviation Found: 1.107e-20 | |
566 | // Expected Error Term: -1.106e-20 | |
567 | // Maximum Relative Change in Control Points: 1.709e-04 | |
568 | // Max Error found at long double precision = 1.446908e-20 | |
569 | static const T Y = 0.5405750274658203125f; | |
570 | static const T P[] = { | |
571 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0029527671653097284033), | |
572 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0141853245895495604051), | |
573 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0104959584626432293901), | |
574 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.00343963795976100077626), | |
575 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.00059065441194877637899), | |
576 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.523435380636174008685e-4), | |
577 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.189896043050331257262e-5), | |
578 | }; | |
579 | static const T Q[] = { | |
580 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), | |
581 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.19352160185285642574), | |
582 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.603256964363454392857), | |
583 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.165411142458540585835), | |
584 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0259729870946203166468), | |
585 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.00221657568292893699158), | |
586 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.804149464190309799804e-4), | |
587 | }; | |
588 | result = Y + tools::evaluate_polynomial(P, T(z - 3.5f)) / tools::evaluate_polynomial(Q, T(z - 3.5f)); | |
92f5a8d4 TL |
589 | T hi, lo; |
590 | int expon; | |
591 | hi = floor(ldexp(frexp(z, &expon), 32)); | |
592 | hi = ldexp(hi, expon - 32); | |
593 | lo = z - hi; | |
594 | T sq = z * z; | |
595 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
596 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
597 | } |
598 | else | |
599 | { | |
600 | // Max Error found at long double precision = 7.961166e-21 | |
601 | // Maximum Deviation Found: 6.677e-21 | |
602 | // Expected Error Term: 6.676e-21 | |
603 | // Maximum Relative Change in Control Points: 2.319e-05 | |
604 | static const T Y = 0.55825519561767578125f; | |
605 | static const T P[] = { | |
606 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.00593438793008050214106), | |
607 | BOOST_MATH_BIG_CONSTANT(T, 64, 0.0280666231009089713937), | |
608 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.141597835204583050043), | |
609 | BOOST_MATH_BIG_CONSTANT(T, 64, -0.978088201154300548842), | |
610 | BOOST_MATH_BIG_CONSTANT(T, 64, -5.47351527796012049443), | |
611 | BOOST_MATH_BIG_CONSTANT(T, 64, -13.8677304660245326627), | |
612 | BOOST_MATH_BIG_CONSTANT(T, 64, -27.1274948720539821722), | |
613 | BOOST_MATH_BIG_CONSTANT(T, 64, -29.2545152747009461519), | |
614 | BOOST_MATH_BIG_CONSTANT(T, 64, -16.8865774499799676937), | |
615 | }; | |
616 | static const T Q[] = { | |
617 | BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), | |
618 | BOOST_MATH_BIG_CONSTANT(T, 64, 4.72948911186645394541), | |
619 | BOOST_MATH_BIG_CONSTANT(T, 64, 23.6750543147695749212), | |
620 | BOOST_MATH_BIG_CONSTANT(T, 64, 60.0021517335693186785), | |
621 | BOOST_MATH_BIG_CONSTANT(T, 64, 131.766251645149522868), | |
622 | BOOST_MATH_BIG_CONSTANT(T, 64, 178.167924971283482513), | |
623 | BOOST_MATH_BIG_CONSTANT(T, 64, 182.499390505915222699), | |
624 | BOOST_MATH_BIG_CONSTANT(T, 64, 104.365251479578577989), | |
625 | BOOST_MATH_BIG_CONSTANT(T, 64, 30.8365511891224291717), | |
626 | }; | |
627 | result = Y + tools::evaluate_polynomial(P, T(1 / z)) / tools::evaluate_polynomial(Q, T(1 / z)); | |
92f5a8d4 TL |
628 | T hi, lo; |
629 | int expon; | |
630 | hi = floor(ldexp(frexp(z, &expon), 32)); | |
631 | hi = ldexp(hi, expon - 32); | |
632 | lo = z - hi; | |
633 | T sq = z * z; | |
634 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
635 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
636 | } |
637 | } | |
638 | else | |
639 | { | |
640 | // | |
641 | // Any value of z larger than 110 will underflow to zero: | |
642 | // | |
643 | result = 0; | |
644 | invert = !invert; | |
645 | } | |
646 | ||
647 | if(invert) | |
648 | { | |
649 | result = 1 - result; | |
650 | } | |
651 | ||
652 | return result; | |
1e59de90 | 653 | } // template <class T, class Lanczos>T erf_imp(T z, bool invert, const Lanczos& l, const std::integral_constant<int, 64>& t) |
7c673cae FG |
654 | |
655 | ||
656 | template <class T, class Policy> | |
1e59de90 | 657 | T erf_imp(T z, bool invert, const Policy& pol, const std::integral_constant<int, 113>& t) |
7c673cae FG |
658 | { |
659 | BOOST_MATH_STD_USING | |
660 | ||
661 | BOOST_MATH_INSTRUMENT_CODE("113-bit precision erf_imp called"); | |
662 | ||
663 | if(z < 0) | |
664 | { | |
665 | if(!invert) | |
666 | return -erf_imp(T(-z), invert, pol, t); | |
667 | else if(z < -0.5) | |
668 | return 2 - erf_imp(T(-z), invert, pol, t); | |
669 | else | |
670 | return 1 + erf_imp(T(-z), false, pol, t); | |
671 | } | |
672 | ||
673 | T result; | |
674 | ||
675 | // | |
676 | // Big bunch of selection statements now to pick which | |
677 | // implementation to use, try to put most likely options | |
678 | // first: | |
679 | // | |
680 | if(z < 0.5) | |
681 | { | |
682 | // | |
683 | // We're going to calculate erf: | |
684 | // | |
685 | if(z == 0) | |
686 | { | |
687 | result = 0; | |
688 | } | |
689 | else if(z < 1e-20) | |
690 | { | |
691 | static const T c = BOOST_MATH_BIG_CONSTANT(T, 113, 0.003379167095512573896158903121545171688); | |
692 | result = z * 1.125 + z * c; | |
693 | } | |
694 | else | |
695 | { | |
696 | // Max Error found at long double precision = 2.342380e-35 | |
697 | // Maximum Deviation Found: 6.124e-36 | |
698 | // Expected Error Term: -6.124e-36 | |
699 | // Maximum Relative Change in Control Points: 3.492e-10 | |
700 | static const T Y = 1.0841522216796875f; | |
701 | static const T P[] = { | |
702 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0442269454158250738961589031215451778), | |
703 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.35549265736002144875335323556961233), | |
704 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.0582179564566667896225454670863270393), | |
705 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.0112694696904802304229950538453123925), | |
706 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.000805730648981801146251825329609079099), | |
707 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.566304966591936566229702842075966273e-4), | |
708 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.169655010425186987820201021510002265e-5), | |
709 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.344448249920445916714548295433198544e-7), | |
710 | }; | |
711 | static const T Q[] = { | |
712 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
713 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.466542092785657604666906909196052522), | |
714 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.100005087012526447295176964142107611), | |
715 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0128341535890117646540050072234142603), | |
716 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00107150448466867929159660677016658186), | |
717 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.586168368028999183607733369248338474e-4), | |
718 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.196230608502104324965623171516808796e-5), | |
719 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.313388521582925207734229967907890146e-7), | |
720 | }; | |
721 | result = z * (Y + tools::evaluate_polynomial(P, T(z * z)) / tools::evaluate_polynomial(Q, T(z * z))); | |
722 | } | |
723 | } | |
724 | else if(invert ? (z < 110) : (z < 8.65f)) | |
725 | { | |
726 | // | |
727 | // We'll be calculating erfc: | |
728 | // | |
729 | invert = !invert; | |
730 | if(z < 1) | |
731 | { | |
732 | // Max Error found at long double precision = 3.246278e-35 | |
733 | // Maximum Deviation Found: 1.388e-35 | |
734 | // Expected Error Term: 1.387e-35 | |
735 | // Maximum Relative Change in Control Points: 6.127e-05 | |
736 | static const T Y = 0.371877193450927734375f; | |
737 | static const T P[] = { | |
738 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.0640320213544647969396032886581290455), | |
739 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.200769874440155895637857443946706731), | |
740 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.378447199873537170666487408805779826), | |
741 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.30521399466465939450398642044975127), | |
742 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.146890026406815277906781824723458196), | |
743 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0464837937749539978247589252732769567), | |
744 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00987895759019540115099100165904822903), | |
745 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00137507575429025512038051025154301132), | |
746 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0001144764551085935580772512359680516), | |
747 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.436544865032836914773944382339900079e-5), | |
748 | }; | |
749 | static const T Q[] = { | |
750 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
751 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.47651182872457465043733800302427977), | |
752 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.78706486002517996428836400245547955), | |
753 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.87295924621659627926365005293130693), | |
754 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.829375825174365625428280908787261065), | |
755 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.251334771307848291593780143950311514), | |
756 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0522110268876176186719436765734722473), | |
757 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00718332151250963182233267040106902368), | |
758 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000595279058621482041084986219276392459), | |
759 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.226988669466501655990637599399326874e-4), | |
760 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.270666232259029102353426738909226413e-10), | |
761 | }; | |
762 | result = Y + tools::evaluate_polynomial(P, T(z - 0.5f)) / tools::evaluate_polynomial(Q, T(z - 0.5f)); | |
92f5a8d4 TL |
763 | T hi, lo; |
764 | int expon; | |
765 | hi = floor(ldexp(frexp(z, &expon), 56)); | |
766 | hi = ldexp(hi, expon - 56); | |
767 | lo = z - hi; | |
768 | T sq = z * z; | |
769 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
770 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
771 | } |
772 | else if(z < 1.5) | |
773 | { | |
774 | // Max Error found at long double precision = 2.215785e-35 | |
775 | // Maximum Deviation Found: 1.539e-35 | |
776 | // Expected Error Term: 1.538e-35 | |
777 | // Maximum Relative Change in Control Points: 6.104e-05 | |
778 | static const T Y = 0.45658016204833984375f; | |
779 | static const T P[] = { | |
780 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.0289965858925328393392496555094848345), | |
781 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0868181194868601184627743162571779226), | |
782 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.169373435121178901746317404936356745), | |
783 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.13350446515949251201104889028133486), | |
784 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0617447837290183627136837688446313313), | |
785 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0185618495228251406703152962489700468), | |
786 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00371949406491883508764162050169531013), | |
787 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000485121708792921297742105775823900772), | |
788 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.376494706741453489892108068231400061e-4), | |
789 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.133166058052466262415271732172490045e-5), | |
790 | }; | |
791 | static const T Q[] = { | |
792 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
793 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.32970330146503867261275580968135126), | |
794 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.46325715420422771961250513514928746), | |
795 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.55307882560757679068505047390857842), | |
796 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.644274289865972449441174485441409076), | |
797 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.182609091063258208068606847453955649), | |
798 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0354171651271241474946129665801606795), | |
799 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00454060370165285246451879969534083997), | |
800 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000349871943711566546821198612518656486), | |
801 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.123749319840299552925421880481085392e-4), | |
802 | }; | |
803 | result = Y + tools::evaluate_polynomial(P, T(z - 1.0f)) / tools::evaluate_polynomial(Q, T(z - 1.0f)); | |
92f5a8d4 TL |
804 | T hi, lo; |
805 | int expon; | |
806 | hi = floor(ldexp(frexp(z, &expon), 56)); | |
807 | hi = ldexp(hi, expon - 56); | |
808 | lo = z - hi; | |
809 | T sq = z * z; | |
810 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
811 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
812 | } |
813 | else if(z < 2.25) | |
814 | { | |
815 | // Maximum Deviation Found: 1.418e-35 | |
816 | // Expected Error Term: 1.418e-35 | |
817 | // Maximum Relative Change in Control Points: 1.316e-04 | |
818 | // Max Error found at long double precision = 1.998462e-35 | |
819 | static const T Y = 0.50250148773193359375f; | |
820 | static const T P[] = { | |
821 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.0201233630504573402185161184151016606), | |
822 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0331864357574860196516686996302305002), | |
823 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0716562720864787193337475444413405461), | |
824 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0545835322082103985114927569724880658), | |
825 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0236692635189696678976549720784989593), | |
826 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00656970902163248872837262539337601845), | |
827 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00120282643299089441390490459256235021), | |
828 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000142123229065182650020762792081622986), | |
829 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.991531438367015135346716277792989347e-5), | |
830 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.312857043762117596999398067153076051e-6), | |
831 | }; | |
832 | static const T Q[] = { | |
833 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
834 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.13506082409097783827103424943508554), | |
835 | BOOST_MATH_BIG_CONSTANT(T, 113, 2.06399257267556230937723190496806215), | |
836 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.18678481279932541314830499880691109), | |
837 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.447733186643051752513538142316799562), | |
838 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.11505680005657879437196953047542148), | |
839 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.020163993632192726170219663831914034), | |
840 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00232708971840141388847728782209730585), | |
841 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000160733201627963528519726484608224112), | |
842 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.507158721790721802724402992033269266e-5), | |
843 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.18647774409821470950544212696270639e-12), | |
844 | }; | |
845 | result = Y + tools::evaluate_polynomial(P, T(z - 1.5f)) / tools::evaluate_polynomial(Q, T(z - 1.5f)); | |
92f5a8d4 TL |
846 | T hi, lo; |
847 | int expon; | |
848 | hi = floor(ldexp(frexp(z, &expon), 56)); | |
849 | hi = ldexp(hi, expon - 56); | |
850 | lo = z - hi; | |
851 | T sq = z * z; | |
852 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
853 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
854 | } |
855 | else if (z < 3) | |
856 | { | |
857 | // Maximum Deviation Found: 3.575e-36 | |
858 | // Expected Error Term: 3.575e-36 | |
859 | // Maximum Relative Change in Control Points: 7.103e-05 | |
860 | // Max Error found at long double precision = 5.794737e-36 | |
861 | static const T Y = 0.52896785736083984375f; | |
862 | static const T P[] = { | |
863 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.00902152521745813634562524098263360074), | |
864 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0145207142776691539346923710537580927), | |
865 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0301681239582193983824211995978678571), | |
866 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0215548540823305814379020678660434461), | |
867 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00864683476267958365678294164340749949), | |
868 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00219693096885585491739823283511049902), | |
869 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000364961639163319762492184502159894371), | |
870 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.388174251026723752769264051548703059e-4), | |
871 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.241918026931789436000532513553594321e-5), | |
872 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.676586625472423508158937481943649258e-7), | |
873 | }; | |
874 | static const T Q[] = { | |
875 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
876 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.93669171363907292305550231764920001), | |
877 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.69468476144051356810672506101377494), | |
878 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.880023580986436640372794392579985511), | |
879 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.299099106711315090710836273697708402), | |
880 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0690593962363545715997445583603382337), | |
881 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0108427016361318921960863149875360222), | |
882 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00111747247208044534520499324234317695), | |
883 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.686843205749767250666787987163701209e-4), | |
884 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.192093541425429248675532015101904262e-5), | |
885 | }; | |
886 | result = Y + tools::evaluate_polynomial(P, T(z - 2.25f)) / tools::evaluate_polynomial(Q, T(z - 2.25f)); | |
92f5a8d4 TL |
887 | T hi, lo; |
888 | int expon; | |
889 | hi = floor(ldexp(frexp(z, &expon), 56)); | |
890 | hi = ldexp(hi, expon - 56); | |
891 | lo = z - hi; | |
892 | T sq = z * z; | |
893 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
894 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
895 | } |
896 | else if(z < 3.5) | |
897 | { | |
898 | // Maximum Deviation Found: 8.126e-37 | |
899 | // Expected Error Term: -8.126e-37 | |
900 | // Maximum Relative Change in Control Points: 1.363e-04 | |
901 | // Max Error found at long double precision = 1.747062e-36 | |
902 | static const T Y = 0.54037380218505859375f; | |
903 | static const T P[] = { | |
904 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.0033703486408887424921155540591370375), | |
905 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0104948043110005245215286678898115811), | |
906 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0148530118504000311502310457390417795), | |
907 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00816693029245443090102738825536188916), | |
908 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00249716579989140882491939681805594585), | |
909 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0004655591010047353023978045800916647), | |
910 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.531129557920045295895085236636025323e-4), | |
911 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.343526765122727069515775194111741049e-5), | |
912 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.971120407556888763695313774578711839e-7), | |
913 | }; | |
914 | static const T Q[] = { | |
915 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
916 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.59911256167540354915906501335919317), | |
917 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.136006830764025173864831382946934), | |
918 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.468565867990030871678574840738423023), | |
919 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.122821824954470343413956476900662236), | |
920 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0209670914950115943338996513330141633), | |
921 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00227845718243186165620199012883547257), | |
922 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000144243326443913171313947613547085553), | |
923 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.407763415954267700941230249989140046e-5), | |
924 | }; | |
925 | result = Y + tools::evaluate_polynomial(P, T(z - 3.0f)) / tools::evaluate_polynomial(Q, T(z - 3.0f)); | |
92f5a8d4 TL |
926 | T hi, lo; |
927 | int expon; | |
928 | hi = floor(ldexp(frexp(z, &expon), 56)); | |
929 | hi = ldexp(hi, expon - 56); | |
930 | lo = z - hi; | |
931 | T sq = z * z; | |
932 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
933 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
934 | } |
935 | else if(z < 5.5) | |
936 | { | |
937 | // Maximum Deviation Found: 5.804e-36 | |
938 | // Expected Error Term: -5.803e-36 | |
939 | // Maximum Relative Change in Control Points: 2.475e-05 | |
940 | // Max Error found at long double precision = 1.349545e-35 | |
941 | static const T Y = 0.55000019073486328125f; | |
942 | static const T P[] = { | |
943 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00118142849742309772151454518093813615), | |
944 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0072201822885703318172366893469382745), | |
945 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0078782276276860110721875733778481505), | |
946 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00418229166204362376187593976656261146), | |
947 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00134198400587769200074194304298642705), | |
948 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000283210387078004063264777611497435572), | |
949 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.405687064094911866569295610914844928e-4), | |
950 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.39348283801568113807887364414008292e-5), | |
951 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.248798540917787001526976889284624449e-6), | |
952 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.929502490223452372919607105387474751e-8), | |
953 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.156161469668275442569286723236274457e-9), | |
954 | }; | |
955 | static const T Q[] = { | |
956 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
957 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.52955245103668419479878456656709381), | |
958 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.06263944820093830054635017117417064), | |
959 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.441684612681607364321013134378316463), | |
960 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.121665258426166960049773715928906382), | |
961 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0232134512374747691424978642874321434), | |
962 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00310778180686296328582860464875562636), | |
963 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000288361770756174705123674838640161693), | |
964 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.177529187194133944622193191942300132e-4), | |
965 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.655068544833064069223029299070876623e-6), | |
966 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.11005507545746069573608988651927452e-7), | |
967 | }; | |
968 | result = Y + tools::evaluate_polynomial(P, T(z - 4.5f)) / tools::evaluate_polynomial(Q, T(z - 4.5f)); | |
92f5a8d4 TL |
969 | T hi, lo; |
970 | int expon; | |
971 | hi = floor(ldexp(frexp(z, &expon), 56)); | |
972 | hi = ldexp(hi, expon - 56); | |
973 | lo = z - hi; | |
974 | T sq = z * z; | |
975 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
976 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
977 | } |
978 | else if(z < 7.5) | |
979 | { | |
980 | // Maximum Deviation Found: 1.007e-36 | |
981 | // Expected Error Term: 1.007e-36 | |
982 | // Maximum Relative Change in Control Points: 1.027e-03 | |
983 | // Max Error found at long double precision = 2.646420e-36 | |
984 | static const T Y = 0.5574436187744140625f; | |
985 | static const T P[] = { | |
986 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000293236907400849056269309713064107674), | |
987 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00225110719535060642692275221961480162), | |
988 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00190984458121502831421717207849429799), | |
989 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000747757733460111743833929141001680706), | |
990 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000170663175280949889583158597373928096), | |
991 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.246441188958013822253071608197514058e-4), | |
992 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.229818000860544644974205957895688106e-5), | |
993 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.134886977703388748488480980637704864e-6), | |
994 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.454764611880548962757125070106650958e-8), | |
995 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.673002744115866600294723141176820155e-10), | |
996 | }; | |
997 | static const T Q[] = { | |
998 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
999 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.12843690320861239631195353379313367), | |
1000 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.569900657061622955362493442186537259), | |
1001 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.169094404206844928112348730277514273), | |
1002 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0324887449084220415058158657252147063), | |
1003 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00419252877436825753042680842608219552), | |
1004 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00036344133176118603523976748563178578), | |
1005 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.204123895931375107397698245752850347e-4), | |
1006 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.674128352521481412232785122943508729e-6), | |
1007 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.997637501418963696542159244436245077e-8), | |
1008 | }; | |
1009 | result = Y + tools::evaluate_polynomial(P, T(z - 6.5f)) / tools::evaluate_polynomial(Q, T(z - 6.5f)); | |
92f5a8d4 TL |
1010 | T hi, lo; |
1011 | int expon; | |
1012 | hi = floor(ldexp(frexp(z, &expon), 56)); | |
1013 | hi = ldexp(hi, expon - 56); | |
1014 | lo = z - hi; | |
1015 | T sq = z * z; | |
1016 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
1017 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
1018 | } |
1019 | else if(z < 11.5) | |
1020 | { | |
1021 | // Maximum Deviation Found: 8.380e-36 | |
1022 | // Expected Error Term: 8.380e-36 | |
1023 | // Maximum Relative Change in Control Points: 2.632e-06 | |
1024 | // Max Error found at long double precision = 9.849522e-36 | |
1025 | static const T Y = 0.56083202362060546875f; | |
1026 | static const T P[] = { | |
1027 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000282420728751494363613829834891390121), | |
1028 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00175387065018002823433704079355125161), | |
1029 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0021344978564889819420775336322920375), | |
1030 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00124151356560137532655039683963075661), | |
1031 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000423600733566948018555157026862139644), | |
1032 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.914030340865175237133613697319509698e-4), | |
1033 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.126999927156823363353809747017945494e-4), | |
1034 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.110610959842869849776179749369376402e-5), | |
1035 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.55075079477173482096725348704634529e-7), | |
1036 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.119735694018906705225870691331543806e-8), | |
1037 | }; | |
1038 | static const T Q[] = { | |
1039 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
1040 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.69889613396167354566098060039549882), | |
1041 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.28824647372749624464956031163282674), | |
1042 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.572297795434934493541628008224078717), | |
1043 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.164157697425571712377043857240773164), | |
1044 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.0315311145224594430281219516531649562), | |
1045 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00405588922155632380812945849777127458), | |
1046 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000336929033691445666232029762868642417), | |
1047 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.164033049810404773469413526427932109e-4), | |
1048 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.356615210500531410114914617294694857e-6), | |
1049 | }; | |
1050 | result = Y + tools::evaluate_polynomial(P, T(z / 2 - 4.75f)) / tools::evaluate_polynomial(Q, T(z / 2 - 4.75f)); | |
92f5a8d4 TL |
1051 | T hi, lo; |
1052 | int expon; | |
1053 | hi = floor(ldexp(frexp(z, &expon), 56)); | |
1054 | hi = ldexp(hi, expon - 56); | |
1055 | lo = z - hi; | |
1056 | T sq = z * z; | |
1057 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
1058 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
1059 | } |
1060 | else | |
1061 | { | |
1062 | // Maximum Deviation Found: 1.132e-35 | |
1063 | // Expected Error Term: -1.132e-35 | |
1064 | // Maximum Relative Change in Control Points: 4.674e-04 | |
1065 | // Max Error found at long double precision = 1.162590e-35 | |
1066 | static const T Y = 0.5632686614990234375f; | |
1067 | static const T P[] = { | |
1068 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.000920922048732849448079451574171836943), | |
1069 | BOOST_MATH_BIG_CONSTANT(T, 113, 0.00321439044532288750501700028748922439), | |
1070 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.250455263029390118657884864261823431), | |
1071 | BOOST_MATH_BIG_CONSTANT(T, 113, -0.906807635364090342031792404764598142), | |
1072 | BOOST_MATH_BIG_CONSTANT(T, 113, -8.92233572835991735876688745989985565), | |
1073 | BOOST_MATH_BIG_CONSTANT(T, 113, -21.7797433494422564811782116907878495), | |
1074 | BOOST_MATH_BIG_CONSTANT(T, 113, -91.1451915251976354349734589601171659), | |
1075 | BOOST_MATH_BIG_CONSTANT(T, 113, -144.1279109655993927069052125017673), | |
1076 | BOOST_MATH_BIG_CONSTANT(T, 113, -313.845076581796338665519022313775589), | |
1077 | BOOST_MATH_BIG_CONSTANT(T, 113, -273.11378811923343424081101235736475), | |
1078 | BOOST_MATH_BIG_CONSTANT(T, 113, -271.651566205951067025696102600443452), | |
1079 | BOOST_MATH_BIG_CONSTANT(T, 113, -60.0530577077238079968843307523245547), | |
1080 | }; | |
1081 | static const T Q[] = { | |
1082 | BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), | |
1083 | BOOST_MATH_BIG_CONSTANT(T, 113, 3.49040448075464744191022350947892036), | |
1084 | BOOST_MATH_BIG_CONSTANT(T, 113, 34.3563592467165971295915749548313227), | |
1085 | BOOST_MATH_BIG_CONSTANT(T, 113, 84.4993232033879023178285731843850461), | |
1086 | BOOST_MATH_BIG_CONSTANT(T, 113, 376.005865281206894120659401340373818), | |
1087 | BOOST_MATH_BIG_CONSTANT(T, 113, 629.95369438888946233003926191755125), | |
1088 | BOOST_MATH_BIG_CONSTANT(T, 113, 1568.35771983533158591604513304269098), | |
1089 | BOOST_MATH_BIG_CONSTANT(T, 113, 1646.02452040831961063640827116581021), | |
1090 | BOOST_MATH_BIG_CONSTANT(T, 113, 2299.96860633240298708910425594484895), | |
1091 | BOOST_MATH_BIG_CONSTANT(T, 113, 1222.73204392037452750381340219906374), | |
1092 | BOOST_MATH_BIG_CONSTANT(T, 113, 799.359797306084372350264298361110448), | |
1093 | BOOST_MATH_BIG_CONSTANT(T, 113, 72.7415265778588087243442792401576737), | |
1094 | }; | |
1095 | result = Y + tools::evaluate_polynomial(P, T(1 / z)) / tools::evaluate_polynomial(Q, T(1 / z)); | |
92f5a8d4 TL |
1096 | T hi, lo; |
1097 | int expon; | |
1098 | hi = floor(ldexp(frexp(z, &expon), 56)); | |
1099 | hi = ldexp(hi, expon - 56); | |
1100 | lo = z - hi; | |
1101 | T sq = z * z; | |
1102 | T err_sqr = ((hi * hi - sq) + 2 * hi * lo) + lo * lo; | |
1103 | result *= exp(-sq) * exp(-err_sqr) / z; | |
7c673cae FG |
1104 | } |
1105 | } | |
1106 | else | |
1107 | { | |
1108 | // | |
1109 | // Any value of z larger than 110 will underflow to zero: | |
1110 | // | |
1111 | result = 0; | |
1112 | invert = !invert; | |
1113 | } | |
1114 | ||
1115 | if(invert) | |
1116 | { | |
1117 | result = 1 - result; | |
1118 | } | |
1119 | ||
1120 | return result; | |
1e59de90 | 1121 | } // template <class T, class Lanczos>T erf_imp(T z, bool invert, const Lanczos& l, const std::integral_constant<int, 113>& t) |
7c673cae FG |
1122 | |
1123 | template <class T, class Policy, class tag> | |
1124 | struct erf_initializer | |
1125 | { | |
1126 | struct init | |
1127 | { | |
1128 | init() | |
1129 | { | |
1130 | do_init(tag()); | |
1131 | } | |
1e59de90 TL |
1132 | static void do_init(const std::integral_constant<int, 0>&){} |
1133 | static void do_init(const std::integral_constant<int, 53>&) | |
7c673cae FG |
1134 | { |
1135 | boost::math::erf(static_cast<T>(1e-12), Policy()); | |
1136 | boost::math::erf(static_cast<T>(0.25), Policy()); | |
1137 | boost::math::erf(static_cast<T>(1.25), Policy()); | |
1138 | boost::math::erf(static_cast<T>(2.25), Policy()); | |
1139 | boost::math::erf(static_cast<T>(4.25), Policy()); | |
1140 | boost::math::erf(static_cast<T>(5.25), Policy()); | |
1141 | } | |
1e59de90 | 1142 | static void do_init(const std::integral_constant<int, 64>&) |
7c673cae FG |
1143 | { |
1144 | boost::math::erf(static_cast<T>(1e-12), Policy()); | |
1145 | boost::math::erf(static_cast<T>(0.25), Policy()); | |
1146 | boost::math::erf(static_cast<T>(1.25), Policy()); | |
1147 | boost::math::erf(static_cast<T>(2.25), Policy()); | |
1148 | boost::math::erf(static_cast<T>(4.25), Policy()); | |
1149 | boost::math::erf(static_cast<T>(5.25), Policy()); | |
1150 | } | |
1e59de90 | 1151 | static void do_init(const std::integral_constant<int, 113>&) |
7c673cae FG |
1152 | { |
1153 | boost::math::erf(static_cast<T>(1e-22), Policy()); | |
1154 | boost::math::erf(static_cast<T>(0.25), Policy()); | |
1155 | boost::math::erf(static_cast<T>(1.25), Policy()); | |
1156 | boost::math::erf(static_cast<T>(2.125), Policy()); | |
1157 | boost::math::erf(static_cast<T>(2.75), Policy()); | |
1158 | boost::math::erf(static_cast<T>(3.25), Policy()); | |
1159 | boost::math::erf(static_cast<T>(5.25), Policy()); | |
1160 | boost::math::erf(static_cast<T>(7.25), Policy()); | |
1161 | boost::math::erf(static_cast<T>(11.25), Policy()); | |
1162 | boost::math::erf(static_cast<T>(12.5), Policy()); | |
1163 | } | |
1164 | void force_instantiate()const{} | |
1165 | }; | |
1166 | static const init initializer; | |
1167 | static void force_instantiate() | |
1168 | { | |
1169 | initializer.force_instantiate(); | |
1170 | } | |
1171 | }; | |
1172 | ||
1173 | template <class T, class Policy, class tag> | |
1174 | const typename erf_initializer<T, Policy, tag>::init erf_initializer<T, Policy, tag>::initializer; | |
1175 | ||
1176 | } // namespace detail | |
1177 | ||
1178 | template <class T, class Policy> | |
1179 | inline typename tools::promote_args<T>::type erf(T z, const Policy& /* pol */) | |
1180 | { | |
1181 | typedef typename tools::promote_args<T>::type result_type; | |
1182 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1183 | typedef typename policies::precision<result_type, Policy>::type precision_type; | |
1184 | typedef typename policies::normalise< | |
1185 | Policy, | |
1186 | policies::promote_float<false>, | |
1187 | policies::promote_double<false>, | |
1188 | policies::discrete_quantile<>, | |
1189 | policies::assert_undefined<> >::type forwarding_policy; | |
1190 | ||
1191 | BOOST_MATH_INSTRUMENT_CODE("result_type = " << typeid(result_type).name()); | |
1192 | BOOST_MATH_INSTRUMENT_CODE("value_type = " << typeid(value_type).name()); | |
1193 | BOOST_MATH_INSTRUMENT_CODE("precision_type = " << typeid(precision_type).name()); | |
1194 | ||
1e59de90 | 1195 | typedef std::integral_constant<int, |
f67539c2 TL |
1196 | precision_type::value <= 0 ? 0 : |
1197 | precision_type::value <= 53 ? 53 : | |
1198 | precision_type::value <= 64 ? 64 : | |
1199 | precision_type::value <= 113 ? 113 : 0 | |
1200 | > tag_type; | |
7c673cae FG |
1201 | |
1202 | BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name()); | |
1203 | ||
1204 | detail::erf_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); // Force constants to be initialized before main | |
1205 | ||
1206 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::erf_imp( | |
1207 | static_cast<value_type>(z), | |
1208 | false, | |
1209 | forwarding_policy(), | |
1210 | tag_type()), "boost::math::erf<%1%>(%1%, %1%)"); | |
1211 | } | |
1212 | ||
1213 | template <class T, class Policy> | |
1214 | inline typename tools::promote_args<T>::type erfc(T z, const Policy& /* pol */) | |
1215 | { | |
1216 | typedef typename tools::promote_args<T>::type result_type; | |
1217 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1218 | typedef typename policies::precision<result_type, Policy>::type precision_type; | |
1219 | typedef typename policies::normalise< | |
1220 | Policy, | |
1221 | policies::promote_float<false>, | |
1222 | policies::promote_double<false>, | |
1223 | policies::discrete_quantile<>, | |
1224 | policies::assert_undefined<> >::type forwarding_policy; | |
1225 | ||
1226 | BOOST_MATH_INSTRUMENT_CODE("result_type = " << typeid(result_type).name()); | |
1227 | BOOST_MATH_INSTRUMENT_CODE("value_type = " << typeid(value_type).name()); | |
1228 | BOOST_MATH_INSTRUMENT_CODE("precision_type = " << typeid(precision_type).name()); | |
1229 | ||
1e59de90 | 1230 | typedef std::integral_constant<int, |
f67539c2 TL |
1231 | precision_type::value <= 0 ? 0 : |
1232 | precision_type::value <= 53 ? 53 : | |
1233 | precision_type::value <= 64 ? 64 : | |
1234 | precision_type::value <= 113 ? 113 : 0 | |
1235 | > tag_type; | |
7c673cae FG |
1236 | |
1237 | BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name()); | |
1238 | ||
1239 | detail::erf_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); // Force constants to be initialized before main | |
1240 | ||
1241 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::erf_imp( | |
1242 | static_cast<value_type>(z), | |
1243 | true, | |
1244 | forwarding_policy(), | |
1245 | tag_type()), "boost::math::erfc<%1%>(%1%, %1%)"); | |
1246 | } | |
1247 | ||
1248 | template <class T> | |
1249 | inline typename tools::promote_args<T>::type erf(T z) | |
1250 | { | |
1251 | return boost::math::erf(z, policies::policy<>()); | |
1252 | } | |
1253 | ||
1254 | template <class T> | |
1255 | inline typename tools::promote_args<T>::type erfc(T z) | |
1256 | { | |
1257 | return boost::math::erfc(z, policies::policy<>()); | |
1258 | } | |
1259 | ||
1260 | } // namespace math | |
1261 | } // namespace boost | |
1262 | ||
1263 | #include <boost/math/special_functions/detail/erf_inv.hpp> | |
1264 | ||
1265 | #endif // BOOST_MATH_SPECIAL_ERF_HPP | |
1266 | ||
1267 | ||
1268 | ||
1269 |