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