7 files changed, 1075 insertions, 3 deletions

diff --git a/libs/math/tools/carlson_ellint_data.cpp b/libs/math/tools/carlson_ellint_data.cpp
index c234b4b17..5bd980733 100644
--- a/libs/math/tools/carlson_ellint_data.cpp
+++ b/libs/math/tools/carlson_ellint_data.cpp
@@ -6,7 +6,8 @@
 
 #include <boost/math/tools/test_data.hpp>
 #include <boost/test/included/prg_exec_monitor.hpp>
-#include <boost/math/special_functions/ellint_3.hpp>
+#include <boost/math/special_functions/ellint_rj.hpp>
+#include <boost/math/special_functions/ellint_rd.hpp>
 #include <fstream>
 #include <boost/math/tools/test_data.hpp>
 #include <boost/random.hpp>
@@ -20,6 +21,377 @@ float truncate_to_float(float const * pf)
    return *pf;
 }
 
+//
+// Archived here is the original implementation of this
+// function by Xiaogang Zhang, we can use this to
+// generate special test cases for the new version:
+//
+template <typename T, typename Policy>
+T ellint_rj_old(T x, T y, T z, T p, const Policy& pol)
+{
+   T value, u, lambda, alpha, beta, sigma, factor, tolerance;
+   T X, Y, Z, P, EA, EB, EC, E2, E3, S1, S2, S3;
+   unsigned long k;
+
+   BOOST_MATH_STD_USING
+      using namespace boost::math;
+
+   static const char* function = "boost::math::ellint_rj<%1%>(%1%,%1%,%1%)";
+
+   if(x < 0)
+   {
+      return policies::raise_domain_error<T>(function,
+         "Argument x must be non-negative, but got x = %1%", x, pol);
+   }
+   if(y < 0)
+   {
+      return policies::raise_domain_error<T>(function,
+         "Argument y must be non-negative, but got y = %1%", y, pol);
+   }
+   if(z < 0)
+   {
+      return policies::raise_domain_error<T>(function,
+         "Argument z must be non-negative, but got z = %1%", z, pol);
+   }
+   if(p == 0)
+   {
+      return policies::raise_domain_error<T>(function,
+         "Argument p must not be zero, but got p = %1%", p, pol);
+   }
+   if(x + y == 0 || y + z == 0 || z + x == 0)
+   {
+      return policies::raise_domain_error<T>(function,
+         "At most one argument can be zero, "
+         "only possible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol);
+   }
+
+   // error scales as the 6th power of tolerance
+   tolerance = pow(T(1) * tools::epsilon<T>() / 3, T(1) / 6);
+
+   // for p < 0, the integral is singular, return Cauchy principal value
+   if(p < 0)
+   {
+      //
+      // We must ensure that (z - y) * (y - x) is positive.
+      // Since the integral is symmetrical in x, y and z
+      // we can just permute the values:
+      //
+      if(x > y)
+         std::swap(x, y);
+      if(y > z)
+         std::swap(y, z);
+      if(x > y)
+         std::swap(x, y);
+
+      T q = -p;
+      T pmy = (z - y) * (y - x) / (y + q);  // p - y
+
+      BOOST_ASSERT(pmy >= 0);
+
+      p = pmy + y;
+      value = ellint_rj_old(x, y, z, p, pol);
+      value *= pmy;
+      value -= 3 * boost::math::ellint_rf(x, y, z, pol);
+      value += 3 * sqrt((x * y * z) / (x * z + p * q)) * boost::math::ellint_rc(x * z + p * q, p * q, pol);
+      value /= (y + q);
+      return value;
+   }
+
+   // duplication
+   sigma = 0;
+   factor = 1;
+   k = 1;
+   do
+   {
+      u = (x + y + z + p + p) / 5;
+      X = (u - x) / u;
+      Y = (u - y) / u;
+      Z = (u - z) / u;
+      P = (u - p) / u;
+
+      if((tools::max)(abs(X), abs(Y), abs(Z), abs(P)) < tolerance)
+         break;
+
+      T sx = sqrt(x);
+      T sy = sqrt(y);
+      T sz = sqrt(z);
+
+      lambda = sy * (sx + sz) + sz * sx;
+      alpha = p * (sx + sy + sz) + sx * sy * sz;
+      alpha *= alpha;
+      beta = p * (p + lambda) * (p + lambda);
+      sigma += factor * boost::math::ellint_rc(alpha, beta, pol);
+      factor /= 4;
+      x = (x + lambda) / 4;
+      y = (y + lambda) / 4;
+      z = (z + lambda) / 4;
+      p = (p + lambda) / 4;
+      ++k;
+   } while(k < policies::get_max_series_iterations<Policy>());
+
+   // Check to see if we gave up too soon:
+   policies::check_series_iterations<T>(function, k, pol);
+
+   // Taylor series expansion to the 5th order
+   EA = X * Y + Y * Z + Z * X;
+   EB = X * Y * Z;
+   EC = P * P;
+   E2 = EA - 3 * EC;
+   E3 = EB + 2 * P * (EA - EC);
+   S1 = 1 + E2 * (E2 * T(9) / 88 - E3 * T(9) / 52 - T(3) / 14);
+   S2 = EB * (T(1) / 6 + P * (T(-6) / 22 + P * T(3) / 26));
+   S3 = P * ((EA - EC) / 3 - P * EA * T(3) / 22);
+   value = 3 * sigma + factor * (S1 + S2 + S3) / (u * sqrt(u));
+
+   return value;
+}
+
+template <typename T, typename Policy>
+T ellint_rd_imp_old(T x, T y, T z, const Policy& pol)
+{
+   T value, u, lambda, sigma, factor, tolerance;
+   T X, Y, Z, EA, EB, EC, ED, EE, S1, S2;
+   unsigned long k;
+
+   BOOST_MATH_STD_USING
+   using namespace boost::math;
+
+   static const char* function = "boost::math::ellint_rd<%1%>(%1%,%1%,%1%)";
+
+   if(x < 0)
+   {
+      return policies::raise_domain_error<T>(function,
+         "Argument x must be >= 0, but got %1%", x, pol);
+   }
+   if(y < 0)
+   {
+      return policies::raise_domain_error<T>(function,
+         "Argument y must be >= 0, but got %1%", y, pol);
+   }
+   if(z <= 0)
+   {
+      return policies::raise_domain_error<T>(function,
+         "Argument z must be > 0, but got %1%", z, pol);
+   }
+   if(x + y == 0)
+   {
+      return policies::raise_domain_error<T>(function,
+         "At most one argument can be zero, but got, x + y = %1%", x + y, pol);
+   }
+
+   // error scales as the 6th power of tolerance
+   tolerance = pow(tools::epsilon<T>() / 3, T(1) / 6);
+
+   // duplication
+   sigma = 0;
+   factor = 1;
+   k = 1;
+   do
+   {
+      u = (x + y + z + z + z) / 5;
+      X = (u - x) / u;
+      Y = (u - y) / u;
+      Z = (u - z) / u;
+      if((tools::max)(abs(X), abs(Y), abs(Z)) < tolerance)
+         break;
+      T sx = sqrt(x);
+      T sy = sqrt(y);
+      T sz = sqrt(z);
+      lambda = sy * (sx + sz) + sz * sx; //sqrt(x * y) + sqrt(y * z) + sqrt(z * x);
+      sigma += factor / (sz * (z + lambda));
+      factor /= 4;
+      x = (x + lambda) / 4;
+      y = (y + lambda) / 4;
+      z = (z + lambda) / 4;
+      ++k;
+   } while(k < policies::get_max_series_iterations<Policy>());
+
+   // Check to see if we gave up too soon:
+   policies::check_series_iterations<T>(function, k, pol);
+
+   // Taylor series expansion to the 5th order
+   EA = X * Y;
+   EB = Z * Z;
+   EC = EA - EB;
+   ED = EA - 6 * EB;
+   EE = ED + EC + EC;
+   S1 = ED * (ED * T(9) / 88 - Z * EE * T(9) / 52 - T(3) / 14);
+   S2 = Z * (EE / 6 + Z * (-EC * T(9) / 22 + Z * EA * T(3) / 26));
+   value = 3 * sigma + factor * (1 + S1 + S2) / (u * sqrt(u));
+
+   return value;
+}
+
+template <typename T, typename Policy>
+T ellint_rf_imp_old(T x, T y, T z, const Policy& pol)
+{
+   T value, X, Y, Z, E2, E3, u, lambda, tolerance;
+   unsigned long k;
+   BOOST_MATH_STD_USING
+   using namespace boost::math;
+   static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)";
+   if(x < 0 || y < 0 || z < 0)
+   {
+      return policies::raise_domain_error<T>(function,
+         "domain error, all arguments must be non-negative, "
+         "only sensible result is %1%.",
+         std::numeric_limits<T>::quiet_NaN(), pol);
+   }
+   if(x + y == 0 || y + z == 0 || z + x == 0)
+   {
+      return policies::raise_domain_error<T>(function,
+         "domain error, at most one argument can be zero, "
+         "only sensible result is %1%.",
+         std::numeric_limits<T>::quiet_NaN(), pol);
+   }
+   // Carlson scales error as the 6th power of tolerance,
+   // but this seems not to work for types larger than
+   // 80-bit reals, this heuristic seems to work OK:
+   if(policies::digits<T, Policy>() > 64)
+   {
+      tolerance = pow(tools::epsilon<T>(), T(1) / 4.25f);
+      BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);
+   }
+   else
+   {
+      tolerance = pow(4 * tools::epsilon<T>(), T(1) / 6);
+      BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);
+   }
+   // duplication
+   k = 1;
+   do
+   {
+      u = (x + y + z) / 3;
+      X = (u - x) / u;
+      Y = (u - y) / u;
+      Z = (u - z) / u;
+      // Termination condition:
+      if((tools::max)(abs(X), abs(Y), abs(Z)) < tolerance)
+         break;
+      T sx = sqrt(x);
+      T sy = sqrt(y);
+      T sz = sqrt(z);
+      lambda = sy * (sx + sz) + sz * sx;
+      x = (x + lambda) / 4;
+      y = (y + lambda) / 4;
+      z = (z + lambda) / 4;
+      ++k;
+   } while(k < policies::get_max_series_iterations<Policy>());
+   // Check to see if we gave up too soon:
+   policies::check_series_iterations<T>(function, k, pol);
+   BOOST_MATH_INSTRUMENT_VARIABLE(k);
+   // Taylor series expansion to the 5th order
+   E2 = X * Y - Z * Z;
+   E3 = X * Y * Z;
+   value = (1 + E2*(E2 / 24 - E3*T(3) / 44 - T(0.1)) + E3 / 14) / sqrt(u);
+   BOOST_MATH_INSTRUMENT_VARIABLE(value);
+   return value;
+}
+
+
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rj_data_4e(mp_t n)
+{
+   mp_t result = ellint_rj_old(n, n, n, n, boost::math::policies::policy<>());
+   return boost::math::make_tuple(n, n, n, result);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> generate_rj_data_3e(mp_t x, mp_t p)
+{
+   mp_t r = ellint_rj_old(x, x, x, p, boost::math::policies::policy<>());
+   return boost::math::make_tuple(x, x, x, p, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> generate_rj_data_2e_1(mp_t x, mp_t y, mp_t p)
+{
+   mp_t r = ellint_rj_old(x, x, y, p, boost::math::policies::policy<>());
+   return boost::math::make_tuple(x, x, y, p, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> generate_rj_data_2e_2(mp_t x, mp_t y, mp_t p)
+{
+   mp_t r = ellint_rj_old(x, y, x, p, boost::math::policies::policy<>());
+   return boost::math::make_tuple(x, y, x, p, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> generate_rj_data_2e_3(mp_t x, mp_t y, mp_t p)
+{
+   mp_t r = ellint_rj_old(y, x, x, p, boost::math::policies::policy<>());
+   return boost::math::make_tuple(y, x, x, p, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> generate_rj_data_2e_4(mp_t x, mp_t y, mp_t p)
+{
+   mp_t r = ellint_rj_old(x, y, p, p, boost::math::policies::policy<>());
+   return boost::math::make_tuple(x, y, p, p, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data_2e_1(mp_t x, mp_t y)
+{
+   mp_t r = ellint_rd_imp_old(x, y, y, boost::math::policies::policy<>());
+   return boost::math::make_tuple(x, y, y, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data_2e_2(mp_t x, mp_t y)
+{
+   mp_t r = ellint_rd_imp_old(x, x, y, boost::math::policies::policy<>());
+   return boost::math::make_tuple(x, x, y, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data_2e_3(mp_t x)
+{
+   mp_t r = ellint_rd_imp_old(mp_t(0), x, x, boost::math::policies::policy<>());
+   return boost::math::make_tuple(0, x, x, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data_3e(mp_t x)
+{
+   mp_t r = ellint_rd_imp_old(x, x, x, boost::math::policies::policy<>());
+   return boost::math::make_tuple(x, x, x, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data_0xy(mp_t x, mp_t y)
+{
+   mp_t r = ellint_rd_imp_old(mp_t(0), x, y, boost::math::policies::policy<>());
+   return boost::math::make_tuple(mp_t(0), x, y, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data_xxx(mp_t x)
+{
+   mp_t r = ellint_rf_imp_old(x, x, x, boost::math::policies::policy<>());
+   return boost::math::make_tuple(x, x, x, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data_xyy(mp_t x, mp_t y)
+{
+   mp_t r = ellint_rf_imp_old(x, y, y, boost::math::policies::policy<>());
+   return boost::math::make_tuple(x, y, y, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data_xxy(mp_t x, mp_t y)
+{
+   mp_t r = ellint_rf_imp_old(x, x, y, boost::math::policies::policy<>());
+   return boost::math::make_tuple(x, x, y, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data_xyx(mp_t x, mp_t y)
+{
+   mp_t r = ellint_rf_imp_old(x, y, x, boost::math::policies::policy<>());
+   return boost::math::make_tuple(x, y, x, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data_0yy(mp_t y)
+{
+   mp_t r = ellint_rf_imp_old(mp_t(0), y, y, boost::math::policies::policy<>());
+   return boost::math::make_tuple(mp_t(0), y, y, r);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data_xy0(mp_t x, mp_t y)
+{
+   mp_t r = ellint_rf_imp_old(x, y, mp_t(0), boost::math::policies::policy<>());
+   return boost::math::make_tuple(x, y, mp_t(0), r);
+}
+
 boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data(mp_t n)
 {
    static boost::mt19937 r;
@@ -99,11 +471,107 @@ boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data(mp_t n)
    return boost::math::make_tuple(xr, yr, zr, result);
 }
 
-int cpp_main(int argc, char*argv [])
+mp_t rg_imp(mp_t x, mp_t y, mp_t z)
+{
+   using std::swap;
+   // If z is zero permute so the call to RD is valid:
+   if(z == 0)
+      swap(x, z);
+   return (z * ellint_rf_imp_old(x, y, z, boost::math::policies::policy<>())
+      - (x - z) * (y - z) * ellint_rd_imp_old(x, y, z, boost::math::policies::policy<>()) / 3
+      + sqrt(x * y / z)) / 2;
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_data(mp_t n)
+{
+   static boost::mt19937 r;
+   boost::uniform_real<float> ur(0, 1);
+   boost::uniform_int<int> ui(-100, 100);
+   float x = ur(r);
+   x = ldexp(x, ui(r));
+   mp_t xr(truncate_to_float(&x));
+   float y = ur(r);
+   y = ldexp(y, ui(r));
+   mp_t yr(truncate_to_float(&y));
+   float z = ur(r);
+   z = ldexp(z, ui(r));
+   mp_t zr(truncate_to_float(&z));
+
+   mp_t result = rg_imp(xr, yr, zr);
+   return boost::math::make_tuple(xr, yr, zr, result);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_xxx(mp_t x)
+{
+   mp_t result = rg_imp(x, x, x);
+   return boost::math::make_tuple(x, x, x, result);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_xyy(mp_t x, mp_t y)
+{
+   mp_t result = rg_imp(x, y, y);
+   return boost::math::make_tuple(x, y, y, result);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_xxy(mp_t x, mp_t y)
+{
+   mp_t result = rg_imp(x, x, y);
+   return boost::math::make_tuple(x, x, y, result);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_xyx(mp_t x, mp_t y)
+{
+   mp_t result = rg_imp(x, y, x);
+   return boost::math::make_tuple(x, y, x, result);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_0xx(mp_t x)
+{
+   mp_t result = rg_imp(mp_t(0), x, x);
+   return boost::math::make_tuple(mp_t(0), x, x, result);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_x0x(mp_t x)
+{
+   mp_t result = rg_imp(x, mp_t(0), x);
+   return boost::math::make_tuple(x, mp_t(0), x, result);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_xx0(mp_t x)
+{
+   mp_t result = rg_imp(x, x, mp_t(0));
+   return boost::math::make_tuple(x, x, mp_t(0), result);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_00x(mp_t x)
+{
+   mp_t result = sqrt(x) / 2;
+   return boost::math::make_tuple(mp_t(0), mp_t(0), x, result);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_0x0(mp_t x)
+{
+   mp_t result = sqrt(x) / 2;
+   return boost::math::make_tuple(mp_t(0), x, mp_t(0), result);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_x00(mp_t x)
+{
+   mp_t result = sqrt(x) / 2;
+   return boost::math::make_tuple(x, mp_t(0), mp_t(0), result);
+}
+
+boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_xy0(mp_t x, mp_t y)
+{
+   mp_t result = rg_imp(x, y, mp_t(0));
+   return boost::math::make_tuple(x, y, mp_t(0), result);
+}
+
+int cpp_main(int argc, char*argv[])
 {
    using namespace boost::math::tools;
 
-   parameter_info<mp_t> arg1, arg2;
+   parameter_info<mp_t> arg1, arg2, arg3;
    test_data<mp_t> data;
 
    bool cont;
@@ -113,6 +581,7 @@ int cpp_main(int argc, char*argv [])
       return 1;
 
    do{
+#if 0
       int count;
       std::cout << "Number of points: ";
       std::cin >> count;
@@ -129,6 +598,20 @@ int cpp_main(int argc, char*argv [])
       std::getline(std::cin, line);
       boost::algorithm::trim(line);
       cont = (line == "y");
+#else
+      get_user_parameter_info(arg1, "x");
+      get_user_parameter_info(arg2, "y");
+      //get_user_parameter_info(arg3, "p");
+      arg1.type |= dummy_param;
+      arg2.type |= dummy_param;
+      //arg3.type |= dummy_param;
+      data.insert(generate_rd_data_0xy, arg1, arg2);
+
+      std::cout << "Any more data [y/n]?";
+      std::getline(std::cin, line);
+      boost::algorithm::trim(line);
+      cont = (line == "y");
+#endif
    }while(cont);
 
    std::cout << "Enter name of test data file [default=ellint_rf_data.ipp]";
diff --git a/libs/math/tools/doc/cstdfloat.qbk b/libs/math/tools/doc/cstdfloat.qbk
new file mode 100644
index 000000000..5d9b2fd78
--- /dev/null
+++ b/libs/math/tools/doc/cstdfloat.qbk
@@ -0,0 +1,249 @@
+[book Standardized Floating-Point typedefs for C and C++
+
+    [quickbook 1.7]
+    [copyright 2014  Christopher Kormanyos, John Maddock, Paul A. Bristow]
+    [license
+        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])
+    ]
+    [authors [Kormanyos, Christopher],  [Maddock, John], [Bristow, Paul A.] ]
+    [last-revision $Date$]
+    [/version 1.8.3]
+]
+
+[template tr1[] [@http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2005/n1836.pdf Technical Report on C++ Library Extensions]]
+[template C99[] [@http://www.open-std.org/JTC1/SC22/WG14/www/docs/n1256.pdf C99 Standard ISO/IEC 9899:1999]]
+
+[def __gsl [@http://www.gnu.org/software/gsl/ GSL-1.9]]
+[def __glibc [@http://www.gnu.org/software/libc/ GNU C Lib]]
+[def __hpc [@http://docs.hp.com/en/B9106-90010/index.html HP-UX C Library]]
+[def __cephes [@http://www.netlib.org/cephes/ Cephes]]
+[def __NTL [@http://www.shoup.net/ntl/ NTL A Library for doing Number Theory]]
+[def __NTL_RR [@http://shoup.net/ntl/doc/RR.txt NTL::RR]]
+[def __NTL_quad_float [@http://shoup.net/ntl/doc/quad_float.txt NTL::quad_float]]
+[def __MPFR [@http://www.mpfr.org/ GNU MPFR library]]
+[def __GMP [@http://gmplib.org/ GNU Multiple Precision Arithmetic Library]]
+[def __multiprecision [@http://www.boost.org/doc/libs/1_53_0_beta1/libs/multiprecision/doc/html/index.html Boost.Multiprecision]]
+[def __cpp_dec_float [@http://www.boost.org/doc/libs/1_53_0_beta1/libs/multiprecision/doc/html/boost_multiprecision/tut/floats/cpp_dec_float.html cpp_dec_float]]
+[def __R [@http://www.r-project.org/ The R Project for Statistical Computing]]
+[def __godfrey [link godfrey Godfrey]]
+[def __pugh [link pugh Pugh]]
+[def __NaN [@http://en.wikipedia.org/wiki/NaN NaN]]
+[def __errno [@http://en.wikipedia.org/wiki/Errno `::errno`]]
+[def __Mathworld [@http://mathworld.wolfram.com Wolfram MathWorld]]
+[def __Mathematica [@http://www.wolfram.com/products/mathematica/index.html Wolfram Mathematica]]
+[def __WolframAlpha [@http://www.wolframalpha.com/ Wolfram Alpha]]
+[def __TOMS748 [@http://portal.acm.org/citation.cfm?id=210111 TOMS Algorithm 748: enclosing zeros of continuous functions]]
+[def __TOMS910 [@http://portal.acm.org/citation.cfm?id=1916469 TOMS Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations]]
+[def __why_complements [link why_complements why complements?]]
+[def __complements [link math_toolkit.stat_tut.overview.complements complements]]
+[def __performance [link perf performance]]
+[def __building [link math_toolkit.building building libraries]]
+[def __e_float [@http://calgo.acm.org/910.zip e_float (TOMS Algorithm 910)]]
+[def __Abramowitz_Stegun M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, NBS  (1964)]
+[def __DMLF [@http://dlmf.nist.gov/ NIST Digital Library of Mathematical Functions]]
+[def __IEEE754  [@http://en.wikipedia.org/wiki/IEEE_floating_point IEEE_floating_point]]
+[def __N3626  [@http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3626.pdf N3626]]
+[def __N1703 [@http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1703.pdf N1703]]
+
+[/ Some composite templates]
+[template super[x]'''<superscript>'''[x]'''</superscript>''']
+[template sub[x]'''<subscript>'''[x]'''</subscript>''']
+[template floor[x]'''&#x230A;'''[x]'''&#x230B;''']
+[template floorlr[x][lfloor][x][rfloor]]
+[template ceil[x] '''&#x2308;'''[x]'''&#x2309;''']
+
+[/template header_file[file] [@../../../../[file] [file]]]
+
+[note A printer-friendly PDF version of this manual is also available.]
+
+[section:overview Overview]
+
+The header `<boost/cstdfloat.hpp>` provides optional standardized
+floating-point `typedef`s having specified widths.
+These are useful for writing portable code because they
+should behave identically on all platforms.
+All `typedef`s are in `namespace boost`.
+
+The `typedef`s include `float16_t, float32_t, float64_t, float128_t`,
+their corresponding least and fast types,
+and the corresponding maximum-width type.
+The `typedef`s are based on underlying built-in types
+such as `float`, `double`, or `long double`, or based on other compiler-specific
+non-standardized types such as `__float128`.
+The underlying types of these typedef's must conform with
+the corresponding specifications of binary16, binary32, binary64,
+and binary128 in __IEEE754 floating-point format
+[@http://en.wikipedia.org/wiki/IEEE_floating_point].
+
+The typedef's are based on __N3626
+proposed for a new C++14 standard header `<cstdfloat>` and
+__N1703 proposed for a new C language standard header `<stdfloat.h>`.
+
+The 128-bit floating-point type, of great interest in scientific and
+numeric programming, is not required in the boost header,
+and may not be supplied for all platforms/compilers, because compiler
+support for a 128-bit floating-point type is not mandated by either
+the C standard or the C++ standard.
+
+The following code uses `<boost/cstdfloat.hpp>` in combination with
+`<boost/math/special_functions.hpp>` to compute a simplified
+version of the Jahnke-Emden-Lambda function. Here, we use
+a floating-point type with exactly 64 bits (i.e., `float64_t`).
+If we were to use, for instance, built-in `double`,
+then there would be no guarantee that the code would
+behave identically on all platforms. With `float64_t` from
+`<boost/cstdfloat.hpp>`, however, this is very likely.
+Using `float64_t`, we know that
+this code is portable and uses a floating-point type
+with approximately 15 decimal digits of precision.
+
+  #include <cmath>
+  #include <boost/cstdfloat.hpp>
+  #include <boost/math/special_functions.hpp>
+
+  boost::float64_t jahnke_emden_lambda(boost::float64_t v, boost::float64_t x)
+  {
+    const boost::float64_t gamma_v_plus_one = boost::math::tgamma(v + 1);
+    const boost::float64_t x_half_pow_v     = std::pow(x / 2, v);
+
+    return gamma_v_plus_one * boost::math::cyl_bessel_j(x, v) / x_half_pow_v;
+  }
+
+See `cstdfloat_test.cpp` for a more detailed test program.
+
+[endsect] [/section:overview Overview]
+
+[section:rationale Rationale]
+
+The implementation of `<boost/cstdfloat.hpp>` is designed to utilize `<float.h>`,
+defined in the 1989 C standard. The preprocessor is used to query certain
+preprocessor definitions in `<float.h>` such as FLT_MAX, DBL_MAX, etc.
+Based on the results of these queries, an attempt is made to automatically
+detect the presence of built-in floating-point types having specified widths.
+An unequivocal test regarding conformance with __IEEE754 (IEC599) based on
+[@ http://en.cppreference.com/w/cpp/types/numeric_limits/is_iec559 `std::numeric_limits<>::is_iec559`]
+is performed with `BOOST_STATIC_ASSERT`.
+
+The header `<boost/cstdfloat.hpp>` makes the standardized floating-point
+`typedef`s safely available in `namespace boost` without placing any names
+in `namespace std`. The intention is to complement rather than compete
+with a potential future C++ Standard Library that may contain these `typedef`s.
+Should some future C++ standard include `<stdfloat.h>` and `<cstdfloat>`,
+then `<boost/cstdfloat.hpp>` will continue to function, but will become redundant
+and may be safely deprecated.
+
+Because `<boost/cstdfloat.hpp>` is a boost header, its name conforms to the
+boost header naming conventions, not the C++ Standard Library header
+naming conventions.
+
+[note
+<boost/cstdfloat.hpp> [*cannot synthesize or create
+a `typedef` if the underlying type is not provided by the compiler].
+For example, if a compiler does not have an underlying floating-point
+type with 128 bits (highly sought-after in scientific and numeric programming),
+then `float128_t` and its corresponding least and fast types are not
+provided by `<boost/cstdfloat.hpp`>.]
+
+[warning
+As an implementation artifact, certain C macro names from `<float.h>`
+may possibly be visible to users of `<boost/cstdfloat.hpp>`.
+Don't rely on using these macros; they are not part of any Boost-specified interface.
+Use `std::numeric_limits<>` for floating-point ranges, etc. instead.]
+
+[endsect] [/section:rationale Rationale]
+
+[section:exact_typdefs Exact-Width Floating-Point `typedef`s]
+
+The `typedef float#_t`, with # replaced by the width, designates a
+floating-point type of exactly # bits. For example `float32_t` denotes
+a single-precision floating-point type with approximately
+7 decimal digits of precision (equivalent to binary32 in __IEEE754).
+
+Floating-point types specified in C and C++ are allowed to have
+implementation-specific widths and formats.
+However, if a platform supports underlying floating-point types
+(conformant with __IEEE754) with widths of 16, 32, 64, 128 bits,
+or any combination thereof,
+then `<boost/cstdfloat.hpp>` does provide the corresponding `typedef`s
+`float16_t, float32_t, float64_t, float128_t,`
+their corresponding least and fast types,
+and the corresponding maximum-width type
+
+The absence of `float128_t` is indicated by the macro `BOOST_NO_FLOAT128_T`.
+
+[endsect] [/section:exact_typdefs Exact-Width Floating-Point `typedef`s]
+
+
+[section:fastest_typdefs Fastest minimum-width floating-point `typedef`s]
+
+The `typedef float_least#_t`, with # replaced by the width, designates a
+floating-point type with a [*width of at least # bits], such that no
+floating-point type with lesser size has at least the specified width.
+Thus, `float_least32_t` denotes the smallest floating-point type with
+a width of at least 32 bits.
+
+Minimum-width floating-point types are provided for all existing
+exact-width floating-point types on a given platform.
+
+For example, if a platfrom supports `float32_t` and `float64_t`,
+then `float_least32_t` and `float_least64_t` will also be supported, etc.
+
+[endsect] [/section:fastest_typdefs Fastest minimum-width floating-point `typedef`s]
+
+[section:fastest_typdefs Fastest minimum-width floating-point `typedef`s]
+
+The typedef `float_fast#_t`, with # replaced by the width, designates
+the [*fastest] floating-point type with a width of at least # bits.
+
+There is no absolute guarantee that these types are the fastest for all purposes.
+In any case, however, they satisfy the precision and width requirements.
+
+Fastest minimum-width floating-point types are provided for all existing
+exact-width floating-point types on a given platform.
+
+For example, if a platform supports `float32_t` and `float64_t`,
+then `float_fast32_t` and `float_fast64_t` will also be supported, etc.
+
+[endsect] [/section:fastest_typdefs Fastest minimum-width floating-point `typedef`s]
+
+[section:greatest_typdefs Greatest-width floating-point typedef]
+
+The `typedef floatmax_t` designates a floating-point type capable of representing
+any value of any floating-point type in a given platform.
+
+The greatest-width typedef is provided for all platforms.
+
+[endsect] [/section:greatest_typdefs Greatest-width floating-point typedef]
+
+[section:macros Floating-Point Constant Macros]
+
+All macros of the type `BOOST_FLOAT16_C, BOOST_FLOAT32_C, BOOST_FLOAT64_C,
+BOOST_FLOAT128_C, BOOST_FLOATMAX_C` are always defined after inclusion of
+`<boost/cstdfloat.hpp>`. These allow floating-point constants of at
+least the specified width to be declared.
+
+For example:
+
+  #include <boost/cstdfloat.hpp>
+
+  // Declare Pythagoras' constant with approximately 7 decimal digits of precision.
+  static const boost::float32_t pi = BOOST_FLOAT32_C(3.1415926536);
+
+  // Declare the Euler-gamma constant with approximately 34 decimal digits of precision.
+  static const boost::float128_t euler = BOOST_FLOAT128_C(0.57721566490153286060651209008240243104216);
+
+[endsect] [/section:macros Floating-Point Constant Macros]
+
+
+[/ cstdfloat.qbk
+  Copyright 2014 Christopher Kormanyos, John Maddock and Paul A. Bristow.
+  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).
+]
+
+
+
+
diff --git a/libs/math/tools/ellint_d2_data.cpp b/libs/math/tools/ellint_d2_data.cpp
new file mode 100644
index 000000000..79263711b
--- /dev/null
+++ b/libs/math/tools/ellint_d2_data.cpp
@@ -0,0 +1,64 @@
+// Copyright John Maddock 2006.
+// 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 <boost/math/tools/test_data.hpp>
+#include <boost/test/included/prg_exec_monitor.hpp>
+#include <boost/math/special_functions/ellint_1.hpp>
+#include <boost/math/special_functions/ellint_2.hpp>
+#include <fstream>
+#include <boost/math/tools/test_data.hpp>
+#include "mp_t.hpp"
+
+using namespace boost::math::tools;
+using namespace boost::math;
+using namespace std;
+
+mp_t ellint_d(mp_t phi, mp_t k)
+{
+   return (boost::math::ellint_1(k, phi) - boost::math::ellint_2(k, phi)) / (k * k);
+}
+
+int cpp_main(int argc, char*argv [])
+{
+   using namespace boost::math::tools;
+
+   parameter_info<mp_t> arg1, arg2;
+   test_data<mp_t> data;
+
+   bool cont;
+   std::string line;
+
+   if(argc < 1)
+      return 1;
+
+   do{
+      if(0 == get_user_parameter_info(arg1, "phi"))
+         return 1;
+      if(0 == get_user_parameter_info(arg2, "k"))
+         return 1;
+
+      mp_t(*fp)(mp_t, mp_t) = &ellint_d;
+      data.insert(fp, arg1, arg2);
+
+      std::cout << "Any more data [y/n]?";
+      std::getline(std::cin, line);
+      boost::algorithm::trim(line);
+      cont = (line == "y");
+   }while(cont);
+
+   std::cout << "Enter name of test data file [default=ellint_d2_data.ipp]";
+   std::getline(std::cin, line);
+   boost::algorithm::trim(line);
+   if(line == "")
+      line = "ellint_d2_data.ipp";
+   std::ofstream ofs(line.c_str());
+   line.erase(line.find('.'));
+   ofs << std::scientific << std::setprecision(40);
+   write_code(ofs, data, line.c_str());
+
+   return 0;
+}
+
+
diff --git a/libs/math/tools/ellint_d_data.cpp b/libs/math/tools/ellint_d_data.cpp
new file mode 100644
index 000000000..f2c30f37c
--- /dev/null
+++ b/libs/math/tools/ellint_d_data.cpp
@@ -0,0 +1,62 @@
+// Copyright John Maddock 2006.
+// 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 <boost/math/tools/test_data.hpp>
+#include <boost/test/included/prg_exec_monitor.hpp>
+#include <boost/math/special_functions/ellint_1.hpp>
+#include <boost/math/special_functions/ellint_2.hpp>
+#include <fstream>
+#include <boost/math/tools/test_data.hpp>
+#include "mp_t.hpp"
+
+using namespace boost::math::tools;
+using namespace boost::math;
+using namespace std;
+
+mp_t ellint_d(mp_t k)
+{
+   return (boost::math::ellint_1(k) - boost::math::ellint_2(k)) / (k * k);
+}
+
+int cpp_main(int argc, char*argv [])
+{
+   using namespace boost::math::tools;
+
+   parameter_info<mp_t> arg1;
+   test_data<mp_t> data;
+
+   bool cont;
+   std::string line;
+
+   if(argc < 1)
+      return 1;
+
+   do{
+      if(0 == get_user_parameter_info(arg1, "k"))
+         return 1;
+
+      mp_t(*fp)(mp_t) = &ellint_d;
+      data.insert(fp, arg1);
+
+      std::cout << "Any more data [y/n]?";
+      std::getline(std::cin, line);
+      boost::algorithm::trim(line);
+      cont = (line == "y");
+   }while(cont);
+
+   std::cout << "Enter name of test data file [default=ellint_d_data.ipp]";
+   std::getline(std::cin, line);
+   boost::algorithm::trim(line);
+   if(line == "")
+      line = "ellint_d_data.ipp";
+   std::ofstream ofs(line.c_str());
+   line.erase(line.find('.'));
+   ofs << std::scientific << std::setprecision(40);
+   write_code(ofs, data, line.c_str());
+
+   return 0;
+}
+
+
diff --git a/libs/math/tools/heuman_lambda_data.cpp b/libs/math/tools/heuman_lambda_data.cpp
new file mode 100644
index 000000000..44aa3803b
--- /dev/null
+++ b/libs/math/tools/heuman_lambda_data.cpp
@@ -0,0 +1,67 @@
+// Copyright John Maddock 2006.
+// 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 <boost/math/tools/test_data.hpp>
+#include <boost/test/included/prg_exec_monitor.hpp>
+#include <boost/math/special_functions/ellint_1.hpp>
+#include <boost/math/special_functions/jacobi_zeta.hpp>
+#include <fstream>
+#include <boost/math/tools/test_data.hpp>
+#include "mp_t.hpp"
+
+using namespace boost::math::tools;
+using namespace boost::math;
+using namespace std;
+
+mp_t heuman_lambda(mp_t phi, mp_t k)
+{
+   mp_t kp = sqrt(1 - k *k);
+   if((k * k < tools::epsilon<float>()) && (fabs(phi) >= constants::half_pi<mp_t>()))
+      throw std::domain_error("");
+   return ellint_1(kp, phi) / ellint_1(kp) + ellint_1(k) * jacobi_zeta(kp, phi) / constants::half_pi<mp_t>();
+}
+
+int cpp_main(int argc, char*argv [])
+{
+   using namespace boost::math::tools;
+
+   parameter_info<mp_t> arg1, arg2;
+   test_data<mp_t> data;
+
+   bool cont;
+   std::string line;
+
+   if(argc < 1)
+      return 1;
+
+   do{
+      if(0 == get_user_parameter_info(arg1, "phi"))
+         return 1;
+      if(0 == get_user_parameter_info(arg2, "k"))
+         return 1;
+
+      mp_t(*fp)(mp_t, mp_t) = &heuman_lambda;
+      data.insert(fp, arg1, arg2);
+
+      std::cout << "Any more data [y/n]?";
+      std::getline(std::cin, line);
+      boost::algorithm::trim(line);
+      cont = (line == "y");
+   }while(cont);
+
+   std::cout << "Enter name of test data file [default=heuman_lambda_data.ipp]";
+   std::getline(std::cin, line);
+   boost::algorithm::trim(line);
+   if(line == "")
+      line = "heuman_lambda_data.ipp";
+   std::ofstream ofs(line.c_str());
+   line.erase(line.find('.'));
+   ofs << std::scientific << std::setprecision(40);
+   write_code(ofs, data, line.c_str());
+
+   return 0;
+}
+
+
diff --git a/libs/math/tools/jacobi_zeta_data.cpp b/libs/math/tools/jacobi_zeta_data.cpp
new file mode 100644
index 000000000..ded652301
--- /dev/null
+++ b/libs/math/tools/jacobi_zeta_data.cpp
@@ -0,0 +1,64 @@
+// Copyright John Maddock 2006.
+// 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 <boost/math/tools/test_data.hpp>
+#include <boost/test/included/prg_exec_monitor.hpp>
+#include <boost/math/special_functions/ellint_1.hpp>
+#include <boost/math/special_functions/ellint_2.hpp>
+#include <fstream>
+#include <boost/math/tools/test_data.hpp>
+#include "mp_t.hpp"
+
+using namespace boost::math::tools;
+using namespace boost::math;
+using namespace std;
+
+mp_t jacobi_zeta(mp_t phi, mp_t k)
+{
+   return ellint_2(k, phi) - ellint_2(k) * ellint_1(k, phi) / ellint_1(k);
+}
+
+int cpp_main(int argc, char*argv [])
+{
+   using namespace boost::math::tools;
+
+   parameter_info<mp_t> arg1, arg2;
+   test_data<mp_t> data;
+
+   bool cont;
+   std::string line;
+
+   if(argc < 1)
+      return 1;
+
+   do{
+      if(0 == get_user_parameter_info(arg1, "phi"))
+         return 1;
+      if(0 == get_user_parameter_info(arg2, "k"))
+         return 1;
+
+      mp_t(*fp)(mp_t, mp_t) = &jacobi_zeta;
+      data.insert(fp, arg1, arg2);
+
+      std::cout << "Any more data [y/n]?";
+      std::getline(std::cin, line);
+      boost::algorithm::trim(line);
+      cont = (line == "y");
+   }while(cont);
+
+   std::cout << "Enter name of test data file [default=jacobi_zeta_data.ipp]";
+   std::getline(std::cin, line);
+   boost::algorithm::trim(line);
+   if(line == "")
+      line = "jacobi_zeta_data.ipp";
+   std::ofstream ofs(line.c_str());
+   line.erase(line.find('.'));
+   ofs << std::scientific << std::setprecision(40);
+   write_code(ofs, data, line.c_str());
+
+   return 0;
+}
+
+
diff --git a/libs/math/tools/trig_data.cpp b/libs/math/tools/trig_data.cpp
new file mode 100644
index 000000000..41233578f
--- /dev/null
+++ b/libs/math/tools/trig_data.cpp
@@ -0,0 +1,83 @@
+//  (C) Copyright John Maddock 2006.
+//  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 <boost/math/constants/constants.hpp>
+#include <fstream>
+#include <boost/math/tools/test_data.hpp>
+#include "mp_t.hpp"
+
+using namespace boost::math::tools;
+using namespace std;
+
+float external_f;
+float force_truncate(const float* f)
+{
+   external_f = *f;
+   return external_f;
+}
+
+float truncate_to_float(mp_t r)
+{
+   float f = boost::math::tools::real_cast<float>(r);
+   return force_truncate(&f);
+}
+
+struct trig_data_generator
+{
+   boost::math::tuple<mp_t, mp_t> operator()(mp_t z)
+   {
+      return boost::math::make_tuple(sin(z * boost::math::constants::pi<mp_t>()), cos(z * boost::math::constants::pi<mp_t>()));
+   }
+};
+
+
+int main(int argc, char*argv [])
+{
+   parameter_info<mp_t> arg1;
+   test_data<mp_t> data;
+
+   bool cont;
+   std::string line;
+
+   std::cout << "Welcome.\n"
+      "This program will generate spot tests for the cos_pi and sin_pi functions:\n";
+
+   do{
+      if(0 == get_user_parameter_info(arg1, "a"))
+         return 1;
+      data.insert(trig_data_generator(), arg1);
+
+      std::cout << "Any more data [y/n]?";
+      std::getline(std::cin, line);
+      boost::algorithm::trim(line);
+      cont = (line == "y");
+   }while(cont);
+
+   std::cout << "Enter name of test data file [default=trig_data.ipp]";
+   std::getline(std::cin, line);
+   boost::algorithm::trim(line);
+   if(line == "")
+      line = "trig_data.ipp";
+   std::ofstream ofs(line.c_str());
+   ofs << std::scientific << std::setprecision(40);
+   write_code(ofs, data, "trig_data");
+   
+   return 0;
+}
+
+/* Output for asymptotic limits:
+
+Erf asymptotic limit for 24 bit numbers is 2.8 after approximately 6 terms.
+Erfc asymptotic limit for 24 bit numbers is 4.12064 after approximately 17 terms.
+Erf asymptotic limit for 53 bit numbers is 4.3 after approximately 11 terms.
+Erfc asymptotic limit for 53 bit numbers is 6.19035 after approximately 29 terms.
+Erf asymptotic limit for 64 bit numbers is 4.8 after approximately 12 terms.
+Erfc asymptotic limit for 64 bit numbers is 7.06004 after approximately 29 terms.
+Erf asymptotic limit for 106 bit numbers is 6.5 after approximately 14 terms.
+Erfc asymptotic limit for 106 bit numbers is 11.6626 after approximately 29 terms.
+Erf asymptotic limit for 113 bit numbers is 6.8 after approximately 14 terms.
+Erfc asymptotic limit for 113 bit numbers is 12.6802 after approximately 29 terms.
+*/
+