]> git.donarmstrong.com Git - rsem.git/blobdiff - boost/math/special_functions/erf.hpp
Updated boost to v1.55.0
[rsem.git] / boost / math / special_functions / erf.hpp
index 5c23f49784d149b31ed637c8d9d6ad0e9ef42d9a..04954078a76118993a49b2256ef2d3ffa8e6999f 100644 (file)
@@ -15,6 +15,7 @@
 #include <boost/math/special_functions/gamma.hpp>
 #include <boost/math/tools/roots.hpp>
 #include <boost/math/policies/error_handling.hpp>
+#include <boost/math/tools/big_constant.hpp>
 
 namespace boost{ namespace math{
 
@@ -134,7 +135,7 @@ T erf_imp(T z, bool invert, const Policy& pol, const Tag& t)
       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::get_epsilon<T, Policy>(), max_iter, 1);
-      policies::check_series_iterations("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol);
+      policies::check_series_iterations<T>("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol);
    }
    else
    {
@@ -178,11 +179,11 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<53>& t)
    if(z < 0)
    {
       if(!invert)
-         return -erf_imp(-z, invert, pol, t);
+         return -erf_imp(T(-z), invert, pol, t);
       else if(z < -0.5)
-         return 2 - erf_imp(-z, invert, pol, t);
+         return 2 - erf_imp(T(-z), invert, pol, t);
       else
-         return 1 + erf_imp(-z, false, pol, t);
+         return 1 + erf_imp(T(-z), false, pol, t);
    }
 
    T result;
@@ -205,7 +206,8 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<53>& t)
          }
          else
          {
-            result = static_cast<T>(z * 1.125f + z * 0.003379167095512573896158903121545171688L);
+            static const T c = BOOST_MATH_BIG_CONSTANT(T, 53, 0.003379167095512573896158903121545171688);
+            result = static_cast<T>(z * 1.125f + z * c);
          }
       }
       else
@@ -217,18 +219,18 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<53>& t)
 
          static const T Y = 1.044948577880859375f;
          static const T P[] = {    
-            0.0834305892146531832907L,
-            -0.338165134459360935041L,
-            -0.0509990735146777432841L,
-            -0.00772758345802133288487L,
-            -0.000322780120964605683831L,
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0834305892146531832907),
+            BOOST_MATH_BIG_CONSTANT(T, 53, -0.338165134459360935041),
+            BOOST_MATH_BIG_CONSTANT(T, 53, -0.0509990735146777432841),
+            BOOST_MATH_BIG_CONSTANT(T, 53, -0.00772758345802133288487),
+            BOOST_MATH_BIG_CONSTANT(T, 53, -0.000322780120964605683831),
          };
          static const T Q[] = {    
-            1L,
-            0.455004033050794024546L,
-            0.0875222600142252549554L,
-            0.00858571925074406212772L,
-            0.000370900071787748000569L,
+            BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.455004033050794024546),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0875222600142252549554),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.00858571925074406212772),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.000370900071787748000569),
          };
          T zz = z * z;
          result = z * (Y + tools::evaluate_polynomial(P, zz) / tools::evaluate_polynomial(Q, zz));
@@ -248,24 +250,30 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<53>& t)
          // Max Error found at double precision =        4.841816e-17
          static const T Y = 0.405935764312744140625f;
          static const T P[] = {    
-            -0.098090592216281240205L,
-            0.178114665841120341155L,
-            0.191003695796775433986L,
-            0.0888900368967884466578L,
-            0.0195049001251218801359L,
-            0.00180424538297014223957L,
+            BOOST_MATH_BIG_CONSTANT(T, 53, -0.098090592216281240205),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.178114665841120341155),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.191003695796775433986),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0888900368967884466578),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0195049001251218801359),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.00180424538297014223957),
          };
          static const T Q[] = {    
-            1L,
-            1.84759070983002217845L,
-            1.42628004845511324508L,
-            0.578052804889902404909L,
-            0.12385097467900864233L,
-            0.0113385233577001411017L,
-            0.337511472483094676155e-5L,
+            BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 1.84759070983002217845),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 1.42628004845511324508),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.578052804889902404909),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.12385097467900864233),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0113385233577001411017),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.337511472483094676155e-5),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 0.5) / tools::evaluate_polynomial(Q, z - 0.5);
+         BOOST_MATH_INSTRUMENT_VARIABLE(Y);
+         BOOST_MATH_INSTRUMENT_VARIABLE(P[0]);
+         BOOST_MATH_INSTRUMENT_VARIABLE(Q[0]);
+         BOOST_MATH_INSTRUMENT_VARIABLE(z);
+         result = Y + tools::evaluate_polynomial(P, T(z - 0.5)) / tools::evaluate_polynomial(Q, T(z - 0.5));
+         BOOST_MATH_INSTRUMENT_VARIABLE(result);
          result *= exp(-z * z) / z;
+         BOOST_MATH_INSTRUMENT_VARIABLE(result);
       }
       else if(z < 2.5f)
       {
@@ -275,22 +283,22 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<53>& t)
          // Maximum Relative Change in Control Points:   9.886e-05
          static const T Y = 0.50672817230224609375f;
          static const T P[] = {    
-            -0.0243500476207698441272L,
-            0.0386540375035707201728L,
-            0.04394818964209516296L,
-            0.0175679436311802092299L,
-            0.00323962406290842133584L,
-            0.000235839115596880717416L,
+            BOOST_MATH_BIG_CONSTANT(T, 53, -0.0243500476207698441272),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0386540375035707201728),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.04394818964209516296),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0175679436311802092299),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.00323962406290842133584),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.000235839115596880717416),
          };
          static const T Q[] = {    
-            1L,
-            1.53991494948552447182L,
-            0.982403709157920235114L,
-            0.325732924782444448493L,
-            0.0563921837420478160373L,
-            0.00410369723978904575884L,
+            BOOST_MATH_BIG_CONSTANT(T, 53, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 1.53991494948552447182),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.982403709157920235114),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.325732924782444448493),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0563921837420478160373),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.00410369723978904575884),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 1.5) / tools::evaluate_polynomial(Q, z - 1.5);
+         result = Y + tools::evaluate_polynomial(P, T(z - 1.5)) / tools::evaluate_polynomial(Q, T(z - 1.5));
          result *= exp(-z * z) / z;
       }
       else if(z < 4.5f)
