diff --git a/doc/quadrature/adaptive_trapezoidal.qbk b/doc/quadrature/adaptive_trapezoidal.qbk
new file mode 100644
index 000000000..548de93f2
--- /dev/null
+++ b/doc/quadrature/adaptive_trapezoidal.qbk
@@ -0,0 +1,61 @@
+[/
+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:adaptive_trapezoidal Adaptive Trapezoidal Quadrature]
+
+[heading Synopsis]
+
+``
+  #include <boost/math/quadrature/adaptive_trapezoidal.hpp>
+``
+
+[heading Description]
+
+The function __adaptive_trapezoidal calculates the integral of a function /f/ using the surprisingly simple trapezoidal rule.
+If we assume only that the integrand is twice continuously differentiable, we can prove that the error of the composite trapezoidal rule
+is [bigo](h^2).
+Hence halving the interval only cuts the error by about a fourth, which in turn implies that we must evaluate the function many times before an acceptable accuracy can be achieved.
+
+However, the trapezoidal rule has an astonishing property: If the integrand is periodic, and we integrate it over a period, then the trapezoidal rule converges exponentially fast.
+This can be seen by examination of the Euler-Maclaurin formula, which relates a definite integral to its trapezoidal sum and error terms proportional to the derivatives of the function at the endpoints.
+If the derivatives at the endpoints are the same or vanish, then the error very nearly vanishes.
+Hence the trapezoidal rule is essentially optimal for periodic integrands.
+
+In it's simplest form, an integration can be performed by the following code
+
+    __auto f = [](double x) { return 1/(5 - 4*cos(x)); };
+    __double I = boost::math::adaptive_trapezoidal(f, 0, two_pi<double>());
+
+Since the routine is adaptive, step sizes are halved continuously until a tolerance is reached.
+In order to control this tolerance, simply call the routine with an additional argument
+
+    __double I = adaptive_trapezoidal(f, 0, two_pi<double>(), 1e-6);
+
+The routine stops when successive estimates of the integral /I1/ and /I0/ differ by less than the tolerance multiplied by the estimated L1 norm of the function.
+A question arises as to what to do when successive estimates never pass below this threshold.
+The stepsize would be halved until it eventually would be flushed to zero, leading to an infinite loop.
+As such, you may pass an optional argument __max_refinements which controls how many times the interval may be halved before giving up.
+By default, this maximum number of refinement steps is 15, leading to ~30,000 function evaluations.
+For certain problems, this is clearly excessive, so specifying a smaller number is reasonable
+
+    __size_t max_refinements = 10;
+    __double I = adaptive_trapezoidal(f, 0, two_pi<double>(), 1e-6, max_refinements);
+
+Finally, you might wonder what the final error estimate actually was: This is achieved by passing a pointer into the last argument
+
+    __double error_estimate;
+    __double I = adaptive_trapezoidal(f, 0, two_pi<double>(), 1e-6, 10, &error_estimate);
+
+If need be, you can query this value, and if it is unacceptably large, use another method.
+
+
+References:
+
+Stoer, Josef, and Roland Bulirsch. Introduction to numerical analysis. Vol. 12. Springer Science & Business Media, 2013.
+
+
+[endsect]
diff --git a/example/adaptive_trapezoidal_example.cpp b/example/adaptive_trapezoidal_example.cpp
new file mode 100644
index 000000000..af833dfd6
--- /dev/null
+++ b/example/adaptive_trapezoidal_example.cpp
@@ -0,0 +1,25 @@
+/*
+ * Copyright Nick Thompson, 2017
+ *
+ * This example shows to to numerically integrate a periodic function using the adaptive_trapezoidal routine provided by boost.
+ */
+
+#include <iostream>
+#include <cmath>
+#include <limits>
+#include <boost/math/quadrature/adaptive_trapezoidal.hpp>
+
+int main()
+{
+    using boost::math::constants::two_pi;
+    using boost::math::constants::third;
+    // This function has an analytic form for its integral over a period: 2pi/3.
+    auto f = [](double x) { return 1/(5 - 4*cos(x)); };
+
+    double Q = boost::math::adaptive_trapezoidal(f, (double) 0, two_pi<double>());
+
+    std::cout << std::setprecision(std::numeric_limits<double>::digits10);
+    std::cout << "The adaptive trapezoidal rule gives the integral of our function as " << Q << "\n";
+    std::cout << "The exact result is                                                 " << two_pi<double>()*third<double>() << "\n";
+
+}
diff --git a/include/boost/math/quadrature/adaptive_trapezoidal.hpp b/include/boost/math/quadrature/adaptive_trapezoidal.hpp
new file mode 100644
index 000000000..42bb41ffe
--- /dev/null
+++ b/include/boost/math/quadrature/adaptive_trapezoidal.hpp
@@ -0,0 +1,90 @@
+/*
+ * 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)
+ *
+ * Use the adaptive trapezoidal rule to estimate the integral of periodic functions over a period,
+ * or to integrate a function whose derivative vanishes at the endpoints.
+ *
+ * If your function does not satisfy these conditions, and instead is simply continuous and bounded
+ * over the whole interval, then this routine will still converge, albeit slowly. However, there
+ * are much more efficient methods in this case, including Romberg, Simpson, and double exponential quadrature.
+ */
+
+#ifndef BOOST_MATH_QUADRATURE_ADAPTIVE_TRAPEZOIDAL_HPP
+#define BOOST_MATH_QUADRATURE_ADAPTIVE_TRAPEZOIDAL_HPP
+
+#include <cmath>
+#include <limits>
+#include <boost/math/constants/constants.hpp>
+
+namespace boost{ namespace math{
+
+template<class F, class Real>
+Real adaptive_trapezoidal(F f, Real a, Real b, Real tol = sqrt(std::numeric_limits<Real>::epsilon()), size_t max_refinements = 15, Real* error_estimate = nullptr)
+{
+    using std::abs;
+    using std::isfinite;
+    using boost::math::constants::half;
+    if(a >= b)
+    {
+        throw std::domain_error("a < b for integration over the region [a, b] is required.\n");
+    }
+    if (!isfinite(a))
+    {
+        throw std::domain_error("Left endpoint of integration must be finite for adaptive trapezoidal integration.\n");
+    }
+    if (!isfinite(b))
+    {
+        throw std::domain_error("Right endpoint of integration must be finite for adaptive trapedzoidal integration.\n");
+    }
+
+    Real ya = f(a);
+    Real yb = f(b);
+    Real h = (b - a)*half<Real>();
+    Real I0 = (ya + yb)*h;
+    Real IL0 = (abs(ya) + abs(yb))*h;
+
+    Real yh = f(a + h);
+    Real I1 = half<Real>()*I0 + yh*h;
+    Real IL1 = half<Real>()*IL0 + abs(yh)*h;
+
+    // The recursion is:
+    // I_k = 1/2 I_{k-1} + 1/2^k \sum_{j=1; j odd, j < 2^k} f(a + j(b-a)/2^k)
+    size_t k = 2;
+    // We want to go through at least 5 levels so we have sampled the function at least 33 times.
+    // Otherwise, we could terminate prematurely and miss essential features.
+    // This is of course possible anyway, but 33 samples seems to be a reasonable compromise.
+    Real error = abs(I0 - I1);
+    while (k < 5 || (k < max_refinements && error > tol*IL1) )
+    {
+        I0 = I1;
+        IL0 = IL1;
+
+        I1 = half<Real>()*I0;
+        IL1 = half<Real>()*IL0;
+        size_t p = 1 << k;
+        h *= half<Real>();
+        Real sum = 0;
+        Real absum = 0;
+        for(size_t j = 1; j < p; j += 2)
+        {
+            Real y = f(a + j*h);
+            sum += y;
+            absum += abs(y);
+        }
+        I1 += sum*h;
+        IL1 += absum*h;
+        ++k;
+        error = abs(I0 - I1);
+        if(error_estimate)
+        {
+            *error_estimate = error;
+        }
+    }
+    return I1;
+}
+
+}}
+#endif
diff --git a/test/test_adaptive_trapezoidal.cpp b/test/test_adaptive_trapezoidal.cpp
new file mode 100644
index 000000000..93fd6cb6a
--- /dev/null
+++ b/test/test_adaptive_trapezoidal.cpp
@@ -0,0 +1,154 @@
+#define BOOST_TEST_MODULE adaptive_trapezoidal
+
+#include <boost/type_index.hpp>
+#include <boost/test/included/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+#include <boost/math/quadrature/adaptive_trapezoidal.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+#ifdef __GNUC__
+#ifndef __clang__
+#include <boost/multiprecision/float128.hpp>
+#endif
+#endif
+
+using boost::multiprecision::cpp_bin_float_50;
+using boost::multiprecision::cpp_bin_float_100;
+
+template<class Real>
+void test_constant()
+{
+    std::cout << "Testing constants are integrated correctly by the adaptive trapezoidal routine on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
+
+    auto f = [](Real x) { return boost::math::constants::half<Real>(); };
+
+    Real Q = boost::math::adaptive_trapezoidal<decltype(f), Real>(f, (Real) 0.0, (Real) 10.0);
+
+    BOOST_CHECK_CLOSE(Q, 5.0, 100*std::numeric_limits<Real>::epsilon());
+}
+
+
+template<class Real>
+void test_rational_periodic()
+{
+    using boost::math::constants::two_pi;
+    using boost::math::constants::third;
+    std::cout << "Testing that rational periodic functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
+
+    auto f = [](Real x) { return 1/(5 - 4*cos(x)); };
+
+    Real tol = 100*std::numeric_limits<Real>::epsilon();
+    Real Q = boost::math::adaptive_trapezoidal(f, (Real) 0.0, two_pi<Real>(), tol);
+
+    BOOST_CHECK_CLOSE(Q, two_pi<Real>()*third<Real>(), 200*tol);
+}
+
+template<class Real>
+void test_bump_function()
+{
+    std::cout << "Testing that bump functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
+    auto f = [](Real x) {
+        if( x>= 1 || x <= -1)
+        {
+            return (Real) 0;
+        }
+        return (Real) exp(-(Real) 1/(1-x*x));
+    };
+    Real tol = 100*std::numeric_limits<Real>::epsilon();
+    Real Q = boost::math::adaptive_trapezoidal(f, (Real) -1, (Real) 1, tol);
+    // This is all the digits Mathematica gave me!
+    BOOST_CHECK_CLOSE(Q, 0.443994, 1e-4);
+}
+
+template<class Real>
+void test_zero_function()
+{
+    std::cout << "Testing that zero functions are integrated correctly by trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
+    auto f = [](Real x) { return (Real) 0;};
+    Real tol = 100*std::numeric_limits<Real>::epsilon();
+    Real Q = boost::math::adaptive_trapezoidal(f, (Real) -1, (Real) 1, tol);
+    BOOST_CHECK_SMALL(Q, 100*tol);
+}
+
+template<class Real>
+void test_sinsq()
+{
+    std::cout << "Testing that sin(x)^2 is integrated correctly by the trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
+    auto f = [](Real x) { return sin(10*x)*sin(10*x); };
+    Real tol = 100*std::numeric_limits<Real>::epsilon();
+    Real Q = boost::math::adaptive_trapezoidal(f, (Real) 0, (Real) boost::math::constants::pi<Real>(), tol);
+    BOOST_CHECK_CLOSE(Q, boost::math::constants::half_pi<Real>(), 100*tol);
+
+}
+
+template<class Real>
+void test_slowly_converging()
+{
+    std::cout << "Testing that non-periodic functions are correctly integrated by the trapezoidal rule, even if slowly, on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
+    // This function is not periodic, so it should not be fast to converge:
+    auto f = [](Real x) { return sqrt(1 - x*x); };
+
+    Real tol = sqrt(sqrt(std::numeric_limits<Real>::epsilon()));
+    Real error_estimate;
+    Real Q = boost::math::adaptive_trapezoidal(f, (Real) 0, (Real) 1, tol, 15, &error_estimate);
+    BOOST_CHECK_CLOSE(Q, boost::math::constants::half_pi<Real>()/2, 1000*tol);
+}
+
+template<class Real>
+void test_rational_sin()
+{
+    using std::pow;
+    using boost::math::constants::two_pi;
+    std::cout << "Testing that a rational sin function is integrated correctly by the trapezoidal rule on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
+    Real a = 5;
+    auto f = [=](Real x) { Real t = a + sin(x); return 1.0/(t*t); };
+
+    Real expected = two_pi<Real>()*a/pow(a*a - 1, 1.5);
+    Real tol = 100*std::numeric_limits<Real>::epsilon();
+    Real Q = boost::math::adaptive_trapezoidal(f, (Real) 0, (Real) boost::math::constants::two_pi<Real>(), tol);
+    BOOST_CHECK_CLOSE(Q, expected, 100*tol);
+}
+
+BOOST_AUTO_TEST_CASE(adaptive_trapezoidal)
+{
+    test_constant<float>();
+    test_constant<double>();
+    test_constant<long double>();
+    test_constant<cpp_bin_float_50>();
+    test_constant<cpp_bin_float_100>();
+
+    test_rational_periodic<float>();
+    test_rational_periodic<double>();
+    test_rational_periodic<long double>();
+    test_rational_periodic<cpp_bin_float_50>();
+    test_rational_periodic<cpp_bin_float_100>();
+
+    test_bump_function<long double>();
+
+    test_zero_function<float>();
+    test_zero_function<double>();
+    test_zero_function<long double>();
+    test_zero_function<cpp_bin_float_50>();
+    test_zero_function<cpp_bin_float_100>();
+
+    test_sinsq<float>();
+    test_sinsq<double>();
+    test_sinsq<long double>();
+    test_sinsq<cpp_bin_float_50>();
+    test_sinsq<cpp_bin_float_100>();
+
+    test_slowly_converging<float>();
+    test_slowly_converging<double>();
+    test_slowly_converging<long double>();
+
+    test_rational_sin<float>();
+    test_rational_sin<double>();
+    test_rational_sin<long double>();
+
+#ifdef __GNUC__
+#ifndef __clang__
+    test_constant<boost::multiprecision::float128>();
+    test_rational_periodic<boost::multiprecision::float128>();
+#endif
+#endif
+}