boost/math/special_functions/erf.hpp
// (C) Copyright John Maddock 2006. // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_MATH_SPECIAL_ERF_HPP #define BOOST_MATH_SPECIAL_ERF_HPP #include <boost/math/special_functions/math_fwd.hpp> #include <boost/math/tools/config.hpp> #include <boost/math/special_functions/gamma.hpp> #include <boost/math/tools/roots.hpp> #include <boost/math/policies/error_handling.hpp> namespace boost{ namespace math{ namespace detail { // // Asymptotic series for large z: // template <class T> struct erf_asympt_series_t { erf_asympt_series_t(T z) : xx(2 * -z * z), tk(1) { BOOST_MATH_STD_USING result = -exp(-z * z) / sqrt(boost::math::constants::pi<T>()); result /= z; } typedef T result_type; T operator()() { BOOST_MATH_STD_USING T r = result; result *= tk / xx; tk += 2; if( fabs(r) < fabs(result)) result = 0; return r; } private: T result; T xx; int tk; }; // // How large z has to be in order to ensure that the series converges: // template <class T> inline float erf_asymptotic_limit_N(const T&) { return (std::numeric_limits<float>::max)(); } inline float erf_asymptotic_limit_N(const mpl::int_<24>&) { return 2.8F; } inline float erf_asymptotic_limit_N(const mpl::int_<53>&) { return 4.3F; } inline float erf_asymptotic_limit_N(const mpl::int_<64>&) { return 4.8F; } inline float erf_asymptotic_limit_N(const mpl::int_<106>&) { return 6.5F; } inline float erf_asymptotic_limit_N(const mpl::int_<113>&) { return 6.8F; } template <class T, class Policy> inline T erf_asymptotic_limit() { typedef typename policies::precision<T, Policy>::type precision_type; typedef typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<24> >, typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<0> >, mpl::int_<0>, mpl::int_<24> >::type, typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<53> >, mpl::int_<53>, typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<64> >, mpl::int_<64>, typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<106> >, mpl::int_<106>, typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<113> >, mpl::int_<113>, mpl::int_<0> >::type >::type >::type >::type >::type tag_type; return erf_asymptotic_limit_N(tag_type()); } template <class T, class Policy, class Tag> T erf_imp(T z, bool invert, const Policy& pol, const Tag& t) { BOOST_MATH_STD_USING BOOST_MATH_INSTRUMENT_CODE("Generic erf_imp called"); if(z < 0) { if(!invert) return -erf_imp(-z, invert, pol, t); else return 1 + erf_imp(-z, false, pol, t); } T result; if(!invert && (z > detail::erf_asymptotic_limit<T, Policy>())) { detail::erf_asympt_series_t<T> s(z); boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); result = boost::math::tools::sum_series(s, policies::digits<T, Policy>(), max_iter, 1); policies::check_series_iterations("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol); } else { T x = z * z; if(x < 0.6) { // Compute P: result = z * exp(-x); result /= sqrt(boost::math::constants::pi<T>()); if(result != 0) result *= 2 * detail::lower_gamma_series(T(0.5f), x, pol); } else if(x < 1.1f) { // Compute Q: invert = !invert; result = tgamma_small_upper_part(T(0.5f), x, pol); result /= sqrt(boost::math::constants::pi<T>()); } else { // Compute Q: invert = !invert; result = z * exp(-x); result /= sqrt(boost::math::constants::pi<T>()); result *= upper_gamma_fraction(T(0.5f), x, policies::digits<T, Policy>()); } } if(invert) result = 1 - result; return result; } template <class T, class Policy> T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<53>& t) { BOOST_MATH_STD_USING BOOST_MATH_INSTRUMENT_CODE("53-bit precision erf_imp called"); if(z < 0) { if(!invert) return -erf_imp(-z, invert, pol, t); else if(z < -0.5) return 2 - erf_imp(-z, invert, pol, t); else return 1 + erf_imp(-z, false, pol, t); } T result; // // Big bunch of selection statements now to pick // which implementation to use, // try to put most likely options first: // if(z < 0.5) { // // We're going to calculate erf: // if(z == 0) { result = T(0); } else if(z < 1e-10) { result = static_cast<T>(z * 1.125f + z * 0.003379167095512573896158903121545171688L); } else { static const T n[7] = { 0.00337916709551257778174L, -0.000147024115786688745475L, -0.37463022236812520164L, 0.0163061594494816999803L, -0.0534354147807331748737L, 0.00161898096813581982844L, -0.0059528010489182840404L }; static const T d[7] = { 1, -0.0435089806536379531594L, 0.442761965043509204727L, -0.017375974533016704678L, 0.0772756490303260060769L, -0.00210552465858669941879L, 0.00544772980263244037286L }; result = static_cast<T>(z * 1.125f + z * tools::evaluate_polynomial(n, z) / tools::evaluate_polynomial(d, z)); } } else if((z < 14) || ((z < 28) && invert)) { // // We'll be calculating erfc: // invert = !invert; T r, b; if(z < 0.75) { // Worst case absolute error found: 8.554649561e-018 static const T n[5] = { -0.0361790390718262468222L, 0.301888464724047222196L, 0.201731143672633894981L, 0.0659353268087389983319L, 0.00721876720062364930761L }; static const T d[6] = { 1, 1.58814245739127341535L, 0.99354580430196422336L, 0.291753007176902027213L, 0.033994791234913855515L, -0.000104234653166533504303L }; static const float f0 = 0.3440242112F; r = tools::evaluate_polynomial(n, z - 0.5) / tools::evaluate_polynomial(d, z - 0.5); b = f0; } else if(z < 1.25) { // Worst case absolute error found: 6.50251514e-018 static const T n[6] = { -0.039787689261113685983L, 0.160309168830518003303L, 0.163049978514596540313L, 0.0710685660158400750009L, 0.01497188097404877543L, 0.00130080628375002584279L }; static const T d[6] = { 1, 1.77564880074171280407L, 1.31438791181040008779L, 0.509359151038517059748L, 0.103958527905812829559L, 0.00901292460643094469406L }; static const float f0 = 0.419990927F; r = tools::evaluate_polynomial(n, z - 0.75) / tools::evaluate_polynomial(d, z - 0.75); b = f0; } else if(z < 2.25) { // Worst case absolute error found: 1.132743504e-017 static const T n[6] = { -0.0300838560557949724172L, 0.0592886319615167248092L, 0.0622294724048409148736L, 0.0248575228109427909578L, 0.00463781847004901844581L, 0.000347305179334822548368L }; static const T d[7] = { 1, 1.57915060645728571344L, 1.03342495188878679417L, 0.35158678814344218974L, 0.062469256580984456783L, 0.00466640448020624599948L, 0.290106403940303572448e-6L }; static const float f0 = 0.4898625016F; r = tools::evaluate_polynomial(n, z - 1.25) / tools::evaluate_polynomial(d, z - 1.25); b = f0; } else if(z < 3.5) { // Worst case absolute error found: 3.446364609e-018 static const T n[6] = { -0.0117907570137227857015L, 0.0162667227692515660221L, 0.0175329212378413544794L, 0.00620897681269247137578L, 0.000986614895094589251706L, 0.601354618401624353425e-4L }; static const T d[6] = { 1, 1.33374851361555383557L, 0.73227756904205983415L, 0.207410266363727673685L, 0.0304034048466731110163L, 0.00185296959991832048613L }; static const float f0 = 0.5317370892F; r = tools::evaluate_polynomial(n, z - 2.25) / tools::evaluate_polynomial(d, z - 2.25); b = f0; } else if(z < 5.5) { // Worst case absolute error found: 1.579588208e-018 static const T n[6] = { -0.00588219091116732271979L, 0.00434428684527812140098L, 0.00466899990542371512895L, 0.00139937567253199794533L, 0.000179205902444982389766L, 0.845033527560949509345e-5L }; static const T d[6] = { 1, 1.07389345953392962127L, 0.470965611895885060643L, 0.105594730223366124873L, 0.0121252833787344059719L, 0.000571755036133730341579L }; static const float f0 = 0.5494099855F; r = tools::evaluate_polynomial(n, z - 3.5) / tools::evaluate_polynomial(d, z - 3.5); b = f0; } else if(z < 9) { // Worst case absolute error found: 1.410768708e-017 static const T n[5] = { -0.00273864253749621265032L, 0.0013089921066773026803L, 0.000775841526778089659703L, 0.000110909476102006410909L, 0.472577590124068298534e-5L }; static const T d[6] = { 1, 0.650694792327863647878L, 0.161126734432670927888L, 0.0180081468446110640846L, 0.000767341359508884026192L, -0.287636719206664167616e-9L }; static const float f0 = 0.5580308437F; r = tools::evaluate_polynomial(n, z - 5.5) / tools::evaluate_polynomial(d, z - 5.5); b = f0; } else if(z < 14) { // Worst case absolute error found: 1.458310511e-018 static const T n[5] = { -0.000995856413171151859346L, 0.000320252910249376187643L, 0.000129085624923151780987L, 0.121577881306587454509e-4L, 0.33293110334156470348e-6L }; static const T d[5] = { 1, 0.428034987547594828954L, 0.0692297359775940896439L, 0.00501515176145997560701L, 0.00013733589151338416322L }; static const float f0 = 0.5617653728F; r = tools::evaluate_polynomial(n, z - 9) / tools::evaluate_polynomial(d, z - 9); b = f0; } else if(z < 21) { // Worst case absolute error found: 1.08182873e-019 static const T n[5] = { -0.000395463268432048215535L, 0.91155953112698182321e-4L, 0.237451641259281193813e-4L, 0.145759953022524466816e-5L, 0.259395907606548998142e-7L }; static const T d[5] = { 1, 0.281604524251560309285L, 0.0298468482900092392397L, 0.00141114575715338885136L, 0.251128951158576064819e-4L }; static const float f0 = 0.5631566644F; r = tools::evaluate_polynomial(n, z - 14) / tools::evaluate_polynomial(d, z - 14); b = f0; } else { // Worst case absolute error found: 7.010370259e-018 static const T n[4] = { -0.000139182098873874523526L, 0.395254617101737287826e-4L, 0.376801239136290345387e-5L, 0.629017242098850415839e-7L }; static const T d[4] = { 1, 0.15077096006891495258L, 0.00756136203065884121997L, 0.000126226197336507576933L }; static const float f0 = 0.5636912584F; r = tools::evaluate_polynomial(n, z - 21) / tools::evaluate_polynomial(d, z - 21); b = f0; } T g = exp(-z * z) / z; result = g * b + g * r; } else { // // Any value of z larger than 28 will underflow to zero: // result = 0; invert = !invert; } if(invert) { result = 1 - result; } return result; } // template <class T, class L>T erf_imp(T z, bool invert, const L& l, const mpl::int_<53>& t) template <class T, class Policy> T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<64>& t) { BOOST_MATH_STD_USING BOOST_MATH_INSTRUMENT_CODE("64-bit precision erf_imp called"); if(z < 0) { if(!invert) return -erf_imp(-z, invert, pol, t); else if(z < -0.5) return 2 - erf_imp(-z, invert, pol, t); else return 1 + erf_imp(-z, false, pol, t); } T result; // // Big bunch of selection statements now to pick which // implementation to use, try to put most likely options // first: // if(z < 0.5) { // // We're going to calculate erf: // if(z == 0) { result = 0; } else if(z < 1e-10) { result = z * 1.125 + z * 0.003379167095512573896158903121545171688L; } else { // Worst case absolute error found: 6.688618532e-21 static const T n[8] = { 0.00337916709551257388990745L, -0.00073695653048167948530905L, -0.374732337392919607868241L, 0.0817442448733587196071743L, -0.0421089319936548595203468L, 0.0070165709512095756344528L, -0.00495091255982435110337458L, 0.000871646599037922480317225L }; static const T d[8] = { 1L, -0.218088218087924645390535L, 0.412542972725442099083918L, -0.0841891147873106755410271L, 0.0655338856400241519690695L, -0.0120019604454941768171266L, 0.00408165558926174048329689L, -0.000615900721557769691924509L }; result = z * 1.125 + z * tools::evaluate_polynomial(n, z) / tools::evaluate_polynomial(d, z); } } else if((z < 110) || ((z < 110) && invert)) // TODO FIXME!!! { // // We'll be calculating erfc: // invert = !invert; T r, b; if(z < 0.75) { // Worst case absolute error found: 5.582813374e-21 static const T n[6] = { -0.0361790390718262471360258L, 0.292251883444882683221149L, 0.281447041797604512774415L, 0.125610208862766947294894L, 0.0274135028268930549240776L, 0.00250839672168065762786937L }; static const T d[6] = { 1L, 1.8545005897903486499845L, 1.43575803037831418074962L, 0.582827658753036572454135L, 0.124810476932949746447682L, 0.0113724176546353285778481L }; static const float f0 = 0.3440242112F; r = tools::evaluate_polynomial(n, z - 0.5) / tools::evaluate_polynomial(d, z - 0.5); b = f0; } else if(z < 1.25) { // Worst case absolute error found: 4.01854729e-21 static const T n[7] = { -0.0397876892611136856954425L, 0.153165212467878293257683L, 0.191260295600936245503129L, 0.10276327061989304213645L, 0.029637090615738836726027L, 0.0046093486780275489468812L, 0.000307607820348680180548455L }; static const T d[7] = { 1L, 1.95520072987627704987886L, 1.64762317199384860109595L, 0.768238607022126250082483L, 0.209793185936509782784315L, 0.0319569316899913392596356L, 0.00213363160895785378615014L }; static const float f0 = 0.419990927F; r = tools::evaluate_polynomial(n, z - 0.75) / tools::evaluate_polynomial(d, z - 0.75); b = f0; } else if(z < 2.25) { // Worst case absolute error found: 2.866005373e-21 static const T n[7] = { -0.0300838560557949717328341L, 0.0538578829844454508530552L, 0.0726211541651914182692959L, 0.0367628469888049348429018L, 0.00964629015572527529605267L, 0.00133453480075291076745275L, 0.778087599782504251917881e-4L }; static const T d[8] = { 1L, 1.75967098147167528287343L, 1.32883571437961120556307L, 0.552528596508757581287907L, 0.133793056941332861912279L, 0.0179509645176280768640766L, 0.00104712440019937356634038L, -0.106640381820357337177643e-7L }; static const float f0 = 0.4898625016F; r = tools::evaluate_polynomial(n, z - 1.25) / tools::evaluate_polynomial(d, z - 1.25); b = f0; } else if(z < 3.5) { // Worst case absolute error found: 1.045355789e-21 static const T n[7] = { -0.0117907570137227847827732L, 0.014262132090538809896674L, 0.0202234435902960820020765L, 0.00930668299990432009042239L, 0.00213357802422065994322516L, 0.00025022987386460102395382L, 0.120534912219588189822126e-4L }; static const T d[7] = { 1L, 1.50376225203620482047419L, 0.965397786204462896346934L, 0.339265230476796681555511L, 0.0689740649541569716897427L, 0.00771060262491768307365526L, 0.000371421101531069302990367L }; static const float f0 = 0.5317370892F; r = tools::evaluate_polynomial(n, z - 2.25) / tools::evaluate_polynomial(d, z - 2.25); b = f0; } else if(z < 5.25) { // Worst case absolute error found: 8.300028706e-22 static const T n[7] = { -0.00546954795538729307482955L, 0.00404190278731707110245394L, 0.0054963369553161170521356L, 0.00212616472603945399437862L, 0.000394984014495083900689956L, 0.365565477064442377259271e-4L, 0.135485897109932323253786e-5L }; static const T d[8] = { 1L, 1.21019697773630784832251L, 0.620914668221143886601045L, 0.173038430661142762569515L, 0.0276550813773432047594539L, 0.00240625974424309709745382L, 0.891811817251336577241006e-4L, -0.465528836283382684461025e-11L }; static const float f0 = 0.5489973426F; r = tools::evaluate_polynomial(n, z - 3.5) / tools::evaluate_polynomial(d, z - 3.5); b = f0; } else if(z < 8) { // Worst case absolute error found: 1.700157534e-21 static const T n[6] = { -0.00270722535905778347999196L, 0.0013187563425029400461378L, 0.00119925933261002333923989L, 0.00027849619811344664248235L, 0.267822988218331849989363e-4L, 0.923043672315028197865066e-6L }; static const T d[7] = { 1L, 0.814632808543141591118279L, 0.268901665856299542168425L, 0.0449877216103041118694989L, 0.00381759663320248459168994L, 0.000131571897888596914350697L, 0.404815359675764138445257e-11L }; static const float f0 = 0.5571740866F; r = tools::evaluate_polynomial(n, z - 5.25) / tools::evaluate_polynomial(d, z - 5.25); b = f0; } else if(z < 11.5) { //Worst case absolute error found: 3.002278011e-22 static const T n[6] = { -0.00109946720691742196814323L, 0.000406425442750422675169153L, 0.000274499489416900707787024L, 0.465293770646659383436343e-4L, 0.320955425395767463401993e-5L, 0.778286018145020892261936e-7L }; static const T d[6] = { 1L, 0.588173710611846046373373L, 0.139363331289409746077541L, 0.0166329340417083678763028L, 0.00100023921310234908642639L, 0.24254837521587225125068e-4L }; static const float f0 = 0.5609807968F; r = tools::evaluate_polynomial(n, z - 8) / tools::evaluate_polynomial(d, z - 8); b = f0; } else if(z < 17) { //Worst case absolute error found: 6.741114695e-21 static const T n[5] = { -0.00056907993601094962855594L, 0.000169498540373762264416984L, 0.518472354581100890120501e-4L, 0.382819312231928859704678e-5L, 0.824989931281894431781794e-7L }; static const T d[6] = { 1L, 0.339637250051139347430323L, 0.043472647870310663055044L, 0.00248549335224637114641629L, 0.535633305337152900549536e-4L, -0.117490944405459578783846e-12L }; static const float f0 = 0.5626493692F; r = tools::evaluate_polynomial(n, z - 11.5) / tools::evaluate_polynomial(d, z - 11.5); b = f0; } else if(z < 24) { // Worst case absolute error found: 7.802346984e-22 static const T n[5] = { -0.000241313599483991337479091L, 0.574224975202501512365975e-4L, 0.115998962927383778460557e-4L, 0.581762134402593739370875e-6L, 0.853971555085673614607418e-8L }; static const T d[5] = { 1L, 0.233044138299687841018015L, 0.0204186940546440312625597L, 0.000797185647564398289151125L, 0.117019281670172327758019e-4L }; static const float f0 = 0.5634598136F; r = tools::evaluate_polynomial(n, z - 17) / tools::evaluate_polynomial(d, z - 17); b = f0; } else if(z < 38) { // Worst case absolute error found: 2.414228989e-22 static const T n[5] = { -0.000146674699277760365803642L, 0.162666552112280519955647e-4L, 0.269116248509165239294897e-5L, 0.979584479468091935086972e-7L, 0.101994647625723465722285e-8L }; static const T d[5] = { 1L, 0.165907812944847226546036L, 0.0103361716191505884359634L, 0.000286593026373868366935721L, 0.298401570840900340874568e-5L }; static const float f0 = 0.5638477802F; r = tools::evaluate_polynomial(n, z - 24) / tools::evaluate_polynomial(d, z - 24); b = f0; } else if(z < 60) { // Worst case absolute error found: 5.896543869e-24 static const T n[5] = { -0.583905797629771786720406e-4L, 0.412510325105496173512992e-5L, 0.431790922420250949096906e-6L, 0.993365155590013193345569e-8L, 0.653480510020104699270084e-10L }; static const T d[5] = { 1L, 0.105077086072039915406159L, 0.00414278428675475620830226L, 0.726338754644523769144108e-4L, 0.477818471047398785369849e-6L }; static const float f0 = 0.5640528202F; r = tools::evaluate_polynomial(n, z - 38) / tools::evaluate_polynomial(d, z - 38); b = f0; } else if(z < 85) { // Worst case absolute error found: 3.080612264e-21 static const T n[4] = { -0.196457797609229579459841e-4L, 0.157243887666800692441195e-5L, 0.543902511192700878690335e-7L, 0.317472492369117710852685e-9L }; static const T d[5] = { 1L, 0.052803989240957632204885L, 0.000926876069151753290378112L, 0.541011723226630257077328e-5L, 0.535093845803642394908747e-15L }; static const float f0 = 0.5641309023F; r = tools::evaluate_polynomial(n, z - 60) / tools::evaluate_polynomial(d, z - 60); b = f0; } else { // Worst case absolute error found: 8.094633491e-22 static const T n[4] = { -0.789224703978722689089794e-5L, 0.622088451660986955124162e-6L, 0.145728445676882396797184e-7L, 0.603715505542715364529243e-10L }; static const T d[4] = { 1L, 0.0375328846356293715248719L, 0.000467919535974625308126054L, 0.193847039275845656900547e-5L }; static const float f0 = 0.5641584396F; r = tools::evaluate_polynomial(n, z - 85) / tools::evaluate_polynomial(d, z - 85); b = f0; } T g = exp(-z * z) / z; result = g * b + g * r; BOOST_MATH_INSTRUMENT_CODE("r = " << r); BOOST_MATH_INSTRUMENT_CODE("b = " << b); BOOST_MATH_INSTRUMENT_CODE("g = " << g); } else { // // Any value of z larger than 28 will underflow to zero: // result = 0; invert = !invert; } if(invert) { result = 1 - result; } return result; } // template <class T, class L>T erf_imp(T z, bool invert, const L& l, const mpl::int_<64>& t) template <class T, class Policy> T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t) { BOOST_MATH_STD_USING BOOST_MATH_INSTRUMENT_CODE("113-bit precision erf_imp called"); if(z < 0) { if(!invert) return -erf_imp(-z, invert, pol, t); else if(z < -0.5) return 2 - erf_imp(-z, invert, pol, t); else return 1 + erf_imp(-z, false, pol, t); } T result; // // Big bunch of selection statements now to pick which // implementation to use, try to put most likely options // first: // if(z < 0.5) { // // We're going to calculate erf: // if(z == 0) { result = 0; } else if(z < 1e-20) { result = z * 1.125 + z * 0.003379167095512573896158903121545171688L; } else { // Worst case absolute error found: 1.928180863e-35 static const T n[13] = { 0.0033791670955125738961589031215451706772L, -0.000356604747854533671135323429762519216044L, -0.374476838669183581687167228866769133591L, 0.0395338132469809122364498388174446488042L, -0.070405473508786506375820161461872523315L, 0.00575264725772369303419496752516485264994L, -0.0122324470706306942925087773122510971344L, 0.000982833333252586078523570049842642796291L, -0.000937806155615159592441487275938040285833L, 0.485407838108763091860415874932955355755e-4L, -0.50171236926234625577876479444632561922e-4L, 0.19406068817888598455243350289053451571e-5L, -0.119351103792049576459000102632508734863e-5L }; static const T d[13] = { 1L, -0.105530368216503232473476334993759958083L, 0.488152943026846232046726653294817930988L, -0.0470361716364117780901924633553851211874L, 0.107663671943702835026199580597519084906L, -0.00919493879447389180633447493128337242362L, 0.0138231121717229362691899919242806829805L, -0.000994048559663865788847688218108232247441L, 0.00109769834527023265969224251892094019735L, -0.600458401801636062015615549258555311545e-4L, 0.51530723974502946291624848874654212384e-4L, -0.164121264470361558910636548509486296153e-5L, 0.112643498977070218963888579607359294396e-5L }; result = z * 1.125 + z * tools::evaluate_rational(n, d, z); } } else if((z < 9) || ((z < 110) && invert)) // TODO FIXME!! { // // We'll be calculating erfc: // invert = !invert; T r, b; if(z < 0.75) { // Worst case absolute error found: 9.46579566e-36 static const T n[10] = { -0.0361790390718262471349157886581290316118L, 0.268773785250238404882137450640472787307L, 0.46350995084084251624649426251701042395L, 0.368375435727102373204587584306335625665L, 0.177618123820303858190236222513516291818L, 0.0566304173556669007529719743050764079095L, 0.0121631149481817424284077180037019529004L, 0.00171397353209314111395429418066990259845L, 0.000144662387395699594624184141956722488753L, 0.559870522050008635715382724858714587198e-5L }; static const T d[10] = { 1L, 2.50344259590701770420935329380375393716L, 2.84905597172139276093882199286535521011L, 1.93691730181297099541395314232750876411L, 0.868059574796050528229446630538462280596L, 0.266360035323208212078527036132085926692L, 0.0560555526482963925944703505114360693216L, 0.0078174400311465420803366235814673576269L, 0.000657067309046405057499687417839930873806L, 0.254293850077789079098316521097979388983e-4L }; static const float f0 = 0.3440242112F; r = tools::evaluate_rational(n, d, z - 0.5); b = f0; } else if(z < 1.25) { // Worst case absolute error found: 1.222145602e-35 static const T n[10] = { -0.03978768926111368569548863384587917014L, 0.136218360681765349252731304877153919181L, 0.252782160406474440925641829129129001834L, 0.198264231106182362320012632943145619752L, 0.0923045825293507328801206570363391760624L, 0.0281157216148097885766639832985410722743L, 0.00573041663561645197870019701493117161792L, 0.000762341440133027349203518836487137709687L, 0.60471020134417423449877859375492618899e-4L, 0.219005333943510376644902615714724932217e-5L }; static const T d[11] = { 1L, 2.38113277319993574121349184069891082204L, 2.57380422881476860215664207822277590181L, 1.65937045609044738941173490190122101824L, 0.704055811320312044285417250966993014161L, 0.20414913933328592198279939394283925451L, 0.0405162285360227740710964820549709038107L, 0.00531638867177288975915820230980317499728L, 0.000419364368135139398723983192742319455284L, 0.151874665979234971229096136924566078234e-4L, 0.807869459506748684117962248796937508011e-11L }; static const float f0 = 0.419990927F; r = tools::evaluate_polynomial(n, z - 0.75) / tools::evaluate_polynomial(d, z - 0.75); b = f0; } else if(z < 2) { // Worst case absolute error found: 5.893842955e-36 static const T n[11] = { -0.0255063683486569102096736247449691465143L, 0.045782379672906795594927072060091308408L, 0.113248439610400562258072020811195716817L, 0.0996016254422112410086711272219455446695L, 0.0508749250027894453228337309651895478017L, 0.0171081937013828309576540212196644542209L, 0.00395354196550210630440706013523069756354L, 0.000629022203390154585475081628606234279007L, 0.664903286194855400689101617763591823345e-4L, 0.423935693893425355108071655059640137208e-5L, 0.124304036910852727351487636048151737214e-6L }; static const T d[11] = { 1L, 2.39207679390801118396945702674440915308L, 2.62237869972649377524874287442154430843L, 1.73645189911172997548091140085423843725L, 0.769812706091926741262865732006953282036L, 0.238986814887891859065369830215615790694L, 0.0526759147274379214313767032352419949829L, 0.00814993801398361741777997755108018659382L, 0.00084829993036396244429607826663252633817L, 0.537276435448416921594616417908288527881e-4L, 0.157537193656690184073389824392669625417e-5L }; static const float f0 = 0.4852850139F; r = tools::evaluate_rational(n, d, z - 1.25); b = f0; } else if(z < 2.75) { // Worst case absolute error found: 4.024770853e-36 static const T n[10] = { -0.0108897177067473013323228381829144739013L, 0.0202210475357865979950082670101965963435L, 0.0403242149078991892599316678797913295452L, 0.0288492313188655084113941326565482276741L, 0.0116982647742533555237861890442866083526L, 0.00301908913020213436098518520436147263177L, 0.000511140165864993121203730804407968689429L, 0.555507897975436549741754647662158917103e-4L, 0.354571088276496768574495922471690102061e-5L, 0.101789333060641482357520298518780163915e-6L }; static const T d[11] = { 1L, 1.98184247277299581801610266575921553552L, 1.77518826058984218945272444617044495028L, 0.943934346956188464279312722940302202684L, 0.328630002685235061519039528479761588138L, 0.0777535542351039388345270792222146705189L, 0.0125143974181120800829065248546370953609L, 0.00132270605931460450441108147393979771563L, 0.834118048375722123506409257130329786209e-4L, 0.239456257167492104073765911366304033453e-5L, 0.197067742893785800814802969598122120027e-13L }; static const float f0 = 0.5216810703F; r = tools::evaluate_polynomial(n, z - 2) / tools::evaluate_polynomial(d, z - 2); b = f0; } else if(z < 3.75) { // Worst case absolute error found: 2.119317982e-36 static const T n[10] = { -0.00669534469911386821762042893742722704704L, 0.00779239067529714323524154862288379264056L, 0.0167670669587476509267036865033136655094L, 0.0113887916348251443051357686146040093464L, 0.00426976750247047946700539147728477144579L, 0.00100469100574064832665606356894416652764L, 0.000153533145320881108157829902752192859922L, 0.149337551064665413462766906201269176262e-4L, 0.846377837919513024118176704010972579138e-6L, 0.214045773545256889299689737489755489478e-7L }; static const T d[11] = { 1L, 1.78724215851193637544287795626580411105L, 1.44052576962222794702178612920219772782L, 0.687639905608366315245841860669607532265L, 0.214374367225822611754443822822738563207L, 0.0452948320968245754796139856381504201504L, 0.00649108394178118005887118777181540680812L, 0.000608904720665003139414993591868256489088L, 0.33959064390570911588709483563995284603e-4L, 0.858811916085393051026834431509997486704e-6L, 0.618878592093890510233502654703683447468e-15L }; static const float f0 = 0.5392661095F; r = tools::evaluate_polynomial(n, z - 2.75) / tools::evaluate_polynomial(d, z - 2.75); b = f0; } else if(z < 5) { // Worst case absolute error found: 3.131502824e-36 static const T n[10] = { -0.00378088626017041998499190989910098718437L, 0.0029008905292996011997575492874095588328L, 0.00662431938391549599757244232386689480515L, 0.00417809740399798845564363621020984935218L, 0.00142019749135652688012034919213168974543L, 0.000299107980170253223293530710056814995102L, 0.405161538841561396150507786831930770839e-4L, 0.346344371670880861875666253626975779945e-5L, 0.171091054330494778613793054233437928605e-6L, 0.373924433717749484258186454458704819755e-8L }; static const T d[10] = { 1L, 1.5810750672399887547849540367499389454L, 1.12479852885403050655678225945856872694L, 0.47277272679268851560291322980574597267L, 0.129444913616967588584693095240544707208L, 0.0239544490709674941887988416318107990646L, 0.00299775294496053944060700963645084591246L, 0.00024478412843088575835960648397300177201L, 0.118424712755145205302405346348931402917e-4L, 0.2588206250858483868392167535345681119e-6L }; static const float f0 = 0.549742341F; r = tools::evaluate_rational(n, d, z - 3.75); b = f0; } else if(z < 6.5) { // Worst case absolute error found: 3.352877573e-35 static const T n[9] = { -0.00210683958249012010206456754123471415706L, 0.00146329021314062287670019911742786780092L, 0.00242029064025351202243048169807220157512L, 0.0011321990764681390160708960047630195582L, 0.000277123780759366982673218537550876769487L, 0.401236501288775561636453586216146028714e-4L, 0.347799938323835778216424009916867086167e-5L, 0.167678812729975154456097184107934455429e-6L, 0.346722083660429051057284107535869165775e-8L }; static const T d[10] = { 1L, 1.22334833521124956366395053444841951468L, 0.661433457507589455018784737495591428263L, 0.206503622658778732280997770028712044451L, 0.0407323027711252752353388616742333806362L, 0.00519969889874877079704615143005539754407L, 0.000419679230772141031030427156828631265963L, 0.195896640254171597965013007459411704085e-4L, 0.405070207572551760879797790899826058473e-6L, -0.949400883467250846930389103621356900319e-17L }; static const float f0 = 0.5556300282F; r = tools::evaluate_polynomial(n, z - 5) / tools::evaluate_polynomial(d, z - 5); b = f0; } else if(z < 8) { // Worst case absolute error found: 2.10254551e-36 static const T n[9] = { -0.00107224589975857477185569028693588970638L, 0.00081159959093417892787006088639848404179L, 0.00105587689576932891666032146026668833287L, 0.000416243954540394829165805666548948771809L, 0.861189599093384016322579078144012057531e-4L, 0.105064862265647286966467927732505059558e-4L, 0.763674245263385902692134637353517251296e-6L, 0.306964079269190247526141442183490066292e-7L, 0.525762928524110511201313708396204710874e-9L }; static const T d[9] = { 1L, 1.03391233015873996503551085347368889767L, 0.471295765635712684194436077437130977978L, 0.123736066749315618886080242926593910851L, 0.0204690897886919138685460664198600282119L, 0.00218521816743913946855947853274936296576L, 0.000147057386621380823003258590658747813774L, 0.570514093706434168568509838021466564264e-5L, 0.977166974697066620826028345712327325748e-7L }; static const float f0 = 0.5588091016F; r = tools::evaluate_rational(n, d, z - 6.5); b = f0; } else if(z < 10) { // Worst case absolute error found: 8.006848023e-37 static const T n[9] = { -0.000764310289345400483607004463638641680792L, 0.000375959928342415987920641866513058424701L, 0.000477297615927227258961911005347799033473L, 0.000165573027183931250708334564452626717999L, 0.296617493157889183515359094859055088531e-4L, 0.310848843632513098624143817615592253324e-5L, 0.192932185180269434271810693046412848027e-6L, 0.658630702075882625552897504189436319318e-8L, 0.952880249971934343233104137698620618851e-10L }; static const T d[9] = { 1L, 0.885953297549629585611202995885653291163L, 0.345435500008639080844920938390800739845L, 0.0774289910213823638414558561872084410517L, 0.0109142290775434200609901181618181718948L, 0.000990815749814504617347317658063197107511L, 0.565801964505889288566864514277149126893e-4L, 0.185846832474795168475989869562411740416e-5L, 0.268875677412705938842038805073576012915e-7L }; static const float f0 = 0.5606456399F; r = tools::evaluate_rational(n, d, z - 8); b = f0; } else if(z < 12.5) { // Worst case absolute error found: 1.920487881e-37 static const T n[9] = { -0.00049563543942917072170112779531688430711L, 0.000181627479476470045833245371263435213396L, 0.000205326116055098869964168775605689851547L, 0.602131211119027087508128340443602490578e-4L, 0.904046610725767796892834734953040318357e-5L, 0.790026495730279360618915285828083295597e-6L, 0.407154634490633661408148126521656550974e-7L, 0.114947371191075676274146563385045083492e-8L, 0.136959786076422777905476122283384578138e-10L }; static const T d[9] = { 1L, 0.738936412939629363226408445024616124644L, 0.239911614287295587917806937612822134282L, 0.0447038494121568822243673246874110377585L, 0.00522920850086320874490830611785659031521L, 0.000393238307986133717620560371515858231357L, 0.18566831172454022627187937328433090172e-4L, 0.503267648562793696253090536560738864711e-6L, 0.599643373938283798258195761814169705225e-8L }; static const float f0 = 0.5619055629F; r = tools::evaluate_rational(n, d, z - 10); b = f0; } else if(z < 15.5) { // Worst case absolute error found: 1.257535398e-36 static const T n[8] = { -0.000310716603972278501158850534578560332617L, 0.00011678083970731888844953949114132723885L, 0.800190190430061077253903477346761222901e-4L, 0.156297703728913451581100601056534652076e-4L, 0.151432948873760306237990776980248434808e-5L, 0.798843820137664551611401320879346402013e-7L, 0.219398460602279142004550137383605410162e-8L, 0.244774638611231694720102050826438123042e-10L }; static const T d[9] = { 1L, 0.536396467095662522242258924420035847546L, 0.12368534422025248132153213926057643819L, 0.0158935327926481775020421628983323726346L, 0.00122923842249710594941031294559763761829L, 0.572250268558063378795115723535491980973e-4L, 0.148480028895780193980507551658637328542e-5L, 0.165653602646537759781637790321962585489e-7L, 0.159564067829807076335030582566349996674e-21L }; static const float f0 = 0.5627119541F; r = tools::evaluate_polynomial(n, z - 12.5) / tools::evaluate_polynomial(d, z - 12.5); b = f0; } else if(z < 20) { // Worst case absolute error found: 8.757869781e-37 static const T n[8] = { -0.000232165799411403604511899928682939571679L, 0.473280617901953093430938763139716601257e-4L, 0.322197983787538821545512424866901113822e-4L, 0.53341606003892328294004958671866936958e-5L, 0.429392533762301420884052643595872389588e-6L, 0.186707830002841677949013638707964861935e-7L, 0.420531779988891521855765212230340499168e-9L, 0.38311153882270641561622347190918577449e-11L }; static const T d[9] = { 1L, 0.440698415779467873589164469370853517393L, 0.0834079018864546179293964148377602724235L, 0.00878845776227345123101968908701592510214L, 0.000556792236817609981676018545344816931862L, 0.212109697877283363015270621439889468055e-4L, 0.449886743173619367191275104721917344569e-6L, 0.409854405005280582884236774230760402868e-8L, 0.185071069100878210011727114952196630136e-23L }; static const float f0 = 0.5632548332F; r = tools::evaluate_polynomial(n, z - 15.5) / tools::evaluate_polynomial(d, z - 15.5); b = f0; } else if(z < 26) { // Worst case absolute error found: 6.571456853e-37 static const T n[8] = { -0.000143129243915019877016909310584833592972L, 0.203555078616578752774553439209122798188e-4L, 0.116664173863253360297276766667046754002e-4L, 0.153607643549360281089355497044679566329e-5L, 0.976486412462035616905934994399021437658e-7L, 0.33416295362298759817564985788402188186e-8L, 0.590812660432887787644458409396105030781e-10L, 0.421471133628743008309458424282421874073e-12L }; static const T d[8] = { 1L, 0.346848503043261151874606241871432055165L, 0.0516236057301875770334953837585483539519L, 0.0042740199483843978391191804633398315544L, 0.000212586135154357046245694533825075007902L, 0.635258188334431428038271842518161892901e-5L, 0.105600415847309067601860633252823051505e-6L, 0.753327238218310630683194279382721367469e-9L }; static const float f0 = 0.5636301041F; r = tools::evaluate_polynomial(n, z - 20) / tools::evaluate_polynomial(d, z - 20); b = f0; } else if(z < 34) { // Worst case absolute error found: 2.529847851e-38 static const T n[8] = { -0.863162280463127451272161303767299107489e-4L, 0.876391266908317792353253474461536821127e-5L, 0.407624004719762912183133597708988715137e-5L, 0.418371930808379615690956857824418800194e-6L, 0.206376792034344913360458422974245946873e-7L, 0.546878311460852031076829190722479684e-9L, 0.74760305098509923829341087347518769626e-11L, 0.411832046806658129073165530819472782663e-13L }; static const T d[9] = { 1L, 0.268714425336129161136892060316889824437L, 0.0309686025544536788982104017649851516554L, 0.00198426609900573235086828549632290797514L, 0.763402107420631006526499294753645914808e-4L, 0.176354119917411263184292422389304506735e-5L, 0.226504369338582253171306523992412653547e-7L, 0.124774448034213281307712525982862926228e-9L, 0.541581693417048102342931921476367282462e-28L }; static const float f0 = 0.5638595223F; r = tools::evaluate_polynomial(n, z - 26) / tools::evaluate_polynomial(d, z - 26); b = f0; } else if(z < 46) { // Worst case absolute error found: 4.880939969e-36 static const T n[7] = { -0.552701065525823960321509114250962372959e-4L, 0.459798238451342341380837226471283129117e-5L, 0.117462487430397988459985818542278619172e-5L, 0.704965053290383355071079647747711714696e-7L, 0.191250482739845326510859812663255140286e-8L, 0.245913419605775857221974833985059356312e-10L, 0.120466123381598216554945834019040289508e-12L }; static const T d[8] = { 1L, 0.175852634004021068883919193145643406305L, 0.0128923775281953424344037205817061020944L, 0.000504385604406829419856756492198778141939L, 0.111061123996047697713589874603922096536e-4L, 0.130499191758882778721577274330215931713e-6L, 0.639279131688964342759306361922751772829e-9L, -0.10503012469906917938721140272475203112e-26L }; static const float f0 = 0.564001143F; r = tools::evaluate_polynomial(n, z - 34) / tools::evaluate_polynomial(d, z - 34); b = f0; } else if(z < 62) { // Worst case absolute error found: 8.634961697e-37 static const T n[7] = { -0.299551937061912926289705712581858190494e-4L, 0.188783643501597286680709990049243153264e-5L, 0.353403900815908094401075506391935032602e-6L, 0.156779149815977342177830075875441013335e-7L, 0.31456307256618424444443489905810774564e-9L, 0.29912544832845265905351204765862702307e-11L, 0.108360038290329929702659864116372774364e-13L }; static const T d[7] = { 1L, 0.13020345740128026085404079010013007005L, 0.00706598549835633149505923515407700416484L, 0.000204578844720147762058725212299415091588L, 0.333280847765611305843836201217690887394e-5L, 0.289666303356425428524772241712503072453e-7L, 0.104933404691983708511908027657434756019e-9L }; static const float f0 = 0.564086318F; r = tools::evaluate_polynomial(n, z - 46) / tools::evaluate_polynomial(d, z - 46); b = f0; } else if(z < 80) { // Worst case absolute error found: 6.10700166e-39 static const T n[7] = { -0.146162619512779884168960178459825270831e-4L, 0.952303834226435420147786300516273108344e-6L, 0.114559890033396819882468040960469980168e-6L, 0.368994353517438072494331750992706039868e-8L, 0.545158660461829568567388818070830588533e-10L, 0.383569328729331398089160922704933035281e-12L, 0.103104113324271568678309317039606225294e-14L }; static const T d[7] = { 1L, 0.0966822058944670599111360985490948186413L, 0.00389546596914826027568119425228340291579L, 0.837237328321088890541798035513749762375e-4L, 0.101236677684940943809478588316884539423e-5L, 0.652985528810044575089151077574382356898e-8L, 0.175523663960756825512727785416137345625e-10L }; static const float f0 = 0.5641308427F; r = tools::evaluate_polynomial(n, z - 62) / tools::evaluate_polynomial(d, z - 62); b = f0; } else { // Worst case absolute error found: 3.409761622e-39 static const T n[7] = { -0.103600733768855845278685109220598569282e-4L, 0.324846376725138276091774803419773168357e-6L, 0.376580042454826796817160322889753111347e-7L, 0.971125540805193703472871419997820785081e-9L, 0.112499588573918670534994853431614458564e-10L, 0.6161310085325929027114924903443564594e-13L, 0.128358125353302055468778305481781957985e-15L }; static const T d[7] = { 1L, 0.0749579802028981336679035635535776767532L, 0.00234137768051079846068630120744888560113L, 0.390095348825177571970348898222511084593e-4L, 0.365628111271063883320224395719085602867e-6L, 0.182790705425846241876560215158605843026e-8L, 0.380806548535012144489621218246876843627e-11L }; static const float f0 = 0.5641558766F; r = tools::evaluate_polynomial(n, z - 80) / tools::evaluate_polynomial(d, z - 80); b = f0; } T g = exp(-z * z) / z; result = g * b + g * r; } else { // // Any value of z larger than 110 will underflow to zero: // result = 0; invert = !invert; } if(invert) { result = 1 - result; } return result; } // template <class T, class L>T erf_imp(T z, bool invert, const L& l, const mpl::int_<113>& t) } // namespace detail template <class T, class Policy> inline typename tools::promote_args<T>::type erf(T z, const Policy& /* pol */) { typedef typename tools::promote_args<T>::type result_type; typedef typename policies::evaluation<result_type, Policy>::type value_type; typedef typename policies::precision<result_type, Policy>::type precision_type; typedef typename policies::normalise< Policy, policies::promote_float<false>, policies::promote_double<false>, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; BOOST_MATH_INSTRUMENT_CODE("result_type = " << typeid(result_type).name()); BOOST_MATH_INSTRUMENT_CODE("value_type = " << typeid(value_type).name()); BOOST_MATH_INSTRUMENT_CODE("precision_type = " << typeid(precision_type).name()); typedef typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<0> >, mpl::int_<0>, typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<53> >, mpl::int_<53>, // double typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<64> >, mpl::int_<64>, // 80-bit long double typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<113> >, mpl::int_<113>, // 128-bit long double mpl::int_<0> // too many bits, use generic version. >::type >::type >::type >::type tag_type; BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name()); return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::erf_imp( static_cast<value_type>(z), false, forwarding_policy(), tag_type()), "boost::math::erf<%1%>(%1%, %1%)"); } template <class T, class Policy> inline typename tools::promote_args<T>::type erfc(T z, const Policy& /* pol */) { typedef typename tools::promote_args<T>::type result_type; typedef typename policies::evaluation<result_type, Policy>::type value_type; typedef typename policies::precision<result_type, Policy>::type precision_type; typedef typename policies::normalise< Policy, policies::promote_float<false>, policies::promote_double<false>, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; BOOST_MATH_INSTRUMENT_CODE("result_type = " << typeid(result_type).name()); BOOST_MATH_INSTRUMENT_CODE("value_type = " << typeid(value_type).name()); BOOST_MATH_INSTRUMENT_CODE("precision_type = " << typeid(precision_type).name()); typedef typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<0> >, mpl::int_<0>, typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<53> >, mpl::int_<53>, // double typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<64> >, mpl::int_<64>, // 80-bit long double typename mpl::if_< mpl::less_equal<precision_type, mpl::int_<113> >, mpl::int_<113>, // 128-bit long double mpl::int_<0> // too many bits, use generic version. >::type >::type >::type >::type tag_type; BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name()); return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::erf_imp( static_cast<value_type>(z), true, forwarding_policy(), tag_type()), "boost::math::erfc<%1%>(%1%, %1%)"); } template <class T> inline typename tools::promote_args<T>::type erf(T z) { return boost::math::erf(z, policies::policy<>()); } template <class T> inline typename tools::promote_args<T>::type erfc(T z) { return boost::math::erfc(z, policies::policy<>()); } } // namespace math } // namespace boost #include <boost/math/special_functions/detail/erf_inv.hpp> #endif // BOOST_MATH_SPECIAL_ERF_HPP