@@ -301,22 +309,22 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<53>& t)
          // Max Error found at double precision =        2.062515e-17
          static const T Y = 0.5405750274658203125f;
          static const T P[] = {    
-            0.00295276716530971662634L,
-            0.0137384425896355332126L,
-            0.00840807615555585383007L,
-            0.00212825620914618649141L,
-            0.000250269961544794627958L,
-            0.113212406648847561139e-4L,
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.00295276716530971662634),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0137384425896355332126),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.00840807615555585383007),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.00212825620914618649141),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.000250269961544794627958),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.113212406648847561139e-4),
          };
          static const T Q[] = {    
-            1L,
-            1.04217814166938418171L,
-            0.442597659481563127003L,
-            0.0958492726301061423444L,
-            0.0105982906484876531489L,
-            0.000479411269521714493907L,
+            BOOST_MATH_BIG_CONSTANT(T, 53, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 1.04217814166938418171),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.442597659481563127003),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0958492726301061423444),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0105982906484876531489),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.000479411269521714493907),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 3.5) / tools::evaluate_polynomial(Q, z - 3.5);
+         result = Y + tools::evaluate_polynomial(P, T(z - 3.5)) / tools::evaluate_polynomial(Q, T(z - 3.5));
          result *= exp(-z * z) / z;
       }
       else
@@ -327,24 +335,24 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<53>& t)
          // Maximum Relative Change in Control Points:   1.357e-05
          static const T Y = 0.5579090118408203125f;
          static const T P[] = {    
-            0.00628057170626964891937L,
-            0.0175389834052493308818L,
-            -0.212652252872804219852L,
-            -0.687717681153649930619L,
-            -2.5518551727311523996L,
-            -3.22729451764143718517L,
-            -2.8175401114513378771L,
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.00628057170626964891937),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0175389834052493308818),
+            BOOST_MATH_BIG_CONSTANT(T, 53, -0.212652252872804219852),
+            BOOST_MATH_BIG_CONSTANT(T, 53, -0.687717681153649930619),
+            BOOST_MATH_BIG_CONSTANT(T, 53, -2.5518551727311523996),
+            BOOST_MATH_BIG_CONSTANT(T, 53, -3.22729451764143718517),
+            BOOST_MATH_BIG_CONSTANT(T, 53, -2.8175401114513378771),
          };
          static const T Q[] = {    
-            1L,
-            2.79257750980575282228L,
-            11.0567237927800161565L,
-            15.930646027911794143L,
-            22.9367376522880577224L,
-            13.5064170191802889145L,
-            5.48409182238641741584L,
+            BOOST_MATH_BIG_CONSTANT(T, 53, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 2.79257750980575282228),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 11.0567237927800161565),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 15.930646027911794143),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 22.9367376522880577224),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 13.5064170191802889145),
+            BOOST_MATH_BIG_CONSTANT(T, 53, 5.48409182238641741584),
          };
-         result = Y + tools::evaluate_polynomial(P, 1 / z) / tools::evaluate_polynomial(Q, 1 / z);
+         result = Y + tools::evaluate_polynomial(P, T(1 / z)) / tools::evaluate_polynomial(Q, T(1 / z));
          result *= exp(-z * z) / z;
       }
    }
@@ -363,7 +371,7 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<53>& t)
    }
 
    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 Lanczos>T erf_imp(T z, bool invert, const Lanczos& l, const mpl::int_<53>& t)
 
 
 template <class T, class Policy>
@@ -376,11 +384,11 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<64>& t)
    if(z < 0)
    {
       if(!invert)
-         return -erf_imp(-z, invert, pol, t);
+         return -erf_imp(T(-z), invert, pol, t);
       else if(z < -0.5)
-         return 2 - erf_imp(-z, invert, pol, t);
+         return 2 - erf_imp(T(-z), invert, pol, t);
       else
-         return 1 + erf_imp(-z, false, pol, t);
+         return 1 + erf_imp(T(-z), false, pol, t);
    }
 
    T result;
@@ -401,7 +409,8 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<64>& t)
       }
       else if(z < 1e-10)
       {
-         result = z * 1.125 + z * 0.003379167095512573896158903121545171688L;
+         static const T c = BOOST_MATH_BIG_CONSTANT(T, 64, 0.003379167095512573896158903121545171688);
+         result = z * 1.125 + z * c;
       }
       else
       {
@@ -411,22 +420,22 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<64>& t)
          // Maximum Relative Change in Control Points:   1.474e-04
          static const T Y = 1.044948577880859375f;
          static const T P[] = {    
-            0.0834305892146531988966L,
-            -0.338097283075565413695L,
-            -0.0509602734406067204596L,
-            -0.00904906346158537794396L,
-            -0.000489468651464798669181L,
-            -0.200305626366151877759e-4L,
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0834305892146531988966),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -0.338097283075565413695),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -0.0509602734406067204596),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -0.00904906346158537794396),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -0.000489468651464798669181),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -0.200305626366151877759e-4),
          };
          static const T Q[] = {    
-            1L,
-            0.455817300515875172439L,
-            0.0916537354356241792007L,
-            0.0102722652675910031202L,
-            0.000650511752687851548735L,
-            0.189532519105655496778e-4L,
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.455817300515875172439),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0916537354356241792007),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0102722652675910031202),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.000650511752687851548735),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.189532519105655496778e-4),
          };
