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
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 {
erf_asympt_series_tboost::math::detail::erf_asympt_series_t41 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
operator ()boost::math::detail::erf_asympt_series_t50 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>
erf_asymptotic_limit_N(const T &)69 inline float erf_asymptotic_limit_N(const T&)
70 {
71 return (std::numeric_limits<float>::max)();
72 }
erf_asymptotic_limit_N(const boost::integral_constant<int,24> &)73 inline float erf_asymptotic_limit_N(const boost::integral_constant<int, 24>&)
74 {
75 return 2.8F;
76 }
erf_asymptotic_limit_N(const boost::integral_constant<int,53> &)77 inline float erf_asymptotic_limit_N(const boost::integral_constant<int, 53>&)
78 {
79 return 4.3F;
80 }
erf_asymptotic_limit_N(const boost::integral_constant<int,64> &)81 inline float erf_asymptotic_limit_N(const boost::integral_constant<int, 64>&)
82 {
83 return 4.8F;
84 }
erf_asymptotic_limit_N(const boost::integral_constant<int,106> &)85 inline float erf_asymptotic_limit_N(const boost::integral_constant<int, 106>&)
86 {
87 return 6.5F;
88 }
erf_asymptotic_limit_N(const boost::integral_constant<int,113> &)89 inline float erf_asymptotic_limit_N(const boost::integral_constant<int, 113>&)
90 {
91 return 6.8F;
92 }
93
94 template <class T, class Policy>
erf_asymptotic_limit()95 inline T erf_asymptotic_limit()
96 {
97 typedef typename policies::precision<T, Policy>::type precision_type;
98 typedef boost::integral_constant<int,
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;
105 return erf_asymptotic_limit_N(tag_type());
106 }
107
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;
erf_series_near_zeroboost::math::detail::erf_series_near_zero115 erf_series_near_zero(const T& z) : term(z), zz(-z * z), k(0) {}
116
operator ()boost::math::detail::erf_series_near_zero117 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);
133 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
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
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);
159 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
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;
166 if(z < 1.3f)
167 {
168 // Compute P:
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);
172 }
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 }
179 else
180 {
181 // Compute Q:
182 invert = !invert;
183 result = z * exp(-x);
184 result /= boost::math::constants::root_pi<T>();
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>
194 T erf_imp(T z, bool invert, const Policy& pol, const boost::integral_constant<int, 53>& t)
195 {
196 BOOST_MATH_STD_USING
197
198 BOOST_MATH_INSTRUMENT_CODE("53-bit precision erf_imp called");
199
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
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 }
263 else if(invert ? (z < 28) : (z < 5.8f))
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));
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;
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));
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;
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));
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;
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;
419 } // template <class T, class Lanczos>T erf_imp(T z, bool invert, const Lanczos& l, const boost::integral_constant<int, 53>& t)
420
421
422 template <class T, class Policy>
423 T erf_imp(T z, bool invert, const Policy& pol, const boost::integral_constant<int, 64>& t)
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 }
486 else if(invert ? (z < 110) : (z < 6.4f))
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));
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;
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));
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;
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));
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;
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));
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;
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;
653 } // template <class T, class Lanczos>T erf_imp(T z, bool invert, const Lanczos& l, const boost::integral_constant<int, 64>& t)
654
655
656 template <class T, class Policy>
657 T erf_imp(T z, bool invert, const Policy& pol, const boost::integral_constant<int, 113>& t)
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));
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;
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));
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;
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));
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;
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));
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;
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));
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;
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));
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;
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));
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;
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));
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;
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));
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;
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;
1121 } // template <class T, class Lanczos>T erf_imp(T z, bool invert, const Lanczos& l, const boost::integral_constant<int, 113>& t)
1122
1123 template <class T, class Policy, class tag>
1124 struct erf_initializer
1125 {
1126 struct init
1127 {
initboost::math::detail::erf_initializer::init1128 init()
1129 {
1130 do_init(tag());
1131 }
do_initboost::math::detail::erf_initializer::init1132 static void do_init(const boost::integral_constant<int, 0>&){}
do_initboost::math::detail::erf_initializer::init1133 static void do_init(const boost::integral_constant<int, 53>&)
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 }
do_initboost::math::detail::erf_initializer::init1142 static void do_init(const boost::integral_constant<int, 64>&)
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 }
do_initboost::math::detail::erf_initializer::init1151 static void do_init(const boost::integral_constant<int, 113>&)
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 }
force_instantiateboost::math::detail::erf_initializer::init1164 void force_instantiate()const{}
1165 };
1166 static const init initializer;
force_instantiateboost::math::detail::erf_initializer1167 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>
erf(T z,const 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
1195 typedef boost::integral_constant<int,
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;
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>
erfc(T z,const 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
1230 typedef boost::integral_constant<int,
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;
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>
erf(T z)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>
erfc(T z)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
1270