Ticket #8621: erf.hpp

File erf.hpp, 53.0 KB (added by John Maddock, 9 years ago)

Debugging version of erf.hpp

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