Adding upstream version 14.2.21.upstream/14.2.21upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

28 files changed, 5080 insertions, 0 deletions

diff --git a/src/boost/libs/multiprecision/example/Jamfile.v2 b/src/boost/libs/multiprecision/example/Jamfile.v2
new file mode 100644
index 00000000..0e09c848
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/Jamfile.v2
@@ -0,0 +1,112 @@
+# \libs\math\example\jamfile.v2
+# Runs multiprecision examples.
+# Copyright 2014 John Maddock
+# Copyright Paul A. Bristow 2014.
+# Copyright Christpher Kormanyos 2014
+
+# Distributed under 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)
+
+# bring in the rules for testing
+import testing ;
+import modules ;
+import path ;
+import ../../config/checks/config : requires ;
+
+local ntl-path = [ modules.peek : NTL_PATH ] ;
+local gmp_path = [ modules.peek : GMP_PATH ] ;
+local mpfr_path = [ modules.peek : MPFR_PATH ] ;
+local mpfi_path = [ modules.peek : MPFI_PATH ] ;
+local tommath_path = [ modules.peek : TOMMATH_PATH ] ;
+
+project
+    : requirements
+     <include>$(gmp_path)
+     <include>$(gmp_path)/mpfr
+     <include>$(gmp_path)/gmpfrxx
+     <include>$(mpfr_path)
+     <include>$(mpfi_path)
+     <include>$(mpfi_path)/src
+     <include>$(tommath_path)
+     <include>../include
+     <include>../../..
+     
+     <toolset>gcc:<cxxflags>-Wno-missing-braces
+      
+      # Assembler error "File too big" caused by lots of C++ templates, for example, math/floating_point_examples.cpp.
+      # Some projects on some toolsets may require
+      #   <toolset>gcc-mingw:<cxxflags>\"-Wa,-mbig-obj\"
+      # See https://digitalkarabela.com/mingw-w64-how-to-fix-file-too-big-too-many-sections/
+      # <toolset>gcc-mingw:<cxxflags>-Wa,-mbig-obj # Some projects may overflow assembler and require equivalent of MSVC /bigobj.
+      # Requires version 2.30 of GNU binutils.
+      # Best applied only to projects that require this, see run math/floating_point_examples.cpp below.
+
+      <toolset>darwin:<cxxflags>-Wno-missing-braces
+      <toolset>acc:<cxxflags>+W2068,2461,2236,4070
+      <toolset>intel:<cxxflags>-Qwd264,239
+      <toolset>msvc:<runtime-link>static
+      <toolset>msvc:<link>static
+      <toolset>msvc:<warnings>all
+      <toolset>msvc:<asynch-exceptions>on
+      <toolset>msvc:<define>_CRT_SECURE_NO_DEPRECATE
+      <toolset>msvc:<define>_SCL_SECURE_NO_DEPRECATE
+      <toolset>msvc:<define>_SCL_SECURE_NO_WARNINGS
+      <toolset>msvc:<define>_CRT_SECURE_NO_WARNINGS
+      <toolset>msvc:<cxxflags>/wd4996
+      <toolset>msvc:<cxxflags>/wd4512
+      <toolset>msvc:<cxxflags>/wd4610
+      <toolset>msvc:<cxxflags>/wd4510
+      <toolset>msvc:<cxxflags>/wd4127
+      <toolset>msvc:<cxxflags>/wd4701
+      <toolset>msvc:<cxxflags>/wd4127
+      <toolset>msvc:<cxxflags>/wd4305
+    ;
+
+lib gmp : : <search>$(gmp_path) ;
+lib mpfr : : <search>$(gmp_path) <search>$(mpfr_path) <search>$(mpfr_path)/build.vc10/lib/Win32/Debug ;
+lib mpfi : : <search>$(gmp_path) <search>$(mpfr_path) <search>$(mpfr_path)/build.vc10/lib/Win32/Debug <search>$(mpfi_path) <search>$(mpfi_path)/src ;
+lib quadmath ;
+
+if $(tommath_path)
+{
+   lib tommath : [ GLOB $(tommath_path) : *.c ] ;
+   TOMMATH = tommath ;
+}
+else
+{
+   lib tommath : : <search>$(tommath_path) ;
+   TOMMATH = tommath ;
+}
+
+lib no_eh_eg_support : ../test/no_eh_test_support.cpp ;
+
+test-suite examples :
+
+      [ run cpp_int_snips.cpp no_eh_eg_support ]
+      [ run cpp_int_import_export.cpp no_eh_eg_support ]
+      [ run cpp_bin_float_import_export.cpp no_eh_eg_support ]
+
+      [ run cpp_dec_float_snips.cpp no_eh_eg_support ]
+
+      [ run cpp_bin_float_snips.cpp no_eh_eg_support ]
+
+      [ run debug_adaptor_snips.cpp no_eh_eg_support ]
+      [ run float128_snips.cpp quadmath no_eh_eg_support : : : [ check-target-builds ../config//has_float128 : : <build>no ] ]
+
+      [ run floating_point_examples.cpp no_eh_eg_support : : : <toolset>gcc-mingw:<cxxflags>-Wa,-mbig-obj ] # See note above.
+      [ run gauss_laguerre_quadrature.cpp no_eh_eg_support : : : release [ requires cxx11_lambdas ] ]
+      [ run hypergeometric_luke_algorithms.cpp no_eh_eg_support ../../chrono/build//boost_chrono ../../system/build//boost_system : : : [ requires cxx11_nullptr ]  ]
+      [ run integer_examples.cpp no_eh_eg_support ]
+      [ run logged_adaptor.cpp no_eh_eg_support mpfi mpfr gmp :  :  : [ check-target-builds ../config//has_mpfi : : <build>no ] ]
+      [ run mixed_integer_arithmetic.cpp no_eh_eg_support ]
+      [ run numeric_limits_snips.cpp no_eh_eg_support /boost//test_exec_monitor : : : [ requires cxx11_numeric_limits ]  [ check-target-builds ../config//has_float128 : <source>quadmath ] ]
+      [ run random_snips.cpp gmp no_eh_eg_support : : : [ requires cxx11_explicit_conversion_operators ] [ check-target-builds ../config//has_gmp : : <build>no ] ]
+      [ run safe_prime.cpp no_eh_eg_support ]
+
+      [ run gmp_snips.cpp gmp no_eh_eg_support : : : [ check-target-builds ../config//has_gmp : : <build>no ] ]
+      [ run mpfi_snips.cpp mpfi mpfr gmp no_eh_eg_support :  :  : [ check-target-builds ../config//has_mpfi : : <build>no ] ]
+      [ run mpfr_snips.cpp mpfr gmp no_eh_eg_support : : : [ check-target-builds ../config//has_mpfr : : <build>no ] ]
+      [ run tommath_snips.cpp $(TOMMATH) no_eh_eg_support : : : [ check-target-builds ../config//has_tommath : : <build>no ] ]
+      [ compile constexpr_float_arithmetic_examples.cpp : [ requires cxx14_constexpr cxx17_if_constexpr ] ]
+
+;
diff --git a/src/boost/libs/multiprecision/example/complex128_examples.cpp b/src/boost/libs/multiprecision/example/complex128_examples.cpp
new file mode 100644
index 00000000..49f57d18
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/complex128_examples.cpp
@@ -0,0 +1,120 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2018 Nick Thompson. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+/*`This example demonstrates the usage of the MPC backend for multiprecision complex numbers.
+In the following, we will show how using MPC backend allows for the same operations as the C++ standard library complex numbers.
+*/
+
+//[complex128_eg
+#include <iostream>
+#include <complex>
+#include <boost/multiprecision/complex128.hpp>
+
+template<class Complex>
+void complex_number_examples()
+{
+    Complex z1{0, 1};
+    std::cout << std::setprecision(std::numeric_limits<typename Complex::value_type>::digits10);
+    std::cout << std::scientific << std::fixed;
+    std::cout << "Print a complex number: " << z1 << std::endl;
+    std::cout << "Square it             : " << z1*z1 << std::endl;
+    std::cout << "Real part             : " << z1.real() << " = " << real(z1) << std::endl;
+    std::cout << "Imaginary part        : " << z1.imag() << " = " << imag(z1) << std::endl;
+    using std::abs;
+    std::cout << "Absolute value        : " << abs(z1) << std::endl;
+    std::cout << "Argument              : " << arg(z1) << std::endl;
+    std::cout << "Norm                  : " << norm(z1) << std::endl;
+    std::cout << "Complex conjugate     : " << conj(z1) << std::endl;
+    std::cout << "Projection onto Riemann sphere: " <<  proj(z1) << std::endl;
+    typename Complex::value_type r = 1;
+    typename Complex::value_type theta = 0.8;
+    using std::polar;
+    std::cout << "Polar coordinates (phase = 0)    : " << polar(r) << std::endl;
+    std::cout << "Polar coordinates (phase !=0)    : " << polar(r, theta) << std::endl;
+
+    std::cout << "\nElementary special functions:\n";
+    using std::exp;
+    std::cout << "exp(z1) = " << exp(z1) << std::endl;
+    using std::log;
+    std::cout << "log(z1) = " << log(z1) << std::endl;
+    using std::log10;
+    std::cout << "log10(z1) = " << log10(z1) << std::endl;
+    using std::pow;
+    std::cout << "pow(z1, z1) = " << pow(z1, z1) << std::endl;
+    using std::sqrt;
+    std::cout << "Take its square root  : " << sqrt(z1) << std::endl;
+    using std::sin;
+    std::cout << "sin(z1) = " << sin(z1) << std::endl;
+    using std::cos;
+    std::cout << "cos(z1) = " << cos(z1) << std::endl;
+    using std::tan;
+    std::cout << "tan(z1) = " << tan(z1) << std::endl;
+    using std::asin;
+    std::cout << "asin(z1) = " << asin(z1) << std::endl;
+    using std::acos;
+    std::cout << "acos(z1) = " << acos(z1) << std::endl;
+    using std::atan;
+    std::cout << "atan(z1) = " << atan(z1) << std::endl;
+    using std::sinh;
+    std::cout << "sinh(z1) = " << sinh(z1) << std::endl;
+    using std::cosh;
+    std::cout << "cosh(z1) = " << cosh(z1) << std::endl;
+    using std::tanh;
+    std::cout << "tanh(z1) = " << tanh(z1) << std::endl;
+    using std::asinh;
+    std::cout << "asinh(z1) = " << asinh(z1) << std::endl;
+    using std::acosh;
+    std::cout << "acosh(z1) = " << acosh(z1) << std::endl;
+    using std::atanh;
+    std::cout << "atanh(z1) = " << atanh(z1) << std::endl;
+}
+
+int main()
+{
+    std::cout << "First, some operations we usually perform with std::complex:\n";
+    complex_number_examples<std::complex<double>>();
+    std::cout << "\nNow the same operations performed using quad precision complex numbers:\n";
+    complex_number_examples<boost::multiprecision::complex128>();
+
+    return 0;
+}
+//]
+
+/*
+
+//[complex128_out
+
+Print a complex number: (0.000000000000000000000000000000000,1.000000000000000000000000000000000)
+Square it             : -1.000000000000000000000000000000000
+Real part             : 0.000000000000000000000000000000000 = 0.000000000000000000000000000000000
+Imaginary part        : 1.000000000000000000000000000000000 = 1.000000000000000000000000000000000
+Absolute value        : 1.000000000000000000000000000000000
+Argument              : 1.570796326794896619231321691639751
+Norm                  : 1.000000000000000000000000000000000
+Complex conjugate     : (0.000000000000000000000000000000000,-1.000000000000000000000000000000000)
+Projection onto Riemann sphere: (0.000000000000000000000000000000000,1.000000000000000000000000000000000)
+Polar coordinates (phase = 0)    : 1.000000000000000000000000000000000
+Polar coordinates (phase !=0)    : (0.696706709347165389063740022772449,0.717356090899522792567167815703377)
+
+Elementary special functions:
+exp(z1) = (0.540302305868139717400936607442977,0.841470984807896506652502321630299)
+log(z1) = (0.000000000000000000000000000000000,1.570796326794896619231321691639751)
+log10(z1) = (0.000000000000000000000000000000000,0.682188176920920673742891812715678)
+pow(z1, z1) = 0.207879576350761908546955619834979
+Take its square root  : (0.707106781186547524400844362104849,0.707106781186547524400844362104849)
+sin(z1) = (0.000000000000000000000000000000000,1.175201193643801456882381850595601)
+cos(z1) = 1.543080634815243778477905620757061
+tan(z1) = (0.000000000000000000000000000000000,0.761594155955764888119458282604794)
+asin(z1) = (0.000000000000000000000000000000000,0.881373587019543025232609324979792)
+acos(z1) = (1.570796326794896619231321691639751,-0.881373587019543025232609324979792)
+atan(z1) = (0.000000000000000000000000000000000,inf)
+sinh(z1) = (0.000000000000000000000000000000000,0.841470984807896506652502321630299)
+cosh(z1) = 0.540302305868139717400936607442977
+tanh(z1) = (0.000000000000000000000000000000000,1.557407724654902230506974807458360)
+asinh(z1) = (0.000000000000000000000000000000000,1.570796326794896619231321691639751)
+acosh(z1) = (0.881373587019543025232609324979792,1.570796326794896619231321691639751)
+atanh(z1) = (0.000000000000000000000000000000000,0.785398163397448309615660845819876)
+//]
+*/
diff --git a/src/boost/libs/multiprecision/example/constexpr_float_arithmetic_examples.cpp b/src/boost/libs/multiprecision/example/constexpr_float_arithmetic_examples.cpp
new file mode 100644
index 00000000..ef2c4579
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/constexpr_float_arithmetic_examples.cpp
@@ -0,0 +1,374 @@
+//  (C) Copyright John Maddock 2019.
+//  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 <iostream>
+#include <boost/math/constants/constants.hpp>
+#ifdef BOOST_HAS_FLOAT128
+#include <boost/multiprecision/float128.hpp>
+#endif
+
+//[constexpr_circle
+
+template <class T>
+inline constexpr T circumference(T radius)
+{
+   return 2 * boost::math::constants::pi<T>() * radius;
+}
+
+template <class T>
+inline constexpr T area(T radius)
+{
+   return boost::math::constants::pi<T>() * radius * radius;
+}
+//]
+
+template <class T, unsigned Order>
+struct const_polynomial
+{
+ public:
+   T data[Order + 1];
+
+ public:
+   constexpr const_polynomial(T val = 0) : data{val} {}
+   constexpr const_polynomial(const std::initializer_list<T>& init) : data{}
+   {
+      if (init.size() > Order + 1)
+         throw std::range_error("Too many initializers in list");
+      for (unsigned i = 0; i < init.size(); ++i)
+         data[i] = init.begin()[i];
+   }
+   constexpr T& operator[](std::size_t N)
+   {
+      return data[N];
+   }
+   constexpr const T& operator[](std::size_t N) const
+   {
+      return data[N];
+   }
+   template <class U>
+   constexpr T operator()(U val)const
+   {
+      T result = data[Order];
+      for (unsigned i = Order; i > 0; --i)
+      {
+         result *= val;
+         result += data[i - 1];
+      }
+      return result;
+   }
+   constexpr const_polynomial<T, Order - 1> derivative() const
+   {
+      const_polynomial<T, Order - 1> result;
+      for (unsigned i = 1; i <= Order; ++i)
+      {
+         result[i - 1] = (*this)[i] * i;
+      }
+      return result;
+   }
+   constexpr const_polynomial operator-()
+   {
+      const_polynomial t(*this);
+      for (unsigned i = 0; i <= Order; ++i)
+         t[i] = -t[i];
+      return t;
+   }
+   template <class U>
+   constexpr const_polynomial& operator*=(U val)
+   {
+      for (unsigned i = 0; i <= Order; ++i)
+         data[i] = data[i] * val;
+      return *this;
+   }
+   template <class U>
+   constexpr const_polynomial& operator/=(U val)
+   {
+      for (unsigned i = 0; i <= Order; ++i)
+         data[i] = data[i] / val;
+      return *this;
+   }
+   template <class U>
+   constexpr const_polynomial& operator+=(U val)
+   {
+      data[0] += val;
+      return *this;
+   }
+   template <class U>
+   constexpr const_polynomial& operator-=(U val)
+   {
+      data[0] -= val;
+      return *this;
+   }
+};
+
+template <class T, unsigned Order1, unsigned Order2>
+inline constexpr const_polynomial<T, (Order1 > Order2 ? Order1 : Order2)> operator+(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b)
+{
+   if
+      constexpr(Order1 > Order2)
+      {
+         const_polynomial<T, Order1> result(a);
+         for (unsigned i = 0; i <= Order2; ++i)
+            result[i] += b[i];
+         return result;
+      }
+   else
+   {
+      const_polynomial<T, Order2> result(b);
+      for (unsigned i = 0; i <= Order1; ++i)
+         result[i] += a[i];
+      return result;
+   }
+}
+template <class T, unsigned Order1, unsigned Order2>
+inline constexpr const_polynomial<T, (Order1 > Order2 ? Order1 : Order2)> operator-(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b)
+{
+   if
+      constexpr(Order1 > Order2)
+      {
+         const_polynomial<T, Order1> result(a);
+         for (unsigned i = 0; i <= Order2; ++i)
+            result[i] -= b[i];
+         return result;
+      }
+   else
+   {
+      const_polynomial<T, Order2> result(b);
+      for (unsigned i = 0; i <= Order1; ++i)
+         result[i] = a[i] - b[i];
+      return result;
+   }
+}
+template <class T, unsigned Order1, unsigned Order2>
+inline constexpr const_polynomial<T, Order1 + Order2> operator*(const const_polynomial<T, Order1>& a, const const_polynomial<T, Order2>& b)
+{
+   const_polynomial<T, Order1 + Order2> result;
+   for (unsigned i = 0; i <= Order1; ++i)
+   {
+      for (unsigned j = 0; j <= Order2; ++j)
+      {
+         result[i + j] += a[i] * b[j];
+      }
+   }
+   return result;
+}
+template <class T, unsigned Order, class U>
+inline constexpr const_polynomial<T, Order> operator*(const const_polynomial<T, Order>& a, const U& b)
+{
+   const_polynomial<T, Order> result(a);
+   for (unsigned i = 0; i <= Order; ++i)
+   {
+      result[i] *= b;
+   }
+   return result;
+}
+template <class U, class T, unsigned Order>
+inline constexpr const_polynomial<T, Order> operator*(const U& b, const const_polynomial<T, Order>& a)
+{
+   const_polynomial<T, Order> result(a);
+   for (unsigned i = 0; i <= Order; ++i)
+   {
+      result[i] *= b;
+   }
+   return result;
+}
+template <class T, unsigned Order, class U>
+inline constexpr const_polynomial<T, Order> operator/(const const_polynomial<T, Order>& a, const U& b)
+{
+   const_polynomial<T, Order> result;
+   for (unsigned i = 0; i <= Order; ++i)
+   {
+      result[i] /= b;
+   }
+   return result;
+}
+
+//[hermite_example
+template <class T, unsigned Order>
+class hermite_polynomial
+{
+   const_polynomial<T, Order> m_data;
+
+ public:
+   constexpr hermite_polynomial() : m_data(hermite_polynomial<T, Order - 1>().data() * const_polynomial<T, 1>{0, 2} - hermite_polynomial<T, Order - 1>().data().derivative())
+   {
+   }
+   constexpr const const_polynomial<T, Order>& data() const
+   {
+      return m_data;
+   }
+   constexpr const T& operator[](std::size_t N)const
+   {
+      return m_data[N];
+   }
+   template <class U>
+   constexpr T operator()(U val)const
+   {
+      return m_data(val);
+   }
+};
+//]
+//[hermite_example2
+template <class T>
+class hermite_polynomial<T, 0>
+{
+   const_polynomial<T, 0> m_data;
+
+ public:
+   constexpr hermite_polynomial() : m_data{1} {}
+   constexpr const const_polynomial<T, 0>& data() const
+   {
+      return m_data;
+   }
+   constexpr const T& operator[](std::size_t N) const
+   {
+      return m_data[N];
+   }
+   template <class U>
+   constexpr T operator()(U val)
+   {
+      return m_data(val);
+   }
+};
+
+template <class T>
+class hermite_polynomial<T, 1>
+{
+   const_polynomial<T, 1> m_data;
+
+ public:
+   constexpr hermite_polynomial() : m_data{0, 2} {}
+   constexpr const const_polynomial<T, 1>& data() const
+   {
+      return m_data;
+   }
+   constexpr const T& operator[](std::size_t N) const
+   {
+      return m_data[N];
+   }
+   template <class U>
+   constexpr T operator()(U val)
+   {
+      return m_data(val);
+   }
+};
+//]
+
+void test_double()
+{
+   constexpr double radius = 2.25;
+   constexpr double c      = circumference(radius);
+   constexpr double a      = area(radius);
+
+   std::cout << "Circumference = " << c << std::endl;
+   std::cout << "Area = " << a << std::endl;
+
+   constexpr const_polynomial<double, 2> pa = {3, 4};
+   constexpr const_polynomial<double, 2> pb = {5, 6};
+   static_assert(pa[0] == 3);
+   static_assert(pa[1] == 4);
+   constexpr auto pc = pa * 2;
+   static_assert(pc[0] == 6);
+   static_assert(pc[1] == 8);
+   constexpr auto pd = 3 * pa;
+   static_assert(pd[0] == 3 * 3);
+   static_assert(pd[1] == 4 * 3);
+   constexpr auto pe = pa + pb;
+   static_assert(pe[0] == 3 + 5);
+   static_assert(pe[1] == 4 + 6);
+   constexpr auto pf = pa - pb;
+   static_assert(pf[0] == 3 - 5);
+   static_assert(pf[1] == 4 - 6);
+   constexpr auto pg = pa * pb;
+   static_assert(pg[0] == 15);
+   static_assert(pg[1] == 38);
+   static_assert(pg[2] == 24);
+
+   constexpr hermite_polynomial<double, 2> h1;
+   static_assert(h1[0] == -2);
+   static_assert(h1[1] == 0);
+   static_assert(h1[2] == 4);
+
+   constexpr hermite_polynomial<double, 3> h3;
+   static_assert(h3[0] == 0);
+   static_assert(h3[1] == -12);
+   static_assert(h3[2] == 0);
+   static_assert(h3[3] == 8);
+
+   constexpr hermite_polynomial<double, 9> h9;
+   static_assert(h9[0] == 0);
+   static_assert(h9[1] == 30240);
+   static_assert(h9[2] == 0);
+   static_assert(h9[3] == -80640);
+   static_assert(h9[4] == 0);
+   static_assert(h9[5] == 48384);
+   static_assert(h9[6] == 0);
+   static_assert(h9[7] == -9216);
+   static_assert(h9[8] == 0);
+   static_assert(h9[9] == 512);
+
+   static_assert(h9(0.5) == 6481);
+
+}
+
+void test_float128()
+{
+#ifdef BOOST_HAS_FLOAT128
+   //[constexpr_circle_usage
+
+   using boost::multiprecision::float128;
+
+   constexpr float128 radius = 2.25;
+   constexpr float128 c      = circumference(radius);
+   constexpr float128 a      = area(radius);
+
+   std::cout << "Circumference = " << c << std::endl;
+   std::cout << "Area = " << a << std::endl;
+
+   //]
+
+   constexpr hermite_polynomial<float128, 2> h1;
+   static_assert(h1[0] == -2);
+   static_assert(h1[1] == 0);
+   static_assert(h1[2] == 4);
+
+   constexpr hermite_polynomial<float128, 3> h3;
+   static_assert(h3[0] == 0);
+   static_assert(h3[1] == -12);
+   static_assert(h3[2] == 0);
+   static_assert(h3[3] == 8);
+
+   //[hermite_example3
+   constexpr hermite_polynomial<float128, 9> h9;
+   //
+   // Verify that the polynomial's coefficients match the known values:
+   //
+   static_assert(h9[0] == 0);
+   static_assert(h9[1] == 30240);
+   static_assert(h9[2] == 0);
+   static_assert(h9[3] == -80640);
+   static_assert(h9[4] == 0);
+   static_assert(h9[5] == 48384);
+   static_assert(h9[6] == 0);
+   static_assert(h9[7] == -9216);
+   static_assert(h9[8] == 0);
+   static_assert(h9[9] == 512);
+   //
+   // Define an abscissa value to evaluate at:
+   //
+   constexpr float128 abscissa(0.5);
+   //
+   // Evaluate H_9(0.5) using all constexpr arithmetic:
+   //
+   static_assert(h9(abscissa) == 6481);
+   //]
+#endif
+}
+
+int main()
+{
+   test_double();
+   test_float128();
+   std::cout << "Done!" << std::endl;
+}
diff --git a/src/boost/libs/multiprecision/example/cpp_bin_float_import_export.cpp b/src/boost/libs/multiprecision/example/cpp_bin_float_import_export.cpp
new file mode 100644
index 00000000..3c9bc143
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/cpp_bin_float_import_export.cpp
@@ -0,0 +1,49 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2015 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <iterator>
+
+// Contains Quickbook snippets in comments.
+
+//[IE2
+
+/*`
+Importing or exporting cpp_bin_float is similar, but we must procede via an intermediate integer:
+*/
+/*=
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <iterator>
+*/
+
+int main()
+{
+   using boost::multiprecision::cpp_bin_float_100;
+   using boost::multiprecision::cpp_int;
+   // Create a cpp_bin_float to import/export:
+   cpp_bin_float_100 f(1);
+   f /= 3;
+   // export into 8-bit unsigned values, most significant bit first:
+   std::vector<unsigned char> v;
+   export_bits(cpp_int(f.backend().bits()), std::back_inserter(v), 8);
+   // Grab the exponent as well:
+   int e = f.backend().exponent();
+   // Import back again, and check for equality, we have to procede via
+   // an intermediate integer:
+   cpp_int i;
+   import_bits(i, v.begin(), v.end());
+   cpp_bin_float_100 g(i);
+   g.backend().exponent() = e;
+   BOOST_ASSERT(f == g);
+}
+
+//]
+
diff --git a/src/boost/libs/multiprecision/example/cpp_bin_float_snips.cpp b/src/boost/libs/multiprecision/example/cpp_bin_float_snips.cpp
new file mode 100644
index 00000000..ac378e05
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/cpp_bin_float_snips.cpp
@@ -0,0 +1,32 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2013 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+//[cpp_bin_float_eg
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/math/special_functions/gamma.hpp>
+#include <iostream>
+
+int main()
+{
+   using namespace boost::multiprecision;
+
+   // Operations at fixed precision and full numeric_limits support:
+   cpp_bin_float_100 b = 2;
+   std::cout << std::numeric_limits<cpp_bin_float_100>::digits << std::endl;
+   std::cout << std::numeric_limits<cpp_bin_float_100>::digits10 << std::endl;
+   // We can use any C++ std lib function, lets print all the digits as well:
+   std::cout << std::setprecision(std::numeric_limits<cpp_bin_float_100>::max_digits10) 
+      << log(b) << std::endl; // print log(2)
+   // We can also use any function from Boost.Math:
+   std::cout << boost::math::tgamma(b) << std::endl;
+   // These even work when the argument is an expression template:
+   std::cout << boost::math::tgamma(b * b) << std::endl;
+   // And since we have an extended exponent range we can generate some really large 
+   // numbers here (4.0238726007709377354370243e+2564):
+   std::cout << boost::math::tgamma(cpp_bin_float_100(1000)) << std::endl;
+   return 0;
+}
+
+//]
diff --git a/src/boost/libs/multiprecision/example/cpp_complex_examples.cpp b/src/boost/libs/multiprecision/example/cpp_complex_examples.cpp
new file mode 100644
index 00000000..7d00cc73
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/cpp_complex_examples.cpp
@@ -0,0 +1,120 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2018 Nick Thompson. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+/*`This example demonstrates the usage of the MPC backend for multiprecision complex numbers.
+In the following, we will show how using MPC backend allows for the same operations as the C++ standard library complex numbers.
+*/
+
+//[cpp_complex_eg
+#include <iostream>
+#include <complex>
+#include <boost/multiprecision/cpp_complex.hpp>
+
+template<class Complex>
+void complex_number_examples()
+{
+    Complex z1{0, 1};
+    std::cout << std::setprecision(std::numeric_limits<typename Complex::value_type>::digits10);
+    std::cout << std::scientific << std::fixed;
+    std::cout << "Print a complex number: " << z1 << std::endl;
+    std::cout << "Square it             : " << z1*z1 << std::endl;
+    std::cout << "Real part             : " << z1.real() << " = " << real(z1) << std::endl;
+    std::cout << "Imaginary part        : " << z1.imag() << " = " << imag(z1) << std::endl;
+    using std::abs;
+    std::cout << "Absolute value        : " << abs(z1) << std::endl;
+    std::cout << "Argument              : " << arg(z1) << std::endl;
+    std::cout << "Norm                  : " << norm(z1) << std::endl;
+    std::cout << "Complex conjugate     : " << conj(z1) << std::endl;
+    std::cout << "Projection onto Riemann sphere: " <<  proj(z1) << std::endl;
+    typename Complex::value_type r = 1;
+    typename Complex::value_type theta = 0.8;
+    using std::polar;
+    std::cout << "Polar coordinates (phase = 0)    : " << polar(r) << std::endl;
+    std::cout << "Polar coordinates (phase !=0)    : " << polar(r, theta) << std::endl;
+
+    std::cout << "\nElementary special functions:\n";
+    using std::exp;
+    std::cout << "exp(z1) = " << exp(z1) << std::endl;
+    using std::log;
+    std::cout << "log(z1) = " << log(z1) << std::endl;
+    using std::log10;
+    std::cout << "log10(z1) = " << log10(z1) << std::endl;
+    using std::pow;
+    std::cout << "pow(z1, z1) = " << pow(z1, z1) << std::endl;
+    using std::sqrt;
+    std::cout << "Take its square root  : " << sqrt(z1) << std::endl;
+    using std::sin;
+    std::cout << "sin(z1) = " << sin(z1) << std::endl;
+    using std::cos;
+    std::cout << "cos(z1) = " << cos(z1) << std::endl;
+    using std::tan;
+    std::cout << "tan(z1) = " << tan(z1) << std::endl;
+    using std::asin;
+    std::cout << "asin(z1) = " << asin(z1) << std::endl;
+    using std::acos;
+    std::cout << "acos(z1) = " << acos(z1) << std::endl;
+    using std::atan;
+    std::cout << "atan(z1) = " << atan(z1) << std::endl;
+    using std::sinh;
+    std::cout << "sinh(z1) = " << sinh(z1) << std::endl;
+    using std::cosh;
+    std::cout << "cosh(z1) = " << cosh(z1) << std::endl;
+    using std::tanh;
+    std::cout << "tanh(z1) = " << tanh(z1) << std::endl;
+    using std::asinh;
+    std::cout << "asinh(z1) = " << asinh(z1) << std::endl;
+    using std::acosh;
+    std::cout << "acosh(z1) = " << acosh(z1) << std::endl;
+    using std::atanh;
+    std::cout << "atanh(z1) = " << atanh(z1) << std::endl;
+}
+
+int main()
+{
+    std::cout << "First, some operations we usually perform with std::complex:\n";
+    complex_number_examples<std::complex<double>>();
+    std::cout << "\nNow the same operations performed using quad precision complex numbers:\n";
+    complex_number_examples<boost::multiprecision::cpp_complex_quad>();
+
+    return 0;
+}
+//]
+
+/*
+
+//[cpp_complex_out
+
+Print a complex number: (0.000000000000000000000000000000000,1.000000000000000000000000000000000)
+Square it             : -1.000000000000000000000000000000000
+Real part             : 0.000000000000000000000000000000000 = 0.000000000000000000000000000000000
+Imaginary part        : 1.000000000000000000000000000000000 = 1.000000000000000000000000000000000
+Absolute value        : 1.000000000000000000000000000000000
+Argument              : 1.570796326794896619231321691639751
+Norm                  : 1.000000000000000000000000000000000
+Complex conjugate     : (0.000000000000000000000000000000000,-1.000000000000000000000000000000000)
+Projection onto Riemann sphere: (0.000000000000000000000000000000000,1.000000000000000000000000000000000)
+Polar coordinates (phase = 0)    : 1.000000000000000000000000000000000
+Polar coordinates (phase !=0)    : (0.696706709347165389063740022772448,0.717356090899522792567167815703377)
+
+Elementary special functions:
+exp(z1) = (0.540302305868139717400936607442977,0.841470984807896506652502321630299)
+log(z1) = (0.000000000000000000000000000000000,1.570796326794896619231321691639751)
+log10(z1) = (0.000000000000000000000000000000000,0.682188176920920673742891812715678)
+pow(z1, z1) = 0.207879576350761908546955619834979
+Take its square root  : (0.707106781186547524400844362104849,0.707106781186547524400844362104849)
+sin(z1) = (0.000000000000000000000000000000000,1.175201193643801456882381850595601)
+cos(z1) = 1.543080634815243778477905620757062
+tan(z1) = (0.000000000000000000000000000000000,0.761594155955764888119458282604794)
+asin(z1) = (0.000000000000000000000000000000000,0.881373587019543025232609324979793)
+acos(z1) = (1.570796326794896619231321691639751,-0.881373587019543025232609324979793)
+atan(z1) = (0.000000000000000000000000000000000,inf)
+sinh(z1) = (0.000000000000000000000000000000000,0.841470984807896506652502321630299)
+cosh(z1) = 0.540302305868139717400936607442977
+tanh(z1) = (0.000000000000000000000000000000000,1.557407724654902230506974807458360)
+asinh(z1) = (0.000000000000000000000000000000000,1.570796326794896619231321691639751)
+acosh(z1) = (0.881373587019543025232609324979792,1.570796326794896619231321691639751)
+atanh(z1) = (0.000000000000000000000000000000000,0.785398163397448309615660845819876)
+//]
+*/
diff --git a/src/boost/libs/multiprecision/example/cpp_dec_float_snips.cpp b/src/boost/libs/multiprecision/example/cpp_dec_float_snips.cpp
new file mode 100644
index 00000000..4bbd89f9
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/cpp_dec_float_snips.cpp
@@ -0,0 +1,33 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2011 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+   //[cpp_dec_float_eg
+#include <boost/multiprecision/cpp_dec_float.hpp>
+#include <boost/math/special_functions/gamma.hpp>
+#include <iostream>
+
+int main()
+{
+   using namespace boost::multiprecision;
+
+   // Operations at fixed precision and full numeric_limits support:
+   cpp_dec_float_100 b = 2;
+   std::cout << std::numeric_limits<cpp_dec_float_100>::digits << std::endl;
+   // Note that digits10 is the same as digits, since we're base 10! :
+   std::cout << std::numeric_limits<cpp_dec_float_100>::digits10 << std::endl;
+   // We can use any C++ std lib function, lets print all the digits as well:
+   std::cout << std::setprecision(std::numeric_limits<cpp_dec_float_100>::max_digits10) 
+      << log(b) << std::endl; // print log(2)
+   // We can also use any function from Boost.Math:
+   std::cout << boost::math::tgamma(b) << std::endl;
+   // These even work when the argument is an expression template:
+   std::cout << boost::math::tgamma(b * b) << std::endl;
+   // And since we have an extended exponent range we can generate some really large 
+   // numbers here (4.0238726007709377354370243e+2564):
+   std::cout << boost::math::tgamma(cpp_dec_float_100(1000)) << std::endl;
+   return 0;
+}
+//]
+
diff --git a/src/boost/libs/multiprecision/example/cpp_int_import_export.cpp b/src/boost/libs/multiprecision/example/cpp_int_import_export.cpp
new file mode 100644
index 00000000..0e907118
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/cpp_int_import_export.cpp
@@ -0,0 +1,43 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2015 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/multiprecision/cpp_int.hpp>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <iterator>
+
+//[IE1
+
+/*`
+In this simple example, we'll import/export the bits of a cpp_int
+to a vector of 8-bit unsigned values:
+*/
+/*=
+#include <boost/multiprecision/cpp_int.hpp>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <iterator>
+*/
+
+int main()
+{
+   using boost::multiprecision::cpp_int;
+   // Create a cpp_int with just a couple of bits set:
+   cpp_int i;
+   bit_set(i, 5000); // set the 5000'th bit
+   bit_set(i, 200);
+   bit_set(i, 50);
+   // export into 8-bit unsigned values, most significant bit first:
+   std::vector<unsigned char> v;
+   export_bits(i, std::back_inserter(v), 8);
+   // import back again, and check for equality:
+   cpp_int j;
+   import_bits(j, v.begin(), v.end());
+   BOOST_ASSERT(i == j);
+}
+
+//]
diff --git a/src/boost/libs/multiprecision/example/cpp_int_snips.cpp b/src/boost/libs/multiprecision/example/cpp_int_snips.cpp
new file mode 100644
index 00000000..4175db3f
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/cpp_int_snips.cpp
@@ -0,0 +1,75 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2011 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/multiprecision/cpp_int.hpp>
+#include <iostream>
+
+void t1()
+{
+//[cpp_int_eg
+//=#include <boost/multiprecision/cpp_int.hpp>
+//=#include <iostream>
+//=
+//=int main()
+//={
+   using namespace boost::multiprecision;
+
+   int128_t v = 1;
+
+   // Do some fixed precision arithmetic:
+   for(unsigned i = 1; i <= 20; ++i)
+      v *= i;
+
+   std::cout << v << std::endl; // prints 2432902008176640000 (i.e. 20!)
+
+   // Repeat at arbitrary precision:
+   cpp_int u = 1;
+   for(unsigned i = 1; i <= 100; ++i)
+      u *= i;
+   
+   // prints 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000 (i.e. 100!)
+   std::cout << u << std::endl;
+
+//=   return 0;
+//=}
+//]
+}
+
+void t3()
+{
+//[cpp_rational_eg
+//=#include <boost/multiprecision/cpp_int.hpp>
+//=#include <iostream>
+//=
+//=int main()
+//={
+   using namespace boost::multiprecision;
+
+   cpp_rational v = 1;
+
+   // Do some arithmetic:
+   for(unsigned i = 1; i <= 1000; ++i)
+      v *= i;
+   v /= 10;
+
+   std::cout << v << std::endl; // prints 1000! / 10
+   std::cout << numerator(v) << std::endl;
+   std::cout << denominator(v) << std::endl;
+
+   cpp_rational w(2, 3);  // component wise constructor
+   std::cout << w << std::endl; // prints 2/3
+//=   return 0;
+//=}
+//]
+}
+
+
+int main()
+{
+   t1();
+   t3();
+   return 0;
+}
+
diff --git a/src/boost/libs/multiprecision/example/debug_adaptor_snips.cpp b/src/boost/libs/multiprecision/example/debug_adaptor_snips.cpp
new file mode 100644
index 00000000..2d0efe9a
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/debug_adaptor_snips.cpp
@@ -0,0 +1,40 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2011 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/multiprecision/cpp_dec_float.hpp>
+#include <boost/multiprecision/debug_adaptor.hpp>
+#include <iostream>
+
+void t1()
+{
+   //[debug_adaptor_eg
+   //=#include <boost/multiprecision/debug_adaptor.hpp>
+   //=#include <boost/multiprecision/cpp_dec_float.hpp>
+
+   using namespace boost::multiprecision;
+
+   typedef number<debug_adaptor<cpp_dec_float<50> > > fp_type;
+
+   fp_type denom = 1;
+   fp_type sum = 1;
+
+   for(unsigned i = 2; i < 50; ++i)
+   {
+      denom *= i;
+      sum += 1 / denom;
+   }
+
+   std::cout << std::setprecision(std::numeric_limits<fp_type>::digits) << sum << std::endl;
+   //]
+}
+
+int main()
+{
+   t1();
+   return 0;
+}
+
+
+
diff --git a/src/boost/libs/multiprecision/example/eigen_example.cpp b/src/boost/libs/multiprecision/example/eigen_example.cpp
new file mode 100644
index 00000000..a70e3fbc
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/eigen_example.cpp
@@ -0,0 +1,52 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2018 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+//[eigen_eg
+#include <iostream>
+#include <boost/multiprecision/cpp_complex.hpp>
+#include <boost/multiprecision/eigen.hpp>
+#include <Eigen/Dense>
+
+int main()
+{
+   using namespace Eigen;
+   typedef boost::multiprecision::cpp_complex_quad complex_type;
+   //
+   // We want to solve Ax = b for x,
+   // define A and b first:
+   //
+   Matrix<complex_type, 2, 2> A, b;
+   A << complex_type(2, 3), complex_type(-1, -2), complex_type(-1, -4), complex_type(3, 6);
+   b << 1, 2, 3, 1;
+   std::cout << "Here is the matrix A:\n" << A << std::endl;
+   std::cout << "Here is the right hand side b:\n" << b << std::endl;
+   //
+   // Solve for x:
+   //
+   Matrix<complex_type, 2, 2> x = A.fullPivHouseholderQr().solve(b);
+   std::cout << "The solution is:\n" << x << std::endl;
+   //
+   // Compute the error in the solution by using the norms of Ax - b and b:
+   //
+   complex_type::value_type relative_error = (A*x - b).norm() / b.norm();
+   std::cout << "The relative error is: " << relative_error << std::endl;
+   return 0;
+}
+//]
+
+/*
+//[eigen_out
+Here is the matrix A:
+(2,3) (-1,-2)
+(-1,-4)   (3,6)
+Here is the right hand side b:
+1 2
+3 1
+The solution is:
+(0.6,-0.6)   (0.7,-0.7)
+(0.64,-0.68) (0.58,-0.46)
+The relative error is: 2.63132e-34
+//]
+*/
diff --git a/src/boost/libs/multiprecision/example/float128_snips.cpp b/src/boost/libs/multiprecision/example/float128_snips.cpp
new file mode 100644
index 00000000..5d059252
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/float128_snips.cpp
@@ -0,0 +1,42 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2013 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+//[float128_eg
+#include <boost/multiprecision/float128.hpp>
+#include <boost/math/special_functions/gamma.hpp>
+#include <iostream>
+
+int main()
+{
+   using namespace boost::multiprecision;
+
+   // Operations at 128-bit precision and full numeric_limits support:
+   float128 b = 2;
+   // There are 113-bits of precision:
+   std::cout << std::numeric_limits<float128>::digits << std::endl;
+   // Or 34 decimal places:
+   std::cout << std::numeric_limits<float128>::digits10 << std::endl;
+   // We can use any C++ std lib function, lets print all the digits as well:
+   std::cout << std::setprecision(std::numeric_limits<float128>::max_digits10) 
+      << log(b) << std::endl; // print log(2) = 0.693147180559945309417232121458176575
+   // We can also use any function from Boost.Math:
+   std::cout << boost::math::tgamma(b) << std::endl;
+   // And since we have an extended exponent range we can generate some really large 
+   // numbers here (4.02387260077093773543702433923004111e+2564):
+   std::cout << boost::math::tgamma(float128(1000)) << std::endl;
+   //
+   // We can declare constants using GCC or Intel's native types, and the Q suffix,
+   // these can be declared constexpr if required:
+   /*<-*/
+#ifndef BOOST_NO_CXX11_CONSTEXPR
+   /*->*/
+   constexpr float128 pi = 3.1415926535897932384626433832795028841971693993751058Q;
+   /*<-*/
+#endif
+   /*->*/
+   return 0;
+}
+//]
+
diff --git a/src/boost/libs/multiprecision/example/floating_point_examples.cpp b/src/boost/libs/multiprecision/example/floating_point_examples.cpp
new file mode 100644
index 00000000..f34e278f
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/floating_point_examples.cpp
@@ -0,0 +1,697 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2012 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/math/constants/constants.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+#include <boost/math/special_functions/gamma.hpp>
+#include <boost/math/special_functions/bessel.hpp>
+#include <iostream>
+#include <iomanip>
+
+#if !defined(BOOST_NO_CXX11_HDR_ARRAY) && !defined(BOOST_NO_CXX11_LAMBDAS) && !(defined(CI_SUPPRESS_KNOWN_ISSUES) && defined(__GNUC__) && defined(_WIN32))
+
+#include <array>
+
+//[AOS1
+
+/*`Generic numeric programming employs templates to use the same code for different
+floating-point types and functions. Consider the area of a circle a of radius r, given by
+
+[:['a = [pi] * r[super 2]]]
+
+The area of a circle can be computed in generic programming using Boost.Math
+for the constant [pi] as shown below:
+
+*/
+
+//=#include <boost/math/constants/constants.hpp>
+
+template<typename T>
+inline T area_of_a_circle(T r)
+{
+   using boost::math::constants::pi;
+   return pi<T>() * r * r;
+}
+
+/*`
+It is possible to use `area_of_a_circle()` with built-in floating-point types as
+well as floating-point types from Boost.Multiprecision. In particular, consider a
+system with 4-byte single-precision float, 8-byte double-precision double and also the
+`cpp_dec_float_50` data type from Boost.Multiprecision with 50 decimal digits
+of precision.
+
+We can compute and print the approximate area of a circle with radius 123/100 for
+`float`, `double` and `cpp_dec_float_50` with the program below.
+
+*/
+
+//]
+
+//[AOS3
+
+/*`In the next example we'll look at calling both standard library and Boost.Math functions from within generic code.
+We'll also show how to cope with template arguments which are expression-templates rather than number types.*/
+
+//]
+
+//[JEL
+
+/*`
+In this example we'll show several implementations of the 
+[@http://mathworld.wolfram.com/LambdaFunction.html Jahnke and Emden Lambda function], 
+each implementation a little more sophisticated than the last.
+
+The Jahnke-Emden Lambda function is defined by the equation:
+
+[:['JahnkeEmden(v, z) = [Gamma](v+1) * J[sub v](z) / (z / 2)[super v]]]
+
+If we were to implement this at double precision using Boost.Math's facilities for the Gamma and Bessel
+function calls it would look like this:
+
+*/
+
+double JEL1(double v, double z)
+{
+   return boost::math::tgamma(v + 1) * boost::math::cyl_bessel_j(v, z) / std::pow(z / 2, v);
+}
+
+/*`
+Calling this function as:
+
+   std::cout << std::scientific << std::setprecision(std::numeric_limits<double>::digits10);
+   std::cout << JEL1(2.5, 0.5) << std::endl;
+
+Yields the output:
+
+[pre 9.822663964796047e-001]
+
+Now let's implement the function again, but this time using the multiprecision type
+`cpp_dec_float_50` as the argument type:
+
+*/   
+
+boost::multiprecision::cpp_dec_float_50 
+   JEL2(boost::multiprecision::cpp_dec_float_50 v, boost::multiprecision::cpp_dec_float_50 z)
+{
+   return boost::math::tgamma(v + 1) * boost::math::cyl_bessel_j(v, z) / boost::multiprecision::pow(z / 2, v);
+}
+
+/*`
+The implementation is almost the same as before, but with one key difference - we can no longer call
+`std::pow`, instead we must call the version inside the `boost::multiprecision` namespace.  In point of
+fact, we could have omitted the namespace prefix on the call to `pow` since the right overload would 
+have been found via [@http://en.wikipedia.org/wiki/Argument-dependent_name_lookup 
+argument dependent lookup] in any case.
+
+Note also that the first argument to `pow` along with the argument to `tgamma` in the above code
+are actually expression templates.  The `pow` and `tgamma` functions will handle these arguments
+just fine.
+
+Here's an example of how the function may be called:
+
+   std::cout << std::scientific << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10);
+   std::cout << JEL2(cpp_dec_float_50(2.5), cpp_dec_float_50(0.5)) << std::endl;
+
+Which outputs:
+
+[pre 9.82266396479604757017335009796882833995903762577173e-01]
+
+Now that we've seen some non-template examples, lets repeat the code again, but this time as a template
+that can be called either with a builtin type (`float`, `double` etc), or with a multiprecision type:
+
+*/
+
+template <class Float>
+Float JEL3(Float v, Float z)
+{
+   using std::pow;
+   return boost::math::tgamma(v + 1) * boost::math::cyl_bessel_j(v, z) / pow(z / 2, v);
+}
+
+/*`
+
+Once again the code is almost the same as before, but the call to `pow` has changed yet again.
+We need the call to resolve to either `std::pow` (when the argument is a builtin type), or
+to `boost::multiprecision::pow` (when the argument is a multiprecision type).  We do that by
+making the call unqualified so that versions of `pow` defined in the same namespace as type
+`Float` are found via argument dependent lookup, while the `using std::pow` directive makes
+the standard library versions visible for builtin floating point types.
+
+Let's call the function with both `double` and multiprecision arguments:
+
+   std::cout << std::scientific << std::setprecision(std::numeric_limits<double>::digits10);
+   std::cout << JEL3(2.5, 0.5) << std::endl;
+   std::cout << std::scientific << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10);
+   std::cout << JEL3(cpp_dec_float_50(2.5), cpp_dec_float_50(0.5)) << std::endl;
+
+Which outputs:
+
+[pre
+9.822663964796047e-001
+9.82266396479604757017335009796882833995903762577173e-01
+]
+
+Unfortunately there is a problem with this version: if we were to call it like this:
+
+   boost::multiprecision::cpp_dec_float_50 v(2), z(0.5);
+   JEL3(v + 0.5, z);
+
+Then we would get a long and inscrutable error message from the compiler: the problem here is that the first
+argument to `JEL3` is not a number type, but an expression template.  We could obviously add a typecast to
+fix the issue:
+
+   JEL(cpp_dec_float_50(v + 0.5), z);
+
+However, if we want the function JEL to be truly reusable, then a better solution might be preferred.
+To achieve this we can borrow some code from Boost.Math which calculates the return type of mixed-argument
+functions, here's how the new code looks now:
+
+*/
+
+template <class Float1, class Float2>
+typename boost::math::tools::promote_args<Float1, Float2>::type 
+   JEL4(Float1 v, Float2 z)
+{
+   using std::pow;
+   return boost::math::tgamma(v + 1) * boost::math::cyl_bessel_j(v, z) / pow(z / 2, v);
+}
+
+/*`
+
+As you can see the two arguments to the function are now separate template types, and
+the return type is computed using the `promote_args` metafunction from Boost.Math.
+
+Now we can call:
+
+   std::cout << std::scientific << std::setprecision(std::numeric_limits<cpp_dec_float_100>::digits10);
+   std::cout << JEL4(cpp_dec_float_100(2) + 0.5, cpp_dec_float_100(0.5)) << std::endl;
+
+And get 100 digits of output:
+
+[pre 9.8226639647960475701733500979688283399590376257717309069410413822165082248153638454147004236848917775e-01]
+
+As a bonus, we can now call the function not just with expression templates, but with other mixed types as well:
+for example `float` and `double` or `int` and `double`, and the correct return type will be computed in each case.
+
+Note that while in this case we didn't have to change the body of the function, in the general case
+any function like this which creates local variables internally would have to use `promote_args`
+to work out what type those variables should be, for example:
+
+   template <class Float1, class Float2>
+   typename boost::math::tools::promote_args<Float1, Float2>::type 
+      JEL5(Float1 v, Float2 z)
+   {
+      using std::pow;
+      typedef typename boost::math::tools::promote_args<Float1, Float2>::type variable_type;
+      variable_type t = pow(z / 2, v);
+      return boost::math::tgamma(v + 1) * boost::math::cyl_bessel_j(v, z) / t;
+   }
+
+*/
+
+//]
+
+//[ND1
+
+/*`
+In this example we'll add even more power to generic numeric programming using not only different
+floating-point types but also function objects as template parameters. Consider
+some well-known central difference rules for numerically computing the first derivative
+of a function ['f[prime](x)] with ['x [isin] [real]]:
+
+[equation floating_point_eg1]
+
+Where the difference terms ['m[sub n]] are given by:
+
+[equation floating_point_eg2]
+
+and ['dx] is the step-size of the derivative.
+
+The third formula in Equation 1 is a three-point central difference rule. It calculates
+the first derivative of ['f[prime](x)] to ['O(dx[super 6])], where ['dx] is the given step-size. 
+For example, if
+the step-size is 0.01 this derivative calculation has about 6 decimal digits of precision - 
+just about right for the 7 decimal digits of single-precision float.
+Let's make a generic template subroutine using this three-point central difference
+rule.  In particular:
+*/
+
+template<typename value_type, typename function_type>
+   value_type derivative(const value_type x, const value_type dx, function_type func)
+{
+   // Compute d/dx[func(*first)] using a three-point
+   // central difference rule of O(dx^6).
+
+   const value_type dx1 = dx;
+   const value_type dx2 = dx1 * 2;
+   const value_type dx3 = dx1 * 3;
+
+   const value_type m1 = (func(x + dx1) - func(x - dx1)) / 2;
+   const value_type m2 = (func(x + dx2) - func(x - dx2)) / 4;
+   const value_type m3 = (func(x + dx3) - func(x - dx3)) / 6;
+
+   const value_type fifteen_m1 = 15 * m1;
+   const value_type six_m2     =  6 * m2;
+   const value_type ten_dx1    = 10 * dx1;
+
+   return ((fifteen_m1 - six_m2) + m3) / ten_dx1;
+}
+
+/*`The `derivative()` template function can be used to compute the first derivative
+of any function to ['O(dx[super 6])]. For example, consider the first derivative of ['sin(x)] evaluated
+at ['x = [pi]/3]. In other words,
+
+[equation floating_point_eg3]
+
+The code below computes the derivative in Equation 3 for float, double and boost's
+multiple-precision type cpp_dec_float_50.
+*/
+
+//]
+
+//[GI1
+
+/*`
+Similar to the generic derivative example, we can calculate integrals in a similar manner:
+*/
+
+template<typename value_type, typename function_type>
+inline value_type integral(const value_type a, 
+                           const value_type b, 
+                           const value_type tol, 
+                           function_type func)
+{
+   unsigned n = 1U;
+
+   value_type h = (b - a);
+   value_type I = (func(a) + func(b)) * (h / 2);
+
+   for(unsigned k = 0U; k < 8U; k++)
+   {
+      h /= 2;
+
+      value_type sum(0);
+      for(unsigned j = 1U; j <= n; j++)
+      {
+         sum += func(a + (value_type((j * 2) - 1) * h));
+      }
+
+      const value_type I0 = I;
+      I = (I / 2) + (h * sum);
+
+      const value_type ratio     = I0 / I;
+      const value_type delta     = ratio - 1;
+      const value_type delta_abs = ((delta < 0) ? -delta : delta);
+
+      if((k > 1U) && (delta_abs < tol))
+      {
+         break;
+      }
+
+      n *= 2U;
+   }
+
+   return I;
+}
+
+/*`
+The following sample program shows how the function can be called, we begin
+by defining a function object, which when integrated should yield the Bessel J
+function:
+*/
+
+template<typename value_type>
+class cyl_bessel_j_integral_rep
+{
+public:
+   cyl_bessel_j_integral_rep(const unsigned N,
+      const value_type& X) : n(N), x(X) { }
+
+   value_type operator()(const value_type& t) const
+   {
+      // pi * Jn(x) = Int_0^pi [cos(x * sin(t) - n*t) dt]
+      return cos(x * sin(t) - (n * t));
+   }
+
+private:
+   const unsigned n;
+   const value_type x;
+};
+
+
+//]
+
+//[POLY
+
+/*`
+In this example we'll look at polynomial evaluation, this is not only an important
+use case, but it's one that `number` performs particularly well at because the
+expression templates ['completely eliminate all temporaries] from a 
+[@http://en.wikipedia.org/wiki/Horner%27s_method Horner polynomial
+evaluation scheme].
+
+The following code evaluates `sin(x)` as a polynomial, accurate to at least 64 decimal places:
+
+*/
+
+using boost::multiprecision::cpp_dec_float;
+typedef boost::multiprecision::number<cpp_dec_float<64> > mp_type;
+
+mp_type mysin(const mp_type& x)
+{
+  // Approximation of sin(x * pi/2) for -1 <= x <= 1, using an order 63 polynomial.
+  static const std::array<mp_type, 32U> coefs =
+  {{
+    mp_type("+1.5707963267948966192313216916397514420985846996875529104874722961539082031431044993140174126711"), //"),
+    mp_type("-0.64596409750624625365575656389794573337969351178927307696134454382929989411386887578263960484"), // ^3
+    mp_type("+0.07969262624616704512050554949047802252091164235106119545663865720995702920146198554317279"), // ^5
+    mp_type("-0.0046817541353186881006854639339534378594950280185010575749538605102665157913157426229824"), // ^7
+    mp_type("+0.00016044118478735982187266087016347332970280754062061156858775174056686380286868007443"), // ^9
+    mp_type("-3.598843235212085340458540018208389404888495232432127661083907575106196374913134E-6"), // ^11
+    mp_type("+5.692172921967926811775255303592184372902829756054598109818158853197797542565E-8"), // ^13
+    mp_type("-6.688035109811467232478226335783138689956270985704278659373558497256423498E-10"), // ^15
+    mp_type("+6.066935731106195667101445665327140070166203261129845646380005577490472E-12"), // ^17
+    mp_type("-4.377065467313742277184271313776319094862897030084226361576452003432E-14"), // ^19
+    mp_type("+2.571422892860473866153865950420487369167895373255729246889168337E-16"), // ^21
+    mp_type("-1.253899540535457665340073300390626396596970180355253776711660E-18"), // ^23
+    mp_type("+5.15645517658028233395375998562329055050964428219501277474E-21"), // ^25
+    mp_type("-1.812399312848887477410034071087545686586497030654642705E-23"), // ^27
+    mp_type("+5.50728578652238583570585513920522536675023562254864E-26"), // ^29
+    mp_type("-1.461148710664467988723468673933026649943084902958E-28"), // ^31
+    mp_type("+3.41405297003316172502972039913417222912445427E-31"), // ^33
+    mp_type("-7.07885550810745570069916712806856538290251E-34"), // ^35
+    mp_type("+1.31128947968267628970845439024155655665E-36"), // ^37
+    mp_type("-2.18318293181145698535113946654065918E-39"), // ^39
+    mp_type("+3.28462680978498856345937578502923E-42"), // ^41
+    mp_type("-4.48753699028101089490067137298E-45"), // ^43
+    mp_type("+5.59219884208696457859353716E-48"), // ^45
+    mp_type("-6.38214503973500471720565E-51"), // ^47
+    mp_type("+6.69528558381794452556E-54"), // ^49
+    mp_type("-6.47841373182350206E-57"), // ^51
+    mp_type("+5.800016389666445E-60"), // ^53
+    mp_type("-4.818507347289E-63"), // ^55
+    mp_type("+3.724683686E-66"), // ^57
+    mp_type("-2.6856479E-69"), // ^59
+    mp_type("+1.81046E-72"), // ^61
+    mp_type("-1.133E-75"), // ^63
+  }};
+
+  const mp_type v = x * 2 / boost::math::constants::pi<mp_type>();
+  const mp_type x2 = (v * v);
+  //
+  // Polynomial evaluation follows, if mp_type allocates memory then
+  // just one such allocation occurs - to initialize the variable "sum" -
+  // and no temporaries are created at all.
+  //
+  const mp_type sum = (((((((((((((((((((((((((((((((     + coefs[31U]
+                                                     * x2 + coefs[30U])
+                                                     * x2 + coefs[29U])
+                                                     * x2 + coefs[28U])
+                                                     * x2 + coefs[27U])
+                                                     * x2 + coefs[26U])
+                                                     * x2 + coefs[25U])
+                                                     * x2 + coefs[24U])
+                                                     * x2 + coefs[23U])
+                                                     * x2 + coefs[22U])
+                                                     * x2 + coefs[21U])
+                                                     * x2 + coefs[20U])
+                                                     * x2 + coefs[19U])
+                                                     * x2 + coefs[18U])
+                                                     * x2 + coefs[17U])
+                                                     * x2 + coefs[16U])
+                                                     * x2 + coefs[15U])
+                                                     * x2 + coefs[14U])
+                                                     * x2 + coefs[13U])
+                                                     * x2 + coefs[12U])
+                                                     * x2 + coefs[11U])
+                                                     * x2 + coefs[10U])
+                                                     * x2 + coefs[9U])
+                                                     * x2 + coefs[8U])
+                                                     * x2 + coefs[7U])
+                                                     * x2 + coefs[6U])
+                                                     * x2 + coefs[5U])
+                                                     * x2 + coefs[4U])
+                                                     * x2 + coefs[3U])
+                                                     * x2 + coefs[2U])
+                                                     * x2 + coefs[1U])
+                                                     * x2 + coefs[0U])
+                                                     * v;
+
+  return sum;
+}
+
+/*`
+Calling the function like so:
+
+   mp_type pid4 = boost::math::constants::pi<mp_type>() / 4;
+   std::cout << std::setprecision(std::numeric_limits< ::mp_type>::digits10) << std::scientific;
+   std::cout << mysin(pid4) << std::endl;
+
+Yields the expected output:
+
+[pre 7.0710678118654752440084436210484903928483593768847403658833986900e-01]
+
+*/
+
+//]
+
+
+int main()
+{
+   using namespace boost::multiprecision;
+   std::cout << std::scientific << std::setprecision(std::numeric_limits<double>::digits10);
+   std::cout << JEL1(2.5, 0.5) << std::endl;
+   std::cout << std::scientific << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10);
+   std::cout << JEL2(cpp_dec_float_50(2.5), cpp_dec_float_50(0.5)) << std::endl;
+   std::cout << std::scientific << std::setprecision(std::numeric_limits<double>::digits10);
+   std::cout << JEL3(2.5, 0.5) << std::endl;
+   std::cout << std::scientific << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10);
+   std::cout << JEL3(cpp_dec_float_50(2.5), cpp_dec_float_50(0.5)) << std::endl;
+   std::cout << std::scientific << std::setprecision(std::numeric_limits<cpp_dec_float_100>::digits10);
+   std::cout << JEL4(cpp_dec_float_100(2) + 0.5, cpp_dec_float_100(0.5)) << std::endl;
+
+   //[AOS2
+
+/*=#include <iostream>
+#include <iomanip>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+
+using boost::multiprecision::cpp_dec_float_50;
+
+int main(int, char**)
+{*/
+   const float r_f(float(123) / 100);
+   const float a_f = area_of_a_circle(r_f);
+
+   const double r_d(double(123) / 100);
+   const double a_d = area_of_a_circle(r_d);
+
+   const cpp_dec_float_50 r_mp(cpp_dec_float_50(123) / 100);
+   const cpp_dec_float_50 a_mp = area_of_a_circle(r_mp);
+
+   // 4.75292
+   std::cout
+      << std::setprecision(std::numeric_limits<float>::digits10)
+      << a_f
+      << std::endl;
+
+   // 4.752915525616
+   std::cout
+      << std::setprecision(std::numeric_limits<double>::digits10)
+      << a_d
+      << std::endl;
+
+   // 4.7529155256159981904701331745635599135018975843146
+   std::cout
+      << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10)
+      << a_mp
+      << std::endl;
+/*=}*/
+
+   //]
+
+   //[ND2
+/*=
+#include <iostream>
+#include <iomanip>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+#include <boost/math/constants/constants.hpp>
+
+
+int main(int, char**)
+{*/
+   using boost::math::constants::pi;
+   using boost::multiprecision::cpp_dec_float_50;
+   //
+   // We'll pass a function pointer for the function object passed to derivative,
+   // the typecast is needed to select the correct overload of std::sin:
+   //
+   const float d_f = derivative(
+      pi<float>() / 3,
+      0.01F,
+      static_cast<float(*)(float)>(std::sin)
+   );
+
+   const double d_d = derivative(
+      pi<double>() / 3,
+      0.001,
+      static_cast<double(*)(double)>(std::sin)
+      );
+   //
+   // In the cpp_dec_float_50 case, the sin function is multiply overloaded
+   // to handle expression templates etc.  As a result it's hard to take its
+   // address without knowing about its implementation details.  We'll use a 
+   // C++11 lambda expression to capture the call.
+   // We also need a typecast on the first argument so we don't accidentally pass
+   // an expression template to a template function:
+   //
+   const cpp_dec_float_50 d_mp = derivative(
+      cpp_dec_float_50(pi<cpp_dec_float_50>() / 3),
+      cpp_dec_float_50(1.0E-9),
+      [](const cpp_dec_float_50& x) -> cpp_dec_float_50
+      {
+         return sin(x);
+      }
+      );
+
+   // 5.000029e-001
+   std::cout
+      << std::setprecision(std::numeric_limits<float>::digits10)
+      << d_f
+      << std::endl;
+
+   // 4.999999999998876e-001
+   std::cout
+      << std::setprecision(std::numeric_limits<double>::digits10)
+      << d_d
+      << std::endl;
+
+   // 4.99999999999999999999999999999999999999999999999999e-01
+   std::cout
+      << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10)
+      << d_mp
+      << std::endl;
+//=}
+
+   /*`
+   The expected value of the derivative is 0.5. This central difference rule in this
+   example is ill-conditioned, meaning it suffers from slight loss of precision. With that
+   in mind, the results agree with the expected value of 0.5.*/
+
+   //]
+
+   //[ND3
+
+   /*`
+   We can take this a step further and use our derivative function to compute
+   a partial derivative.  For example if we take the incomplete gamma function
+   ['P(a, z)], and take the derivative with respect to /z/ at /(2,2)/ then we
+   can calculate the result as shown below, for good measure we'll compare with
+   the "correct" result obtained from a call to ['gamma_p_derivative], the results
+   agree to approximately 44 digits:
+   */
+
+   cpp_dec_float_50 gd = derivative(
+      cpp_dec_float_50(2),
+      cpp_dec_float_50(1.0E-9),
+      [](const cpp_dec_float_50& x) ->cpp_dec_float_50
+      {
+         return boost::math::gamma_p(2, x);
+      }
+   );
+   // 2.70670566473225383787998989944968806815263091819151e-01
+   std::cout
+      << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10)
+      << gd
+      << std::endl;
+   // 2.70670566473225383787998989944968806815253190143120e-01
+   std::cout << boost::math::gamma_p_derivative(cpp_dec_float_50(2), cpp_dec_float_50(2)) << std::endl;
+   //]
+
+   //[GI2
+
+   /* The function can now be called as follows: */
+/*=int main(int, char**)
+{*/
+   using boost::math::constants::pi;
+   typedef boost::multiprecision::cpp_dec_float_50 mp_type;
+
+   const float j2_f =
+      integral(0.0F,
+      pi<float>(),
+      0.01F,
+      cyl_bessel_j_integral_rep<float>(2U, 1.23F)) / pi<float>();
+
+   const double j2_d =
+      integral(0.0,
+      pi<double>(),
+      0.0001,
+      cyl_bessel_j_integral_rep<double>(2U, 1.23)) / pi<double>();
+
+   const mp_type j2_mp =
+      integral(mp_type(0),
+      pi<mp_type>(),
+      mp_type(1.0E-20),
+      cyl_bessel_j_integral_rep<mp_type>(2U, mp_type(123) / 100)) / pi<mp_type>();
+
+   // 0.166369
+   std::cout
+      << std::setprecision(std::numeric_limits<float>::digits10)
+      << j2_f
+      << std::endl;
+
+   // 0.166369383786814
+   std::cout
+      << std::setprecision(std::numeric_limits<double>::digits10)
+      << j2_d
+      << std::endl;
+
+   // 0.16636938378681407351267852431513159437103348245333
+   std::cout
+      << std::setprecision(std::numeric_limits<mp_type>::digits10)
+      << j2_mp
+      << std::endl;
+
+   //
+   // Print true value for comparison:
+   // 0.166369383786814073512678524315131594371033482453329
+   std::cout << boost::math::cyl_bessel_j(2, mp_type(123) / 100) << std::endl;
+//=}
+
+   //]
+
+   std::cout << std::setprecision(std::numeric_limits< ::mp_type>::digits10) << std::scientific;
+   std::cout << mysin(boost::math::constants::pi< ::mp_type>() / 4) << std::endl;
+   std::cout << boost::multiprecision::sin(boost::math::constants::pi< ::mp_type>() / 4) << std::endl;
+
+   return 0;
+}
+
+/*
+
+Program output:
+
+9.822663964796047e-001
+9.82266396479604757017335009796882833995903762577173e-01
+9.822663964796047e-001
+9.82266396479604757017335009796882833995903762577173e-01
+9.8226639647960475701733500979688283399590376257717309069410413822165082248153638454147004236848917775e-01
+4.752916e+000
+4.752915525615998e+000
+4.75291552561599819047013317456355991350189758431460e+00
+5.000029e-001
+4.999999999998876e-001
+4.99999999999999999999999999999999999999999999999999e-01
+2.70670566473225383787998989944968806815263091819151e-01
+2.70670566473225383787998989944968806815253190143120e-01
+7.0710678118654752440084436210484903928483593768847403658833986900e-01
+7.0710678118654752440084436210484903928483593768847403658833986900e-01
+*/
+
+#else
+
+int main() { return 0; }
+
+#endif
diff --git a/src/boost/libs/multiprecision/example/gauss_laguerre_quadrature.cpp b/src/boost/libs/multiprecision/example/gauss_laguerre_quadrature.cpp
new file mode 100644
index 00000000..1ac329a4
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/gauss_laguerre_quadrature.cpp
@@ -0,0 +1,521 @@
+// Copyright Nick Thompson, 2017
+// Copyright John Maddock 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 <cmath>
+#include <cstdint>
+#include <functional>
+#include <iomanip>
+#include <iostream>
+#include <numeric>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/special_functions/cbrt.hpp>
+#include <boost/math/special_functions/factorials.hpp>
+#include <boost/math/special_functions/gamma.hpp>
+#include <boost/math/tools/roots.hpp>
+#include <boost/noncopyable.hpp>
+
+#define CPP_BIN_FLOAT  1
+#define CPP_DEC_FLOAT  2
+#define CPP_MPFR_FLOAT 3
+
+//#define MP_TYPE CPP_BIN_FLOAT
+#define MP_TYPE CPP_DEC_FLOAT
+//#define MP_TYPE CPP_MPFR_FLOAT
+
+namespace
+{
+  struct digits_characteristics
+  {
+    static const int digits10     = 300;
+    static const int guard_digits =   6;
+  };
+}
+
+#if (MP_TYPE == CPP_BIN_FLOAT)
+  #include <boost/multiprecision/cpp_bin_float.hpp>
+  namespace mp = boost::multiprecision;
+  typedef mp::number<mp::cpp_bin_float<digits_characteristics::digits10 + digits_characteristics::guard_digits>, mp::et_off> mp_type;
+#elif (MP_TYPE == CPP_DEC_FLOAT)
+  #include <boost/multiprecision/cpp_dec_float.hpp>
+  namespace mp = boost::multiprecision;
+  typedef mp::number<mp::cpp_dec_float<digits_characteristics::digits10 + digits_characteristics::guard_digits>, mp::et_off> mp_type;
+#elif (MP_TYPE == CPP_MPFR_FLOAT)
+  #include <boost/multiprecision/mpfr.hpp>
+  namespace mp = boost::multiprecision;
+  typedef mp::number<mp::mpfr_float_backend<digits_characteristics::digits10 + digits_characteristics::guard_digits>, mp::et_off> mp_type;
+#else
+#error MP_TYPE is undefined
+#endif
+
+template<typename T>
+class laguerre_function_object
+{
+public:
+  laguerre_function_object(const int n, const T a) : order(n),
+                                                     alpha(a),
+                                                     p1   (0),
+                                                     d2   (0) { }
+
+  laguerre_function_object(const laguerre_function_object& other) : order(other.order),
+                                                                    alpha(other.alpha),
+                                                                    p1   (other.p1),
+                                                                    d2   (other.d2) { }
+
+  ~laguerre_function_object() { }
+
+  T operator()(const T& x) const
+  {
+    // Calculate (via forward recursion):
+    // * the value of the Laguerre function L(n, alpha, x), called (p2),
+    // * the value of the derivative of the Laguerre function (d2),
+    // * and the value of the corresponding Laguerre function of
+    //   previous order (p1).
+
+    // Return the value of the function (p2) in order to be used as a
+    // function object with Boost.Math root-finding. Store the values
+    // of the Laguerre function derivative (d2) and the Laguerre function
+    // of previous order (p1) in class members for later use.
+
+      p1 = T(0);
+    T p2 = T(1);
+      d2 = T(0);
+
+    T j_plus_alpha(alpha);
+    T two_j_plus_one_plus_alpha_minus_x(1 + alpha - x);
+
+    int j;
+
+    const T my_two(2);
+
+    for(j = 0; j < order; ++j)
+    {
+      const T p0(p1);
+
+      // Set the value of the previous Laguerre function.
+      p1 = p2;
+
+      // Use a recurrence relation to compute the value of the Laguerre function.
+      p2 = ((two_j_plus_one_plus_alpha_minus_x * p1) - (j_plus_alpha * p0)) / (j + 1);
+
+      ++j_plus_alpha;
+      two_j_plus_one_plus_alpha_minus_x += my_two;
+    }
+
+    // Set the value of the derivative of the Laguerre function.
+    d2 = ((p2 * j) - (j_plus_alpha * p1)) / x;
+
+    // Return the value of the Laguerre function.
+    return p2;
+  }
+
+  const T& previous  () const { return p1; }
+  const T& derivative() const { return d2; }
+
+  static bool root_tolerance(const T& a, const T& b)
+  {
+    using std::abs;
+
+    // The relative tolerance here is: ((a - b) * 2) / (a + b).
+    return (abs((a - b) * 2) < ((a + b) * boost::math::tools::epsilon<T>()));
+  }
+
+private:
+  const int order;
+  const T   alpha;
+  mutable T p1;
+  mutable T d2;
+
+  laguerre_function_object();
+
+  const laguerre_function_object& operator=(const laguerre_function_object&);
+};
+
+template<typename T>
+class guass_laguerre_abscissas_and_weights : private boost::noncopyable
+{
+public:
+  guass_laguerre_abscissas_and_weights(const int n, const T a) : order(n),
+                                                                 alpha(a),
+                                                                 valid(true),
+                                                                 xi   (),
+                                                                 wi   ()
+  {
+    if(alpha < -20.0F)
+    {
+      // TBD: If we ever boostify this, throw a range error here.
+      // If so, then also document it in the docs.
+      std::cout << "Range error: the order of the Laguerre function must exceed -20.0." << std::endl;
+    }
+    else
+    {
+      calculate();
+    }
+  }
+
+  virtual ~guass_laguerre_abscissas_and_weights() { }
+
+  const std::vector<T>& abscissas() const { return xi; }
+  const std::vector<T>& weights  () const { return wi; }
+
+  bool get_valid() const { return valid; }
+
+private:
+  const int order;
+  const T   alpha;
+  bool      valid;
+
+  std::vector<T> xi;
+  std::vector<T> wi;
+
+  void calculate()
+  {
+    using std::abs;
+
+    std::cout << "finding approximate roots..." << std::endl;
+
+    std::vector<boost::math::tuple<T, T> > root_estimates;
+
+    root_estimates.reserve(static_cast<typename std::vector<boost::math::tuple<T, T> >::size_type>(order));
+
+    const laguerre_function_object<T> laguerre_object(order, alpha);
+
+    // Set the initial values of the step size and the running step
+    // to be used for finding the estimate of the first root.
+    T step_size  = 0.01F;
+    T step       = step_size;
+
+    T first_laguerre_root = 0.0F;
+
+    bool first_laguerre_root_has_been_found = true;
+
+    if(alpha < -1.0F)
+    {
+      // Iteratively step through the Laguerre function using a
+      // small step-size in order to find a rough estimate of
+      // the first zero.
+
+      bool this_laguerre_value_is_negative = (laguerre_object(mp_type(0)) < 0);
+
+      static const int j_max = 10000;
+
+      int j;
+
+      for(j = 0; (j < j_max) && (this_laguerre_value_is_negative != (laguerre_object(step) < 0)); ++j)
+      {
+        // Increment the step size until the sign of the Laguerre function
+        // switches. This indicates a zero-crossing, signalling the next root.
+        step += step_size;
+      }
+
+      if(j >= j_max)
+      {
+        first_laguerre_root_has_been_found = false;
+      }
+      else
+      {
+        // We have found the first zero-crossing. Put a loose bracket around
+        // the root using a window. Here, we know that the first root lies
+        // between (x - step_size) < root < x.
+
+        // Before storing the approximate root, perform a couple of
+        // bisection steps in order to tighten up the root bracket.
+        boost::uintmax_t a_couple_of_iterations = 3U;
+        const std::pair<T, T>
+          first_laguerre_root = boost::math::tools::bisect(laguerre_object,
+                                                           step - step_size,
+                                                           step,
+                                                           laguerre_function_object<T>::root_tolerance,
+                                                           a_couple_of_iterations);
+
+        static_cast<void>(a_couple_of_iterations);
+      }
+    }
+    else
+    {
+      // Calculate an estimate of the 1st root of a generalized Laguerre
+      // function using either a Taylor series or an expansion in Bessel
+      // function zeros. The Bessel function zeros expansion is from Tricomi.
+
+      // Here, we obtain an estimate of the first zero of J_alpha(x).
+
+      T j_alpha_m1;
+
+      if(alpha < 1.4F)
+      {
+        // For small alpha, use a short series obtained from Mathematica(R).
+        // Series[BesselJZero[v, 1], {v, 0, 3}]
+        // N[%, 12]
+        j_alpha_m1 = (((          0.09748661784476F
+                        * alpha - 0.17549359276115F)
+                        * alpha + 1.54288974259931F)
+                        * alpha + 2.40482555769577F);
+      }
+      else
+      {
+        // For larger alpha, use the first line of Eqs. 10.21.40 in the NIST Handbook.
+        const T alpha_pow_third(boost::math::cbrt(alpha));
+        const T alpha_pow_minus_two_thirds(T(1) / (alpha_pow_third * alpha_pow_third));
+
+        j_alpha_m1 = alpha * (((((                             + 0.043F
+                                  * alpha_pow_minus_two_thirds - 0.0908F)
+                                  * alpha_pow_minus_two_thirds - 0.00397F)
+                                  * alpha_pow_minus_two_thirds + 1.033150F)
+                                  * alpha_pow_minus_two_thirds + 1.8557571F)
+                                  * alpha_pow_minus_two_thirds + 1.0F);
+      }
+
+      const T vf             = ((order * 4.0F) + (alpha * 2.0F) + 2.0F);
+      const T vf2            = vf * vf;
+      const T j_alpha_m1_sqr = j_alpha_m1 * j_alpha_m1;
+
+      first_laguerre_root = (j_alpha_m1_sqr * (-0.6666666666667F + ((0.6666666666667F * alpha) * alpha) + (0.3333333333333F * j_alpha_m1_sqr) + vf2)) / (vf2 * vf);
+    }
+
+    if(first_laguerre_root_has_been_found)
+    {
+      bool this_laguerre_value_is_negative = (laguerre_object(mp_type(0)) < 0);
+
+      // Re-set the initial value of the step-size based on the
+      // estimate of the first root.
+      step_size = first_laguerre_root / 2;
+      step      = step_size;
+
+      // Step through the Laguerre function using a step-size
+      // of dynamic width in order to find the zero crossings
+      // of the Laguerre function, providing rough estimates
+      // of the roots. Refine the brackets with a few bisection
+      // steps, and store the results as bracketed root estimates.
+
+      while(static_cast<int>(root_estimates.size()) < order)
+      {
+        // Increment the step size until the sign of the Laguerre function
+        // switches. This indicates a zero-crossing, signalling the next root.
+        step += step_size;
+
+        if(this_laguerre_value_is_negative != (laguerre_object(step) < 0))
+        {
+          // We have found the next zero-crossing.
+
+          // Change the running sign of the Laguerre function.
+          this_laguerre_value_is_negative = (!this_laguerre_value_is_negative);
+
+          // We have found the first zero-crossing. Put a loose bracket around
+          // the root using a window. Here, we know that the first root lies
+          // between (x - step_size) < root < x.
+
+          // Before storing the approximate root, perform a couple of
+          // bisection steps in order to tighten up the root bracket.
+          boost::uintmax_t a_couple_of_iterations = 3U;
+          const std::pair<T, T>
+            root_estimate_bracket = boost::math::tools::bisect(laguerre_object,
+                                                               step - step_size,
+                                                               step,
+                                                               laguerre_function_object<T>::root_tolerance,
+                                                               a_couple_of_iterations);
+
+          static_cast<void>(a_couple_of_iterations);
+
+          // Store the refined root estimate as a bracketed range in a tuple.
+          root_estimates.push_back(boost::math::tuple<T, T>(root_estimate_bracket.first,
+                                                            root_estimate_bracket.second));
+
+          if(root_estimates.size() >= static_cast<std::size_t>(2U))
+          {
+            // Determine the next step size. This is based on the distance between
+            // the previous two roots, whereby the estimates of the previous roots
+            // are computed by taking the average of the lower and upper range of
+            // the root-estimate bracket.
+
+            const T r0 = (  boost::math::get<0>(*(root_estimates.rbegin() + 1U))
+                          + boost::math::get<1>(*(root_estimates.rbegin() + 1U))) / 2;
+
+            const T r1 = (  boost::math::get<0>(*root_estimates.rbegin())
+                          + boost::math::get<1>(*root_estimates.rbegin())) / 2;
+
+            const T distance_between_previous_roots = r1 - r0;
+
+            step_size = distance_between_previous_roots / 3;
+          }
+        }
+      }
+
+      const T norm_g =
+        ((alpha == 0) ? T(-1)
+                      : -boost::math::tgamma(alpha + order) / boost::math::factorial<T>(order - 1));
+
+      xi.reserve(root_estimates.size());
+      wi.reserve(root_estimates.size());
+
+      // Calculate the abscissas and weights to full precision.
+      for(std::size_t i = static_cast<std::size_t>(0U); i < root_estimates.size(); ++i)
+      {
+        std::cout << "calculating abscissa and weight for index: " << i << std::endl;
+
+        // Calculate the abscissas using iterative root-finding.
+
+        // Select the maximum allowed iterations, being at least 20.
+        // The determination of the maximum allowed iterations is
+        // based on the number of decimal digits in the numerical
+        // type T.
+        const int my_digits10 = static_cast<int>(static_cast<float>(boost::math::tools::digits<T>()) * 0.301F);
+        const boost::uintmax_t number_of_iterations_allowed = (std::max)(20, my_digits10 / 2);
+
+        boost::uintmax_t number_of_iterations_used = number_of_iterations_allowed;
+
+        // Perform the root-finding using ACM TOMS 748 from Boost.Math.
+        const std::pair<T, T>
+          laguerre_root_bracket = boost::math::tools::toms748_solve(laguerre_object,
+                                                                    boost::math::get<0>(root_estimates[i]),
+                                                                    boost::math::get<1>(root_estimates[i]),
+                                                                    laguerre_function_object<T>::root_tolerance,
+                                                                    number_of_iterations_used);
+
+        // Based on the result of *each* root-finding operation, re-assess
+        // the validity of the Guass-Laguerre abscissas and weights object.
+        valid &= (number_of_iterations_used < number_of_iterations_allowed);
+
+        // Compute the Laguerre root as the average of the values from
+        // the solved root bracket.
+        const T laguerre_root = (  laguerre_root_bracket.first
+                                 + laguerre_root_bracket.second) / 2;
+
+        // Calculate the weight for this Laguerre root. Here, we calculate
+        // the derivative of the Laguerre function and the value of the
+        // previous Laguerre function on the x-axis at the value of this
+        // Laguerre root.
+        static_cast<void>(laguerre_object(laguerre_root));
+
+        // Store the abscissa and weight for this index.
+        xi.push_back(laguerre_root);
+        wi.push_back(norm_g / ((laguerre_object.derivative() * order) * laguerre_object.previous()));
+      }
+    }
+  }
+};
+
+namespace
+{
+  template<typename T>
+  struct gauss_laguerre_ai
+  {
+  public:
+    gauss_laguerre_ai(const T X) : x(X)
+    {
+      using std::exp;
+      using std::sqrt;
+
+      zeta = ((sqrt(x) * x) * 2) / 3;
+
+      const T zeta_times_48_pow_sixth = sqrt(boost::math::cbrt(zeta * 48));
+
+      factor = 1 / ((sqrt(boost::math::constants::pi<T>()) * zeta_times_48_pow_sixth) * (exp(zeta) * gamma_of_five_sixths()));
+    }
+
+    gauss_laguerre_ai(const gauss_laguerre_ai& other) : x     (other.x),
+                                                        zeta  (other.zeta),
+                                                        factor(other.factor) { }
+
+    T operator()(const T& t) const
+    {
+      using std::sqrt;
+
+      return factor / sqrt(boost::math::cbrt(2 + (t / zeta)));
+    }
+
+  private:
+    const T x;
+    T zeta;
+    T factor;
+
+    static const T& gamma_of_five_sixths()
+    {
+      static const T value = boost::math::tgamma(T(5) / 6);
+
+      return value;
+    }
+
+    const gauss_laguerre_ai& operator=(const gauss_laguerre_ai&);
+  };
+
+  template<typename T>
+  T gauss_laguerre_airy_ai(const T x)
+  {
+    static const float digits_factor  = static_cast<float>(std::numeric_limits<mp_type>::digits10) / 300.0F;
+    static const int   laguerre_order = static_cast<int>(600.0F * digits_factor);
+
+    static const guass_laguerre_abscissas_and_weights<T> abscissas_and_weights(laguerre_order, -T(1) / 6);
+
+    T airy_ai_result;
+
+    if(abscissas_and_weights.get_valid())
+    {
+      const gauss_laguerre_ai<T> this_gauss_laguerre_ai(x);
+
+      airy_ai_result =
+        std::inner_product(abscissas_and_weights.abscissas().begin(),
+                           abscissas_and_weights.abscissas().end(),
+                           abscissas_and_weights.weights().begin(),
+                           T(0),
+                           std::plus<T>(),
+                           [&this_gauss_laguerre_ai](const T& this_abscissa, const T& this_weight) -> T
+                           {
+                             return this_gauss_laguerre_ai(this_abscissa) * this_weight;
+                           });
+    }
+    else
+    {
+      // TBD: Consider an error message.
+      airy_ai_result = T(0);
+    }
+
+    return airy_ai_result;
+  }
+}
+
+int main()
+{
+  // Use Gauss-Laguerre integration to compute airy_ai(120 / 7).
+
+  // 9 digits
+  // 3.89904210e-22
+
+  // 10 digits
+  // 3.899042098e-22
+
+  // 50 digits.
+  // 3.8990420982303275013276114626640705170145070824318e-22
+
+  // 100 digits.
+  // 3.899042098230327501327611462664070517014507082431797677146153303523108862015228
+  // 864136051942933142648e-22
+
+  // 200 digits.
+  // 3.899042098230327501327611462664070517014507082431797677146153303523108862015228
+  // 86413605194293314264788265460938200890998546786740097437064263800719644346113699
+  // 77010905030516409847054404055843899790277e-22
+
+  // 300 digits.
+  // 3.899042098230327501327611462664070517014507082431797677146153303523108862015228
+  // 86413605194293314264788265460938200890998546786740097437064263800719644346113699
+  // 77010905030516409847054404055843899790277083960877617919088116211775232728792242
+  // 9346416823281460245814808276654088201413901972239996130752528e-22
+
+  // 500 digits.
+  // 3.899042098230327501327611462664070517014507082431797677146153303523108862015228
+  // 86413605194293314264788265460938200890998546786740097437064263800719644346113699
+  // 77010905030516409847054404055843899790277083960877617919088116211775232728792242
+  // 93464168232814602458148082766540882014139019722399961307525276722937464859521685
+  // 42826483602153339361960948844649799257455597165900957281659632186012043089610827
+  // 78871305322190941528281744734605934497977375094921646511687434038062987482900167
+  // 45127557400365419545e-22
+
+  // Mathematica(R) or Wolfram's Alpha:
+  // N[AiryAi[120 / 7], 300]
+  std::cout << std::setprecision(digits_characteristics::digits10)
+            << gauss_laguerre_airy_ai(mp_type(120) / 7)
+            << std::endl;
+}
diff --git a/src/boost/libs/multiprecision/example/gmp_snips.cpp b/src/boost/libs/multiprecision/example/gmp_snips.cpp
new file mode 100644
index 00000000..f41e77bf
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/gmp_snips.cpp
@@ -0,0 +1,117 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2011 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/multiprecision/gmp.hpp>
+#include <boost/math/special_functions/gamma.hpp>
+#include <iostream>
+
+void t1()
+{
+   //[mpz_eg
+//=#include <boost/multiprecision/gmp.hpp>
+//=#include <iostream>
+//=
+//=int main()
+//={
+   using namespace boost::multiprecision;
+
+   mpz_int v = 1;
+
+   // Do some arithmetic:
+   for(unsigned i = 1; i <= 1000; ++i)
+      v *= i;
+
+   std::cout << v << std::endl; // prints 1000!
+
+   // Access the underlying representation:
+   mpz_t z;
+   mpz_init(z);
+   mpz_set(z, v.backend().data());
+   mpz_clear(z);
+//=   return 0;
+//=}
+//]
+}
+
+void t2()
+{
+//[mpf_eg
+//=#include <boost/multiprecision/gmp.hpp>
+//=#include <boost/math/special_functions/gamma.hpp>
+//=#include <iostream>
+//=
+//=int main()
+//={
+   using namespace boost::multiprecision;
+
+   // Operations at variable precision and limited standard library support:
+   mpf_float a = 2;
+   mpf_float::default_precision(1000);
+   std::cout << mpf_float::default_precision() << std::endl;
+   std::cout << sqrt(a) << std::endl; // print root-2
+
+   // Operations at fixed precision and full standard library support:
+   mpf_float_100 b = 2;
+   std::cout << std::numeric_limits<mpf_float_100>::digits << std::endl;
+   // We can use any C++ std lib function:
+   std::cout << log(b) << std::endl; // print log(2)
+   // We can also use any function from Boost.Math:
+   std::cout << boost::math::tgamma(b) << std::endl;
+   // These even work when the argument is an expression template:
+   std::cout << boost::math::tgamma(b * b) << std::endl;
+
+   // Access the underlying representation:
+   mpf_t f;
+   mpf_init(f);
+   mpf_set(f, a.backend().data());
+   mpf_clear(f);
+//=   return 0;
+//=}
+//]
+}
+
+void t3()
+{
+//[mpq_eg
+//=#include <boost/multiprecision/gmp.hpp>
+//=#include <boost/multiprecision/gmp.hpp>
+//=#include <iostream>
+//=
+//=int main()
+//={
+   using namespace boost::multiprecision;
+
+   mpq_rational v = 1;
+
+   // Do some arithmetic:
+   for(unsigned i = 1; i <= 1000; ++i)
+      v *= i;
+   v /= 10;
+
+   std::cout << v << std::endl; // prints 1000! / 10
+   std::cout << numerator(v) << std::endl;
+   std::cout << denominator(v) << std::endl;
+
+   mpq_rational w(2, 3);  // component wise constructor
+   std::cout << w << std::endl; // prints 2/3
+
+   // Access the underlying data:
+   mpq_t q;
+   mpq_init(q);
+   mpq_set(q, v.backend().data());
+   mpq_clear(q);
+//=   return 0;
+//=}
+//]
+}
+
+int main()
+{
+   t1();
+   t2();
+   t3();
+   return 0;
+}
+
diff --git a/src/boost/libs/multiprecision/example/hashing_examples.cpp b/src/boost/libs/multiprecision/example/hashing_examples.cpp
new file mode 100644
index 00000000..a0c3fb50
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/hashing_examples.cpp
@@ -0,0 +1,77 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2012 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/multiprecision/cpp_int.hpp>
+#include <boost/random.hpp>
+#include <boost/functional/hash.hpp>
+#include <unordered_set>
+#include <city.h>
+
+//[hash1
+
+/*`
+All of the types in this library support hashing via boost::hash or std::hash.
+That means we can use multiprecision types directly in hashed containers such as std::unordered_set:
+*/
+//]
+
+void t1()
+{
+   //[hash2
+   using namespace boost::multiprecision;
+   using namespace boost::random;
+
+   mt19937 mt;
+   uniform_int_distribution<uint256_t> ui;
+
+   std::unordered_set<uint256_t> set;
+   // Put 1000 random values into the container:
+   for(unsigned i = 0; i < 1000; ++i)
+      set.insert(ui(mt));
+
+   //]
+}
+
+//[hash3
+
+/*` 
+Or we can define our own hash function, for example in this case based on
+Google's CityHash:
+*/
+
+struct cityhash
+{
+   std::size_t operator()(const boost::multiprecision::uint256_t& val)const
+   {
+      // create a hash from all the limbs of the argument, this function is probably x64 specific,
+      // and requires that we access the internals of the data type:
+      std::size_t result = CityHash64(reinterpret_cast<const char*>(val.backend().limbs()), val.backend().size() * sizeof(val.backend().limbs()[0]));
+      // modify the returned hash based on sign:
+      return val < 0 ? ~result : result;
+   }
+};
+
+//]
+
+void t2()
+{
+//[hash4
+
+/*`As before insert some values into a container, this time using our custom hasher:*/
+
+   std::unordered_set<uint256_t, cityhash> set2;
+   for(unsigned i = 0; i < 1000; ++i)
+      set2.insert(ui(mt));
+
+//]
+}
+
+int main()
+{
+   t1();
+   t2();
+   return 0;
+}
+
diff --git a/src/boost/libs/multiprecision/example/hypergeometric_luke_algorithms.cpp b/src/boost/libs/multiprecision/example/hypergeometric_luke_algorithms.cpp
new file mode 100644
index 00000000..4fb4e13b
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/hypergeometric_luke_algorithms.cpp
@@ -0,0 +1,801 @@
+
+///////////////////////////////////////////////////////////////////////////////
+// Copyright Christopher Kormanyos 2013 - 2014.
+// Copyright John Maddock 2013.
+// Distributed under 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)
+//
+// This work is based on an earlier work:
+// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
+// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
+//
+
+#include <algorithm>
+#include <cstdint>
+#include <deque>
+#include <functional>
+#include <iostream>
+#include <limits>
+#include <numeric>
+#include <vector>
+#include <boost/math/constants/constants.hpp>
+#include <boost/noncopyable.hpp>
+
+//#define USE_CPP_BIN_FLOAT
+#define USE_CPP_DEC_FLOAT
+//#define USE_MPFR
+
+#if !defined(DIGIT_COUNT)
+#define DIGIT_COUNT 100
+#endif
+
+#if !defined(BOOST_NO_CXX11_HDR_CHRONO)
+  #include <chrono>
+  #define STD_CHRONO std::chrono
+#else
+  #include <boost/chrono.hpp>
+  #define STD_CHRONO boost::chrono
+#endif
+
+#if defined(USE_CPP_BIN_FLOAT)
+  #include <boost/multiprecision/cpp_bin_float.hpp>
+  typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<DIGIT_COUNT + 10> > mp_type;
+#elif defined(USE_CPP_DEC_FLOAT)
+  #include <boost/multiprecision/cpp_dec_float.hpp>
+  typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<DIGIT_COUNT + 10> > mp_type;
+#elif defined(USE_MPFR)
+  #include <boost/multiprecision/mpfr.hpp>
+  typedef boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<DIGIT_COUNT + 10> > mp_type;
+#else
+  #error no multiprecision floating type is defined
+#endif
+
+template <class clock_type>
+struct stopwatch
+{
+public:
+  typedef typename clock_type::duration duration_type;
+
+  stopwatch() : m_start(clock_type::now()) { }
+
+  stopwatch(const stopwatch& other) : m_start(other.m_start) { }
+
+  stopwatch& operator=(const stopwatch& other)
+  {
+    m_start = other.m_start;
+    return *this;
+  }
+
+  ~stopwatch() { }
+
+  duration_type elapsed() const
+  {
+    return (clock_type::now() - m_start);
+  }
+
+  void reset()
+  {
+    m_start = clock_type::now();
+  }
+
+private:
+  typename clock_type::time_point m_start;
+};
+
+namespace my_math
+{
+  template<class T> T chebyshev_t(const std::int32_t n, const T& x);
+
+  template<class T> T chebyshev_t(const std::uint32_t n, const T& x, std::vector<T>* vp);
+
+  template<class T> bool isneg(const T& x) { return (x < T(0)); }
+
+  template<class T> const T& zero() { static const T value_zero(0); return value_zero; }
+  template<class T> const T& one () { static const T value_one (1); return value_one; }
+  template<class T> const T& two () { static const T value_two (2); return value_two; }
+}
+
+namespace orthogonal_polynomial_series
+{
+  template<typename T> static inline T orthogonal_polynomial_template(const T& x, const std::uint32_t n, std::vector<T>* const vp = static_cast<std::vector<T>*>(0u))
+  {
+    // Compute the value of an orthogonal chebyshev polinomial.
+    // Use stable upward recursion.
+
+    if(vp != nullptr)
+    {
+      vp->clear();
+      vp->reserve(static_cast<std::size_t>(n + 1u));
+    }
+
+    T y0 = my_math::one<T>();
+
+    if(vp != nullptr) { vp->push_back(y0); }
+
+    if(n == static_cast<std::uint32_t>(0u))
+    {
+      return y0;
+    }
+
+    T y1 = x;
+
+    if(vp != nullptr) { vp->push_back(y1); }
+
+    if(n == static_cast<std::uint32_t>(1u))
+    {
+      return y1;
+    }
+
+    T a = my_math::two <T>();
+    T b = my_math::zero<T>();
+    T c = my_math::one <T>();
+
+    T yk;
+
+    // Calculate higher orders using the recurrence relation.
+    // The direction of stability is upward recursion.
+    for(std::int32_t k = static_cast<std::int32_t>(2); k <= static_cast<std::int32_t>(n); ++k)
+    {
+      yk = (((a * x) + b) * y1) - (c * y0);
+
+      y0 = y1;
+      y1 = yk;
+
+      if(vp != nullptr) { vp->push_back(yk); }
+    }
+
+    return yk;
+  }
+}
+
+template<class T> T my_math::chebyshev_t(const std::int32_t n, const T& x)
+{
+  if(my_math::isneg(x))
+  {
+    const bool b_negate = ((n % static_cast<std::int32_t>(2)) != static_cast<std::int32_t>(0));
+
+    const T y = chebyshev_t(n, -x);
+
+    return (!b_negate ? y : -y);
+  }
+
+  if(n < static_cast<std::int32_t>(0))
+  {
+    const std::int32_t nn = static_cast<std::int32_t>(-n);
+
+    return chebyshev_t(nn, x);
+  }
+  else
+  {
+    return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::uint32_t>(n));
+  }
+}
+
+template<class T> T my_math::chebyshev_t(const std::uint32_t n, const T& x, std::vector<T>* const vp) { return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::int32_t>(n),  vp); }
+
+namespace util
+{
+  template <class T> float digit_scale()
+  {
+    const int d = ((std::max)(std::numeric_limits<T>::digits10, 15));
+    return static_cast<float>(d) / 300.0F;
+  }
+}
+
+namespace examples
+{
+  namespace nr_006
+  {
+    template<typename T> class hypergeometric_pfq_base : private boost::noncopyable
+    {
+    public:
+      virtual ~hypergeometric_pfq_base() { }
+
+      virtual void ccoef() const = 0;
+
+      virtual T series() const
+      {
+        using my_math::chebyshev_t;
+
+        // Compute the Chebyshev coefficients.
+        // Get the values of the shifted Chebyshev polynomials.
+        std::vector<T> chebyshev_t_shifted_values;
+        const T z_shifted = ((Z / W) * static_cast<std::int32_t>(2)) - static_cast<std::int32_t>(1);
+
+        chebyshev_t(static_cast<std::uint32_t>(C.size()),
+                    z_shifted,
+                    &chebyshev_t_shifted_values);
+
+        // Luke: C     ---------- COMPUTE SCALE FACTOR                       ----------
+        // Luke: C
+        // Luke: C     ---------- SCALE THE COEFFICIENTS                     ----------
+        // Luke: C
+
+        // The coefficient scaling is preformed after the Chebyshev summation,
+        // and it is carried out with a single division operation.
+        bool b_neg = false;
+
+        const T scale = std::accumulate(C.begin(),
+                                        C.end(),
+                                        T(0),
+                                        [&b_neg](T scale_sum, const T& ck) -> T
+                                        {
+                                          ((!b_neg) ? (scale_sum += ck) : (scale_sum -= ck));
+                                          b_neg = (!b_neg);
+                                          return scale_sum;
+                                        });
+
+        // Compute the result of the series expansion using unscaled coefficients.
+        const T sum = std::inner_product(C.begin(),
+                                         C.end(),
+                                         chebyshev_t_shifted_values.begin(),
+                                         T(0));
+
+        // Return the properly scaled result.
+        return sum / scale;
+      }
+
+    protected:
+      const   T             Z;
+      const   T             W;
+      mutable std::deque<T> C;
+
+      hypergeometric_pfq_base(const T& z,
+                              const T& w) : Z(z),
+                                            W(w),
+                                            C(0u) { }
+
+      virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale<T>() * 500.0F); }
+    };
+
+    template<typename T> class ccoef4_hypergeometric_0f1 : public hypergeometric_pfq_base<T>
+    {
+    public:
+      ccoef4_hypergeometric_0f1(const T& c,
+                                const T& z,
+                                const T& w) : hypergeometric_pfq_base<T>(z, w),
+                                              CP(c) { }
+
+      virtual ~ccoef4_hypergeometric_0f1() { }
+
+      virtual void ccoef() const
+      {
+        // See Luke 1977 page 80.
+        const std::int32_t N1 = static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(1));
+        const std::int32_t N2 = static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(2));
+
+        // Luke: C     ---------- START COMPUTING COEFFICIENTS USING         ----------
+        // Luke: C     ---------- BACKWARD RECURRENCE SCHEME                 ----------
+        // Luke: C
+        T A3(0);
+        T A2(0);
+        T A1(boost::math::tools::root_epsilon<T>());
+
+        hypergeometric_pfq_base<T>::C.resize(1u, A1);
+
+        std::int32_t X1 = N2;
+
+        T C1 = T(1) - CP;
+
+        const T Z1 = T(4) / hypergeometric_pfq_base<T>::W;
+
+        for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
+        {
+          const T DIVFAC = T(1) / X1;
+
+          --X1;
+
+          // The terms have been slightly re-arranged resulting in lower complexity.
+          // Parentheses have been added to avoid reliance on operator precedence.
+          const T term =   (A2 - ((A3 * DIVFAC) * X1))
+                         + ((A2 * X1) * ((1 + (C1 + X1)) * Z1))
+                         + ((A1 * X1) * ((DIVFAC - (C1 * Z1)) + (X1 * Z1)));
+
+          hypergeometric_pfq_base<T>::C.push_front(term);
+
+          A3 = A2;
+          A2 = A1;
+          A1 = hypergeometric_pfq_base<T>::C.front();
+        }
+
+        hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
+      }
+
+    private:
+      const T CP;
+    };
+
+    template<typename T> class ccoef1_hypergeometric_1f0 : public hypergeometric_pfq_base<T>
+    {
+    public:
+      ccoef1_hypergeometric_1f0(const T& a,
+                                const T& z,
+                                const T& w) : hypergeometric_pfq_base<T>(z, w),
+                                              AP(a) { }
+
+      virtual ~ccoef1_hypergeometric_1f0() { }
+
+      virtual void ccoef() const
+      {
+        // See Luke 1977 page 67.
+        const std::int32_t N1 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1));
+        const std::int32_t N2 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(2));
+
+        // Luke: C     ---------- START COMPUTING COEFFICIENTS USING         ----------
+        // Luke: C     ---------- BACKWARD RECURRENCE SCHEME                 ----------
+        // Luke: C
+        T A2(0);
+        T A1(boost::math::tools::root_epsilon<T>());
+
+        hypergeometric_pfq_base<T>::C.resize(1u, A1);
+
+        std::int32_t X1 = N2;
+
+        T V1 = T(1) - AP;
+
+        // Here, we have corrected what appears to be an error in Luke's code.
+
+        // Luke's original code listing has:
+        //  AFAC = 2 + FOUR/W
+        // But it appears as though the correct form is:
+        //  AFAC = 2 - FOUR/W.
+
+        const T AFAC = 2 - (T(4) / hypergeometric_pfq_base<T>::W);
+
+        for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
+        {
+          --X1;
+
+          // The terms have been slightly re-arranged resulting in lower complexity.
+          // Parentheses have been added to avoid reliance on operator precedence.
+          const T term = -(((X1 * AFAC) * A1) + ((X1 + V1) * A2)) / (X1 - V1);
+
+          hypergeometric_pfq_base<T>::C.push_front(term);
+
+          A2 = A1;
+          A1 = hypergeometric_pfq_base<T>::C.front();
+        }
+
+        hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
+      }
+
+    private:
+      const T AP;
+
+      virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale<T>() * 1600.0F); }
+    };
+
+    template<typename T> class ccoef3_hypergeometric_1f1 : public hypergeometric_pfq_base<T>
+    {
+    public:
+      ccoef3_hypergeometric_1f1(const T& a,
+                                const T& c,
+                                const T& z,
+                                const T& w) : hypergeometric_pfq_base<T>(z, w),
+                                              AP(a),
+                                              CP(c) { }
+
+      virtual ~ccoef3_hypergeometric_1f1() { }
+
+      virtual void ccoef() const
+      {
+        // See Luke 1977 page 74.
+        const std::int32_t N1 = static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(1));
+        const std::int32_t N2 = static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(2));
+
+        // Luke: C     ---------- START COMPUTING COEFFICIENTS USING         ----------
+        // Luke: C     ---------- BACKWARD RECURRENCE SCHEME                 ----------
+        // Luke: C
+        T A3(0);
+        T A2(0);
+        T A1(boost::math::tools::root_epsilon<T>());
+
+        hypergeometric_pfq_base<T>::C.resize(1u, A1);
+
+        std::int32_t X  = N1;
+        std::int32_t X1 = N2;
+
+        T XA  =  X + AP;
+        T X3A = (X + 3) - AP;
+
+        const T Z1 = T(4) / hypergeometric_pfq_base<T>::W;
+
+        for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
+        {
+          --X;
+          --X1;
+          --XA;
+          --X3A;
+
+          const T X3A_over_X2 = X3A / static_cast<std::int32_t>(X + 2);
+
+          // The terms have been slightly re-arranged resulting in lower complexity.
+          // Parentheses have been added to avoid reliance on operator precedence.
+          const T PART1 =  A1 * (((X + CP) * Z1) - X3A_over_X2);
+          const T PART2 =  A2 * (Z1 * ((X + 3) - CP) + (XA / X1));
+          const T PART3 =  A3 * X3A_over_X2;
+
+          const T term = (((PART1 + PART2) + PART3) * X1) / XA;
+
+          hypergeometric_pfq_base<T>::C.push_front(term);
+
+          A3 = A2;
+          A2 = A1;
+          A1 = hypergeometric_pfq_base<T>::C.front();
+        }
+
+        hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
+      }
+
+    private:
+      const T AP;
+      const T CP;
+    };
+
+    template<typename T> class ccoef6_hypergeometric_1f2 : public hypergeometric_pfq_base<T>
+    {
+    public:
+      ccoef6_hypergeometric_1f2(const T& a,
+                                const T& b,
+                                const T& c,
+                                const T& z,
+                                const T& w) : hypergeometric_pfq_base<T>(z, w),
+                                              AP(a),
+                                              BP(b),
+                                              CP(c) { }
+
+      virtual ~ccoef6_hypergeometric_1f2() { }
+
+      virtual void ccoef() const
+      {
+        // See Luke 1977 page 85.
+        const std::int32_t N1 = static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(1));
+
+        // Luke: C     ---------- START COMPUTING COEFFICIENTS USING         ----------
+        // Luke: C     ---------- BACKWARD RECURRENCE SCHEME                 ----------
+        // Luke: C
+        T A4(0);
+        T A3(0);
+        T A2(0);
+        T A1(boost::math::tools::root_epsilon<T>());
+
+        hypergeometric_pfq_base<T>::C.resize(1u, A1);
+
+        std::int32_t X  = N1;
+        T            PP = X + AP;
+
+        const T Z1 = T(4) / hypergeometric_pfq_base<T>::W;
+
+        for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
+        {
+          --X;
+          --PP;
+
+          const std::int32_t TWO_X    = static_cast<std::int32_t>(X * 2);
+          const std::int32_t X_PLUS_1 = static_cast<std::int32_t>(X + 1);
+          const std::int32_t X_PLUS_3 = static_cast<std::int32_t>(X + 3);
+          const std::int32_t X_PLUS_4 = static_cast<std::int32_t>(X + 4);
+
+          const T QQ = T(TWO_X + 3) / static_cast<std::int32_t>(TWO_X + static_cast<std::int32_t>(5));
+          const T SS = (X + BP) * (X + CP);
+
+          // The terms have been slightly re-arranged resulting in lower complexity.
+          // Parentheses have been added to avoid reliance on operator precedence.
+          const T PART1 =   A1 * (((PP - (QQ * (PP + 1))) * 2) + (SS * Z1));
+          const T PART2 =  (A2 * (X + 2)) * ((((TWO_X + 1) * PP) / X_PLUS_1) - ((QQ * 4) * (PP + 1)) + (((TWO_X + 3) * (PP + 2)) / X_PLUS_3) + ((Z1 * 2) * (SS - (QQ * (X_PLUS_1 + BP)) * (X_PLUS_1 + CP))));
+          const T PART3 =   A3 * ((((X_PLUS_3 - AP) - (QQ * (X_PLUS_4 - AP))) * 2) + (((QQ * Z1) * (X_PLUS_4 - BP)) * (X_PLUS_4 - CP)));
+          const T PART4 = ((A4 * QQ) * (X_PLUS_4 - AP)) / X_PLUS_3;
+
+          const T term = (((PART1 - PART2) + (PART3 - PART4)) * X_PLUS_1) / PP;
+
+          hypergeometric_pfq_base<T>::C.push_front(term);
+
+          A4 = A3;
+          A3 = A2;
+          A2 = A1;
+          A1 = hypergeometric_pfq_base<T>::C.front();
+        }
+
+        hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
+      }
+
+    private:
+      const T AP;
+      const T BP;
+      const T CP;
+    };
+
+    template<typename T> class ccoef2_hypergeometric_2f1 : public hypergeometric_pfq_base<T>
+    {
+    public:
+      ccoef2_hypergeometric_2f1(const T& a,
+                                const T& b,
+                                const T& c,
+                                const T& z,
+                                const T& w) : hypergeometric_pfq_base<T>(z, w),
+                                              AP(a),
+                                              BP(b),
+                                              CP(c) { }
+
+      virtual ~ccoef2_hypergeometric_2f1() { }
+
+      virtual void ccoef() const
+      {
+        // See Luke 1977 page 59.
+        const std::int32_t N1 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1));
+        const std::int32_t N2 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(2));
+
+        // Luke: C     ---------- START COMPUTING COEFFICIENTS USING         ----------
+        // Luke: C     ---------- BACKWARD RECURRENCE SCHEME                 ----------
+        // Luke: C
+        T A3(0);
+        T A2(0);
+        T A1(boost::math::tools::root_epsilon<T>());
+
+        hypergeometric_pfq_base<T>::C.resize(1u, A1);
+
+        std::int32_t X  = N1;
+        std::int32_t X1 = N2;
+        std::int32_t X3 = static_cast<std::int32_t>((X * 2) + 3);
+
+        T X3A = (X + 3) - AP;
+        T X3B = (X + 3) - BP;
+
+        const T Z1 = T(4) / hypergeometric_pfq_base<T>::W;
+
+        for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
+        {
+          --X;
+          --X1;
+          --X3A;
+          --X3B;
+          X3 -= 2;
+
+          const std::int32_t X_PLUS_2 = static_cast<std::int32_t>(X + 2);
+
+          const T XAB = T(1) / ((X + AP) * (X + BP));
+
+          // The terms have been slightly re-arranged resulting in lower complexity.
+          // Parentheses have been added to avoid reliance on operator precedence.
+          const T PART1 = (A1 * X1) * (2 - (((AP + X1) * (BP + X1)) * ((T(X3) / X_PLUS_2) * XAB)) + ((CP + X) * (XAB * Z1)));
+          const T PART2 = (A2 * XAB) * ((X3A * X3B) - (X3 * ((X3A + X3B) - 1)) + (((3 - CP) + X) * (X1 * Z1)));
+          const T PART3 = (A3 * X1) * (X3A / X_PLUS_2) * (X3B * XAB);
+
+          const T term = (PART1 + PART2) - PART3;
+
+          hypergeometric_pfq_base<T>::C.push_front(term);
+
+          A3 = A2;
+          A2 = A1;
+          A1 = hypergeometric_pfq_base<T>::C.front();
+        }
+
+        hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
+      }
+
+    private:
+      const T AP;
+      const T BP;
+      const T CP;
+
+      virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale<T>() * 1600.0F); }
+    };
+
+    template<class T> T luke_ccoef4_hypergeometric_0f1(const T& a, const T& x);
+    template<class T> T luke_ccoef1_hypergeometric_1f0(const T& a, const T& x);
+    template<class T> T luke_ccoef3_hypergeometric_1f1(const T& a, const T& b, const T& x);
+    template<class T> T luke_ccoef6_hypergeometric_1f2(const T& a, const T& b, const T& c, const T& x);
+    template<class T> T luke_ccoef2_hypergeometric_2f1(const T& a, const T& b, const T& c, const T& x);
+  }
+}
+
+template<class T>
+T examples::nr_006::luke_ccoef4_hypergeometric_0f1(const T& a, const T& x)
+{
+  const ccoef4_hypergeometric_0f1<T> hypergeometric_0f1_object(a, x, T(-20));
+
+  hypergeometric_0f1_object.ccoef();
+
+  return hypergeometric_0f1_object.series();
+}
+
+template<class T>
+T examples::nr_006::luke_ccoef1_hypergeometric_1f0(const T& a, const T& x)
+{
+  const ccoef1_hypergeometric_1f0<T> hypergeometric_1f0_object(a, x, T(-20));
+
+  hypergeometric_1f0_object.ccoef();
+
+  return hypergeometric_1f0_object.series();
+}
+
+template<class T>
+T examples::nr_006::luke_ccoef3_hypergeometric_1f1(const T& a, const T& b, const T& x)
+{
+  const ccoef3_hypergeometric_1f1<T> hypergeometric_1f1_object(a, b, x, T(-20));
+
+  hypergeometric_1f1_object.ccoef();
+
+  return hypergeometric_1f1_object.series();
+}
+
+template<class T>
+T examples::nr_006::luke_ccoef6_hypergeometric_1f2(const T& a, const T& b, const T& c, const T& x)
+{
+  const ccoef6_hypergeometric_1f2<T> hypergeometric_1f2_object(a, b, c, x, T(-20));
+
+  hypergeometric_1f2_object.ccoef();
+
+  return hypergeometric_1f2_object.series();
+}
+
+template<class T>
+T examples::nr_006::luke_ccoef2_hypergeometric_2f1(const T& a, const T& b, const T& c, const T& x)
+{
+  const ccoef2_hypergeometric_2f1<T> hypergeometric_2f1_object(a, b, c, x, T(-20));
+
+  hypergeometric_2f1_object.ccoef();
+
+  return hypergeometric_2f1_object.series();
+}
+
+int main()
+{
+  stopwatch<STD_CHRONO::high_resolution_clock> my_stopwatch;
+  float total_time = 0.0F;
+
+  std::vector<mp_type> hypergeometric_0f1_results(20U);
+  std::vector<mp_type> hypergeometric_1f0_results(20U);
+  std::vector<mp_type> hypergeometric_1f1_results(20U);
+  std::vector<mp_type> hypergeometric_2f1_results(20U);
+  std::vector<mp_type> hypergeometric_1f2_results(20U);
+
+  const mp_type a(mp_type(3) / 7);
+  const mp_type b(mp_type(2) / 3);
+  const mp_type c(mp_type(1) / 4);
+
+  std::int_least16_t i;
+
+  std::cout << "test hypergeometric_0f1." << std::endl;
+  i = 1U;
+  my_stopwatch.reset();
+
+  // Generate a table of values of Hypergeometric0F1.
+  // Compare with the Mathematica command:
+  // Table[N[HypergeometricPFQ[{}, {3/7}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
+  std::for_each(hypergeometric_0f1_results.begin(),
+                hypergeometric_0f1_results.end(),
+                [&i, &a](mp_type& new_value)
+                {
+                  const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
+
+                  new_value = examples::nr_006::luke_ccoef4_hypergeometric_0f1(a, x);
+
+                  ++i;
+                });
+
+  total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count();
+
+  // Print the values of Hypergeometric0F1.
+  std::for_each(hypergeometric_0f1_results.begin(),
+                hypergeometric_0f1_results.end(),
+                [](const mp_type& h)
+                {
+                  std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
+                });
+
+  std::cout << "test hypergeometric_1f0." << std::endl;
+  i = 1U;
+  my_stopwatch.reset();
+
+  // Generate a table of values of Hypergeometric1F0.
+  // Compare with the Mathematica command:
+  // Table[N[HypergeometricPFQ[{3/7}, {}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
+  std::for_each(hypergeometric_1f0_results.begin(),
+                hypergeometric_1f0_results.end(),
+                [&i, &a](mp_type& new_value)
+                {
+                  const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
+
+                  new_value = examples::nr_006::luke_ccoef1_hypergeometric_1f0(a, x);
+
+                  ++i;
+                });
+
+  total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count();
+
+  // Print the values of Hypergeometric1F0.
+  std::for_each(hypergeometric_1f0_results.begin(),
+                hypergeometric_1f0_results.end(),
+                [](const mp_type& h)
+                {
+                  std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
+                });
+
+  std::cout << "test hypergeometric_1f1." << std::endl;
+  i = 1U;
+  my_stopwatch.reset();
+
+  // Generate a table of values of Hypergeometric1F1.
+  // Compare with the Mathematica command:
+  // Table[N[HypergeometricPFQ[{3/7}, {2/3}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
+  std::for_each(hypergeometric_1f1_results.begin(),
+                hypergeometric_1f1_results.end(),
+                [&i, &a, &b](mp_type& new_value)
+                {
+                  const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
+
+                  new_value = examples::nr_006::luke_ccoef3_hypergeometric_1f1(a, b, x);
+
+                  ++i;
+                });
+
+  total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count();
+
+  // Print the values of Hypergeometric1F1.
+  std::for_each(hypergeometric_1f1_results.begin(),
+                hypergeometric_1f1_results.end(),
+                [](const mp_type& h)
+                {
+                  std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
+                });
+
+  std::cout << "test hypergeometric_1f2." << std::endl;
+  i = 1U;
+  my_stopwatch.reset();
+
+  // Generate a table of values of Hypergeometric1F2.
+  // Compare with the Mathematica command:
+  // Table[N[HypergeometricPFQ[{3/7}, {2/3, 1/4}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
+  std::for_each(hypergeometric_1f2_results.begin(),
+                hypergeometric_1f2_results.end(),
+                [&i, &a, &b, &c](mp_type& new_value)
+                {
+                  const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
+
+                  new_value = examples::nr_006::luke_ccoef6_hypergeometric_1f2(a, b, c, x);
+
+                  ++i;
+                });
+
+  total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count();
+
+  // Print the values of Hypergeometric1F2.
+  std::for_each(hypergeometric_1f2_results.begin(),
+                hypergeometric_1f2_results.end(),
+                [](const mp_type& h)
+                {
+                  std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
+                });
+
+  std::cout << "test hypergeometric_2f1." << std::endl;
+  i = 1U;
+  my_stopwatch.reset();
+
+  // Generate a table of values of Hypergeometric2F1.
+  // Compare with the Mathematica command:
+  // Table[N[HypergeometricPFQ[{3/7, 2/3}, {1/4}, -(i * EulerGamma)], 100], {i, 1, 20, 1}]
+  std::for_each(hypergeometric_2f1_results.begin(),
+                hypergeometric_2f1_results.end(),
+                [&i, &a, &b, &c](mp_type& new_value)
+                {
+                  const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
+
+                  new_value = examples::nr_006::luke_ccoef2_hypergeometric_2f1(a, b, c, x);
+
+                  ++i;
+                });
+
+  total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count();
+
+  // Print the values of Hypergeometric2F1.
+  std::for_each(hypergeometric_2f1_results.begin(),
+                hypergeometric_2f1_results.end(),
+                [](const mp_type& h)
+                {
+                  std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
+                });
+
+  std::cout << "Total execution time = " << std::setprecision(3) << total_time << "s" << std::endl;
+}
diff --git a/src/boost/libs/multiprecision/example/integer_examples.cpp b/src/boost/libs/multiprecision/example/integer_examples.cpp
new file mode 100644
index 00000000..5c2f8147
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/integer_examples.cpp
@@ -0,0 +1,232 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2012 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/multiprecision/cpp_int.hpp>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+
+// Includes Quickbook code snippets as comments.
+
+//[FAC1
+
+/*`
+In this simple example, we'll write a routine to print out all of the factorials
+which will fit into a 128-bit integer.  At the end of the routine we do some
+fancy iostream formatting of the results:
+*/
+/*=
+#include <boost/multiprecision/cpp_int.hpp>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+*/
+
+void print_factorials()
+{
+   using boost::multiprecision::cpp_int;
+   //
+   // Print all the factorials that will fit inside a 128-bit integer.
+   //
+   // Begin by building a big table of factorials, once we know just how
+   // large the largest is, we'll be able to "pretty format" the results.
+   //
+   // Calculate the largest number that will fit inside 128 bits, we could
+   // also have used numeric_limits<int128_t>::max() for this value:
+   cpp_int limit = (cpp_int(1) << 128) - 1;
+   //
+   // Our table of values:
+   std::vector<cpp_int> results;
+   //
+   // Initial values:
+   unsigned i = 1;
+   cpp_int factorial = 1;
+   //
+   // Cycle through the factorials till we reach the limit:
+   while(factorial < limit)
+   {
+      results.push_back(factorial);
+      ++i;
+      factorial *= i;
+   }
+   //
+   // Lets see how many digits the largest factorial was:
+   unsigned digits = results.back().str().size();
+   //
+   // Now print them out, using right justification, while we're at it
+   // we'll indicate the limit of each integer type, so begin by defining
+   // the limits for 16, 32, 64 etc bit integers:
+   cpp_int limits[] = {
+      (cpp_int(1) << 16) - 1,
+      (cpp_int(1) << 32) - 1,
+      (cpp_int(1) << 64) - 1,
+      (cpp_int(1) << 128) - 1,
+   };
+   std::string bit_counts[] = { "16", "32", "64", "128" };
+   unsigned current_limit = 0;
+   for(unsigned j = 0; j < results.size(); ++j)
+   {
+      if(limits[current_limit] < results[j])
+      {
+         std::string message = "Limit of " + bit_counts[current_limit] + " bit integers";
+         std::cout << std::setfill('.') << std::setw(digits+1) << std::right << message << std::setfill(' ') << std::endl;
+         ++current_limit;
+      }
+      std::cout << std::setw(digits + 1) << std::right << results[j] << std::endl;
+   }
+}
+
+/*`
+The output from this routine is:
+[pre
+                                       1
+                                       2
+                                       6
+                                      24
+                                     120
+                                     720
+                                    5040
+                                   40320
+................Limit of 16 bit integers
+                                  362880
+                                 3628800
+                                39916800
+                               479001600
+................Limit of 32 bit integers
+                              6227020800
+                             87178291200
+                           1307674368000
+                          20922789888000
+                         355687428096000
+                        6402373705728000
+                      121645100408832000
+                     2432902008176640000
+................Limit of 64 bit integers
+                    51090942171709440000
+                  1124000727777607680000
+                 25852016738884976640000
+                620448401733239439360000
+              15511210043330985984000000
+             403291461126605635584000000
+           10888869450418352160768000000
+          304888344611713860501504000000
+         8841761993739701954543616000000
+       265252859812191058636308480000000
+      8222838654177922817725562880000000
+    263130836933693530167218012160000000
+   8683317618811886495518194401280000000
+ 295232799039604140847618609643520000000
+]
+*/
+
+//]
+
+//[BITOPS
+
+/*`
+In this example we'll show how individual bits within an integer may be manipulated,
+we'll start with an often needed calculation of ['2[super n] - 1], which we could obviously
+implement like this:
+*/
+
+using boost::multiprecision::cpp_int;
+
+cpp_int b1(unsigned n)
+{
+   cpp_int r(1);
+   return (r << n) - 1;
+}
+
+/*`
+Calling:
+
+   std::cout << std::hex << std::showbase << b1(200) << std::endl;
+
+Yields as expected:
+
+[pre 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF]
+
+However, we could equally just set the n'th bit in the result, like this:
+*/
+
+cpp_int b2(unsigned n)
+{
+   cpp_int r(0);
+   return --bit_set(r, n);
+}
+
+/*`
+Note how the `bit_set` function sets the specified bit in its argument and then returns a reference to the result -
+which we can then simply decrement.  The result from a call to `b2` is the same as that to `b1`.
+
+We can equally test bits, so for example the n'th bit of the result returned from `b2` shouldn't be set
+unless we increment it first:
+
+   BOOST_ASSERT(!bit_test(b1(200), 200));     // OK
+   BOOST_ASSERT(bit_test(++b1(200), 200));    // OK
+
+And of course if we flip the n'th bit after increment, then we should get back to zero:
+
+   BOOST_ASSERT(!bit_flip(++b1(200), 200));   // OK
+*/
+
+//]
+
+int main()
+{
+   print_factorials();
+
+   std::cout << std::hex << std::showbase << b1(200) << std::endl;
+   std::cout << std::hex << std::showbase << b2(200) << std::endl;
+   BOOST_ASSERT(!bit_test(b1(200), 200));  // OK
+   BOOST_ASSERT(bit_test(++b1(200), 200));    // OK
+   BOOST_ASSERT(!bit_flip(++b1(200), 200));   // OK
+   return 0;
+}
+
+/*
+
+Program output:
+
+                                       1
+                                       2
+                                       6
+                                      24
+                                     120
+                                     720
+                                    5040
+                                   40320
+................Limit of 16 bit integers
+                                  362880
+                                 3628800
+                                39916800
+                               479001600
+................Limit of 32 bit integers
+                              6227020800
+                             87178291200
+                           1307674368000
+                          20922789888000
+                         355687428096000
+                        6402373705728000
+                      121645100408832000
+                     2432902008176640000
+................Limit of 64 bit integers
+                    51090942171709440000
+                  1124000727777607680000
+                 25852016738884976640000
+                620448401733239439360000
+              15511210043330985984000000
+             403291461126605635584000000
+           10888869450418352160768000000
+          304888344611713860501504000000
+         8841761993739701954543616000000
+       265252859812191058636308480000000
+      8222838654177922817725562880000000
+    263130836933693530167218012160000000
+   8683317618811886495518194401280000000
+ 295232799039604140847618609643520000000
+ 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
+ 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
+ */
diff --git a/src/boost/libs/multiprecision/example/logged_adaptor.cpp b/src/boost/libs/multiprecision/example/logged_adaptor.cpp
new file mode 100644
index 00000000..f204cf39
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/logged_adaptor.cpp
@@ -0,0 +1,119 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2013 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+//[logged_adaptor
+
+#include <boost/multiprecision/mpfi.hpp>
+#include <boost/multiprecision/logged_adaptor.hpp>
+#include <iostream>
+#include <iomanip>
+//
+// Begin by overloading log_postfix_event so we can capture each arithmetic event as it happens:
+//
+namespace boost{ namespace multiprecision{
+
+template <unsigned D>
+inline void log_postfix_event(const mpfi_float_backend<D>& val, const char* event_description)
+{
+   // Print out the (relative) diameter of the interval:
+   using namespace boost::multiprecision;
+   number<mpfr_float_backend<D> > diam;
+   mpfi_diam(diam.backend().data(), val.data());
+   std::cout << "Diameter was " << diam << " after operation: " << event_description << std::endl;
+}
+template <unsigned D, class T>
+inline void log_postfix_event(const mpfi_float_backend<D>&, const T&, const char* event_description)
+{
+   // This version is never called in this example.
+}
+
+}}
+
+
+int main()
+{
+   using namespace boost::multiprecision;
+   typedef number<logged_adaptor<mpfi_float_backend<17> > > logged_type;
+   //
+   // Test case deliberately introduces cancellation error, relative size of interval
+   // gradually gets larger after each operation:
+   //
+   logged_type a = 1;
+   a /= 10;
+
+   for(unsigned i = 0; i < 13; ++i)
+   {
+      logged_type b = a * 9;
+      b /= 10;
+      a -= b;
+   }
+   std::cout << "Final value was: " << a << std::endl;
+   return 0;
+}
+
+//]
+
+/*
+//[logged_adaptor_output
+
+Diameter was nan after operation: Default construct
+Diameter was 0 after operation: Assignment from arithmetic type
+Diameter was 4.33681e-18 after operation: /=
+Diameter was nan after operation: Default construct
+Diameter was 7.70988e-18 after operation: *
+Diameter was 9.63735e-18 after operation: /=
+Diameter was 1.30104e-16 after operation: -=
+Diameter was nan after operation: Default construct
+Diameter was 1.30104e-16 after operation: *
+Diameter was 1.38537e-16 after operation: /=
+Diameter was 2.54788e-15 after operation: -=
+Diameter was nan after operation: Default construct
+Diameter was 2.54788e-15 after operation: *
+Diameter was 2.54863e-15 after operation: /=
+Diameter was 4.84164e-14 after operation: -=
+Diameter was nan after operation: Default construct
+Diameter was 4.84164e-14 after operation: *
+Diameter was 4.84221e-14 after operation: /=
+Diameter was 9.19962e-13 after operation: -=
+Diameter was nan after operation: Default construct
+Diameter was 9.19962e-13 after operation: *
+Diameter was 9.19966e-13 after operation: /=
+Diameter was 1.74793e-11 after operation: -=
+Diameter was nan after operation: Default construct
+Diameter was 1.74793e-11 after operation: *
+Diameter was 1.74793e-11 after operation: /=
+Diameter was 3.32107e-10 after operation: -=
+Diameter was nan after operation: Default construct
+Diameter was 3.32107e-10 after operation: *
+Diameter was 3.32107e-10 after operation: /=
+Diameter was 6.31003e-09 after operation: -=
+Diameter was nan after operation: Default construct
+Diameter was 6.31003e-09 after operation: *
+Diameter was 6.31003e-09 after operation: /=
+Diameter was 1.19891e-07 after operation: -=
+Diameter was nan after operation: Default construct
+Diameter was 1.19891e-07 after operation: *
+Diameter was 1.19891e-07 after operation: /=
+Diameter was 2.27792e-06 after operation: -=
+Diameter was nan after operation: Default construct
+Diameter was 2.27792e-06 after operation: *
+Diameter was 2.27792e-06 after operation: /=
+Diameter was 4.32805e-05 after operation: -=
+Diameter was nan after operation: Default construct
+Diameter was 4.32805e-05 after operation: *
+Diameter was 4.32805e-05 after operation: /=
+Diameter was 0.00082233 after operation: -=
+Diameter was nan after operation: Default construct
+Diameter was 0.00082233 after operation: *
+Diameter was 0.00082233 after operation: /=
+Diameter was 0.0156243 after operation: -=
+Diameter was nan after operation: Default construct
+Diameter was 0.0156243 after operation: *
+Diameter was 0.0156243 after operation: /=
+Diameter was 0.296861 after operation: -=
+Final value was: {8.51569e-15,1.14843e-14}
+
+//]
+*/
diff --git a/src/boost/libs/multiprecision/example/mixed_integer_arithmetic.cpp b/src/boost/libs/multiprecision/example/mixed_integer_arithmetic.cpp
new file mode 100644
index 00000000..c1ab5b5c
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/mixed_integer_arithmetic.cpp
@@ -0,0 +1,55 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2012 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+//
+// Compare arithmetic results using fixed_int to GMP results.
+//
+
+#ifdef _MSC_VER
+#  define _SCL_SECURE_NO_WARNINGS
+#endif
+
+//[mixed_eg
+#include <boost/multiprecision/cpp_int.hpp>
+
+int main()
+{
+   using namespace boost::multiprecision;
+
+   boost::uint64_t i = (std::numeric_limits<boost::uint64_t>::max)();
+   boost::uint64_t j = 1;
+
+   uint128_t ui128;
+   uint256_t ui256;
+   //
+   // Start by performing arithmetic on 64-bit integers to yield 128-bit results:
+   //
+   std::cout << std::hex << std::showbase << i << std::endl;
+   std::cout << std::hex << std::showbase << add(ui128, i, j) << std::endl;
+   std::cout << std::hex << std::showbase << multiply(ui128, i, i) << std::endl;
+   //
+   // The try squaring a 128-bit integer to yield a 256-bit result:
+   //
+   ui128 = (std::numeric_limits<uint128_t>::max)();
+   std::cout << std::hex << std::showbase << multiply(ui256, ui128, ui128) << std::endl;
+
+   return 0;
+}
+//]
+
+/* 
+
+Program output:
+
+//[mixed_output
+
+0xffffffffffffffff
+0x10000000000000000
+0xFFFFFFFFFFFFFFFE0000000000000001
+0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFE00000000000000000000000000000001
+
+//]
+*/
+
diff --git a/src/boost/libs/multiprecision/example/mpc_examples.cpp b/src/boost/libs/multiprecision/example/mpc_examples.cpp
new file mode 100644
index 00000000..30f9a0bf
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/mpc_examples.cpp
@@ -0,0 +1,121 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2018 Nick Thompson. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+/*`This example demonstrates the usage of the MPC backend for multiprecision complex numbers.
+In the following, we will show how using MPC backend allows for the same operations as the C++ standard library complex numbers.
+*/
+
+//[mpc_eg
+#include <iostream>
+#include <complex>
+#include <boost/multiprecision/mpc.hpp>
+
+template<class Complex>
+void complex_number_examples()
+{
+    Complex z1{0, 1};
+    std::cout << std::setprecision(std::numeric_limits<typename Complex::value_type>::digits10);
+    std::cout << std::scientific << std::fixed;
+    std::cout << "Print a complex number: " << z1 << std::endl;
+    std::cout << "Square it             : " << z1*z1 << std::endl;
+    std::cout << "Real part             : " << z1.real() << " = " << real(z1) << std::endl;
+    std::cout << "Imaginary part        : " << z1.imag() << " = " << imag(z1) << std::endl;
+    using std::abs;
+    std::cout << "Absolute value        : " << abs(z1) << std::endl;
+    std::cout << "Argument              : " << arg(z1) << std::endl;
+    std::cout << "Norm                  : " << norm(z1) << std::endl;
+    std::cout << "Complex conjugate     : " << conj(z1) << std::endl;
+    std::cout << "Projection onto Riemann sphere: " <<  proj(z1) << std::endl;
+    typename Complex::value_type r = 1;
+    typename Complex::value_type theta = 0.8;
+    using std::polar;
+    std::cout << "Polar coordinates (phase = 0)    : " << polar(r) << std::endl;
+    std::cout << "Polar coordinates (phase !=0)    : " << polar(r, theta) << std::endl;
+
+    std::cout << "\nElementary special functions:\n";
+    using std::exp;
+    std::cout << "exp(z1) = " << exp(z1) << std::endl;
+    using std::log;
+    std::cout << "log(z1) = " << log(z1) << std::endl;
+    using std::log10;
+    std::cout << "log10(z1) = " << log10(z1) << std::endl;
+    using std::pow;
+    std::cout << "pow(z1, z1) = " << pow(z1, z1) << std::endl;
+    using std::sqrt;
+    std::cout << "Take its square root  : " << sqrt(z1) << std::endl;
+    using std::sin;
+    std::cout << "sin(z1) = " << sin(z1) << std::endl;
+    using std::cos;
+    std::cout << "cos(z1) = " << cos(z1) << std::endl;
+    using std::tan;
+    std::cout << "tan(z1) = " << tan(z1) << std::endl;
+    using std::asin;
+    std::cout << "asin(z1) = " << asin(z1) << std::endl;
+    using std::acos;
+    std::cout << "acos(z1) = " << acos(z1) << std::endl;
+    using std::atan;
+    std::cout << "atan(z1) = " << atan(z1) << std::endl;
+    using std::sinh;
+    std::cout << "sinh(z1) = " << sinh(z1) << std::endl;
+    using std::cosh;
+    std::cout << "cosh(z1) = " << cosh(z1) << std::endl;
+    using std::tanh;
+    std::cout << "tanh(z1) = " << tanh(z1) << std::endl;
+    using std::asinh;
+    std::cout << "asinh(z1) = " << asinh(z1) << std::endl;
+    using std::acosh;
+    std::cout << "acosh(z1) = " << acosh(z1) << std::endl;
+    using std::atanh;
+    std::cout << "atanh(z1) = " << atanh(z1) << std::endl;
+}
+
+int main()
+{
+    std::cout << "First, some operations we usually perform with std::complex:\n";
+    complex_number_examples<std::complex<double>>();
+    std::cout << "\nNow the same operations performed using the MPC backend:\n";
+    complex_number_examples<boost::multiprecision::mpc_complex_50>();
+
+    return 0;
+}
+//]
+
+/*
+
+//[mpc_out
+
+Print a complex number: (0.00000000000000000000000000000000000000000000000000,1.00000000000000000000000000000000000000000000000000)
+Square it             : -1.00000000000000000000000000000000000000000000000000
+Real part             : 0.00000000000000000000000000000000000000000000000000 = 0.00000000000000000000000000000000000000000000000000
+Imaginary part        : 1.00000000000000000000000000000000000000000000000000 = 1.00000000000000000000000000000000000000000000000000
+Absolute value        : 1.00000000000000000000000000000000000000000000000000
+Argument              : 1.57079632679489661923132169163975144209858469968755
+Norm                  : 1.00000000000000000000000000000000000000000000000000
+Complex conjugate     : (0.00000000000000000000000000000000000000000000000000,-1.00000000000000000000000000000000000000000000000000)
+Projection onto Riemann sphere: (0.00000000000000000000000000000000000000000000000000,1.00000000000000000000000000000000000000000000000000)
+Polar coordinates (phase = 0)    : 1.00000000000000000000000000000000000000000000000000
+Polar coordinates (phase !=0)    : (0.69670670934716538906374002277244853473117519431538,0.71735609089952279256716781570337728075604730751255)
+
+Elementary special functions:
+exp(z1) = (0.54030230586813971740093660744297660373231042061792,0.84147098480789650665250232163029899962256306079837)
+log(z1) = (0.00000000000000000000000000000000000000000000000000,1.57079632679489661923132169163975144209858469968755)
+log10(z1) = (0.00000000000000000000000000000000000000000000000000,0.68218817692092067374289181271567788510506374186196)
+pow(z1, z1) = 0.20787957635076190854695561983497877003387784163177
+Take its square root  : (0.70710678118654752440084436210484903928483593768847,0.70710678118654752440084436210484903928483593768847)
+sin(z1) = (0.00000000000000000000000000000000000000000000000000,1.17520119364380145688238185059560081515571798133410)
+cos(z1) = 1.54308063481524377847790562075706168260152911236587
+tan(z1) = (0.00000000000000000000000000000000000000000000000000,0.76159415595576488811945828260479359041276859725794)
+asin(z1) = (0.00000000000000000000000000000000000000000000000000,0.88137358701954302523260932497979230902816032826163)
+acos(z1) = (1.57079632679489661923132169163975144209858469968755,-0.88137358701954302523260932497979230902816032826163)
+atan(z1) = (0.00000000000000000000000000000000000000000000000000,inf)
+sinh(z1) = (0.00000000000000000000000000000000000000000000000000,0.84147098480789650665250232163029899962256306079837)
+cosh(z1) = 0.54030230586813971740093660744297660373231042061792
+tanh(z1) = (0.00000000000000000000000000000000000000000000000000,1.55740772465490223050697480745836017308725077238152)
+asinh(z1) = (0.00000000000000000000000000000000000000000000000000,1.57079632679489661923132169163975144209858469968755)
+acosh(z1) = (0.88137358701954302523260932497979230902816032826163,1.57079632679489661923132169163975144209858469968755)
+atanh(z1) = (0.00000000000000000000000000000000000000000000000000,0.78539816339744830961566084581987572104929234984378)
+
+//]
+*/
diff --git a/src/boost/libs/multiprecision/example/mpfi_snips.cpp b/src/boost/libs/multiprecision/example/mpfi_snips.cpp
new file mode 100644
index 00000000..b36e22b5
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/mpfi_snips.cpp
@@ -0,0 +1,42 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2011 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+//[mpfi_eg
+#include <boost/multiprecision/mpfi.hpp>
+#include <boost/math/special_functions/gamma.hpp>
+#include <iostream>
+
+int main()
+{
+   using namespace boost::multiprecision;
+
+   // Operations at variable precision and no numeric_limits support:
+   mpfi_float a = 2;
+   mpfi_float::default_precision(1000);
+   std::cout << mpfi_float::default_precision() << std::endl;
+   std::cout << sqrt(a) << std::endl; // print root-2
+
+   // Operations at fixed precision and full numeric_limits support:
+   mpfi_float_100 b = 2;
+   std::cout << std::numeric_limits<mpfi_float_100>::digits << std::endl;
+   // We can use any C++ std lib function:
+   std::cout << log(b) << std::endl; // print log(2)
+
+   // Access the underlying data:
+   mpfi_t r;
+   mpfi_init(r);
+   mpfi_set(r, b.backend().data());
+
+   // Construct some explicit intervals and perform set operations:
+   mpfi_float_50 i1(1, 2), i2(1.5, 2.5);
+   std::cout << intersect(i1, i2) << std::endl;
+   std::cout << hull(i1, i2) << std::endl;
+   std::cout << overlap(i1, i2) << std::endl;
+   std::cout << subset(i1, i2) << std::endl;
+   mpfi_clear(r);
+   return 0;
+}
+//]
+
diff --git a/src/boost/libs/multiprecision/example/mpfr_precision.cpp b/src/boost/libs/multiprecision/example/mpfr_precision.cpp
new file mode 100644
index 00000000..eb559cbe
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/mpfr_precision.cpp
@@ -0,0 +1,252 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2018 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+//[mpfr_variable
+
+/*`
+This example illustrates the use of variable-precision arithmetic with
+the `mpfr_float` number type.  We'll calculate the median of the
+beta distribution to an absurdly high precision and compare the
+accuracy and times taken for various methods.  That is, we want
+to calculate the value of `x` for which ['I[sub x](a, b) = 0.5].
+
+Ultimately we'll use Newtons method and set the precision of
+mpfr_float to have just enough digits at each iteration.
+
+The full source of the this program is in [@../../example/mpfr_precision.cpp]
+
+We'll skip over the #includes and using declations, and go straight to
+some support code, first off a simple stopwatch for performance measurement:
+
+*/
+
+//=template <class clock_type>
+//=struct stopwatch { /*details \*/ };
+
+/*`
+We'll use `stopwatch<std::chono::high_resolution_clock>` as our performance measuring device.
+
+We also have a small utility class for controlling the current precision of mpfr_float:
+
+   struct scoped_precision
+   {
+      unsigned p;
+      scoped_precision(unsigned new_p) : p(mpfr_float::default_precision())
+      {
+         mpfr_float::default_precision(new_p);
+      }
+      ~scoped_precision()
+      {
+         mpfr_float::default_precision(p);
+      }
+   };
+
+*/
+//<-
+#include <boost/multiprecision/mpfr.hpp>
+#include <boost/math/special_functions/beta.hpp>
+#include <boost/math/special_functions/relative_difference.hpp>
+#include <iostream>
+#include <chrono>
+
+using boost::multiprecision::mpfr_float;
+using boost::math::ibeta_inv;
+using namespace boost::math::policies;
+
+template <class clock_type>
+struct stopwatch
+{
+public:
+   typedef typename clock_type::duration duration_type;
+
+   stopwatch() : m_start(clock_type::now()) { }
+
+   stopwatch(const stopwatch& other) : m_start(other.m_start) { }
+
+   stopwatch& operator=(const stopwatch& other)
+   {
+      m_start = other.m_start;
+      return *this;
+   }
+
+   ~stopwatch() { }
+
+   float elapsed() const
+   {
+      return float(std::chrono::nanoseconds((clock_type::now() - m_start)).count()) / 1e9f;
+   }
+
+   void reset()
+   {
+      m_start = clock_type::now();
+   }
+
+private:
+   typename clock_type::time_point m_start;
+};
+
+struct scoped_precision
+{
+   unsigned p;
+   scoped_precision(unsigned new_p) : p(mpfr_float::default_precision())
+   {
+      mpfr_float::default_precision(new_p);
+   }
+   ~scoped_precision()
+   {
+      mpfr_float::default_precision(p);
+   }
+};
+//->
+
+/*`
+We'll begin with a reference method that simply calls the Boost.Math function `ibeta_inv` and uses the
+full working precision of the arguments throughout.  Our reference function takes 3 arguments:
+
+* The 2 parameters `a` and `b` of the beta distribution, and
+* The number of decimal digits precision to achieve in the result.
+
+We begin by setting the default working precision to that requested, and then, since we don't know where
+our arguments `a` and `b` have been or what precision they have, we make a copy of them - note that since
+copying also copies the precision as well as the value, we have to set the precision expicitly with a
+second argument to the copy.  Then we can simply return the result of `ibeta_inv`:
+*/
+mpfr_float beta_distribution_median_method_1(mpfr_float const& a_, mpfr_float const& b_, unsigned digits10)
+{
+   scoped_precision sp(digits10);
+   mpfr_float half(0.5), a(a_, digits10), b(b_, digits10);
+   return ibeta_inv(a, b, half);
+}
+/*`
+You be wondering why we needed to change the precision of our variables `a` and `b` as well as setting the default -
+there are in fact two ways in which this can go wrong if we don't do that:
+
+* The variables have too much precision - this will cause all arithmetic operations involving those types to be
+promoted to the higher precision wasting precious calculation time.
+* The variables have too little precision - this will cause expressions involving only those variables to be
+calculated at the lower precision - for example if we calculate `exp(a)` internally, this will be evaluated at
+the precision of `a`, and not the current default.
+
+Since our reference method carries out all calculations at the full precision requested, an obvious refinement
+would be to calculate a first approximation to `double` precision and then to use Newton steps to refine it further.
+
+Our function begins the same as before: set the new default precision and then make copies of our arguments
+at the correct precision.  We then call `ibeta_inv` with all double precision arguments, promote the result
+to an `mpfr_float` and perform Newton steps to obtain the result.  Note that our termination condition is somewhat
+cude: we simply assume that we have approximately 14 digits correct from the double-precision approximation and
+that the precision doubles with each step.  We also cheat, and use an internal Boost.Math function that calculates
+['I[sub x](a, b)] and it's derivative in one go:
+
+*/
+mpfr_float beta_distribution_median_method_2(mpfr_float const& a_, mpfr_float const& b_, unsigned digits10)
+{
+   scoped_precision sp(digits10);
+   mpfr_float half(0.5), a(a_, digits10), b(b_, digits10);
+   mpfr_float guess = ibeta_inv((double)a, (double)b, 0.5);
+   unsigned current_digits = 14;
+   mpfr_float f, f1;
+   while (current_digits < digits10)
+   {
+      f = boost::math::detail::ibeta_imp(a, b, guess, boost::math::policies::policy<>(), false, true, &f1) - half;
+      guess -= f / f1;
+      current_digits *= 2;
+   }
+   return guess;
+}
+/*`
+Before we refine the method further, it might be wise to take stock and see how method's 1 and 2 compare.
+We'll ask them both for 1500 digit precision, and compare against the value produced by `ibeta_inv` at 1700 digits.
+Here's what the results look like:
+
+[pre
+Method 1 time = 0.611647
+Relative error: 2.99991e-1501
+Method 2 time = 0.646746
+Relative error: 7.55843e-1501
+]
+
+Clearly they are both equally accurate, but Method 1 is actually faster and our plan for improved performance
+hasn't actually worked.  It turns out that we're not actually comparing like with like, because `ibeta_inv` uses
+Halley iteration internally which churns out more digits of precision rather more rapidly than Newton iteration.
+So the time we save by refining an initial `double` approximation, then loose it again by taking more iterations
+to get to the result.
+
+Time for a more refined approach.  It follows the same form as Method 2, but now we set the working precision
+within the Newton iteration loop, to just enough digits to cover the expected precision at each step.  That means
+we also create new copies of our arguments at the correct precision within the loop, and likewise change the precision
+of the current `guess` each time through:
+
+*/
+
+mpfr_float beta_distribution_median_method_3(mpfr_float const& a_, mpfr_float const& b_, unsigned digits10)
+{
+   mpfr_float guess = ibeta_inv((double)a_, (double)b_, 0.5);
+   unsigned current_digits = 14;
+   mpfr_float f(0, current_digits), f1(0, current_digits), delta(1);
+   while (current_digits < digits10)
+   {
+      current_digits *= 2;
+      scoped_precision sp((std::min)(current_digits, digits10));
+      mpfr_float a(a_, mpfr_float::default_precision()), b(b_, mpfr_float::default_precision());
+      guess.precision(mpfr_float::default_precision());
+      f = boost::math::detail::ibeta_imp(a, b, guess, boost::math::policies::policy<>(), false, true, &f1) - 0.5f;
+      guess -= f / f1;
+   }
+   return guess;
+}
+
+/*`
+The new performance results look much more promising:
+
+[pre
+Method 1 time = 0.591244
+Relative error: 2.99991e-1501
+Method 2 time = 0.622679
+Relative error: 7.55843e-1501
+Method 3 time = 0.143393
+Relative error: 4.03898e-1501
+]
+
+This time we're 4x faster than `ibeta_inv`, and no doubt that could be improved a little more by carefully
+optimising the number of iterations and the method (Halley vs Newton) taken.
+
+Finally, here's the driver code for the above methods:
+
+*/
+
+int main()
+{
+   try {
+      mpfr_float a(10), b(20);
+
+      mpfr_float true_value = beta_distribution_median_method_1(a, b, 1700);
+
+      stopwatch<std::chrono::high_resolution_clock> my_stopwatch;
+
+      mpfr_float v1 = beta_distribution_median_method_1(a, b, 1500);
+      float hp_time = my_stopwatch.elapsed();
+      std::cout << "Method 1 time = " << hp_time << std::endl;
+      std::cout << "Relative error: " << boost::math::relative_difference(v1, true_value) << std::endl;
+
+      my_stopwatch.reset();
+      mpfr_float v2 = beta_distribution_median_method_2(a, b, 1500);
+      hp_time = my_stopwatch.elapsed();
+      std::cout << "Method 2 time = " << hp_time << std::endl;
+      std::cout << "Relative error: " << boost::math::relative_difference(v2, true_value) << std::endl;
+
+      my_stopwatch.reset();
+      mpfr_float v3 = beta_distribution_median_method_3(a, b, 1500);
+      hp_time = my_stopwatch.elapsed();
+      std::cout << "Method 3 time = " << hp_time << std::endl;
+      std::cout << "Relative error: " << boost::math::relative_difference(v3, true_value) << std::endl;
+   }
+   catch (const std::exception& e)
+   {
+      std::cout << "Found exception with message: " << e.what() << std::endl;
+   }
+   return 0;
+}
+//]
+
diff --git a/src/boost/libs/multiprecision/example/mpfr_snips.cpp b/src/boost/libs/multiprecision/example/mpfr_snips.cpp
new file mode 100644
index 00000000..e0ac6aeb
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/mpfr_snips.cpp
@@ -0,0 +1,40 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2011 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+//[mpfr_eg
+#include <boost/multiprecision/mpfr.hpp>
+#include <boost/math/special_functions/gamma.hpp>
+#include <iostream>
+
+int main()
+{
+   using namespace boost::multiprecision;
+
+   // Operations at variable precision and no numeric_limits support:
+   mpfr_float a = 2;
+   mpfr_float::default_precision(1000);
+   std::cout << mpfr_float::default_precision() << std::endl;
+   std::cout << sqrt(a) << std::endl; // print root-2
+
+   // Operations at fixed precision and full numeric_limits support:
+   mpfr_float_100 b = 2;
+   std::cout << std::numeric_limits<mpfr_float_100>::digits << std::endl;
+   // We can use any C++ std lib function:
+   std::cout << log(b) << std::endl; // print log(2)
+   // We can also use any function from Boost.Math:
+   std::cout << boost::math::tgamma(b) << std::endl;
+   // These even work when the argument is an expression template:
+   std::cout << boost::math::tgamma(b * b) << std::endl;
+
+   // Access the underlying data:
+   mpfr_t r;
+   mpfr_init(r);
+   mpfr_set(r, b.backend().data(), GMP_RNDN);
+   mpfr_clear(r);
+   return 0;
+}
+//]
+
+
diff --git a/src/boost/libs/multiprecision/example/numeric_limits_snips.cpp b/src/boost/libs/multiprecision/example/numeric_limits_snips.cpp
new file mode 100644
index 00000000..e114f0a2
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/numeric_limits_snips.cpp
@@ -0,0 +1,470 @@
+// Copyright Paul A. Bristow 2013
+// Copyright John Maddock 2013
+// Copyright Christopher Kormanyos
+
+// 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)
+
+// Examples of numeric_limits usage as snippets for multiprecision documentation.
+
+// Includes text as Quickbook comments.
+
+#include <iostream>
+#include <iomanip>
+#include <string>
+#include <sstream>
+#include <limits> // numeric_limits
+#include <iomanip>
+#include <locale>
+#include <boost/assert.hpp>
+
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/special_functions/nonfinite_num_facets.hpp>
+
+#include <boost/math/special_functions/factorials.hpp>
+#include <boost/math/special_functions/next.hpp>
+#include <boost/math/tools/precision.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp> // is decimal.
+#include <boost/multiprecision/cpp_bin_float.hpp> // is binary.
+
+#define BOOST_TEST_MAIN
+#include <boost/test/unit_test.hpp> // Boost.Test
+#include <boost/test/floating_point_comparison.hpp>
+
+static long double const log10Two = 0.30102999566398119521373889472449L; // log10(2.)
+
+template <typename T>
+int max_digits10()
+{
+   int significand_digits = std::numeric_limits<T>::digits;
+  // BOOST_CONSTEXPR_OR_CONST int significand_digits = std::numeric_limits<T>::digits;
+   return static_cast<int>(ceil(1 + significand_digits * log10Two));
+} // template <typename T> int max_digits10()
+
+// Used to test max_digits10<>() function below.
+//#define BOOST_NO_CXX11_NUMERIC_LIMITS
+
+BOOST_AUTO_TEST_CASE(test_numeric_limits_snips)
+{
+#if !(defined(CI_SUPPRESS_KNOWN_ISSUES) && defined(BOOST_MSVC) && (BOOST_MSVC == 1600))
+  try
+  {
+
+// Example of portable way to get `std::numeric_limits<T>::max_digits10`.
+//[max_digits10_1
+
+/*`For example, to be portable (including obselete platforms) for type `T` where `T` may be:
+ `float`, `double`, `long double`, `128-bit quad type`, `cpp_bin_float_50` ...
+*/
+
+  typedef float T;
+
+#if defined BOOST_NO_CXX11_NUMERIC_LIMITS
+   // No max_digits10 implemented.
+    std::cout.precision(max_digits10<T>());
+#else
+  #if(_MSC_VER <= 1600)
+   //  Wrong value for std::numeric_limits<float>::max_digits10.
+    std::cout.precision(max_digits10<T>());
+  #else // Use the C++11 max_digits10.
+     std::cout.precision(std::numeric_limits<T>::max_digits10);
+  #endif
+#endif
+
+  std::cout << "std::cout.precision(max_digits10) = " << std::cout.precision() << std::endl; // 9
+
+  double x = 1.2345678901234567889;
+
+  std::cout << "x = " << x << std::endl; //
+
+/*`which should output:
+
+  std::cout.precision(max_digits10) = 9
+  x = 1.23456789
+*/
+
+//] [/max_digits10_1]
+
+  {
+//[max_digits10_2
+
+  double write = 2./3; // Any arbitrary value that cannot be represented exactly.
+  double read = 0;
+  std::stringstream s;
+  s.precision(std::numeric_limits<double>::digits10); // or `float64_t` for 64-bit IEE754 double.
+  s << write;
+  s >> read;
+  if(read != write)
+  {
+    std::cout <<  std::setprecision(std::numeric_limits<double>::digits10)
+      << read << " != " << write << std::endl;
+  }
+
+//] [/max_digits10_2]
+  // 0.666666666666667 != 0.666666666666667
+  }
+
+  {
+//[max_digits10_3
+
+  double pi = boost::math::double_constants::pi;
+  std::cout.precision(std::numeric_limits<double>::max_digits10);
+  std::cout << pi << std::endl; // 3.1415926535897931
+
+//] [/max_digits10_3]
+  }
+  {
+//[max_digits10_4
+/*`and similarly for a much higher precision type:
+*/
+
+  using namespace boost::multiprecision;
+
+  typedef number<cpp_dec_float<50> > cpp_dec_float_50; // 50 decimal digits.
+
+  using boost::multiprecision::cpp_dec_float_50;
+
+  cpp_dec_float_50 pi = boost::math::constants::pi<cpp_dec_float_50>();
+  std::cout.precision(std::numeric_limits<cpp_dec_float_50>::max_digits10);
+  std::cout << pi << std::endl;
+  // 3.141592653589793238462643383279502884197169399375105820974944592307816406
+//] [/max_digits10_4]
+  }
+
+  {
+//[max_digits10_5
+
+  for (int i = 2; i < 15; i++)
+  {
+    std::cout << std::setw(std::numeric_limits<int>::max_digits10)
+      << boost::math::factorial<double>(i)  << std::endl;
+  }
+
+//] [/max_digits10_5]
+  }
+
+  }
+  catch(std::exception ex)
+  {
+    std::cout << "Caught Exception " << ex.what() << std::endl;
+  }
+
+  {
+//[max_digits10_6
+
+  typedef double T;
+
+  bool denorm = std::numeric_limits<T>::denorm_min() < (std::numeric_limits<T>::min)();
+  BOOST_ASSERT(denorm);
+
+//] [/max_digits10_6]
+  }
+
+  {
+    unsigned char c = 255;
+    std::cout << "char c = " << (int)c << std::endl;
+  }
+
+  {
+//[digits10_1
+    std::cout
+      << std::setw(std::numeric_limits<short>::digits10 +1 +1) // digits10+1, and +1 for sign.
+      << std::showpos << (std::numeric_limits<short>::max)() // +32767
+      << std::endl
+      << std::setw(std::numeric_limits<short>::digits10 +1 +1)
+      << (std::numeric_limits<short>::min)() << std::endl;   // -32767
+//] [/digits10_1]
+  }
+
+  {
+//[digits10_2
+    std::cout
+      << std::setw(std::numeric_limits<unsigned short>::digits10 +1 +1) // digits10+1, and +1 for sign.
+      << std::showpos << (std::numeric_limits<unsigned short>::max)() //  65535
+      << std::endl
+      << std::setw(std::numeric_limits<unsigned short>::digits10 +1 +1) // digits10+1, and +1 for sign.
+      << (std::numeric_limits<unsigned short>::min)() << std::endl;   //      0
+//] [/digits10_2]
+  }
+
+  std::cout <<std::noshowpos << std::endl;
+
+  {
+//[digits10_3
+  std::cout.precision(std::numeric_limits<double>::max_digits10);
+  double d =  1e15;
+  double dp1 = d+1;
+  std::cout << d << "\n" << dp1 << std::endl;
+  // 1000000000000000
+  // 1000000000000001
+  std::cout <<  dp1 - d << std::endl; // 1
+//] [/digits10_3]
+  }
+
+  {
+//[digits10_4
+  std::cout.precision(std::numeric_limits<double>::max_digits10);
+  double d =  1e16;
+  double dp1 = d+1;
+  std::cout << d << "\n" << dp1 << std::endl;
+  // 10000000000000000
+  // 10000000000000000
+    std::cout << dp1 - d << std::endl; // 0 !!!
+//] [/digits10_4]
+  }
+
+  {
+//[epsilon_1
+  std::cout.precision(std::numeric_limits<double>::max_digits10);
+  double d = 1.;
+  double eps = std::numeric_limits<double>::epsilon();
+  double dpeps = d+eps;
+  std::cout << std::showpoint // Ensure all trailing zeros are shown.
+    << d << "\n"           // 1.0000000000000000
+    << dpeps << std::endl; // 2.2204460492503131e-016
+  std::cout << dpeps - d   // 1.0000000000000002
+    << std::endl;
+//] [epsilon_1]
+  }
+
+  {
+//[epsilon_2
+  double one = 1.;
+  double nad = boost::math::float_next(one);
+  std::cout << nad << "\n"  //  1.0000000000000002
+    << nad - one // 2.2204460492503131e-016
+    << std::endl;
+//] [epsilon_2]
+  }
+  {
+//[epsilon_3
+  std::cout.precision(std::numeric_limits<double>::max_digits10);
+  double d = 1.;
+  double eps = std::numeric_limits<double>::epsilon();
+  double dpeps = d + eps/2;
+
+  std::cout << std::showpoint // Ensure all trailing zeros are shown.
+    << dpeps << "\n"       // 1.0000000000000000
+    << eps/2 << std::endl; // 1.1102230246251565e-016
+  std::cout << dpeps - d   // 0.00000000000000000
+    << std::endl;
+//] [epsilon_3]
+  }
+
+  {
+    typedef double RealType;
+//[epsilon_4
+/*`A tolerance might be defined using this version of epsilon thus:
+*/
+    RealType tolerance = boost::math::tools::epsilon<RealType>() * 2;
+//] [epsilon_4]
+  }
+
+  {
+//[digits10_5
+    -(std::numeric_limits<double>::max)() == std::numeric_limits<double>::lowest();
+//] [/digits10_5]
+//  warning C4553: '==': result of expression not used; did you intend '='? is spurious.
+  }
+
+  {
+//[denorm_min_1
+  std::cout.precision(std::numeric_limits<double>::max_digits10);
+  if (std::numeric_limits<double>::has_denorm == std::denorm_present)
+  {
+    double d = std::numeric_limits<double>::denorm_min();
+
+      std::cout << d << std::endl; //  4.9406564584124654e-324
+
+      int exponent;
+
+      double significand = frexp(d, &exponent);
+      std::cout << "exponent = " << std::hex << exponent << std::endl; //  fffffbcf
+      std::cout << "significand = " << std::hex << significand << std::endl; // 0.50000000000000000
+  }
+  else
+  {
+    std::cout << "No denormalization. " << std::endl;
+  }
+//] [denorm_min_1]
+  }
+
+  {
+//[round_error_1
+    double round_err = std::numeric_limits<double>::epsilon() // 2.2204460492503131e-016
+                     * std::numeric_limits<double>::round_error(); // 1/2
+    std::cout << round_err << std::endl; // 1.1102230246251565e-016
+//] [/round_error_1]
+  }
+
+  {
+    typedef double T;
+//[tolerance_1
+/*`For example, if we want a tolerance that might suit about 9 arithmetical operations,
+say sqrt(9) = 3,  we could define:
+*/
+
+    T tolerance =  3 * std::numeric_limits<T>::epsilon();
+
+/*`This is very widely used in Boost.Math testing
+with Boost.Test's macro `BOOST_CHECK_CLOSE_FRACTION`
+*/
+
+    T expected = 1.0;
+    T calculated = 1.0 + std::numeric_limits<T>::epsilon();
+
+    BOOST_CHECK_CLOSE_FRACTION(expected, calculated, tolerance);
+
+//] [/tolerance_1]
+  }
+
+#if !(defined(CI_SUPPRESS_KNOWN_ISSUES) && defined(__GNUC__) && defined(_WIN32))
+  {
+//[tolerance_2
+
+  using boost::multiprecision::number;
+  using boost::multiprecision::cpp_dec_float;
+  using boost::multiprecision::et_off;
+
+  typedef number<cpp_dec_float<50>, et_off > cpp_dec_float_50; // 50 decimal digits.
+/*`[note that Boost.Test does not yet allow floating-point comparisons with expression templates on,
+so the default expression template parameter has been replaced by `et_off`.]
+*/
+
+  cpp_dec_float_50 tolerance =  3 * std::numeric_limits<cpp_dec_float_50>::epsilon();
+  cpp_dec_float_50 expected = boost::math::constants::two_pi<cpp_dec_float_50>();
+  cpp_dec_float_50 calculated = 2 * boost::math::constants::pi<cpp_dec_float_50>();
+
+  BOOST_CHECK_CLOSE_FRACTION(expected, calculated, tolerance);
+
+//] [/tolerance_2]
+  }
+
+  {
+//[tolerance_3
+
+  using boost::multiprecision::cpp_bin_float_quad;
+
+  cpp_bin_float_quad tolerance =  3 * std::numeric_limits<cpp_bin_float_quad>::epsilon();
+  cpp_bin_float_quad expected = boost::math::constants::two_pi<cpp_bin_float_quad>();
+  cpp_bin_float_quad calculated = 2 * boost::math::constants::pi<cpp_bin_float_quad>();
+
+  BOOST_CHECK_CLOSE_FRACTION(expected, calculated, tolerance);
+
+//] [/tolerance_3]
+  }
+
+  {
+//[tolerance_4
+
+  using boost::multiprecision::cpp_bin_float_oct;
+
+  cpp_bin_float_oct tolerance =  3 * std::numeric_limits<cpp_bin_float_oct>::epsilon();
+  cpp_bin_float_oct expected = boost::math::constants::two_pi<cpp_bin_float_oct>();
+  cpp_bin_float_oct calculated = 2 * boost::math::constants::pi<cpp_bin_float_oct>();
+
+  BOOST_CHECK_CLOSE_FRACTION(expected, calculated, tolerance);
+
+//] [/tolerance_4]
+  }
+
+  {
+//[nan_1]
+
+/*`NaN can be used with binary multiprecision types like `cpp_bin_float_quad`:
+*/
+  using boost::multiprecision::cpp_bin_float_quad;
+
+  if (std::numeric_limits<cpp_bin_float_quad>::has_quiet_NaN == true)
+  {
+    cpp_bin_float_quad tolerance =  3 * std::numeric_limits<cpp_bin_float_quad>::epsilon();
+
+    cpp_bin_float_quad NaN =  std::numeric_limits<cpp_bin_float_quad>::quiet_NaN();
+    std::cout << "cpp_bin_float_quad NaN is "  << NaN << std::endl; //   cpp_bin_float_quad NaN is nan
+
+    cpp_bin_float_quad expected = NaN;
+    cpp_bin_float_quad calculated = 2 * NaN;
+    // Comparisons of NaN's always fail:
+    bool b = expected == calculated;
+    std::cout << b << std::endl;
+    BOOST_CHECK_NE(expected, expected);
+    BOOST_CHECK_NE(expected, calculated);
+  }
+  else
+  {
+    std::cout << "Type " << typeid(cpp_bin_float_quad).name() << " does not have NaNs!" << std::endl;
+  }
+
+//]  [/nan_1]
+  }
+
+  {
+//[facet_1]
+
+/*`
+See [@boost:/libs/math/example/nonfinite_facet_sstream.cpp]
+and we also need
+
+  #include <boost/math/special_functions/nonfinite_num_facets.hpp>
+
+Then we can equally well use a multiprecision type cpp_bin_float_quad:
+
+*/
+  using boost::multiprecision::cpp_bin_float_quad;
+
+  typedef cpp_bin_float_quad T;
+
+  using boost::math::nonfinite_num_put;
+  using boost::math::nonfinite_num_get;
+  {
+    std::locale old_locale;
+    std::locale tmp_locale(old_locale, new nonfinite_num_put<char>);
+    std::locale new_locale(tmp_locale, new nonfinite_num_get<char>);
+    std::stringstream ss;
+    ss.imbue(new_locale);
+    T inf = std::numeric_limits<T>::infinity();
+    ss << inf; // Write out.
+   BOOST_ASSERT(ss.str() == "inf");
+    T r;
+    ss >> r; // Read back in.
+    BOOST_ASSERT(inf == r); // Confirms that the floating-point values really are identical.
+    std::cout << "infinity output was " << ss.str() << std::endl;
+    std::cout << "infinity input was " << r << std::endl;
+  }
+
+/*`
+``
+  infinity output was inf
+  infinity input was inf
+``
+Similarly we can do the same with NaN (except that we cannot use `assert` (because any comparisons with NaN always return false).
+*/
+  {
+    std::locale old_locale;
+    std::locale tmp_locale(old_locale, new nonfinite_num_put<char>);
+    std::locale new_locale(tmp_locale, new nonfinite_num_get<char>);
+    std::stringstream ss;
+    ss.imbue(new_locale);
+    T n;
+    T NaN = std::numeric_limits<T>::quiet_NaN();
+    ss << NaN; // Write out.
+    BOOST_ASSERT(ss.str() == "nan");
+    std::cout << "NaN output was " << ss.str() << std::endl;
+    ss >> n; // Read back in.
+    std::cout << "NaN input was " << n << std::endl;
+  }
+/*`
+``
+  NaN output was nan
+  NaN input was nan
+``
+*/
+//]  [/facet_1]
+  }
+
+#endif
+#endif
+} // BOOST_AUTO_TEST_CASE(test_numeric_limits_snips)
+
diff --git a/src/boost/libs/multiprecision/example/random_snips.cpp b/src/boost/libs/multiprecision/example/random_snips.cpp
new file mode 100644
index 00000000..e295564b
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/random_snips.cpp
@@ -0,0 +1,315 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2011 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/multiprecision/cpp_int.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/random.hpp>
+#include <boost/scoped_ptr.hpp>
+#include <iostream>
+#include <iomanip>
+
+void t1()
+{
+//[random_eg1
+//=#include <boost/multiprecision/cpp_int.hpp>
+//=#include <boost/random.hpp>
+//=
+//=int main()
+//={
+   using namespace boost::multiprecision;
+   using namespace boost::random;
+
+   //
+   // Declare our random number generator type, the underlying generator
+   // is the Mersenne twister mt19937 engine, and we'll generate 256 bit
+   // random values, independent_bits_engine will make multiple calls
+   // to the underlying engine until we have the requested number of bits:
+   //
+   typedef independent_bits_engine<mt19937, 256, cpp_int> generator_type;
+   generator_type gen;
+   //
+   // Generate some values:
+   //
+   std::cout << std::hex << std::showbase;
+   for(unsigned i = 0; i < 10; ++i)
+      std::cout << gen() << std::endl;
+   //
+   // Alternatively if we wish to generate random values in a fixed-precision
+   // type, then we must use an unsigned type in order to adhere to the
+   // conceptual requirements of the generator:
+   //
+   typedef independent_bits_engine<mt19937, 512, uint512_t> generator512_type;
+   generator512_type gen512;
+   //
+   // Generate some 1024-bit unsigned values:
+   //
+   std::cout << std::hex << std::showbase;
+   for(unsigned i = 0; i < 10; ++i)
+      std::cout << gen512() << std::endl;
+   //=   return 0;
+//=}
+//]
+}
+
+//
+// Output from t1() is:
+//[random_eg1_out
+/*`[pre
+0xD091BB5C22AE9EF6E7E1FAEED5C31F792082352CF807B7DFE9D300053895AFE1
+0xA1E24BBA4EE4092B18F868638C16A625474BA8C43039CD1A8C006D5FFE2D7810
+0xF51F2AE7FF1816E4F702EF59F7BADAFA285954A1B9D09511F878C4B3FB2A0137
+0xF508E4AA1C1FE6527C419418CC50AA59CCDF2E5C4C0A1F3B2452A9DC01397D8D
+0x6BF88C311CCA797AEA6DA4AEA3C78807CACE1969E0E0D4ADF5A14BAB80F00988
+0xA7DE9F4CCC450CBA0924668F5C7DC380D96089C53640AC4CEF1A2E6DAE6D9426
+0xADC1965B6613BA46C1FB41C2BD9B0ECDBE3DEDFC7989C8EE6468FD6E6C0DF032
+0xA7CD66342C826D8B2BD2E4124D4A2DBEB4BF6FA7CC1A89590826328251097330
+0x46E46CB0DF577EC20BD1E364262C556418DDA0C9FE7B45D9D2CE21C9D268409A
+0xB1E049E1200BFA47512D6E73C3851EEEF341C0817D973E4808D17554A9E20D28
+0xD091BB5C22AE9EF6E7E1FAEED5C31F792082352CF807B7DFE9D300053895AFE1A1E24BBA4EE4092B18F868638C16A625474BA8C43039CD1A8C006D5FFE2D7810
+0xF51F2AE7FF1816E4F702EF59F7BADAFA285954A1B9D09511F878C4B3FB2A0137F508E4AA1C1FE6527C419418CC50AA59CCDF2E5C4C0A1F3B2452A9DC01397D8D
+0x6BF88C311CCA797AEA6DA4AEA3C78807CACE1969E0E0D4ADF5A14BAB80F00988A7DE9F4CCC450CBA0924668F5C7DC380D96089C53640AC4CEF1A2E6DAE6D9426
+0xADC1965B6613BA46C1FB41C2BD9B0ECDBE3DEDFC7989C8EE6468FD6E6C0DF032A7CD66342C826D8B2BD2E4124D4A2DBEB4BF6FA7CC1A89590826328251097330
+0x46E46CB0DF577EC20BD1E364262C556418DDA0C9FE7B45D9D2CE21C9D268409AB1E049E1200BFA47512D6E73C3851EEEF341C0817D973E4808D17554A9E20D28
+0x70518CE6203AC30361ADD0AB35D0430CC3F8E8920D1C8509CB92388E095436BF2FD6E20868A29AF97D61330B753EC6FC7211EFEA7CD15133A574C4FFCB41F198
+0xB598EEF6EBBE7347C1332568CEBA5A7046A99459B4AD9F11AE00FEAA00B8B573A7B480B6B5F0B06C29A0EC27A4DAA0101E76A1C574BE91337F94C950C61F6ED6
+0xF5B1C7A192E195F8572384D4E0732C8895D41B68CEE496C3394BBD52048CD47CC05309BED23D2D63414DE9C5D2229F23818666A3F0A8B109B2F6B12769A48341
+0xE4123C566C548C8FF5941F6194B993AA8C1651342876763C237CE42EC300D11B263821CA3AEB820241EC0F84CF4AC36DD7393EE6FD0FC06A4118A30A551B54A4
+0xD074F86F4CC1C54A3E57A70303774CDAEDE43895379CE62759988939E8490DDC325410E1D9352F6A4047080AF47C081D9DB51A85C765D71F79297527FCCA2773
+]
+*/
+//]
+
+
+void t2()
+{
+   std::cout << std::dec;
+//[random_eg2
+//=#include <boost/multiprecision/cpp_int.hpp>
+//=#include <boost/random.hpp>
+//=
+//=int main()
+//={
+   using namespace boost::multiprecision;
+   using namespace boost::random;
+
+   //
+   // Generate integers in a given range using uniform_int,
+   // the underlying generator is invoked multiple times
+   // to generate enough bits:
+   //
+   mt19937 mt;
+   uniform_int_distribution<cpp_int> ui(-(cpp_int(1) << 256), cpp_int(1) << 256);
+   //
+   // Generate the numbers:
+   //
+   for(unsigned i = 0; i < 10; ++i)
+      std::cout << ui(mt) << std::endl;
+
+//=   return 0;
+//=}
+//]
+}
+//[random_eg2_out
+/*`
+Program output is
+
+[pre
+25593993629538149833210527544371584707508847463356155903670894544241785158492
+12721121657520147247744796431842326146296294180809160027132416389225539366745
+106034929479008809862776424170460808190085984129117168803272987114325199071833
+86048861429530654936263414134573980939351899046345384016090167510299251354700
+-23473382144925885755951447143660880642389842563343761080591177733698450031250
+76840269649240973945508128641415259490679375154523618053296924666747244530145
+21638369166612496703991271955994563624044383325105383029306009417224944272131
+18829152205014764576551421737727569993966577957447887116062495161081023584880
+101521572847669971701030312596819435590097618913255156117898217707115132658117
+-97490271301923067621481012355971422109456300816856752380346627103308328292057
+]
+*/
+//]
+
+void t3()
+{
+//[random_eg3
+//=#include <boost/multiprecision/cpp_bin_float.hpp>
+//=#include <boost/random.hpp>
+//=
+//=int main()
+//={
+   using namespace boost::multiprecision;
+   using namespace boost::random;
+
+   mt19937 gen;
+   //
+   // Generate the values:
+   //
+   std::cout << std::setprecision(50);
+   for(unsigned i = 0; i < 20; ++i)
+      std::cout << generate_canonical<cpp_bin_float_50, std::numeric_limits<cpp_bin_float_50>::digits>(gen) << std::endl;
+//=   return 0;
+//=}
+//]
+}
+
+//[random_eg3_out
+/*`
+Which produces the following output:
+
+[pre
+0.96886777112423135248554451482797431507115448261086
+0.54722059636785192454525760726084778627750790023546
+0.99646132554800874317788284808573062871409279729804
+0.98110969177693891782396443737643892769773768718591
+0.29702944955795083040856753579705872634075574515969
+0.63976335709815275010379796044374742646738557798647
+0.79792861516022605265555700991255998690336456180995
+0.68135953856026596523755400091345037778580909233387
+0.47475868061723477935404326837783394169122045199915
+0.30191312687731969398296589840622989141067852863748
+0.87242882006730022427155209451091472382531795659709
+0.82190326480741096300318873712966555706035846579562
+0.49058903962146072778707295967429263659897501512813
+0.2102090745190061764133345429475530760261103345204
+0.4087311609617603484960794513055502599728804206333
+0.79397497154919267900450180642484943996546102712187
+0.70577425166871982574205252142383800792823003687121
+0.64396095652194035523385641523010248768636064728226
+0.5737546665965914620678634509134819579811035412969
+0.017773895576552474810236796736785695789752666554273
+]
+*/
+//]
+
+void t4()
+{
+   std::cout << std::endl;
+//[random_eg4
+//=#include <boost/multiprecision/cpp_bin_float.hpp>
+//=#include <boost/multiprecision/cpp_int.hpp>
+//=#include <boost/random.hpp>
+//=
+//=int main()
+//={
+   using namespace boost::multiprecision;
+   using namespace boost::random;
+   //
+   // Generate some distruted values:
+   //
+   uniform_real_distribution<cpp_bin_float_50> ur(-20, 20);
+   gamma_distribution<cpp_bin_float_50> gd(20);
+   independent_bits_engine<mt19937, std::numeric_limits<cpp_bin_float_50>::digits, cpp_int> gen;
+   //
+   // Generate some values:
+   //
+   std::cout << std::setprecision(50);
+   for(unsigned i = 0; i < 20; ++i)
+      std::cout << ur(gen) << std::endl;
+   for(unsigned i = 0; i < 20; ++i)
+      std::cout << gd(gen) << std::endl;
+//=   return 0;
+//=}
+//]
+}
+
+//[random_eg4_out
+/*`
+Which produces the following output:
+
+[pre
+-18.576837157065858312137736538355805944098004018928
+4.5605477000094480453928920098152026546185388161216
+-1.7611402252150150370944527411235180945558276280598
+-2.471338289511354190492328039842914272146783953149
+-7.4131520453411321647183692139916357315276121488316
+-9.192739117661751364518299455475684051782402347659
+7.0126880787149555595443325648941661436898526919013
+2.8554749162054097111723076181877881960039268668423
+14.390501287552165467965587841551705310012046701036
+-8.9747073123748752412086051960748002945548570524149
+-8.1305063133718605220959174700954037986278348616362
+9.5496899464463627949564295930962040525540578754312
+-15.309681742947663333436391348699943078942921692008
+2.0454914298189175280771944784358385982869708951824
+-10.069253024538932382193363493367304983742246396276
+13.449212808583153116670057807764145176004060370818
+-6.0065092542772507561228141992257782449634820245355
+15.00971466974838379824678369267201922989930663822
+16.158514812070905438581736305533045434508525979205
+-2.1531361299576399413547008719541457739794964378093
+19.398278792113040046930806838893737245011219380822
+12.965216582396067073600685365545292876001524716225
+19.561779374349650983983836397553672788578622096947
+15.982213641588944604037715576313848977716540941271
+23.96044616946856385664151481695038833903083043492
+21.054716943622792848187523422423642819628010070375
+18.596078774135209530930707331338838805575875990091
+19.539530839287848627426769425090194390388333335812
+17.176133236359396942946640290935498641489373354297
+16.228802394876800099035133760539461530246286999827
+23.63807160907473465631049083277558060813997674519
+12.838499607321990428122225501321564153572478845401
+16.878362445712403300584931374939967549572637230102
+20.646246409377134464856282996941395597420615529803
+16.602429236226052406561338766554127142762673418695
+21.680007865714197450495711030406314524681744024329
+21.038948660115771777833205901845639760348321521616
+30.494499676527802078320016654058105593076348727966
+18.704734464995637480940828829962787676146589788572
+22.502216997171061548799304902323434654678156658236
+]
+*/
+//]
+
+void t5()
+{
+//[random_eg5
+//=#include <boost/multiprecision/cpp_bin_float.hpp>
+//=#include <boost/random.hpp>
+//=#include <boost/scoped_ptr.hpp>
+//=
+//=int main()
+//={
+   using namespace boost::multiprecision;
+   using namespace boost::random;
+   //
+   // Generate some multiprecision values, note that the generator is so large
+   // that we have to allocate it on the heap, otherwise we may run out of
+   // stack space!  We could avoid this by using a floating point type which
+   // allocates it's internal storage on the heap - cpp_bin_float will do
+   // this with the correct template parameters, as will the GMP or MPFR
+   // based reals.
+   //
+   typedef lagged_fibonacci_01_engine<cpp_bin_float_50, 48, 44497, 21034 > big_fib_gen;
+   boost::scoped_ptr<big_fib_gen> pgen(new big_fib_gen);
+   //
+   // Generate some values:
+   //
+   std::cout << std::setprecision(50);
+   for(unsigned i = 0; i < 20; ++i)
+      std::cout << (*pgen)() << std::endl;
+   //
+   // try again with a ranlux generator, this is not quite so large
+   // so we can use the heap this time:
+   //
+   typedef subtract_with_carry_01_engine<cpp_bin_float_50, std::numeric_limits<cpp_bin_float_50>::digits - 5, 10, 24 > ranlux_big_base_01;
+   typedef discard_block_engine< ranlux_big_base_01, 389, 24 > big_ranlux;
+   big_ranlux rg;
+   for(unsigned i = 0; i < 20; ++i)
+      std::cout << rg() << std::endl;
+//=   return 0;
+//=}
+//]
+}
+
+int main()
+{
+   t1();
+   t2();
+   t3();
+   t4();
+   t5();
+   return 0;
+}
+
diff --git a/src/boost/libs/multiprecision/example/safe_prime.cpp b/src/boost/libs/multiprecision/example/safe_prime.cpp
new file mode 100644
index 00000000..dbe40ef2
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/safe_prime.cpp
@@ -0,0 +1,45 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2012 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+//[safe_prime
+
+#include <boost/multiprecision/cpp_int.hpp>
+#include <boost/multiprecision/miller_rabin.hpp>
+#include <iostream>
+#include <iomanip>
+
+int main()
+{
+   using namespace boost::random;
+   using namespace boost::multiprecision;
+
+   typedef cpp_int int_type;
+   mt11213b base_gen(clock());
+   independent_bits_engine<mt11213b, 256, int_type> gen(base_gen);
+   //
+   // We must use a different generator for the tests and number generation, otherwise
+   // we get false positives.
+   //
+   mt19937 gen2(clock());
+
+   for(unsigned i = 0; i < 100000; ++i)
+   {
+      int_type n = gen();
+      if(miller_rabin_test(n, 25, gen2))
+      {
+         // Value n is probably prime, see if (n-1)/2 is also prime:
+         std::cout << "We have a probable prime with value: " << std::hex << std::showbase << n << std::endl;
+         if(miller_rabin_test((n-1)/2, 25, gen2))
+         {
+            std::cout << "We have a safe prime with value: " << std::hex << std::showbase << n << std::endl;
+            return 0;
+         }
+      }
+   }
+   std::cout << "Ooops, no safe primes were found - probably a bad choice of seed values!" << std::endl;
+   return 0;
+}
+
+//]
diff --git a/src/boost/libs/multiprecision/example/tommath_snips.cpp b/src/boost/libs/multiprecision/example/tommath_snips.cpp
new file mode 100644
index 00000000..f3b817bd
--- /dev/null
+++ b/src/boost/libs/multiprecision/example/tommath_snips.cpp
@@ -0,0 +1,84 @@
+///////////////////////////////////////////////////////////////
+//  Copyright 2011 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/multiprecision/tommath.hpp>
+#include <iostream>
+
+void t1()
+{
+//[tommath_eg
+//=#include <boost/multiprecision/tommath.hpp>
+//=#include <iostream>
+//=
+//=int main()
+//={
+   boost::multiprecision::tom_int v = 1;
+
+   // Do some arithmetic:
+   for(unsigned i = 1; i <= 1000; ++i)
+      v *= i;
+
+   std::cout << v << std::endl; // prints 1000!
+   std::cout << std::hex << v << std::endl; // prints 1000! in hex format
+
+   try{
+      std::cout << std::hex << -v << std::endl; // Ooops! can't print a negative value in hex format!
+   }
+   catch(const std::runtime_error& e)
+   {
+      std::cout << e.what() << std::endl;
+   }
+
+   try{
+      // v is not a 2's complement type, bitwise operations are only supported
+      // on positive values:
+      v = -v & 2;
+   }
+   catch(const std::runtime_error& e)
+   {
+      std::cout << e.what() << std::endl;
+   }
+
+//=   return 0;
+//=}
+//]
+}
+
+void t3()
+{
+//[mp_rat_eg
+//=#include <boost/multiprecision/tommath.hpp>
+//=#include <iostream>
+//=
+//=int main()
+//={
+   using namespace boost::multiprecision;
+
+   tom_rational v = 1;
+
+   // Do some arithmetic:
+   for(unsigned i = 1; i <= 1000; ++i)
+      v *= i;
+   v /= 10;
+
+   std::cout << v << std::endl; // prints 1000! / 10
+   std::cout << numerator(v) << std::endl;
+   std::cout << denominator(v) << std::endl;
+
+   tom_rational w(2, 3); // Component wise constructor
+   std::cout << w << std::endl; // prints 2/3
+
+//=   return 0;
+//=}
+//]
+}
+
+int main()
+{
+   t1();
+   t3();
+   return 0;
+}
+