diff --git a/doc/math.qbk b/doc/math.qbk
index 6232de51f..fba3de2df 100644
--- a/doc/math.qbk
+++ b/doc/math.qbk
@@ -633,6 +633,7 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [section:quadrature Quadrature]
 [include quadrature/trapezoidal.qbk]
 [include quadrature/double_exponential.qbk]
+[include quadrature/gauss_kronrod.qbk]
 [endsect]
 [endmathpart] [/mathpart roots Root Finding Algorithms]
 
diff --git a/doc/quadrature/gauss_kronrod.qbk b/doc/quadrature/gauss_kronrod.qbk
new file mode 100644
index 000000000..9fc692101
--- /dev/null
+++ b/doc/quadrature/gauss_kronrod.qbk
@@ -0,0 +1,43 @@
+[/
+Copyright (c) 2017 Nick Thompson
+Use, modification and distribution are subject to the
+Boost Software License, Version 1.0. (See accompanying file
+LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+]
+
+[section:gauss_kronrod Gauss-Kronrod Quadrature]
+
+Gauss-Kronrod quadrature is an extension of Gaussian quadrature which provides an a posteriori error estimate for the integral.
+
+The idea behind Gaussian quadrature is to choose /n/ nodes and weights in such a way that polynomials of order /2n-1/ are integrated exactly.
+However, integration of polynomials is trivial, so it is rarely done via numerical methods.
+Instead, transcendental and numerically defined functions are integrated via Gaussian quadrature, and the defining problem becomes how to estimate the remainder.
+Gaussian quadrature alone (without some form on interval splitting) cannot answer this question.
+
+It is possible to compute a Gaussian quadrature of order /n/ and another of order (say) /2n+1/, and use the difference as an error estimate.
+However, this is not optimal, as the zeros of the Legendre polynomials (nodes of the Gaussian quadrature) never the same for different orders, so /3n+1/ function evaluations must be performed.
+Kronrod considered the problem of how to interleave nodes into a Gaussian quadrature in such a way that all previous function evaluations can be reused, while increasing the order of polynomials that can be integrated exactly.
+This allows an a posteriori error estimate to be provided while still preserving exponential convergence.
+Kronrod discovered that by adding /n+1/ nodes (computed from the zeros of the Legendre-Stieltjes polynomials) to a Gaussian quadrature of order /n/, he could integrate polynomials of order /3n+1/.
+
+[section:caveats Caveats]
+
+Although this code allows a user to use any floating point type, the nodes and weights have only been computed to 100 decimal digits.
+We think this should be acceptable for most use cases; if more digits are needed we recommend the user switch to `tanh_sinh` quadrature,
+which is truly arbitrary precision.
+
+For essentially all analytic integrands bounded on the domain, the error estimates provided by the routine are woefully pessimistic.
+In fact, in this case the error is very nearly equal to the error of the Gaussian quadrature formula of order /n/,
+whereas the Kronrod extension converges exponentially beyond the Gaussian estimate.
+Very sophisticated method exist to estimate the error, but all require the integrand to lie in a particular function space.
+A more sophisticated a posteriori error estimate for an element of a particular function space is left to the user.
+
+
+[section:gk_refes References]
+
+* Kronrod, Aleksandr Semenovish (1965), ['Nodes and weights of quadrature formulas. Sixteen-place tables], New York: Consultants Bureau
+* Dirk P. Laurie, ['Calculation of Gauss-Kronrod Quadrature Rules], Mathematics of Computation, Volume 66, Number 219, 1997
+* Gonnet, Pedro, ['A Review of Error Estimation in Adaptive Quadrature], https://arxiv.org/pdf/1003.4629.pdf
+
+[endsect]
+[endsect]
diff --git a/example/legendre_stieltjes_example.cpp b/example/legendre_stieltjes_example.cpp
index ab92e68c4..00ceba277 100644
--- a/example/legendre_stieltjes_example.cpp
+++ b/example/legendre_stieltjes_example.cpp
@@ -8,6 +8,7 @@
 #include <string>
 #include <boost/math/constants/constants.hpp>
 #include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp>
 #include <boost/math/special_functions/legendre_stieltjes.hpp>
 
 using boost::math::legendre_p;
@@ -16,83 +17,94 @@ using boost::math::legendre_p_prime;
 using boost::math::legendre_stieltjes;
 using boost::multiprecision::cpp_bin_float_quad;
 using boost::multiprecision::cpp_bin_float_100;
+using boost::multiprecision::cpp_dec_float_100;
 
 template<class Real>
-void gauss_konrod_rule(size_t order)
+void gauss_kronrod_rule(size_t order)
 {
     std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
     std::cout << std::fixed;
     auto gauss_nodes = boost::math::legendre_p_zeros<Real>(order);
     auto E = legendre_stieltjes<Real>(order + 1);
     std::vector<Real> gauss_weights(gauss_nodes.size(), std::numeric_limits<Real>::quiet_NaN());
-    std::vector<Real> gauss_konrod_weights(gauss_nodes.size(), std::numeric_limits<Real>::quiet_NaN());
+    std::vector<Real> gauss_kronrod_weights(gauss_nodes.size(), std::numeric_limits<Real>::quiet_NaN());
     for (size_t i = 0; i < gauss_nodes.size(); ++i)
     {
         Real node = gauss_nodes[i];
         Real lp = legendre_p_prime<Real>(order, node);
         gauss_weights[i] = 2/( (1-node*node)*lp*lp);
         // P_n(x) = (2n)!/(2^n (n!)^2) pi_n(x), where pi_n is the monic Legendre polynomial.
-        gauss_konrod_weights[i] = gauss_weights[i] + static_cast<Real>(2)/(static_cast<Real>(order+1)*legendre_p_prime(order, node)*E(node));
+        gauss_kronrod_weights[i] = gauss_weights[i] + static_cast<Real>(2)/(static_cast<Real>(order+1)*legendre_p_prime(order, node)*E(node));
     }
 
-    std::cout << "Gauss Nodes:\n";
+    std::cout << "static const std::vector<Real> gauss_nodes {\n";
     for (auto const & node : gauss_nodes)
     {
-        std::cout << node << "\n";
+        std::cout << "    boost::lexical_cast<Real>(\"" << node << "\"),\n";
     }
+    std::cout << "};\n\n";
 
-    std::cout << "Gauss Weights:\n";
+    std::cout << "static const std::vector<Real> gauss_weights {\n";
     for (auto const & weight : gauss_weights)
     {
-        std::cout << weight << "\n";
+        std::cout << "    boost::lexical_cast<Real>(\"" << weight << "\"),\n";
+    }
+    std::cout << "};\n\n";
+
+    std::cout << "static const std::vector<Real> gauss_kronrod_weights {\n";
+    for (auto const & weight : gauss_kronrod_weights)
+    {
+        std::cout << "    boost::lexical_cast<Real>(\"" << weight << "\"),\n";
+    }
+    std::cout << "};\n\n";
+
+    auto kronrod_nodes = E.zeros();
+    std::vector<Real> kronrod_weights(kronrod_nodes.size());
+    for (size_t i = 0; i < kronrod_weights.size(); ++i)
+    {
+        Real node = kronrod_nodes[i];
+        kronrod_weights[i] = static_cast<Real>(2)/(static_cast<Real>(order+1)*legendre_p(order, node)*E.prime(node));
     }
 
-    std::cout << "Gauss-Konrod weights: \n";
-    for (auto const & w : gauss_konrod_weights)
+    std::cout << "static const std::vector<Real> kronrod_nodes {\n";
+    for (auto node : kronrod_nodes)
     {
-        std::cout << w << "\n";
+        std::cout << "    boost::lexical_cast<Real>(\"" << node << "\"),\n";
     }
+    std::cout << "};\n\n";
 
-    auto konrod_nodes = E.zeros();
-    std::vector<Real> konrod_weights(konrod_nodes.size());
-    for (size_t i = 0; i < konrod_weights.size(); ++i)
+    std::cout << "static const std::vector<Real> kronrod_weights {\n";
+    for (auto const & weight : kronrod_weights)
     {
-        Real node = konrod_nodes[i];
-        konrod_weights[i] = static_cast<Real>(2)/(static_cast<Real>(order+1)*legendre_p(order, node)*E.prime(node));
-    }
-
-    std::cout << "Konrod nodes:\n";
-    for (auto node : konrod_nodes)
-    {
-        std::cout << node << "\n";
-    }
-
-    std::cout << "Konrod weights: \n";
-    for (auto const & w : gauss_konrod_weights)
-    {
-        std::cout << w << "\n";
+        std::cout << "    boost::lexical_cast<Real>(\"" << weight << "\"),\n";
     }
+    std::cout << "};\n\n";
 
 }
 
 int main()
 {
-    std::cout << "Gauss-Konrod 7-15 Rule:\n";
-    gauss_konrod_rule<cpp_bin_float_100>(7);
+    //std::cout << "Gauss-Kronrod 7-15 Rule:\n";
+    //gauss_kronrod_rule<cpp_dec_float_100>(7);
 
-    std::cout << "\n\nGauss-Konrod 10-21 Rule:\n";
-    gauss_konrod_rule<cpp_bin_float_100>(10);
+    //std::cout << "\n\nGauss-Kronrod 10-21 Rule:\n";
+    //gauss_kronrod_rule<cpp_dec_float_100>(10);
 
-    std::cout << "\n\nGauss-Konrod 15-31 Rule:\n";
-    gauss_konrod_rule<cpp_bin_float_100>(15);
+    std::cout << "\n\nGauss-Kronrod 15-31 Rule:\n";
+    gauss_kronrod_rule<cpp_dec_float_100>(15);
+    /*
+    std::cout << "\n\nGauss-Kronrod 20-41 Rule:\n";
+    gauss_kronrod_rule<cpp_dec_float_100>(20);
 
-    std::cout << "\n\nGauss-Konrod 20-41 Rule:\n";
-    gauss_konrod_rule<cpp_bin_float_100>(20);
+    std::cout << "\n\nGauss-Kronrod 25-51 Rule:\n";
+    gauss_kronrod_rule<cpp_dec_float_100>(25);
 
-    std::cout << "\n\nGauss-Konrod 25-51 Rule:\n";
-    gauss_konrod_rule<cpp_bin_float_100>(25);
+    std::cout << "\n\nGauss-Kronrod 30-61 Rule:\n";
+    gauss_kronrod_rule<cpp_dec_float_100>(30);
 
-    std::cout << "\n\nGauss-Konrod 30-61 Rule:\n";
-    gauss_konrod_rule<cpp_bin_float_100>(30);
+    std::cout << "\n\nGauss-Kronrod 35-71 Rule:\n";
+    gauss_kronrod_rule<cpp_dec_float_100>(35);
 
+    std::cout << "\n\nGauss-Kronrod 40-81 Rule:\n";
+    gauss_kronrod_rule<cpp_dec_float_100>(40);*/
 }
diff --git a/include/boost/math/quadrature/detail/gauss_kronrod_detail.hpp b/include/boost/math/quadrature/detail/gauss_kronrod_detail.hpp
new file mode 100644
index 000000000..18de59779
--- /dev/null
+++ b/include/boost/math/quadrature/detail/gauss_kronrod_detail.hpp
@@ -0,0 +1,142 @@
+// Copyright Nick Thompson 2017.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <vector>
+
+namespace boost { namespace math { namespace quadrature { namespace detail {
+
+// These nodes and weights can be verified at:
+// https://www.advanpix.com/2011/11/07/gauss-kronrod-quadrature-nodes-weights/
+template<class Real, class F>
+Real g7_k15(const F& f, Real* error = nullptr)
+{
+    static const std::vector<Real> gauss_nodes{
+        boost::lexical_cast<Real>("0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
+        boost::lexical_cast<Real>("0.4058451513773971669066064120769614633473820140993701263870432517946638132261256553283126897277465878"),
+        boost::lexical_cast<Real>("0.7415311855993944398638647732807884070741476471413902601199553519674298746721805137928268323668632471"),
+        boost::lexical_cast<Real>("0.9491079123427585245261896840478512624007709376706177835487691039130633303548401408057307700279257241")
+    };
+
+    static const std::vector<Real> gauss_weights {
+        boost::lexical_cast<Real>("0.4179591836734693877551020408163265306122448979591836734693877551020408163265306122448979591836734694"),
+        boost::lexical_cast<Real>("0.3818300505051189449503697754889751338783650835338627347510834510307055464341297083486846593440448015"),
+        boost::lexical_cast<Real>("0.2797053914892766679014677714237795824869250652265987645370140326936188104305626768132409429011976188"),
+        boost::lexical_cast<Real>("0.1294849661688696932706114326790820183285874022599466639772086387246552349720423087156254181629208451")
+    };
+
+    static const std::vector<Real> gauss_kronrod_weights {
+        boost::lexical_cast<Real>("0.2094821410847278280129991748917142636977620802237043167129980065613751513232564861681690821167594910"),
+        boost::lexical_cast<Real>("0.1903505780647854099132564024210136828260780754553583558854408803674405807241021267960596460510637759"),
+        boost::lexical_cast<Real>("0.1406532597155259187451895905102379203998897572479985755617454689331270809309095040809737912241555591"),
+        boost::lexical_cast<Real>("0.0630920926299785532907006631892042866650711572115507071136055451469839974779648749281991702645044420"),
+    };
+
+    static const std::vector<Real> kronrod_nodes {
+        boost::lexical_cast<Real>("0.2077849550078984676006894037732449134797844071451706497138457346198669384494352022691034322718369853"),
+        boost::lexical_cast<Real>("0.5860872354676911302941448382587295984367807506043609513049928931988037360744440746451167449893594210"),
+        boost::lexical_cast<Real>("0.8648644233597690727897127886409262012109723070740881486014577127670677081325957210358584785960459054"),
+        boost::lexical_cast<Real>("0.9914553711208126392068546975263285166420443383703347012910874135724417393465340723592450350962684176")
+    };
+
+    static const std::vector<Real> kronrod_weights {
+        boost::lexical_cast<Real>("0.2044329400752988924141619992346490847165176041807183574244709531204546769854659887934837429200934755"),
+        boost::lexical_cast<Real>("0.1690047266392679028265834265985502841062449003029442414973400675569568092161902911293670240385535991"),
+        boost::lexical_cast<Real>("0.1047900103222501838398763225415180174437566542138306118933906513396374632157628952416757162750931133"),
+        boost::lexical_cast<Real>("0.0229353220105292249637320080589695919935608112757469922675074302547118157879760759461563681681562895"),
+    };
+
+    Real f0 = f(0);
+    Real I0 = f0*gauss_weights[0];
+    Real I1 = f0*gauss_kronrod_weights[0];
+    for (size_t i = 1; i < gauss_nodes.size(); ++i) {
+        Real yp = f(gauss_nodes[i]);
+        Real ym = f(-gauss_nodes[i]);
+        I0 += gauss_weights[i]*(yp + ym);
+        I1 += gauss_kronrod_weights[i]*(yp + ym);
+    }
+
+    for(size_t i = 0; i < kronrod_weights.size(); ++i) {
+        I1 += kronrod_weights[i]*(f(kronrod_nodes[i]) + f(-kronrod_nodes[i]));
+    }
+
+    if (error) {
+        using std::abs;
+        *error = abs(I1 - I0);
+    }
+    return I1;
+}
+
+
+template<class Real, class F>
+Real g10_k21(const F& f, Real* error = nullptr)
+{
+
+    static const std::vector<Real> gauss_nodes {
+        boost::lexical_cast<Real>("0.1488743389816312108848260011297199846175648594206916957079892535159036173556685213711776297994636912"),
+        boost::lexical_cast<Real>("0.4333953941292471907992659431657841622000718376562464965027015131437669890777035012251027579501177212"),
+        boost::lexical_cast<Real>("0.6794095682990244062343273651148735757692947118348094676648171889525585753950749246150785735704803795"),
+        boost::lexical_cast<Real>("0.8650633666889845107320966884234930485275430149653304525219597318453747551380555613567907289460457707"),
+        boost::lexical_cast<Real>("0.9739065285171717200779640120844520534282699466923821192312120666965952032346361596257235649562685563"),
+    };
+
+    static const std::vector<Real> gauss_weights {
+        boost::lexical_cast<Real>("0.2955242247147528701738929946513383294210467170268536013543080297559959382171523292703565957937542167"),
+        boost::lexical_cast<Real>("0.2692667193099963550912269215694693528597599384608837958005632762421534323191792767642266367092527608"),
+        boost::lexical_cast<Real>("0.2190863625159820439955349342281631924587718705226770898809565436351999106529512812426839931772021928"),
+        boost::lexical_cast<Real>("0.1494513491505805931457763396576973324025566396694273678354772687532386547266300109459472646347319519"),
+        boost::lexical_cast<Real>("0.0666713443086881375935688098933317928578648343201581451286948816134120640840871017767855096850588778"),
+    };
+
+    static const std::vector<Real> gauss_kronrod_weights {
+        boost::lexical_cast<Real>("0.1477391049013384913748415159720680455237316254852066045181919543988599301673569640573270395918288225"),
+        boost::lexical_cast<Real>("0.1347092173114733259280540017717068327609919130085597140663666849132029140012128203667695315948827177"),
+        boost::lexical_cast<Real>("0.1093871588022976418992105903258049602718132998343452200781967582982655037289143216868389943267455384"),
+        boost::lexical_cast<Real>("0.0750396748109199527670431409161900093952193820009100881736970480484304043428584951788138087306465541"),
+        boost::lexical_cast<Real>("0.0325581623079647274788189724593897606173889398456626095715375042327141218201654986923816076053846265"),
+    };
+
+    static const std::vector<Real> kronrod_nodes {
+        boost::lexical_cast<Real>("0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
+        boost::lexical_cast<Real>("0.2943928627014601981311266031038655661626866251569579186488822917272461116633273788844552317826823736"),
+        boost::lexical_cast<Real>("0.5627571346686046833390000992726941408430138819419669588603462145877926635321632754971208785416999242"),
+        boost::lexical_cast<Real>("0.7808177265864168970637175783450423771634075202981571797469485999950560798276142065452697723423899624"),
+        boost::lexical_cast<Real>("0.9301574913557082260012071800595083462251679099819392423034940686682841598309167305501119457285100788"),
+        boost::lexical_cast<Real>("0.9956571630258080807355272806890028479212605872194789243633791611175702304677486735715232599691207672"),
+    };
+
+    static const std::vector<Real> kronrod_weights {
+        boost::lexical_cast<Real>("0.1494455540029169056649364683898212037452363166874728038356085187369896447851184192572103070568954026"),
+        boost::lexical_cast<Real>("0.1427759385770600807970942731387170608859790565319055556074100474397077044990934002781113170628375643"),
+        boost::lexical_cast<Real>("0.1234919762620658510779581098310741595123003495286483276446799412097405423897545468968153862236373823"),
+        boost::lexical_cast<Real>("0.0931254545836976055350654650833663443900188288807600319700850387601777356722007752374141230616158275"),
+        boost::lexical_cast<Real>("0.0547558965743519960313813002445801763737211140583335575244326158047840989278189753251163015690032981"),
+        boost::lexical_cast<Real>("0.0116946388673718742780643960621920483962173324819318889275981475256222220580649926518067367049699673"),
+    };
+
+    Real I0 = 0;
+    Real I1 = 0;
+    for (size_t i = 0; i < gauss_nodes.size(); ++i) {
+        Real yp = f(gauss_nodes[i]);
+        Real ym = f(-gauss_nodes[i]);
+        I0 += gauss_weights[i]*(yp + ym);
+        I1 += gauss_kronrod_weights[i]*(yp + ym);
+    }
+
+    I1 += f(0)*kronrod_weights[0];
+    for(size_t i = 1; i < kronrod_weights.size(); ++i) {
+        I1 += kronrod_weights[i]*(f(kronrod_nodes[i]) + f(-kronrod_nodes[i]));
+    }
+
+    if (error) {
+        using std::abs;
+        *error = abs(I1 - I0);
+    }
+    return I1;
+}
+
+
+
+}}}}
diff --git a/test/gauss_kronrod_test.cpp b/test/gauss_kronrod_test.cpp
new file mode 100644
index 000000000..d2b36576e
--- /dev/null
+++ b/test/gauss_kronrod_test.cpp
@@ -0,0 +1,60 @@
+// Copyright Nick Thompson 2017.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE GaussKronrod
+
+#include <boost/test/unit_test.hpp>
+#include <boost/math/constants/constants.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/math/quadrature/detail/gauss_kronrod_detail.hpp>
+
+using boost::math::constants::half_pi;
+using boost::math::quadrature::detail::g7_k15;
+using boost::math::quadrature::detail::g10_k21;
+using boost::multiprecision::cpp_bin_float_quad;
+
+
+template<class Real>
+void test_g7_k15()
+{
+    std::cout << std::setprecision(std::numeric_limits<Real>::digits10) << std::fixed;
+    auto f = [](Real x) { return 1/(1+x*x); };
+    Real error;
+    Real Q = g7_k15<Real, decltype(f)>(f, &error);
+    std::cout << Q << std::endl;
+    std::cout << "Q - pi/2 = " << Q - half_pi<Real>() << std::endl;
+    std::cout << "Error estimate: " << error << std::endl;
+}
+
+template<class Real>
+void test_g10_k21()
+{
+    std::cout << std::setprecision(std::numeric_limits<Real>::digits10) << std::fixed;
+    auto f = [](Real x) { return 1/(1+x*x); };
+    Real error;
+    Real Q = g10_k21<Real, decltype(f)>(f, &error);
+    std::cout << Q << std::endl;
+    std::cout << "Q - pi/2 = " << Q - half_pi<Real>() << std::endl;
+    std::cout << "Error estimate: " << error << std::endl;
+}
+
+
+
+BOOST_AUTO_TEST_CASE(GaussKronrod)
+{
+    test_g7_k15<double>();
+    test_g7_k15<long double>();
+    test_g7_k15<cpp_bin_float_quad>();
+    test_g7_k15<boost::multiprecision::cpp_bin_float_100>();
+
+    test_g10_k21<double>();
+    test_g10_k21<long double>();
+    test_g10_k21<cpp_bin_float_quad>();
+    test_g10_k21<boost::multiprecision::cpp_bin_float_100>();
+
+}