-         result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
+         result = z * (Y + tools::evaluate_polynomial(P, T(z * z)) / tools::evaluate_polynomial(Q, T(z * z)));
       }
    }
    else if(invert ? (z < 110) : (z < 6.4f))
@@ -443,25 +452,25 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<64>& t)
          // Maximum Relative Change in Control Points:   5.110e-03
          static const T Y = 0.405935764312744140625f;
          static const T P[] = {    
-            -0.0980905922162812031672L,
-            0.159989089922969141329L,
-            0.222359821619935712378L,
-            0.127303921703577362312L,
-            0.0384057530342762400273L,
-            0.00628431160851156719325L,
-            0.000441266654514391746428L,
-            0.266689068336295642561e-7L,
+            BOOST_MATH_BIG_CONSTANT(T, 64, -0.0980905922162812031672),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.159989089922969141329),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.222359821619935712378),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.127303921703577362312),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0384057530342762400273),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00628431160851156719325),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.000441266654514391746428),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.266689068336295642561e-7),
          };
          static const T Q[] = {    
-            1L,
-            2.03237474985469469291L,
-            1.78355454954969405222L,
-            0.867940326293760578231L,
-            0.248025606990021698392L,
-            0.0396649631833002269861L,
-            0.00279220237309449026796L,
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 2.03237474985469469291),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.78355454954969405222),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.867940326293760578231),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.248025606990021698392),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0396649631833002269861),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00279220237309449026796),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 0.5f) / tools::evaluate_polynomial(Q, z - 0.5f);
+         result = Y + tools::evaluate_polynomial(P, T(z - 0.5f)) / tools::evaluate_polynomial(Q, T(z - 0.5f));
          result *= exp(-z * z) / z;
       }
       else if(z < 2.5)
@@ -472,24 +481,24 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<64>& t)
          // Maximum Relative Change in Control Points:   1.793e-04
          static const T Y = 0.50672817230224609375f;
          static const T P[] = {    
-            -0.024350047620769840217L,
-            0.0343522687935671451309L,
-            0.0505420824305544949541L,
-            0.0257479325917757388209L,
-            0.00669349844190354356118L,
-            0.00090807914416099524444L,
-            0.515917266698050027934e-4L,
+            BOOST_MATH_BIG_CONSTANT(T, 64, -0.024350047620769840217),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0343522687935671451309),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0505420824305544949541),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0257479325917757388209),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00669349844190354356118),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00090807914416099524444),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.515917266698050027934e-4),
          };
          static const T Q[] = {    
-            1L,
-            1.71657861671930336344L,
-            1.26409634824280366218L,
-            0.512371437838969015941L,
-            0.120902623051120950935L,
-            0.0158027197831887485261L,
-            0.000897871370778031611439L,
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.71657861671930336344),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.26409634824280366218),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.512371437838969015941),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.120902623051120950935),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0158027197831887485261),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.000897871370778031611439),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 1.5f) / tools::evaluate_polynomial(Q, z - 1.5f);
+         result = Y + tools::evaluate_polynomial(P, T(z - 1.5f)) / tools::evaluate_polynomial(Q, T(z - 1.5f));
          result *= exp(-z * z) / z;
       }
       else if(z < 4.5)
@@ -500,24 +509,24 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<64>& t)
          // Max Error found at long double precision =   1.446908e-20
          static const T Y  = 0.5405750274658203125f;
          static const T P[] = {    
-            0.0029527671653097284033L,
-            0.0141853245895495604051L,
-            0.0104959584626432293901L,
-            0.00343963795976100077626L,
-            0.00059065441194877637899L,
-            0.523435380636174008685e-4L,
-            0.189896043050331257262e-5L,
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0029527671653097284033),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0141853245895495604051),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0104959584626432293901),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00343963795976100077626),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00059065441194877637899),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.523435380636174008685e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.189896043050331257262e-5),
          };
          static const T Q[] = {    
-            1L,
-            1.19352160185285642574L,
-            0.603256964363454392857L,
-            0.165411142458540585835L,
-            0.0259729870946203166468L,
-            0.00221657568292893699158L,
-            0.804149464190309799804e-4L,
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1.19352160185285642574),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.603256964363454392857),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.165411142458540585835),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0259729870946203166468),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00221657568292893699158),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.804149464190309799804e-4),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 3.5f) / tools::evaluate_polynomial(Q, z - 3.5f);
+         result = Y + tools::evaluate_polynomial(P, T(z - 3.5f)) / tools::evaluate_polynomial(Q, T(z - 3.5f));
          result *= exp(-z * z) / z;
       }
       else
