ceph/src/boost/boost/math/special_functions/erf.hpp

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