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

source code of boost/libs/math/include/boost/math/special_functions/erf.hpp