@@ -528,28 +537,28 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<64>& t)
          // Maximum Relative Change in Control Points:   2.319e-05
          static const T Y = 0.55825519561767578125f;
          static const T P[] = {    
-            0.00593438793008050214106L,
-            0.0280666231009089713937L,
-            -0.141597835204583050043L,
-            -0.978088201154300548842L,
-            -5.47351527796012049443L,
-            -13.8677304660245326627L,
-            -27.1274948720539821722L,
-            -29.2545152747009461519L,
-            -16.8865774499799676937L,
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.00593438793008050214106),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 0.0280666231009089713937),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -0.141597835204583050043),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -0.978088201154300548842),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -5.47351527796012049443),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -13.8677304660245326627),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -27.1274948720539821722),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -29.2545152747009461519),
+            BOOST_MATH_BIG_CONSTANT(T, 64, -16.8865774499799676937),
          };
          static const T Q[] = {    
-            1L,
-            4.72948911186645394541L,
-            23.6750543147695749212L,
-            60.0021517335693186785L,
-            131.766251645149522868L,
-            178.167924971283482513L,
-            182.499390505915222699L,
-            104.365251479578577989L,
-            30.8365511891224291717L,
+            BOOST_MATH_BIG_CONSTANT(T, 64, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 4.72948911186645394541),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 23.6750543147695749212),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 60.0021517335693186785),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 131.766251645149522868),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 178.167924971283482513),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 182.499390505915222699),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 104.365251479578577989),
+            BOOST_MATH_BIG_CONSTANT(T, 64, 30.8365511891224291717),
          };
-         result = Y + tools::evaluate_polynomial(P, 1 / z) / tools::evaluate_polynomial(Q, 1 / z);
+         result = Y + tools::evaluate_polynomial(P, T(1 / z)) / tools::evaluate_polynomial(Q, T(1 / z));
          result *= exp(-z * z) / z;
       }
    }
@@ -568,7 +577,7 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<64>& t)
    }
 
    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 Lanczos>T erf_imp(T z, bool invert, const Lanczos& l, const mpl::int_<64>& t)
 
 
 template <class T, class Policy>
@@ -581,11 +590,11 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
    if(z < 0)
    {
       if(!invert)
-         return -erf_imp(-z, invert, pol, t);
+         return -erf_imp(T(-z), invert, pol, t);
       else if(z < -0.5)
-         return 2 - erf_imp(-z, invert, pol, t);
+         return 2 - erf_imp(T(-z), invert, pol, t);
       else
-         return 1 + erf_imp(-z, false, pol, t);
+         return 1 + erf_imp(T(-z), false, pol, t);
    }
 
    T result;
@@ -606,7 +615,8 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
       }
       else if(z < 1e-20)
       {
-         result = z * 1.125 + z * 0.003379167095512573896158903121545171688L;
+         static const T c = BOOST_MATH_BIG_CONSTANT(T, 113, 0.003379167095512573896158903121545171688);
+         result = z * 1.125 + z * c;
       }
       else
       {
@@ -616,26 +626,26 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
          // Maximum Relative Change in Control Points:   3.492e-10
          static const T Y = 1.0841522216796875f;
          static const T P[] = {    
-            0.0442269454158250738961589031215451778L,
-            -0.35549265736002144875335323556961233L,
-            -0.0582179564566667896225454670863270393L,
-            -0.0112694696904802304229950538453123925L,
-            -0.000805730648981801146251825329609079099L,
-            -0.566304966591936566229702842075966273e-4L,
-            -0.169655010425186987820201021510002265e-5L,
-            -0.344448249920445916714548295433198544e-7L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0442269454158250738961589031215451778),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.35549265736002144875335323556961233),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.0582179564566667896225454670863270393),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.0112694696904802304229950538453123925),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.000805730648981801146251825329609079099),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.566304966591936566229702842075966273e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.169655010425186987820201021510002265e-5),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.344448249920445916714548295433198544e-7),
          };
          static const T Q[] = {    
-            1L,
-            0.466542092785657604666906909196052522L,
-            0.100005087012526447295176964142107611L,
-            0.0128341535890117646540050072234142603L,
-            0.00107150448466867929159660677016658186L,
-            0.586168368028999183607733369248338474e-4L,
-            0.196230608502104324965623171516808796e-5L,
-            0.313388521582925207734229967907890146e-7L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.466542092785657604666906909196052522),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.100005087012526447295176964142107611),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0128341535890117646540050072234142603),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00107150448466867929159660677016658186),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.586168368028999183607733369248338474e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.196230608502104324965623171516808796e-5),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.313388521582925207734229967907890146e-7),
          };
-         result = z * (Y + tools::evaluate_polynomial(P, z * z) / tools::evaluate_polynomial(Q, z * z));
+         result = z * (Y + tools::evaluate_polynomial(P, T(z * z)) / tools::evaluate_polynomial(Q, T(z * z)));
       }
    }
    else if(invert ? (z < 110) : (z < 8.65f))
@@ -652,31 +662,31 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
          // Maximum Relative Change in Control Points:   6.127e-05
          static const T Y = 0.371877193450927734375f;
          static const T P[] = {    
-            -0.0640320213544647969396032886581290455L,
-            0.200769874440155895637857443946706731L,
-            0.378447199873537170666487408805779826L,
-            0.30521399466465939450398642044975127L,
-            0.146890026406815277906781824723458196L,
-            0.0464837937749539978247589252732769567L,
-            0.00987895759019540115099100165904822903L,
-            0.00137507575429025512038051025154301132L,
-            0.0001144764551085935580772512359680516L,
-            0.436544865032836914773944382339900079e-5L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.0640320213544647969396032886581290455),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.200769874440155895637857443946706731),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.378447199873537170666487408805779826),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.30521399466465939450398642044975127),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.146890026406815277906781824723458196),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0464837937749539978247589252732769567),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00987895759019540115099100165904822903),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00137507575429025512038051025154301132),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0001144764551085935580772512359680516),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.436544865032836914773944382339900079e-5),
          };
          static const T Q[] = {    
-            1L,
-            2.47651182872457465043733800302427977L,
-            2.78706486002517996428836400245547955L,
-            1.87295924621659627926365005293130693L,
-            0.829375825174365625428280908787261065L,
-            0.251334771307848291593780143950311514L,
-            0.0522110268876176186719436765734722473L,
-            0.00718332151250963182233267040106902368L,
-            0.000595279058621482041084986219276392459L,
-            0.226988669466501655990637599399326874e-4L,
-            0.270666232259029102353426738909226413e-10L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.47651182872457465043733800302427977),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.78706486002517996428836400245547955),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.87295924621659627926365005293130693),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.829375825174365625428280908787261065),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.251334771307848291593780143950311514),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0522110268876176186719436765734722473),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00718332151250963182233267040106902368),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000595279058621482041084986219276392459),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.226988669466501655990637599399326874e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.270666232259029102353426738909226413e-10),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 0.5f) / tools::evaluate_polynomial(Q, z - 0.5f);
+         result = Y + tools::evaluate_polynomial(P, T(z - 0.5f)) / tools::evaluate_polynomial(Q, T(z - 0.5f));
          result *= exp(-z * z) / z;
       }
       else if(z < 1.5)
