]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // (C) Copyright John Maddock 2006. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #ifndef BOOST_MATH_SPECIAL_ERF_HPP | |
7 | #define BOOST_MATH_SPECIAL_ERF_HPP | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #endif | |
12 | ||
13 | #include <boost/math/special_functions/math_fwd.hpp> | |
14 | #include <boost/math/tools/config.hpp> | |
15 | #include <boost/math/special_functions/gamma.hpp> | |
16 | #include <boost/math/tools/roots.hpp> | |
17 | #include <boost/math/policies/error_handling.hpp> | |
18 | #include <boost/math/tools/big_constant.hpp> | |
19 | ||
20 | namespace boost{ namespace math{ | |
21 | ||
22 | namespace detail | |
23 | { | |
24 | ||
25 | // | |
26 | // Asymptotic series for large z: | |
27 | // | |
28 | template <class T> | |
29 | struct 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 | } | |
50 | private: | |
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 | // | |
58 | template <class T> | |
59 | inline float erf_asymptotic_limit_N(const T&) | |
60 | { | |
61 | return (std::numeric_limits<float>::max)(); | |
62 | } | |
63 | inline float erf_asymptotic_limit_N(const mpl::int_<24>&) | |
64 | { | |
65 | return 2.8F; | |
66 | } | |
67 | inline float erf_asymptotic_limit_N(const mpl::int_<53>&) | |
68 | { | |
69 | return 4.3F; | |
70 | } | |
71 | inline float erf_asymptotic_limit_N(const mpl::int_<64>&) | |
72 | { | |
73 | return 4.8F; | |
74 | } | |
75 | inline float erf_asymptotic_limit_N(const mpl::int_<106>&) | |
76 | { | |
77 | return 6.5F; | |
78 | } | |
79 | inline float erf_asymptotic_limit_N(const mpl::int_<113>&) | |
80 | { | |
81 | return 6.8F; | |
82 | } | |
83 | ||
84 | template <class T, class Policy> | |
85 | inline 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 | ||
116 | template <class T, class Policy, class Tag> | |
117 | T 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 | ||
172 | template <class T, class Policy> | |
173 | T 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 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.0), | |
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 | BOOST_MATH_BIG_CONSTANT(T, 53, 1.0), | |
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.0), | |
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.0), | |
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.0), | |
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 | ||
377 | template <class T, class Policy> | |
378 | T 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.0), | |
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.0), | |
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.0), | |
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.0), | |
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.0), | |
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 | ||
583 | template <class T, class Policy> | |
584 | T 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.0), | |
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.0), | |
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.0), | |
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.0), | |
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.0), | |
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.0), | |
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.0), | |
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.0), | |
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.0), | |
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.0), | |
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 | ||
987 | template <class T, class Policy, class tag> | |
988 | struct 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 | ||
1037 | template <class T, class Policy, class tag> | |
1038 | const typename erf_initializer<T, Policy, tag>::init erf_initializer<T, Policy, tag>::initializer; | |
1039 | ||
1040 | } // namespace detail | |
1041 | ||
1042 | template <class T, class Policy> | |
1043 | inline 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 | ||
1088 | template <class T, class Policy> | |
1089 | inline 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 | ||
1134 | template <class T> | |
1135 | inline typename tools::promote_args<T>::type erf(T z) | |
1136 | { | |
1137 | return boost::math::erf(z, policies::policy<>()); | |
1138 | } | |
1139 | ||
1140 | template <class T> | |
1141 | inline 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 |