@@ -687,30 +697,30 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
          // Maximum Relative Change in Control Points:   6.104e-05
          static const T Y = 0.45658016204833984375f;
          static const T P[] = {    
-            -0.0289965858925328393392496555094848345L,
-            0.0868181194868601184627743162571779226L,
-            0.169373435121178901746317404936356745L,
-            0.13350446515949251201104889028133486L,
-            0.0617447837290183627136837688446313313L,
-            0.0185618495228251406703152962489700468L,
-            0.00371949406491883508764162050169531013L,
-            0.000485121708792921297742105775823900772L,
-            0.376494706741453489892108068231400061e-4L,
-            0.133166058052466262415271732172490045e-5L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.0289965858925328393392496555094848345),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0868181194868601184627743162571779226),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.169373435121178901746317404936356745),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.13350446515949251201104889028133486),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0617447837290183627136837688446313313),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0185618495228251406703152962489700468),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00371949406491883508764162050169531013),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000485121708792921297742105775823900772),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.376494706741453489892108068231400061e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.133166058052466262415271732172490045e-5),
          };
          static const T Q[] = {    
-            1L,
-            2.32970330146503867261275580968135126L,
-            2.46325715420422771961250513514928746L,
-            1.55307882560757679068505047390857842L,
-            0.644274289865972449441174485441409076L,
-            0.182609091063258208068606847453955649L,
-            0.0354171651271241474946129665801606795L,
-            0.00454060370165285246451879969534083997L,
-            0.000349871943711566546821198612518656486L,
-            0.123749319840299552925421880481085392e-4L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.32970330146503867261275580968135126),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.46325715420422771961250513514928746),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.55307882560757679068505047390857842),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.644274289865972449441174485441409076),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.182609091063258208068606847453955649),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0354171651271241474946129665801606795),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00454060370165285246451879969534083997),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000349871943711566546821198612518656486),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.123749319840299552925421880481085392e-4),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 1.0f) / tools::evaluate_polynomial(Q, z - 1.0f);
+         result = Y + tools::evaluate_polynomial(P, T(z - 1.0f)) / tools::evaluate_polynomial(Q, T(z - 1.0f));
          result *= exp(-z * z) / z;
       }
       else if(z < 2.25)
@@ -721,31 +731,31 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
          // Max Error found at long double precision =   1.998462e-35
          static const T Y = 0.50250148773193359375f;
          static const T P[] = {    
-            -0.0201233630504573402185161184151016606L,
-            0.0331864357574860196516686996302305002L,
-            0.0716562720864787193337475444413405461L,
-            0.0545835322082103985114927569724880658L,
-            0.0236692635189696678976549720784989593L,
-            0.00656970902163248872837262539337601845L,
-            0.00120282643299089441390490459256235021L,
-            0.000142123229065182650020762792081622986L,
-            0.991531438367015135346716277792989347e-5L,
-            0.312857043762117596999398067153076051e-6L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.0201233630504573402185161184151016606),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0331864357574860196516686996302305002),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0716562720864787193337475444413405461),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0545835322082103985114927569724880658),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0236692635189696678976549720784989593),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00656970902163248872837262539337601845),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00120282643299089441390490459256235021),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000142123229065182650020762792081622986),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.991531438367015135346716277792989347e-5),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.312857043762117596999398067153076051e-6),
          };
          static const T Q[] = {    
-            1L,
-            2.13506082409097783827103424943508554L,
-            2.06399257267556230937723190496806215L,
-            1.18678481279932541314830499880691109L,
-            0.447733186643051752513538142316799562L,
-            0.11505680005657879437196953047542148L,
-            0.020163993632192726170219663831914034L,
-            0.00232708971840141388847728782209730585L,
-            0.000160733201627963528519726484608224112L,
-            0.507158721790721802724402992033269266e-5L,
-            0.18647774409821470950544212696270639e-12L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.13506082409097783827103424943508554),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2.06399257267556230937723190496806215),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.18678481279932541314830499880691109),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.447733186643051752513538142316799562),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.11505680005657879437196953047542148),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.020163993632192726170219663831914034),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00232708971840141388847728782209730585),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000160733201627963528519726484608224112),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.507158721790721802724402992033269266e-5),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.18647774409821470950544212696270639e-12),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 1.5f) / tools::evaluate_polynomial(Q, z - 1.5f);
+         result = Y + tools::evaluate_polynomial(P, T(z - 1.5f)) / tools::evaluate_polynomial(Q, T(z - 1.5f));
          result *= exp(-z * z) / z;
       }
       else if (z < 3)
@@ -756,30 +766,30 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
          // Max Error found at long double precision =   5.794737e-36
          static const T Y = 0.52896785736083984375f;
          static const T P[] = {    
-            -0.00902152521745813634562524098263360074L,
-            0.0145207142776691539346923710537580927L,
-            0.0301681239582193983824211995978678571L,
-            0.0215548540823305814379020678660434461L,
-            0.00864683476267958365678294164340749949L,
-            0.00219693096885585491739823283511049902L,
-            0.000364961639163319762492184502159894371L,
-            0.388174251026723752769264051548703059e-4L,
-            0.241918026931789436000532513553594321e-5L,
-            0.676586625472423508158937481943649258e-7L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.00902152521745813634562524098263360074),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0145207142776691539346923710537580927),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0301681239582193983824211995978678571),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0215548540823305814379020678660434461),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00864683476267958365678294164340749949),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00219693096885585491739823283511049902),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000364961639163319762492184502159894371),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.388174251026723752769264051548703059e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.241918026931789436000532513553594321e-5),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.676586625472423508158937481943649258e-7),
          };
          static const T Q[] = {    
-            1L,
-            1.93669171363907292305550231764920001L,
-            1.69468476144051356810672506101377494L,
-            0.880023580986436640372794392579985511L,
-            0.299099106711315090710836273697708402L,
-            0.0690593962363545715997445583603382337L,
-            0.0108427016361318921960863149875360222L,
-            0.00111747247208044534520499324234317695L,
-            0.686843205749767250666787987163701209e-4L,
-            0.192093541425429248675532015101904262e-5L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.93669171363907292305550231764920001),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.69468476144051356810672506101377494),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.880023580986436640372794392579985511),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.299099106711315090710836273697708402),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0690593962363545715997445583603382337),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0108427016361318921960863149875360222),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00111747247208044534520499324234317695),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.686843205749767250666787987163701209e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.192093541425429248675532015101904262e-5),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 2.25f) / tools::evaluate_polynomial(Q, z - 2.25f);
+         result = Y + tools::evaluate_polynomial(P, T(z - 2.25f)) / tools::evaluate_polynomial(Q, T(z - 2.25f));
          result *= exp(-z * z) / z;
       }
       else if(z < 3.5)
@@ -790,28 +800,28 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
          // Max Error found at long double precision =   1.747062e-36
          static const T Y = 0.54037380218505859375f;
          static const T P[] = {    
-            -0.0033703486408887424921155540591370375L,
-            0.0104948043110005245215286678898115811L,
-            0.0148530118504000311502310457390417795L,
-            0.00816693029245443090102738825536188916L,
-            0.00249716579989140882491939681805594585L,
-            0.0004655591010047353023978045800916647L,
-            0.531129557920045295895085236636025323e-4L,
-            0.343526765122727069515775194111741049e-5L,
-            0.971120407556888763695313774578711839e-7L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.0033703486408887424921155540591370375),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0104948043110005245215286678898115811),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0148530118504000311502310457390417795),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00816693029245443090102738825536188916),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00249716579989140882491939681805594585),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0004655591010047353023978045800916647),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.531129557920045295895085236636025323e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.343526765122727069515775194111741049e-5),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.971120407556888763695313774578711839e-7),
          };
          static const T Q[] = {    
-            1L,
-            1.59911256167540354915906501335919317L,
-            1.136006830764025173864831382946934L,
-            0.468565867990030871678574840738423023L,
-            0.122821824954470343413956476900662236L,
-            0.0209670914950115943338996513330141633L,
-            0.00227845718243186165620199012883547257L,
-            0.000144243326443913171313947613547085553L,
-            0.407763415954267700941230249989140046e-5L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.59911256167540354915906501335919317),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.136006830764025173864831382946934),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.468565867990030871678574840738423023),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.122821824954470343413956476900662236),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0209670914950115943338996513330141633),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00227845718243186165620199012883547257),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000144243326443913171313947613547085553),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.407763415954267700941230249989140046e-5),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 3.0f) / tools::evaluate_polynomial(Q, z - 3.0f);
+         result = Y + tools::evaluate_polynomial(P, T(z - 3.0f)) / tools::evaluate_polynomial(Q, T(z - 3.0f));
          result *= exp(-z * z) / z;
       }
       else if(z < 5.5)
@@ -822,32 +832,32 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
          // Max Error found at long double precision =   1.349545e-35
          static const T Y = 0.55000019073486328125f;
          static const T P[] = {    
-            0.00118142849742309772151454518093813615L,
-            0.0072201822885703318172366893469382745L,
-            0.0078782276276860110721875733778481505L,
-            0.00418229166204362376187593976656261146L,
-            0.00134198400587769200074194304298642705L,
-            0.000283210387078004063264777611497435572L,
-            0.405687064094911866569295610914844928e-4L,
-            0.39348283801568113807887364414008292e-5L,
-            0.248798540917787001526976889284624449e-6L,
-            0.929502490223452372919607105387474751e-8L,
-            0.156161469668275442569286723236274457e-9L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00118142849742309772151454518093813615),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0072201822885703318172366893469382745),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0078782276276860110721875733778481505),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00418229166204362376187593976656261146),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00134198400587769200074194304298642705),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000283210387078004063264777611497435572),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.405687064094911866569295610914844928e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.39348283801568113807887364414008292e-5),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.248798540917787001526976889284624449e-6),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.929502490223452372919607105387474751e-8),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.156161469668275442569286723236274457e-9),
          };
          static const T Q[] = {    
-            1L,
-            1.52955245103668419479878456656709381L,
-            1.06263944820093830054635017117417064L,
-            0.441684612681607364321013134378316463L,
-            0.121665258426166960049773715928906382L,
-            0.0232134512374747691424978642874321434L,
-            0.00310778180686296328582860464875562636L,
-            0.000288361770756174705123674838640161693L,
-            0.177529187194133944622193191942300132e-4L,
-            0.655068544833064069223029299070876623e-6L,
-            0.11005507545746069573608988651927452e-7L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.52955245103668419479878456656709381),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.06263944820093830054635017117417064),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.441684612681607364321013134378316463),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.121665258426166960049773715928906382),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0232134512374747691424978642874321434),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00310778180686296328582860464875562636),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000288361770756174705123674838640161693),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.177529187194133944622193191942300132e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.655068544833064069223029299070876623e-6),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.11005507545746069573608988651927452e-7),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 4.5f) / tools::evaluate_polynomial(Q, z - 4.5f);
+         result = Y + tools::evaluate_polynomial(P, T(z - 4.5f)) / tools::evaluate_polynomial(Q, T(z - 4.5f));
          result *= exp(-z * z) / z;
       }
       else if(z < 7.5)
@@ -858,30 +868,30 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
          // Max Error found at long double precision =   2.646420e-36
          static const T Y = 0.5574436187744140625f;
          static const T P[] = {    
-            0.000293236907400849056269309713064107674L,
-            0.00225110719535060642692275221961480162L,
-            0.00190984458121502831421717207849429799L,
-            0.000747757733460111743833929141001680706L,
-            0.000170663175280949889583158597373928096L,
-            0.246441188958013822253071608197514058e-4L,
-            0.229818000860544644974205957895688106e-5L,
-            0.134886977703388748488480980637704864e-6L,
-            0.454764611880548962757125070106650958e-8L,
-            0.673002744115866600294723141176820155e-10L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000293236907400849056269309713064107674),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00225110719535060642692275221961480162),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00190984458121502831421717207849429799),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000747757733460111743833929141001680706),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000170663175280949889583158597373928096),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.246441188958013822253071608197514058e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.229818000860544644974205957895688106e-5),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.134886977703388748488480980637704864e-6),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.454764611880548962757125070106650958e-8),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.673002744115866600294723141176820155e-10),
          };
          static const T Q[] = {    
-            1L,
-            1.12843690320861239631195353379313367L,
-            0.569900657061622955362493442186537259L,
-            0.169094404206844928112348730277514273L,
-            0.0324887449084220415058158657252147063L,
-            0.00419252877436825753042680842608219552L,
-            0.00036344133176118603523976748563178578L,
-            0.204123895931375107397698245752850347e-4L,
-            0.674128352521481412232785122943508729e-6L,
-            0.997637501418963696542159244436245077e-8L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.12843690320861239631195353379313367),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.569900657061622955362493442186537259),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.169094404206844928112348730277514273),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0324887449084220415058158657252147063),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00419252877436825753042680842608219552),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00036344133176118603523976748563178578),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.204123895931375107397698245752850347e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.674128352521481412232785122943508729e-6),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.997637501418963696542159244436245077e-8),
          };
-         result = Y + tools::evaluate_polynomial(P, z - 6.5f) / tools::evaluate_polynomial(Q, z - 6.5f);
+         result = Y + tools::evaluate_polynomial(P, T(z - 6.5f)) / tools::evaluate_polynomial(Q, T(z - 6.5f));
          result *= exp(-z * z) / z;
       }
       else if(z < 11.5)
@@ -892,30 +902,30 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
          // Max Error found at long double precision =   9.849522e-36
          static const T Y = 0.56083202362060546875f;
          static const T P[] = {    
-            0.000282420728751494363613829834891390121L,
-            0.00175387065018002823433704079355125161L,
-            0.0021344978564889819420775336322920375L,
-            0.00124151356560137532655039683963075661L,
-            0.000423600733566948018555157026862139644L,
-            0.914030340865175237133613697319509698e-4L,
-            0.126999927156823363353809747017945494e-4L,
-            0.110610959842869849776179749369376402e-5L,
-            0.55075079477173482096725348704634529e-7L,
-            0.119735694018906705225870691331543806e-8L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000282420728751494363613829834891390121),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00175387065018002823433704079355125161),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0021344978564889819420775336322920375),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00124151356560137532655039683963075661),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000423600733566948018555157026862139644),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.914030340865175237133613697319509698e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.126999927156823363353809747017945494e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.110610959842869849776179749369376402e-5),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.55075079477173482096725348704634529e-7),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.119735694018906705225870691331543806e-8),
          };
          static const T Q[] = {    
-            1L,
-            1.69889613396167354566098060039549882L,
-            1.28824647372749624464956031163282674L,
-            0.572297795434934493541628008224078717L,
-            0.164157697425571712377043857240773164L,
-            0.0315311145224594430281219516531649562L,
-            0.00405588922155632380812945849777127458L,
-            0.000336929033691445666232029762868642417L,
-            0.164033049810404773469413526427932109e-4L,
-            0.356615210500531410114914617294694857e-6L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.69889613396167354566098060039549882),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1.28824647372749624464956031163282674),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.572297795434934493541628008224078717),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.164157697425571712377043857240773164),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.0315311145224594430281219516531649562),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00405588922155632380812945849777127458),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000336929033691445666232029762868642417),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.164033049810404773469413526427932109e-4),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.356615210500531410114914617294694857e-6),
          };
-         result = Y + tools::evaluate_polynomial(P, z / 2 - 4.75f) / tools::evaluate_polynomial(Q, z / 2 - 4.75f);
+         result = Y + tools::evaluate_polynomial(P, T(z / 2 - 4.75f)) / tools::evaluate_polynomial(Q, T(z / 2 - 4.75f));
          result *= exp(-z * z) / z;
       }
       else
@@ -926,34 +936,34 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
          // Max Error found at long double precision =   1.162590e-35
          static const T Y = 0.5632686614990234375f;
          static const T P[] = {    
-            0.000920922048732849448079451574171836943L,
-            0.00321439044532288750501700028748922439L,
-            -0.250455263029390118657884864261823431L,
-            -0.906807635364090342031792404764598142L,
-            -8.92233572835991735876688745989985565L,
-            -21.7797433494422564811782116907878495L,
-            -91.1451915251976354349734589601171659L,
-            -144.1279109655993927069052125017673L,
-            -313.845076581796338665519022313775589L,
-            -273.11378811923343424081101235736475L,
-            -271.651566205951067025696102600443452L,
-            -60.0530577077238079968843307523245547L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.000920922048732849448079451574171836943),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 0.00321439044532288750501700028748922439),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.250455263029390118657884864261823431),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -0.906807635364090342031792404764598142),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -8.92233572835991735876688745989985565),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -21.7797433494422564811782116907878495),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -91.1451915251976354349734589601171659),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -144.1279109655993927069052125017673),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -313.845076581796338665519022313775589),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -273.11378811923343424081101235736475),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -271.651566205951067025696102600443452),
+            BOOST_MATH_BIG_CONSTANT(T, 113, -60.0530577077238079968843307523245547),
          };
          static const T Q[] = {    
-            1L,
-            3.49040448075464744191022350947892036L,
-            34.3563592467165971295915749548313227L,
-            84.4993232033879023178285731843850461L,
-            376.005865281206894120659401340373818L,
-            629.95369438888946233003926191755125L,
-            1568.35771983533158591604513304269098L,
-            1646.02452040831961063640827116581021L,
-            2299.96860633240298708910425594484895L,
-            1222.73204392037452750381340219906374L,
-            799.359797306084372350264298361110448L,
-            72.7415265778588087243442792401576737L,
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 3.49040448075464744191022350947892036),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 34.3563592467165971295915749548313227),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 84.4993232033879023178285731843850461),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 376.005865281206894120659401340373818),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 629.95369438888946233003926191755125),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1568.35771983533158591604513304269098),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1646.02452040831961063640827116581021),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 2299.96860633240298708910425594484895),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 1222.73204392037452750381340219906374),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 799.359797306084372350264298361110448),
+            BOOST_MATH_BIG_CONSTANT(T, 113, 72.7415265778588087243442792401576737),
          };
-         result = Y + tools::evaluate_polynomial(P, 1 / z) / tools::evaluate_polynomial(Q, 1 / z);
+         result = Y + tools::evaluate_polynomial(P, T(1 / z)) / tools::evaluate_polynomial(Q, T(1 / z));
          result *= exp(-z * z) / z;
       }
    }
@@ -972,7 +982,60 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t)
    }
 
    return result;
-} // template <class T, class L>T erf_imp(T z, bool invert, const L& l, const mpl::int_<113>& t)
+} // template <class T, class Lanczos>T erf_imp(T z, bool invert, const Lanczos& l, const mpl::int_<113>& t)
+
+template <class T, class Policy, class tag>
+struct erf_initializer
+{
+   struct init
+   {
+      init()
+      {
+         do_init(tag());
+      }
+      static void do_init(const mpl::int_<0>&){}
+      static void do_init(const mpl::int_<53>&)
+      {
+         boost::math::erf(static_cast<T>(1e-12), Policy());
+         boost::math::erf(static_cast<T>(0.25), Policy());
+         boost::math::erf(static_cast<T>(1.25), Policy());
+         boost::math::erf(static_cast<T>(2.25), Policy());
+         boost::math::erf(static_cast<T>(4.25), Policy());
+         boost::math::erf(static_cast<T>(5.25), Policy());
+      }
+      static void do_init(const mpl::int_<64>&)
+      {
+         boost::math::erf(static_cast<T>(1e-12), Policy());
+         boost::math::erf(static_cast<T>(0.25), Policy());
+         boost::math::erf(static_cast<T>(1.25), Policy());
+         boost::math::erf(static_cast<T>(2.25), Policy());
+         boost::math::erf(static_cast<T>(4.25), Policy());
+         boost::math::erf(static_cast<T>(5.25), Policy());
+      }
+      static void do_init(const mpl::int_<113>&)
+      {
+         boost::math::erf(static_cast<T>(1e-22), Policy());
+         boost::math::erf(static_cast<T>(0.25), Policy());
+         boost::math::erf(static_cast<T>(1.25), Policy());
+         boost::math::erf(static_cast<T>(2.125), Policy());
+         boost::math::erf(static_cast<T>(2.75), Policy());
+         boost::math::erf(static_cast<T>(3.25), Policy());
+         boost::math::erf(static_cast<T>(5.25), Policy());
+         boost::math::erf(static_cast<T>(7.25), Policy());
+         boost::math::erf(static_cast<T>(11.25), Policy());
+         boost::math::erf(static_cast<T>(12.5), Policy());
+      }
+      void force_instantiate()const{}
+   };
+   static const init initializer;
+   static void force_instantiate()
+   {
+      initializer.force_instantiate();
+   }
+};
+
+template <class T, class Policy, class tag>
+const typename erf_initializer<T, Policy, tag>::init erf_initializer<T, Policy, tag>::initializer;
 
 } // namespace detail
 
@@ -1013,6 +1076,8 @@ inline typename tools::promote_args<T>::type erf(T z, const Policy& /* pol */)
 
    BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name());
 
+   detail::erf_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); // Force constants to be initialized before main
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::erf_imp(
       static_cast<value_type>(z),
       false,
@@ -1057,6 +1122,8 @@ inline typename tools::promote_args<T>::type erfc(T z, const Policy& /* pol */)
 
    BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name());
 
+   detail::erf_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); // Force constants to be initialized before main
+
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::erf_imp(
       static_cast<value_type>(z),
       true,