Adding upstream version 14.2.21.upstream/14.2.21upstream

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

6 files changed, 2177 insertions, 0 deletions

diff --git a/src/boost/libs/math/include_private/boost/math/constants/generate.hpp b/src/boost/libs/math/include_private/boost/math/constants/generate.hpp
new file mode 100644
index 00000000..2c29be71
--- /dev/null
+++ b/src/boost/libs/math/include_private/boost/math/constants/generate.hpp
@@ -0,0 +1,76 @@
+//  Copyright John Maddock 2010.
+//  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)
+
+#ifndef BOOST_MATH_CONSTANTS_GENERATE_INCLUDED
+#define BOOST_MATH_CONSTANTS_GENERATE_INCLUDED
+
+#include <boost/math/constants/constants.hpp>
+#include <boost/regex.hpp>
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+
+#ifdef USE_MPFR
+#include <boost/math/bindings/mpfr.hpp>
+#elif defined(USE_MPREAL)
+#include <boost/math/bindings/mpreal.hpp>
+#elif defined(USE_CPP_FLOAT)
+#include <boost/multiprecision/cpp_dec_float.hpp>
+#else
+#include <boost/math/bindings/rr.hpp>
+#endif
+
+namespace boost{ namespace math{ namespace constants{ 
+
+#ifdef USE_MPFR
+typedef mpfr_class generator_type;
+#elif defined(USE_MPREAL)
+typedef mpfr::mpreal generator_type;
+#elif defined(USE_CPP_FLOAT)
+typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<500> > generator_type;
+#else
+typedef ntl::RR generator_type;
+#endif
+
+inline void print_constant(const char* name, generator_type(*f)(const mpl::int_<0>&))
+{
+#ifdef USE_MPFR
+   mpfr_class::set_dprec(((200 + 1) * 1000L) / 301L);
+#elif defined(USE_MPREAL)
+   mpfr::mpreal::set_default_prec(((200 + 1) * 1000L) / 301L);
+#elif defined(USE_CPP_FLOAT)
+   // Nothing to do, precision is already set.
+#else
+   ntl::RR::SetPrecision(((200 + 1) * 1000L) / 301L);
+   ntl::RR::SetOutputPrecision(102);
+#endif
+   generator_type value = f(boost::mpl::int_<0>());
+   std::stringstream os;
+   os << std::setprecision(110) << std::scientific;
+   os << value;
+   std::string s = os.str();
+   static const regex e("([+-]?\\d+(?:\\.\\d{0,36})?)(\\d*)(?:e([+-]?\\d+))?");
+   smatch what;
+   if(regex_match(s, what, e))
+   {
+      std::cout << 
+         "BOOST_DEFINE_MATH_CONSTANT(" << name << ", " 
+         << what[1] << "e" << (what[3].length() ? what[3].str() : std::string("0")) << ", " 
+         << "\"" << what[1] << what[2] << "e" << (what[3].length() ? what[3].str() : std::string("0")) 
+         << "\");" << std::endl;
+   }
+   else
+   {
+      std::cout << "Format of numeric constant was not recognised!!" << std::endl;
+   }
+}
+
+#define BOOST_CONSTANTS_GENERATE(name) \
+   boost::math::constants::print_constant(#name, \
+   & boost::math::constants::detail::BOOST_JOIN(constant_, name)<boost::math::constants::generator_type>::get)
+
+}}} // namespaces
+
+#endif // BOOST_MATH_CONSTANTS_GENERATE_INCLUDED
diff --git a/src/boost/libs/math/include_private/boost/math/tools/iteration_logger.hpp b/src/boost/libs/math/include_private/boost/math/tools/iteration_logger.hpp
new file mode 100644
index 00000000..c272a17d
--- /dev/null
+++ b/src/boost/libs/math/include_private/boost/math/tools/iteration_logger.hpp
@@ -0,0 +1,77 @@
+//  Copyright John Maddock 2014.
+//  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)
+
+#ifndef BOOST_MATH_TOOLS_ITERATION_LOGGER
+#define BOOST_MATH_TOOLS_ITERATION_LOGGER
+
+#include <map>
+#include <string>
+#include <utility>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <boost/current_function.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+struct logger
+{
+   std::map<std::pair<std::string, std::string>, unsigned long long> data;
+
+   ~logger()
+   {
+      for(std::map<std::pair<std::string, std::string>, unsigned long long>::const_iterator i = data.begin(); i != data.end(); ++i)
+      {
+         std::cout << "Logging iteration data:\n  file: " << i->first.first
+            << "\n  function: " << i->first.second
+            << "\n  count: " << i->second << std::endl;
+         //
+         // Read in existing data:
+         //
+         std::map<std::string, unsigned long long> file_data;
+         std::string filename = i->first.first + ".logfile";
+         std::ifstream is(filename.c_str());
+         while(is.good())
+         {
+            std::string f;
+            unsigned long long c;
+            is >> std::ws;
+            std::getline(is, f);
+            is >> c;
+            if(f.size())
+               file_data[f] = c;
+         }
+         is.close();
+         if(file_data.find(i->first.second) != file_data.end())
+            file_data[i->first.second] += i->second;
+         else
+            file_data[i->first.second] = i->second;
+         //
+         // Write it out again:
+         //
+         std::ofstream os(filename.c_str());
+         for(std::map<std::string, unsigned long long>::const_iterator j = file_data.begin(); j != file_data.end(); ++j)
+            os << j->first << "\n    " << j->second << std::endl;
+         os.close();
+      }
+   }
+};
+
+inline void log_iterations(const char* file, const char* function, unsigned long long count)
+{
+   static logger l;
+   std::pair<std::string, std::string> key(file, function);
+   if(l.data.find(key) == l.data.end())
+      l.data[key] = 0;
+   l.data[key] += count;
+}
+
+#define BOOST_MATH_LOG_COUNT(count) boost::math::detail::log_iterations(__FILE__, BOOST_CURRENT_FUNCTION, count);
+
+
+}}} // namespaces
+
+#endif // BOOST_MATH_TOOLS_ITERATION_LOGGER
+
diff --git a/src/boost/libs/math/include_private/boost/math/tools/remez.hpp b/src/boost/libs/math/include_private/boost/math/tools/remez.hpp
new file mode 100644
index 00000000..44a3412f
--- /dev/null
+++ b/src/boost/libs/math/include_private/boost/math/tools/remez.hpp
@@ -0,0 +1,667 @@
+//  (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)
+
+#ifndef BOOST_MATH_TOOLS_REMEZ_HPP
+#define BOOST_MATH_TOOLS_REMEZ_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <boost/math/tools/solve.hpp>
+#include <boost/math/tools/minima.hpp>
+#include <boost/math/tools/roots.hpp>
+#include <boost/math/tools/polynomial.hpp>
+#include <boost/function/function1.hpp>
+#include <boost/scoped_array.hpp>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/policies/policy.hpp>
+
+namespace boost{ namespace math{ namespace tools{
+
+namespace detail{
+
+//
+// The error function: the difference between F(x) and
+// the current approximation.  This is the function
+// for which we must find the extema.
+//
+template <class T>
+struct remez_error_function
+{
+   typedef boost::function1<T, T const &> function_type;
+public:
+   remez_error_function(
+      function_type f_, 
+      const polynomial<T>& n, 
+      const polynomial<T>& d, 
+      bool rel_err)
+         : f(f_), numerator(n), denominator(d), rel_error(rel_err) {}
+
+   T operator()(const T& z)const
+   {
+      T y = f(z);
+      T abs = y - (numerator.evaluate(z) / denominator.evaluate(z));
+      T err;
+      if(rel_error)
+      {
+         if(y != 0)
+            err = abs / fabs(y);
+         else if(0 == abs)
+         {
+            // we must be at a root, or it's not recoverable:
+            BOOST_ASSERT(0 == abs);
+            err = 0;
+         }
+         else
+         {
+            // We have a divide by zero!
+            // Lets assume that f(x) is zero as a result of
+            // internal cancellation error that occurs as a result
+            // of shifting a root at point z to the origin so that
+            // the approximation can be "pinned" to pass through
+            // the origin: in that case it really
+            // won't matter what our approximation calculates here
+            // as long as it's a small number, return the absolute error:
+            err = abs;
+         }
+      }
+      else
+         err = abs;
+      return err;
+   }
+private:
+   function_type f;
+   polynomial<T> numerator;
+   polynomial<T> denominator;
+   bool rel_error;
+};
+//
+// This function adapts the error function so that it's minima
+// are the extema of the error function.  We can find the minima
+// with standard techniques.
+//
+template <class T>
+struct remez_max_error_function
+{
+   remez_max_error_function(const remez_error_function<T>& f)
+      : func(f) {}
+
+   T operator()(const T& x)
+   {
+      BOOST_MATH_STD_USING
+      return -fabs(func(x));
+   }
+private:
+   remez_error_function<T> func;
+};
+
+} // detail
+
+template <class T>
+class remez_minimax
+{
+public:
+   typedef boost::function1<T, T const &> function_type;
+   typedef boost::numeric::ublas::vector<T> vector_type;
+   typedef boost::numeric::ublas::matrix<T> matrix_type;
+
+   remez_minimax(function_type f, unsigned oN, unsigned oD, T a, T b, bool pin = true, bool rel_err = false, int sk = 0, int bits = 0);
+   remez_minimax(function_type f, unsigned oN, unsigned oD, T a, T b, bool pin, bool rel_err, int sk, int bits, const vector_type& points);
+
+   void reset(unsigned oN, unsigned oD, T a, T b, bool pin = true, bool rel_err = false, int sk = 0, int bits = 0);
+   void reset(unsigned oN, unsigned oD, T a, T b, bool pin, bool rel_err, int sk, int bits, const vector_type& points);
+
+   void set_brake(int b)
+   {
+      BOOST_ASSERT(b < 100);
+      BOOST_ASSERT(b >= 0);
+      m_brake = b;
+   }
+
+   T iterate();
+
+   polynomial<T> denominator()const;
+   polynomial<T> numerator()const;
+
+   vector_type const& chebyshev_points()const
+   {
+      return control_points;
+   }
+
+   vector_type const& zero_points()const
+   {
+      return zeros;
+   }
+
+   T error_term()const
+   {
+      return solution[solution.size() - 1];
+   }
+   T max_error()const
+   {
+      return m_max_error;
+   }
+   T max_change()const
+   {
+      return m_max_change;
+   }
+   void rotate()
+   {
+      --orderN;
+      ++orderD;
+   }
+   void rescale(T a, T b)
+   {
+      T scale = (b - a) / (max - min);
+      for(unsigned i = 0; i < control_points.size(); ++i)
+      {
+         control_points[i] = (control_points[i] - min) * scale + a;
+      }
+      min = a;
+      max = b;
+   }
+private:
+
+   void init_chebyshev();
+
+   function_type func;            // The function to approximate.
+   vector_type control_points;    // Current control points to be used for the next iteration.
+   vector_type solution;          // Solution from the last iteration contains all unknowns including the error term.
+   vector_type zeros;             // Location of points of zero error from last iteration, plus the two end points.
+   vector_type maxima;            // Location of maxima of the error function, actually contains the control points used for the last iteration.
+   T m_max_error;                 // Maximum error found in last approximation.
+   T m_max_change;                // Maximum change in location of control points after last iteration.
+   unsigned orderN;               // Order of the numerator polynomial.
+   unsigned orderD;               // Order of the denominator polynomial.
+   T min, max;                    // End points of the range to optimise over.
+   bool rel_error;                // If true optimise for relative not absolute error.
+   bool pinned;                   // If true the approximation is "pinned" to go through the origin.
+   unsigned unknowns;             // Total number of unknowns.
+   int m_precision;               // Number of bits precision to which the zeros and maxima are found.
+   T m_max_change_history[2];     // Past history of changes to control points.
+   int m_brake;                     // amount to break by in percentage points.
+   int m_skew;                      // amount to skew starting points by in percentage points: -100-100
+};
+
+#ifndef BRAKE
+#define BRAKE 0
+#endif
+#ifndef SKEW
+#define SKEW 0
+#endif
+
+template <class T>
+void remez_minimax<T>::init_chebyshev()
+{
+   BOOST_MATH_STD_USING
+   //
+   // Fill in the zeros:
+   //
+   unsigned terms = pinned ? orderD + orderN : orderD + orderN + 1;
+
+   for(unsigned i = 0; i < terms; ++i)
+   {
+      T cheb = cos((2 * terms - 1 - 2 * i) * constants::pi<T>() / (2 * terms));
+      cheb += 1;
+      cheb /= 2;
+      if(m_skew != 0)
+      {
+         T p = static_cast<T>(200 + m_skew) / 200;
+         cheb = pow(cheb, p);
+      }
+      cheb *= (max - min);
+      cheb += min;
+      zeros[i+1] = cheb;
+   }
+   zeros[0] = min;
+   zeros[unknowns] = max;
+   // perform a regular interpolation fit:
+   matrix_type A(terms, terms);
+   vector_type b(terms);
+   // fill in the y values:
+   for(unsigned i = 0; i < b.size(); ++i)
+   {
+      b[i] = func(zeros[i+1]);
+   }
+   // fill in powers of x evaluated at each of the control points:
+   unsigned offsetN = pinned ? 0 : 1;
+   unsigned offsetD = offsetN + orderN;
+   unsigned maxorder = (std::max)(orderN, orderD);
+   for(unsigned i = 0; i < b.size(); ++i)
+   {
+      T x0 = zeros[i+1];
+      T x = x0;
+      if(!pinned)
+         A(i, 0) = 1;
+      for(unsigned j = 0; j < maxorder; ++j)
+      {
+         if(j < orderN)
+            A(i, j + offsetN) = x;
+         if(j < orderD)
+         {
+            A(i, j + offsetD) = -x * b[i];
+         }
+         x *= x0;
+      }
+   }
+   //
+   // Now go ahead and solve the expression to get our solution:
+   //
+   vector_type l_solution = boost::math::tools::solve(A, b);
+   // need to add a "fake" error term:
+   l_solution.resize(unknowns);
+   l_solution[unknowns-1] = 0;
+   solution = l_solution;
+   //
+   // Now find all the extrema of the error function:
+   //
+   detail::remez_error_function<T> Err(func, this->numerator(), this->denominator(), rel_error);
+   detail::remez_max_error_function<T> Ex(Err);
+   m_max_error = 0;
+   //int max_err_location = 0;
+   for(unsigned i = 0; i < unknowns; ++i)
+   {
+      std::pair<T, T> r = brent_find_minima(Ex, zeros[i], zeros[i+1], m_precision);
+      maxima[i] = r.first;
+      T rel_err = fabs(r.second);
+      if(rel_err > m_max_error)
+      {
+         m_max_error = fabs(r.second);
+         //max_err_location = i;
+      }
+   }
+   control_points = maxima;
+}
+
+template <class T>
+void remez_minimax<T>::reset(
+         unsigned oN, 
+         unsigned oD, 
+         T a, 
+         T b, 
+         bool pin, 
+         bool rel_err, 
+         int sk,
+         int bits)
+{
+   control_points = vector_type(oN + oD + (pin ? 1 : 2));
+   solution = control_points;
+   zeros = vector_type(oN + oD + (pin ? 2 : 3));
+   maxima = control_points;
+   orderN = oN;
+   orderD = oD;
+   rel_error = rel_err;
+   pinned = pin;
+   m_skew = sk;
+   min = a;
+   max = b;
+   m_max_error = 0;
+   unknowns = orderN + orderD + (pinned ? 1 : 2);
+   // guess our initial control points:
+   control_points[0] = min;
+   control_points[unknowns - 1] = max;
+   T interval = (max - min) / (unknowns - 1);
+   T spot = min + interval;
+   for(unsigned i = 1; i < control_points.size(); ++i)
+   {
+      control_points[i] = spot;
+      spot += interval;
+   }
+   solution[unknowns - 1] = 0;
+   m_max_error = 0;
+   if(bits == 0)
+   {
+      // don't bother about more than float precision:
+      m_precision = (std::min)(24, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2);
+   }
+   else
+   {
+      // can't be more accurate than half the bits of T:
+      m_precision = (std::min)(bits, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2);
+   }
+   m_max_change_history[0] = m_max_change_history[1] = 1;
+   init_chebyshev();
+   // do one iteration whatever:
+   //iterate();
+}
+
+template <class T>
+inline remez_minimax<T>::remez_minimax(
+         typename remez_minimax<T>::function_type f, 
+         unsigned oN, 
+         unsigned oD, 
+         T a, 
+         T b, 
+         bool pin, 
+         bool rel_err, 
+         int sk,
+         int bits)
+   : func(f) 
+{
+   m_brake = 0;
+   reset(oN, oD, a, b, pin, rel_err, sk, bits);
+}
+
+template <class T>
+void remez_minimax<T>::reset(
+         unsigned oN, 
+         unsigned oD, 
+         T a, 
+         T b, 
+         bool pin, 
+         bool rel_err, 
+         int sk,
+         int bits,
+         const vector_type& points)
+{
+   control_points = vector_type(oN + oD + (pin ? 1 : 2));
+   solution = control_points;
+   zeros = vector_type(oN + oD + (pin ? 2 : 3));
+   maxima = control_points;
+   orderN = oN;
+   orderD = oD;
+   rel_error = rel_err;
+   pinned = pin;
+   m_skew = sk;
+   min = a;
+   max = b;
+   m_max_error = 0;
+   unknowns = orderN + orderD + (pinned ? 1 : 2);
+   control_points = points;
+   solution[unknowns - 1] = 0;
+   m_max_error = 0;
+   if(bits == 0)
+   {
+      // don't bother about more than float precision:
+      m_precision = (std::min)(24, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2);
+   }
+   else
+   {
+      // can't be more accurate than half the bits of T:
+      m_precision = (std::min)(bits, (boost::math::policies::digits<T, boost::math::policies::policy<> >() / 2) - 2);
+   }
+   m_max_change_history[0] = m_max_change_history[1] = 1;
+   // do one iteration whatever:
+   //iterate();
+}
+
+template <class T>
+inline remez_minimax<T>::remez_minimax(
+         typename remez_minimax<T>::function_type f, 
+         unsigned oN, 
+         unsigned oD, 
+         T a, 
+         T b, 
+         bool pin, 
+         bool rel_err, 
+         int sk,
+         int bits,
+         const vector_type& points)
+   : func(f)
+{
+   m_brake = 0;
+   reset(oN, oD, a, b, pin, rel_err, sk, bits, points);
+}
+
+template <class T>
+T remez_minimax<T>::iterate()
+{
+   BOOST_MATH_STD_USING
+   matrix_type A(unknowns, unknowns);
+   vector_type b(unknowns);
+
+   // fill in evaluation of f(x) at each of the control points:
+   for(unsigned i = 0; i < b.size(); ++i)
+   {
+      // take care that none of our control points are at the origin:
+      if(pinned && (control_points[i] == 0))
+      {
+         if(i)
+            control_points[i] = control_points[i-1] / 3;
+         else
+            control_points[i] = control_points[i+1] / 3;
+      }
+      b[i] = func(control_points[i]);
+   }
+
+   T err_err;
+   unsigned convergence_count = 0;
+   do{
+      // fill in powers of x evaluated at each of the control points:
+      int sign = 1;
+      unsigned offsetN = pinned ? 0 : 1;
+      unsigned offsetD = offsetN + orderN;
+      unsigned maxorder = (std::max)(orderN, orderD);
+      T Elast = solution[unknowns - 1];
+
+      for(unsigned i = 0; i < b.size(); ++i)
+      {
+         T x0 = control_points[i];
+         T x = x0;
+         if(!pinned)
+            A(i, 0) = 1;
+         for(unsigned j = 0; j < maxorder; ++j)
+         {
+            if(j < orderN)
+               A(i, j + offsetN) = x;
+            if(j < orderD)
+            {
+               T mult = rel_error ? T(b[i] - sign * fabs(b[i]) * Elast): T(b[i] - sign * Elast);
+               A(i, j + offsetD) = -x * mult;
+            }
+            x *= x0;
+         }
+         // The last variable to be solved for is the error term, 
+         // sign changes with each control point:
+         T E = rel_error ? T(sign * fabs(b[i])) : T(sign);
+         A(i, unknowns - 1) = E;
+         sign = -sign;
+      }
+
+   #ifdef BOOST_MATH_INSTRUMENT
+      for(unsigned i = 0; i < b.size(); ++i)
+         std::cout << b[i] << " ";
+      std::cout << "\n\n";
+      for(unsigned i = 0; i < b.size(); ++i)
+      {
+         for(unsigned j = 0; j < b.size(); ++ j)
+            std::cout << A(i, j) << " ";
+         std::cout << "\n";
+      }
+      std::cout << std::endl;
+   #endif
+      //
+      // Now go ahead and solve the expression to get our solution:
+      //
+      solution = boost::math::tools::solve(A, b);
+
+      err_err = (Elast != 0) ? T(fabs((fabs(solution[unknowns-1]) - fabs(Elast)) / fabs(Elast))) : T(1);
+   }while(orderD && (convergence_count++ < 80) && (err_err > 0.001));
+
+   //
+   // Perform a sanity check to verify that the solution to the equations
+   // is not so much in error as to be useless.  The matrix inversion can
+   // be very close to singular, so this can be a real problem.
+   //
+   vector_type sanity = prod(A, solution);
+   for(unsigned i = 0; i < b.size(); ++i)
+   {
+      T err = fabs((b[i] - sanity[i]) / fabs(b[i]));
+      if(err > sqrt(epsilon<T>()))
+      {
+         std::cerr << "Sanity check failed: more than half the digits in the found solution are in error." << std::endl;
+      }
+   }
+
+   //
+   // Next comes another sanity check, we want to verify that all the control
+   // points do actually alternate in sign, in practice we may have 
+   // additional roots in the error function that cause this to fail.
+   // Failure here is always fatal: even though this code attempts to correct
+   // the problem it usually only postpones the inevitable.
+   //
+   polynomial<T> num, denom;
+   num = this->numerator();
+   denom = this->denominator();
+   T e1 = b[0] - num.evaluate(control_points[0]) / denom.evaluate(control_points[0]);
+#ifdef BOOST_MATH_INSTRUMENT
+   std::cout << e1;
+#endif
+   for(unsigned i = 1; i < b.size(); ++i)
+   {
+      T e2 = b[i] - num.evaluate(control_points[i]) / denom.evaluate(control_points[i]);
+#ifdef BOOST_MATH_INSTRUMENT
+      std::cout << " " << e2;
+#endif
+      if(e2 * e1 > 0)
+      {
+         std::cerr << std::flush << "Basic sanity check failed: Error term does not alternate in sign, non-recoverable error may follow..." << std::endl;
+         T perturbation = 0.05;
+         do{
+            T point = control_points[i] * (1 - perturbation) + control_points[i-1] * perturbation;
+            e2 = func(point) - num.evaluate(point) / denom.evaluate(point);
+            if(e2 * e1 < 0)
+            {
+               control_points[i] = point;
+               break;
+            }
+            perturbation += 0.05;
+         }while(perturbation < 0.8);
+
+         if((e2 * e1 > 0) && (i + 1 < b.size()))
+         {
+            perturbation = 0.05;
+            do{
+               T point = control_points[i] * (1 - perturbation) + control_points[i+1] * perturbation;
+               e2 = func(point) - num.evaluate(point) / denom.evaluate(point);
+               if(e2 * e1 < 0)
+               {
+                  control_points[i] = point;
+                  break;
+               }
+               perturbation += 0.05;
+            }while(perturbation < 0.8);
+         }
+
+      }
+      e1 = e2;
+   }
+
+#ifdef BOOST_MATH_INSTRUMENT
+   for(unsigned i = 0; i < solution.size(); ++i)
+      std::cout << solution[i] << " ";
+   std::cout << std::endl << this->numerator() << std::endl;
+   std::cout << this->denominator() << std::endl;
+   std::cout << std::endl;
+#endif
+
+   //
+   // The next step is to find all the intervals in which our maxima
+   // lie:
+   //
+   detail::remez_error_function<T> Err(func, this->numerator(), this->denominator(), rel_error);
+   zeros[0] = min;
+   zeros[unknowns] = max;
+   for(unsigned i = 1; i < control_points.size(); ++i)
+   {
+      eps_tolerance<T> tol(m_precision);
+      boost::uintmax_t max_iter = 1000;
+      std::pair<T, T> p = toms748_solve(
+         Err, 
+         control_points[i-1], 
+         control_points[i], 
+         tol, 
+         max_iter);
+      zeros[i] = (p.first + p.second) / 2;
+      //zeros[i] = bisect(Err, control_points[i-1], control_points[i], m_precision);
+   }
+   //
+   // Now find all the extrema of the error function:
+   //
+   detail::remez_max_error_function<T> Ex(Err);
+   m_max_error = 0;
+   //int max_err_location = 0;
+   for(unsigned i = 0; i < unknowns; ++i)
+   {
+      std::pair<T, T> r = brent_find_minima(Ex, zeros[i], zeros[i+1], m_precision);
+      maxima[i] = r.first;
+      T rel_err = fabs(r.second);
+      if(rel_err > m_max_error)
+      {
+         m_max_error = fabs(r.second);
+         //max_err_location = i;
+      }
+   }
+   //
+   // Almost done now! we just need to set our control points
+   // to the extrema, and calculate how much each point has changed
+   // (this will be our termination condition):
+   //
+   swap(control_points, maxima);
+   m_max_change = 0;
+   //int max_change_location = 0;
+   for(unsigned i = 0; i < unknowns; ++i)
+   {
+      control_points[i] = (control_points[i] * (100 - m_brake) + maxima[i] * m_brake) / 100;
+      T change = fabs((control_points[i] - maxima[i]) / control_points[i]);
+#if 0
+      if(change > m_max_change_history[1])
+      {
+         // divergence!!! try capping the change:
+         std::cerr << "Possible divergent step, change will be capped!!" << std::endl;
+         change = m_max_change_history[1];
+         if(control_points[i] < maxima[i])
+            control_points[i] = maxima[i] - change * maxima[i];
+         else
+            control_points[i] = maxima[i] + change * maxima[i];
+      }
+#endif
+      if(change > m_max_change)
+      {
+         m_max_change = change;
+         //max_change_location = i;
+      }
+   }
+   //
+   // store max change information:
+   //
+   m_max_change_history[0] = m_max_change_history[1];
+   m_max_change_history[1] = fabs(m_max_change);
+
+   return m_max_change;
+}
+
+template <class T>
+polynomial<T> remez_minimax<T>::numerator()const
+{
+   boost::scoped_array<T> a(new T[orderN + 1]);
+   if(pinned)
+      a[0] = 0;
+   unsigned terms = pinned ? orderN : orderN + 1;
+   for(unsigned i = 0; i < terms; ++i)
+      a[pinned ? i+1 : i] = solution[i];
+   return boost::math::tools::polynomial<T>(&a[0], orderN);
+}
+
+template <class T>
+polynomial<T> remez_minimax<T>::denominator()const
+{
+   unsigned terms = orderD + 1;
+   unsigned offsetD = pinned ? orderN : (orderN + 1);
+   boost::scoped_array<T> a(new T[terms]);
+   a[0] = 1;
+   for(unsigned i = 0; i < orderD; ++i)
+      a[i+1] = solution[i + offsetD];
+   return boost::math::tools::polynomial<T>(&a[0], orderD);
+}
+
+
+}}} // namespaces
+
+#endif // BOOST_MATH_TOOLS_REMEZ_HPP
+
+
+
diff --git a/src/boost/libs/math/include_private/boost/math/tools/solve.hpp b/src/boost/libs/math/include_private/boost/math/tools/solve.hpp
new file mode 100644
index 00000000..78acde26
--- /dev/null
+++ b/src/boost/libs/math/include_private/boost/math/tools/solve.hpp
@@ -0,0 +1,79 @@
+//  (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)
+
+#ifndef BOOST_MATH_TOOLS_SOLVE_HPP
+#define BOOST_MATH_TOOLS_SOLVE_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <boost/config.hpp>
+#include <boost/assert.hpp>
+
+#ifdef BOOST_MSVC
+#pragma warning(push)
+#pragma warning(disable:4996 4267 4244)
+#endif
+
+#include <boost/numeric/ublas/lu.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+
+#ifdef BOOST_MSVC
+#pragma warning(pop)
+#endif
+
+namespace boost{ namespace math{ namespace tools{
+
+//
+// Find x such that Ax = b
+//
+// Caution: this uses undocumented, and untested ublas code,
+// however short of writing our own LU-decompostion code
+// it's the only game in town.
+//
+template <class T>
+boost::numeric::ublas::vector<T> solve(
+          const boost::numeric::ublas::matrix<T>& A_,
+          const boost::numeric::ublas::vector<T>& b_)
+{
+   //BOOST_ASSERT(A_.size() == b_.size());
+
+   boost::numeric::ublas::matrix<T> A(A_);
+   boost::numeric::ublas::vector<T> b(b_);
+   boost::numeric::ublas::permutation_matrix<> piv(b.size());
+   lu_factorize(A, piv);
+   lu_substitute(A, piv, b);
+   //
+   // iterate to reduce error:
+   //
+   boost::numeric::ublas::vector<T> delta(b.size());
+   for(unsigned k = 0; k < 1; ++k)
+   {
+      noalias(delta) = prod(A_, b);
+      delta -= b_;
+      lu_substitute(A, piv, delta);
+      b -= delta;
+
+      T max_error = 0;
+
+      for(unsigned i = 0; i < delta.size(); ++i)
+      {
+         T err = fabs(delta[i] / b[i]);
+         if(err > max_error)
+            max_error = err;
+      }
+      //std::cout << "Max change in LU error correction: " << max_error << std::endl;
+   }
+
+   return b;
+}
+
+}}} // namespaces
+
+#endif // BOOST_MATH_TOOLS_SOLVE_HPP
+
+
diff --git a/src/boost/libs/math/include_private/boost/math/tools/test.hpp b/src/boost/libs/math/include_private/boost/math/tools/test.hpp
new file mode 100644
index 00000000..14b7e2b4
--- /dev/null
+++ b/src/boost/libs/math/include_private/boost/math/tools/test.hpp
@@ -0,0 +1,311 @@
+//  (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)
+
+#ifndef BOOST_MATH_TOOLS_TEST_HPP
+#define BOOST_MATH_TOOLS_TEST_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <boost/math/tools/config.hpp>
+#include <boost/math/tools/stats.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/special_functions/relative_difference.hpp>
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/test/test_tools.hpp>
+#include <stdexcept>
+#include <iostream>
+#include <iomanip>
+
+namespace boost{ namespace math{ namespace tools{
+
+template <class T>
+struct test_result
+{
+private:
+   boost::math::tools::stats<T> stat;   // Statistics for the test.
+   unsigned worst_case;                 // Index of the worst case test.
+public:
+   test_result() { worst_case = 0; }
+   void set_worst(int i){ worst_case = i; }
+   void add(const T& point){ stat.add(point); }
+   // accessors:
+   unsigned worst()const{ return worst_case; }
+   T min BOOST_PREVENT_MACRO_SUBSTITUTION()const{ return (stat.min)(); }
+   T max BOOST_PREVENT_MACRO_SUBSTITUTION()const{ return (stat.max)(); }
+   T total()const{ return stat.total(); }
+   T mean()const{ return stat.mean(); }
+   boost::uintmax_t count()const{ return stat.count(); }
+   T variance()const{ return stat.variance(); }
+   T variance1()const{ return stat.variance1(); }
+   T rms()const{ return stat.rms(); }
+
+   test_result& operator+=(const test_result& t)
+   {
+      if((t.stat.max)() > (stat.max)())
+         worst_case = t.worst_case;
+      stat += t.stat;
+      return *this;
+   }
+};
+
+template <class T>
+struct calculate_result_type
+{
+   typedef typename T::value_type row_type;
+   typedef typename row_type::value_type value_type;
+};
+
+template <class T>
+T relative_error(T a, T b)
+{
+   return boost::math::relative_difference(a, b);
+}
+
+
+template <class T>
+void set_output_precision(T, std::ostream& os)
+{
+#ifdef BOOST_MSVC
+#pragma warning(push)
+#pragma warning(disable:4127)
+#endif
+   if(std::numeric_limits<T>::digits10)
+   {
+      os << std::setprecision(std::numeric_limits<T>::digits10 + 2);
+   }
+#ifdef BOOST_MSVC
+#pragma warning(pop)
+#endif
+}
+
+template <class Seq>
+void print_row(const Seq& row, std::ostream& os = std::cout)
+{
+   try {
+      set_output_precision(row[0], os);
+      for (unsigned i = 0; i < row.size(); ++i)
+      {
+         if (i)
+            os << ", ";
+         os << row[i];
+      }
+      os << std::endl;
+   }
+   catch (const std::exception&) {}
+}
+
+//
+// Function test accepts an matrix of input values (probably a 2D boost::array)
+// and calls two functors for each row in the array - one calculates a value
+// to test, and one extracts the expected value from the array (or possibly
+// calculates it at high precision).  The two functors are usually simple lambda
+// expressions.
+//
+template <class A, class F1, class F2>
+test_result<typename calculate_result_type<A>::value_type> test(const A& a, F1 test_func, F2 expect_func)
+{
+   typedef typename A::value_type         row_type;
+   typedef typename row_type::value_type  value_type;
+
+   test_result<value_type> result;
+
+   for(unsigned i = 0; i < a.size(); ++i)
+   {
+      const row_type& row = a[i];
+      value_type point;
+#ifndef BOOST_NO_EXCEPTIONS
+      try
+      {
+#endif
+         point = test_func(row);
+#ifndef BOOST_NO_EXCEPTIONS
+      }
+      catch(const std::underflow_error&)
+      {
+         point = 0;
+      }
+      catch(const std::overflow_error&)
+      {
+         point = std::numeric_limits<value_type>::has_infinity ? 
+            std::numeric_limits<value_type>::infinity()
+            : tools::max_value<value_type>();
+      }
+      catch(const std::exception& e)
+      {
+         std::cerr << e.what() << std::endl;
+         print_row(row, std::cerr);
+         BOOST_ERROR("Unexpected exception.");
+         // so we don't get further errors:
+         point = expect_func(row);
+      }
+#endif
+      value_type expected = expect_func(row);
+      value_type err = relative_error(point, expected);
+#ifdef BOOST_INSTRUMENT
+      if(err != 0)
+      {
+         std::cout << row[0] << " " << err;
+         if(std::numeric_limits<value_type>::is_specialized)
+         {
+            std::cout << " (" << err / std::numeric_limits<value_type>::epsilon() << "eps)";
+         }
+         std::cout << std::endl;
+      }
+#endif
+      if(!(boost::math::isfinite)(point) && (boost::math::isfinite)(expected))
+      {
+         std::cerr << "CAUTION: Found non-finite result, when a finite value was expected at entry " << i << "\n";
+         std::cerr << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl;
+         print_row(row, std::cerr);
+         BOOST_ERROR("Unexpected non-finite result");
+      }
+      if(err > 0.5)
+      {
+         std::cerr << "CAUTION: Gross error found at entry " << i << ".\n";
+         std::cerr << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl;
+         print_row(row, std::cerr);
+         BOOST_ERROR("Gross error");
+      }
+      result.add(err);
+      if((result.max)() == err)
+         result.set_worst(i);
+   }
+   return result;
+}
+
+template <class Real, class A, class F1, class F2>
+test_result<Real> test_hetero(const A& a, F1 test_func, F2 expect_func)
+{
+   typedef typename A::value_type         row_type;
+   typedef Real                          value_type;
+
+   test_result<value_type> result;
+
+   for(unsigned i = 0; i < a.size(); ++i)
+   {
+      const row_type& row = a[i];
+      value_type point;
+#ifndef BOOST_NO_EXCEPTIONS
+      try
+      {
+#endif
+         point = test_func(row);
+#ifndef BOOST_NO_EXCEPTIONS
+      }
+      catch(const std::underflow_error&)
+      {
+         point = 0;
+      }
+      catch(const std::overflow_error&)
+      {
+         point = std::numeric_limits<value_type>::has_infinity ? 
+            std::numeric_limits<value_type>::infinity()
+            : tools::max_value<value_type>();
+      }
+      catch(const std::exception& e)
+      {
+         std::cerr << "Unexpected exception at entry: " << i << "\n";
+         std::cerr << e.what() << std::endl;
+         print_row(row, std::cerr);
+         BOOST_ERROR("Unexpected exception.");
+         // so we don't get further errors:
+         point = expect_func(row);
+      }
+#endif
+      value_type expected = expect_func(row);
+      value_type err = relative_error(point, expected);
+#ifdef BOOST_INSTRUMENT
+      if(err != 0)
+      {
+         std::cout << row[0] << " " << err;
+         if(std::numeric_limits<value_type>::is_specialized)
+         {
+            std::cout << " (" << err / std::numeric_limits<value_type>::epsilon() << "eps)";
+         }
+         std::cout << std::endl;
+      }
+#endif
+      if(!(boost::math::isfinite)(point) && (boost::math::isfinite)(expected))
+      {
+         std::cerr << "CAUTION: Found non-finite result, when a finite value was expected at entry " << i << "\n";
+         std::cerr << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl;
+         print_row(row, std::cerr);
+         BOOST_ERROR("Unexpected non-finite result");
+      }
+      if(err > 0.5)
+      {
+         std::cerr << "CAUTION: Gross error found at entry " << i << ".\n";
+         std::cerr << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl;
+         print_row(row, std::cerr);
+         BOOST_ERROR("Gross error");
+      }
+      result.add(err);
+      if((result.max)() == err)
+         result.set_worst(i);
+   }
+   return result;
+}
+
+template <class Val, class Exception>
+void test_check_throw(Val v, Exception e)
+{
+   BOOST_CHECK(errno);
+   errno = 0;
+}
+
+template <class Val>
+void test_check_throw(Val val, std::domain_error const* e)
+{
+   BOOST_CHECK(errno == EDOM);
+   errno = 0;
+   if(std::numeric_limits<Val>::has_quiet_NaN)
+   {
+      BOOST_CHECK((boost::math::isnan)(val));
+   }
+}
+
+template <class Val>
+void test_check_throw(Val v, std::overflow_error const* e)
+{
+   BOOST_CHECK(errno == ERANGE);
+   errno = 0;
+   BOOST_CHECK((v >= boost::math::tools::max_value<Val>()) || (v <= -boost::math::tools::max_value<Val>()));
+}
+
+template <class Val>
+void test_check_throw(Val v, boost::math::rounding_error const* e)
+{
+   BOOST_CHECK(errno == ERANGE);
+   errno = 0;
+   if(std::numeric_limits<Val>::is_specialized && std::numeric_limits<Val>::is_integer)
+   {
+      BOOST_CHECK((v == (std::numeric_limits<Val>::max)()) || (v == (std::numeric_limits<Val>::min)()));
+   }
+   else
+   {
+      BOOST_CHECK((v == boost::math::tools::max_value<Val>()) || (v == -boost::math::tools::max_value<Val>()));
+   }
+}
+
+} // namespace tools
+} // namespace math
+} // namespace boost
+
+
+  //
+  // exception-free testing support, ideally we'd only define this in our tests,
+  // but to keep things simple we really need it somewhere that's always included:
+  //
+#ifdef BOOST_NO_EXCEPTIONS
+#  define BOOST_MATH_CHECK_THROW(x, ExceptionType) boost::math::tools::test_check_throw(x, static_cast<ExceptionType const*>(0));
+#else
+#  define BOOST_MATH_CHECK_THROW(x, y) BOOST_CHECK_THROW(x, y)
+#endif
+
+#endif
+
+
diff --git a/src/boost/libs/math/include_private/boost/math/tools/test_data.hpp b/src/boost/libs/math/include_private/boost/math/tools/test_data.hpp
new file mode 100644
index 00000000..c2b37583
--- /dev/null
+++ b/src/boost/libs/math/include_private/boost/math/tools/test_data.hpp
@@ -0,0 +1,967 @@
+//  (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)
+
+#ifndef BOOST_MATH_TOOLS_TEST_DATA_HPP
+#define BOOST_MATH_TOOLS_TEST_DATA_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <boost/math/tools/config.hpp>
+#include <boost/assert.hpp>
+#ifdef BOOST_MSVC
+#  pragma warning(push)
+#  pragma warning(disable: 4127 4701 4512)
+#  pragma warning(disable: 4130) // '==' : logical operation on address of string constant.
+#endif
+#include <boost/algorithm/string/trim.hpp>
+#include <boost/lexical_cast.hpp>
+#ifdef BOOST_MSVC
+#pragma warning(pop)
+#endif
+#include <boost/type_traits/is_floating_point.hpp>
+#include <boost/type_traits/is_convertible.hpp>
+#include <boost/type_traits/integral_constant.hpp>
+#ifndef BOOST_NO_CXX11_HDR_RANDOM
+#include <random>
+namespace random_ns = std;
+#else
+#include <boost/random.hpp>
+namespace random_ns = boost::random;
+#endif
+#include <boost/math/tools/tuple.hpp>
+#include <boost/math/tools/real_cast.hpp>
+
+#include <set>
+#include <vector>
+#include <iostream>
+
+#ifdef BOOST_MSVC
+#  pragma warning(push)
+#  pragma warning(disable: 4130) // '==' : logical operation on address of string constant.
+// Used as a warning with BOOST_ASSERT
+#endif
+
+namespace boost{ namespace math{ namespace tools{
+
+enum parameter_type
+{
+   random_in_range = 0,
+   periodic_in_range = 1,
+   power_series = 2,
+   single_value = 3,
+   plus_minus_value = 4,
+   dummy_param = 0x80
+};
+
+parameter_type operator | (parameter_type a, parameter_type b)
+{
+   return static_cast<parameter_type>((int)a|(int)b);
+}
+parameter_type& operator |= (parameter_type& a, parameter_type b)
+{
+   a = static_cast<parameter_type>(a|b);
+   return a;
+}
+
+//
+// If type == random_in_range then
+// z1 and r2 are the endpoints of the half open range and n1 is the number of points.
+//
+// If type == periodic_in_range then
+// z1 and r2 are the endpoints of the half open range and n1 is the number of points.
+//
+// If type == power_series then
+// n1 and n2 are the endpoints of the exponents (closed range) and z1 is the basis.
+//
+// If type == single_value then z1 contains the single value to add.
+//
+// If type == plus_minus_value then test at +-z1
+//
+// If type & dummy_param then this data is ignored and not stored in the output, it
+// is passed to the generator function however which can do with it as it sees fit.
+//
+template <class T>
+struct parameter_info
+{
+   parameter_type type;
+   T z1, z2;
+   int n1, n2;
+};
+
+template <class T>
+inline parameter_info<T> make_random_param(T start_range, T end_range, int n_points)
+{
+   parameter_info<T> result = { random_in_range, start_range, end_range, n_points, 0 };
+   return result;
+}
+
+template <class T>
+inline parameter_info<T> make_periodic_param(T start_range, T end_range, int n_points)
+{
+   parameter_info<T> result = { periodic_in_range, start_range, end_range, n_points, 0 };
+   return result;
+}
+
+template <class T>
+inline parameter_info<T> make_power_param(T basis, int start_exponent, int end_exponent)
+{
+   parameter_info<T> result = { power_series, basis, 0, start_exponent, end_exponent };
+   return result;
+}
+
+template <class T>
+inline parameter_info<T> make_single_param(T val)
+{
+   parameter_info<T> result = { single_value, val };
+   return result;
+}
+
+template <class T>
+inline parameter_info<T> make_plus_minus_param(T val)
+{
+   parameter_info<T> result = { plus_minus_value, val };
+   return result;
+}
+
+namespace detail{
+
+template <class Seq, class Item, int N>
+inline void unpack_and_append_tuple(Seq&,
+                                    const Item&,
+                                    const boost::integral_constant<int, N>&,
+                                    const boost::false_type&)
+{
+   // termimation condition nothing to do here
+}
+
+template <class Seq, class Item, int N>
+inline void unpack_and_append_tuple(Seq& s,
+                                    const Item& data,
+                                    const boost::integral_constant<int, N>&,
+                                    const boost::true_type&)
+{
+   // extract the N'th element, append, and recurse:
+   typedef typename Seq::value_type value_type;
+   value_type val = boost::math::get<N>(data);
+   s.push_back(val);
+
+   typedef boost::integral_constant<int, N+1> next_value;
+   typedef boost::integral_constant<bool, (boost::math::tuple_size<Item>::value > N+1)> terminate;
+
+   unpack_and_append_tuple(s, data, next_value(), terminate());
+}
+
+template <class Seq, class Item>
+inline void unpack_and_append(Seq& s, const Item& data, const boost::true_type&)
+{
+   s.push_back(data);
+}
+
+template <class Seq, class Item>
+inline void unpack_and_append(Seq& s, const Item& data, const boost::false_type&)
+{
+   // Item had better be a tuple-like type or we've had it!!!!
+   typedef boost::integral_constant<int, 0> next_value;
+   typedef boost::integral_constant<bool, (boost::math::tuple_size<Item>::value > 0)> terminate;
+
+   unpack_and_append_tuple(s, data, next_value(), terminate());
+}
+
+template <class Seq, class Item>
+inline void unpack_and_append(Seq& s, const Item& data)
+{
+   typedef typename Seq::value_type value_type;
+   unpack_and_append(s, data, ::boost::is_convertible<Item, value_type>());
+}
+
+} // detail
+
+template <class T>
+class test_data
+{
+public:
+   typedef std::vector<T> row_type;
+   typedef row_type value_type;
+private:
+   typedef std::set<row_type> container_type;
+public:
+   typedef typename container_type::reference reference;
+   typedef typename container_type::const_reference const_reference;
+   typedef typename container_type::iterator iterator;
+   typedef typename container_type::const_iterator const_iterator;
+   typedef typename container_type::difference_type difference_type;
+   typedef typename container_type::size_type size_type;
+
+   // creation:
+   test_data(){}
+   template <class F>
+   test_data(F func, const parameter_info<T>& arg1)
+   {
+      insert(func, arg1);
+   }
+
+   // insertion:
+   template <class F>
+   test_data& insert(F func, const parameter_info<T>& arg1)
+   {
+      // generate data for single argument functor F
+
+      typedef typename std::set<T>::const_iterator it_type;
+
+      std::set<T> points;
+      create_test_points(points, arg1);
+      it_type a = points.begin();
+      it_type b = points.end();
+      row_type row;
+      while(a != b)
+      {
+         if((arg1.type & dummy_param) == 0)
+            row.push_back(*a);
+#ifndef BOOST_NO_EXCEPTIONS
+         try{
+#endif
+            // domain_error exceptions from func are swallowed
+            // and this data point is ignored:
+            boost::math::tools::detail::unpack_and_append(row, func(*a));
+            m_data.insert(row);
+#ifndef BOOST_NO_EXCEPTIONS
+         }
+         catch(const std::domain_error&){}
+#endif
+         row.clear();
+         ++a;
+      }
+      return *this;
+   }
+
+   template <class F>
+   test_data& insert(F func, const parameter_info<T>& arg1, const parameter_info<T>& arg2)
+   {
+      // generate data for 2-argument functor F
+
+      typedef typename std::set<T>::const_iterator it_type;
+
+      std::set<T> points1, points2;
+      create_test_points(points1, arg1);
+      create_test_points(points2, arg2);
+      it_type a = points1.begin();
+      it_type b = points1.end();
+      row_type row;
+      while(a != b)
+      {
+         it_type c = points2.begin();
+         it_type d = points2.end();
+         while(c != d)
+         {
+            if((arg1.type & dummy_param) == 0)
+               row.push_back(*a);
+            if((arg2.type & dummy_param) == 0)
+               row.push_back(*c);
+#ifndef BOOST_NO_EXCEPTIONS
+            try{
+#endif
+               // domain_error exceptions from func are swallowed
+               // and this data point is ignored:
+               detail::unpack_and_append(row, func(*a, *c));
+               m_data.insert(row);
+#ifndef BOOST_NO_EXCEPTIONS
+            }
+            catch(const std::domain_error&){}
+#endif
+            row.clear();
+            ++c;
+         }
+         ++a;
+      }
+      return *this;
+   }
+
+   template <class F>
+   test_data& insert(F func, const parameter_info<T>& arg1, const parameter_info<T>& arg2, const parameter_info<T>& arg3)
+   {
+      // generate data for 3-argument functor F
+
+      typedef typename std::set<T>::const_iterator it_type;
+
+      std::set<T> points1, points2, points3;
+      create_test_points(points1, arg1);
+      create_test_points(points2, arg2);
+      create_test_points(points3, arg3);
+      it_type a = points1.begin();
+      it_type b = points1.end();
+      row_type row;
+      while(a != b)
+      {
+         it_type c = points2.begin();
+         it_type d = points2.end();
+         while(c != d)
+         {
+            it_type e = points3.begin();
+            it_type f = points3.end();
+            while(e != f)
+            {
+               if((arg1.type & dummy_param) == 0)
+                  row.push_back(*a);
+               if((arg2.type & dummy_param) == 0)
+                  row.push_back(*c);
+               if((arg3.type & dummy_param) == 0)
+                  row.push_back(*e);
+#ifndef BOOST_NO_EXCEPTIONS
+               try{
+#endif
+                  // domain_error exceptions from func are swallowed
+                  // and this data point is ignored:
+                  detail::unpack_and_append(row, func(*a, *c, *e));
+                  m_data.insert(row);
+#ifndef BOOST_NO_EXCEPTIONS
+               }
+               catch(const std::domain_error&){}
+#endif
+               row.clear();
+               ++e;
+            }
+            ++c;
+         }
+         ++a;
+      }
+      return *this;
+   }
+
+   template <class F>
+   test_data& insert(F func, const parameter_info<T>& arg1, const parameter_info<T>& arg2, const parameter_info<T>& arg3, const parameter_info<T>& arg4)
+   {
+      // generate data for 4-argument functor F
+
+      typedef typename std::set<T>::const_iterator it_type;
+
+      std::set<T> points1, points2, points3, points4;
+      create_test_points(points1, arg1);
+      create_test_points(points2, arg2);
+      create_test_points(points3, arg3);
+      create_test_points(points4, arg4);
+      it_type a = points1.begin();
+      it_type b = points1.end();
+      row_type row;
+      while(a != b)
+      {
+         it_type c = points2.begin();
+         it_type d = points2.end();
+         while(c != d)
+         {
+            it_type e = points3.begin();
+            it_type f = points3.end();
+            while(e != f)
+            {
+               it_type g = points4.begin();
+               it_type h = points4.end();
+               while (g != h)
+               {
+                  if ((arg1.type & dummy_param) == 0)
+                     row.push_back(*a);
+                  if ((arg2.type & dummy_param) == 0)
+                     row.push_back(*c);
+                  if ((arg3.type & dummy_param) == 0)
+                     row.push_back(*e);
+                  if ((arg4.type & dummy_param) == 0)
+                     row.push_back(*g);
+#ifndef BOOST_NO_EXCEPTIONS
+                  try {
+#endif
+                     // domain_error exceptions from func are swallowed
+                     // and this data point is ignored:
+                     detail::unpack_and_append(row, func(*a, *c, *e, *g));
+                     m_data.insert(row);
+#ifndef BOOST_NO_EXCEPTIONS
+                  }
+                  catch (const std::domain_error&) {}
+#endif
+                  row.clear();
+                  ++g;
+               }
+               ++e;
+            }
+            ++c;
+         }
+         ++a;
+      }
+      return *this;
+   }
+
+   template <class F>
+   test_data& insert(F func, const parameter_info<T>& arg1, const parameter_info<T>& arg2, const parameter_info<T>& arg3, const parameter_info<T>& arg4, const parameter_info<T>& arg5)
+   {
+      // generate data for 5-argument functor F
+
+      typedef typename std::set<T>::const_iterator it_type;
+
+      std::set<T> points1, points2, points3, points4, points5;
+      create_test_points(points1, arg1);
+      create_test_points(points2, arg2);
+      create_test_points(points3, arg3);
+      create_test_points(points4, arg4);
+      create_test_points(points5, arg5);
+      it_type a = points1.begin();
+      it_type b = points1.end();
+      row_type row;
+      while(a != b)
+      {
+         it_type c = points2.begin();
+         it_type d = points2.end();
+         while(c != d)
+         {
+            it_type e = points3.begin();
+            it_type f = points3.end();
+            while(e != f)
+            {
+               it_type g = points4.begin();
+               it_type h = points4.end();
+               while (g != h)
+               {
+                  it_type i = points5.begin();
+                  it_type j = points5.end();
+                  while (i != j)
+                  {
+                     if ((arg1.type & dummy_param) == 0)
+                        row.push_back(*a);
+                     if ((arg2.type & dummy_param) == 0)
+                        row.push_back(*c);
+                     if ((arg3.type & dummy_param) == 0)
+                        row.push_back(*e);
+                     if ((arg4.type & dummy_param) == 0)
+                        row.push_back(*g);
+                     if ((arg5.type & dummy_param) == 0)
+                        row.push_back(*i);
+#ifndef BOOST_NO_EXCEPTIONS
+                     try {
+#endif
+                        // domain_error exceptions from func are swallowed
+                        // and this data point is ignored:
+                        detail::unpack_and_append(row, func(*a, *c, *e, *g, *i));
+                        m_data.insert(row);
+#ifndef BOOST_NO_EXCEPTIONS
+                     }
+                     catch (const std::domain_error&) {}
+#endif
+                     row.clear();
+                     ++i;
+                  }
+                  ++g;
+               }
+               ++e;
+            }
+            ++c;
+         }
+         ++a;
+      }
+      return *this;
+   }
+
+   void clear(){ m_data.clear(); }
+
+   // access:
+   iterator begin() { return m_data.begin(); }
+   iterator end() { return m_data.end(); }
+   const_iterator begin()const { return m_data.begin(); }
+   const_iterator end()const { return m_data.end(); }
+   bool operator==(const test_data& d)const{ return m_data == d.m_data; }
+   bool operator!=(const test_data& d)const{ return m_data != d.m_data; }
+   void swap(test_data& other){ m_data.swap(other.m_data); }
+   size_type size()const{ return m_data.size(); }
+   size_type max_size()const{ return m_data.max_size(); }
+   bool empty()const{ return m_data.empty(); }
+
+   bool operator < (const test_data& dat)const{ return m_data < dat.m_data; }
+   bool operator <= (const test_data& dat)const{ return m_data <= dat.m_data; }
+   bool operator > (const test_data& dat)const{ return m_data > dat.m_data; }
+   bool operator >= (const test_data& dat)const{ return m_data >= dat.m_data; }
+
+private:
+   void create_test_points(std::set<T>& points, const parameter_info<T>& arg1);
+   std::set<row_type> m_data;
+
+   static float extern_val;
+   static float truncate_to_float(float const * pf);
+   static float truncate_to_float(float c){ return truncate_to_float(&c); }
+};
+
+//
+// This code exists to bemuse the compiler's optimizer and force a
+// truncation to float-precision only:
+//
+template <class T>
+inline float test_data<T>::truncate_to_float(float const * pf)
+{
+   BOOST_MATH_STD_USING
+   int expon;
+   float f = floor(ldexp(frexp(*pf, &expon), 22));
+   f = ldexp(f, expon - 22);
+   return f;
+
+   //extern_val = *pf;
+   //return *pf;
+}
+
+template <class T>
+float test_data<T>::extern_val = 0;
+
+template <class T>
+void test_data<T>::create_test_points(std::set<T>& points, const parameter_info<T>& arg1)
+{
+   BOOST_MATH_STD_USING
+   //
+   // Generate a set of test points as requested, try and generate points
+   // at only float precision: otherwise when testing float versions of functions
+   // there will be a rounding error in our input values which throws off the results
+   // (Garbage in garbage out etc).
+   //
+   switch(arg1.type & 0x7F)
+   {
+   case random_in_range:
+      {
+         BOOST_ASSERT(arg1.z1 < arg1.z2);
+         BOOST_ASSERT(arg1.n1 > 0);
+         typedef float random_type;
+
+         random_ns::mt19937 rnd;
+         random_ns::uniform_real_distribution<random_type> ur_a(real_cast<random_type>(arg1.z1), real_cast<random_type>(arg1.z2));
+
+         for(int i = 0; i < arg1.n1; ++i)
+         {
+            random_type r = ur_a(rnd);
+            points.insert(truncate_to_float(r));
+         }
+     }
+      break;
+   case periodic_in_range:
+      {
+         BOOST_ASSERT(arg1.z1 < arg1.z2);
+         BOOST_ASSERT(arg1.n1 > 0);
+         float interval = real_cast<float>((arg1.z2 - arg1.z1) / arg1.n1);
+         T val = arg1.z1;
+         while(val < arg1.z2)
+         {
+            points.insert(truncate_to_float(real_cast<float>(val)));
+            val += interval;
+         }
+      }
+      break;
+   case power_series:
+      {
+         BOOST_ASSERT(arg1.n1 < arg1.n2);
+
+         typedef float random_type;
+         typedef typename boost::mpl::if_<
+            ::boost::is_floating_point<T>,
+            T, long double>::type power_type;
+
+         random_ns::mt19937 rnd;
+         random_ns::uniform_real_distribution<random_type> ur_a(1.0, 2.0);
+
+         for(int power = arg1.n1; power <= arg1.n2; ++power)
+         {
+            random_type r = ur_a(rnd);
+            power_type p = ldexp(static_cast<power_type>(r), power);
+            points.insert(truncate_to_float(real_cast<float>(arg1.z1 + p)));
+         }
+      }
+      break;
+   case single_value:
+   {
+      points.insert(truncate_to_float(real_cast<float>(arg1.z1)));
+      break;
+   }
+   case plus_minus_value:
+   {
+      points.insert(truncate_to_float(real_cast<float>(arg1.z1)));
+      points.insert(truncate_to_float(-real_cast<float>(arg1.z1)));
+      break;
+   }
+   default:
+      BOOST_ASSERT(0 == "Invalid parameter_info object");
+      // Assert will fail if get here.
+      // Triggers warning 4130) // '==' : logical operation on address of string constant.
+   }
+}
+
+//
+// Prompt a user for information on a parameter range:
+//
+template <class T>
+bool get_user_parameter_info(parameter_info<T>& info, const char* param_name)
+{
+#ifdef BOOST_MSVC
+#  pragma warning(push)
+#  pragma warning(disable: 4127)
+#endif
+   std::string line;
+   do{
+      std::cout << "What kind of distribution do you require for parameter " << param_name << "?\n"
+         "Choices are:\n"
+         "  r     Random values in a half open range\n"
+         "  p     Evenly spaced periodic values in a half open range\n"
+         "  e     Exponential power series at a particular point: a + 2^b for some range of b\n"
+         "[Default=r]";
+
+      std::getline(std::cin, line);
+      boost::algorithm::trim(line);
+
+      if(line == "r")
+      {
+         info.type = random_in_range;
+         break;
+      }
+      else if(line == "p")
+      {
+         info.type = periodic_in_range;
+         break;
+      }
+      else if(line == "e")
+      {
+         info.type = power_series;
+         break;
+      }
+      else if(line == "")
+      {
+         info.type = random_in_range;
+         break;
+      }
+      //
+      // Ooops, not a valid input....
+      //
+      std::cout << "Sorry don't recognise \"" << line << "\" as a valid input\n"
+         "do you want to try again [y/n]?";
+      std::getline(std::cin, line);
+      boost::algorithm::trim(line);
+      if(line == "n")
+         return false;
+      else if(line == "y")
+         continue;
+      std::cout << "Sorry don't recognise that either, giving up...\n\n";
+      return false;
+   }while(true);
+
+   switch(info.type & ~dummy_param)
+   {
+   case random_in_range:
+   case periodic_in_range:
+      // get start and end points of range:
+      do{
+         std::cout << "Data will be in the half open range a <= x < b,\n"
+            "enter value for the start point fo the range [default=0]:";
+         std::getline(std::cin, line);
+         boost::algorithm::trim(line);
+         if(line == "")
+         {
+            info.z1 = 0;
+            break;
+         }
+#ifndef BOOST_NO_EXCEPTIONS
+         try{
+#endif
+            info.z1 = boost::lexical_cast<T>(line);
+            break;
+#ifndef BOOST_NO_EXCEPTIONS
+         }
+         catch(const boost::bad_lexical_cast&)
+         {
+            std::cout << "Sorry, that was not valid input, try again [y/n]?";
+            std::getline(std::cin, line);
+            boost::algorithm::trim(line);
+            if(line == "y")
+               continue;
+            if(line == "n")
+               return false;
+            std::cout << "Sorry don't recognise that either, giving up...\n\n";
+            return false;
+         }
+#endif
+      }while(true);
+      do{
+         std::cout << "Enter value for the end point fo the range [default=1]:";
+         std::getline(std::cin, line);
+         boost::algorithm::trim(line);
+         if(line == "")
+         {
+            info.z2 = 1;
+         }
+         else
+         {
+#ifndef BOOST_NO_EXCEPTIONS
+            try
+            {
+#endif
+               info.z2 = boost::lexical_cast<T>(line);
+#ifndef BOOST_NO_EXCEPTIONS
+            }
+            catch(const boost::bad_lexical_cast&)
+            {
+               std::cout << "Sorry, that was not valid input, try again [y/n]?";
+               std::getline(std::cin, line);
+               boost::algorithm::trim(line);
+               if(line == "y")
+                  continue;
+               if(line == "n")
+                  return false;
+               std::cout << "Sorry don't recognise that either, giving up...\n\n";
+               return false;
+            }
+#endif
+         }
+         if(info.z1 >= info.z2)
+         {
+            std::cout << "The end point of the range was <= the start point\n"
+               "try a different value for the endpoint [y/n]?";
+            std::getline(std::cin, line);
+            boost::algorithm::trim(line);
+            if(line == "y")
+               continue;
+            if(line == "n")
+               return false;
+            std::cout << "Sorry don't recognise that either, giving up...\n\n";
+            return false;
+         }
+         break;
+      }while(true);
+      do{
+         // get the number of points:
+         std::cout << "How many data points do you want?";
+         std::getline(std::cin, line);
+         boost::algorithm::trim(line);
+#ifndef BOOST_NO_EXCEPTIONS
+         try{
+#endif
+            info.n1 = boost::lexical_cast<int>(line);
+            info.n2 = 0;
+            if(info.n1 <= 0)
+            {
+               std::cout << "The number of points should be > 0\n"
+                  "try again [y/n]?";
+               std::getline(std::cin, line);
+               boost::algorithm::trim(line);
+               if(line == "y")
+                  continue;
+               if(line == "n")
+                  return false;
+               std::cout << "Sorry don't recognise that either, giving up...\n\n";
+               return false;
+            }
+            break;
+#ifndef BOOST_NO_EXCEPTIONS
+         }
+         catch(const boost::bad_lexical_cast&)
+         {
+            std::cout << "Sorry, that was not valid input, try again [y/n]?";
+            std::getline(std::cin, line);
+            boost::algorithm::trim(line);
+            if(line == "y")
+               continue;
+            if(line == "n")
+               return false;
+            std::cout << "Sorry don't recognise that either, giving up...\n\n";
+            return false;
+         }
+#endif
+      }while(true);
+      break;
+   case power_series:
+      // get start and end points of range:
+      info.z2 = 0;
+      do{
+         std::cout << "Data will be in the form a + r*2^b\n"
+            "for random value r,\n"
+            "enter value for the point a [default=0]:";
+         std::getline(std::cin, line);
+         boost::algorithm::trim(line);
+         if(line == "")
+         {
+            info.z1 = 0;
+            break;
+         }
+#ifndef BOOST_NO_EXCEPTIONS
+         try{
+#endif
+            info.z1 = boost::lexical_cast<T>(line);
+            break;
+#ifndef BOOST_NO_EXCEPTIONS
+         }
+         catch(const boost::bad_lexical_cast&)
+         {
+            std::cout << "Sorry, that was not valid input, try again [y/n]?";
+            std::getline(std::cin, line);
+            boost::algorithm::trim(line);
+            if(line == "y")
+               continue;
+            if(line == "n")
+               return false;
+            std::cout << "Sorry don't recognise that either, giving up...\n\n";
+            return false;
+         }
+#endif
+      }while(true);
+
+      do{
+         std::cout << "Data will be in the form a + r*2^b\n"
+            "for random value r,\n"
+            "enter value for the starting exponent b:";
+         std::getline(std::cin, line);
+         boost::algorithm::trim(line);
+#ifndef BOOST_NO_EXCEPTIONS
+         try{
+#endif
+            info.n1 = boost::lexical_cast<int>(line);
+            break;
+#ifndef BOOST_NO_EXCEPTIONS
+         }
+         catch(const boost::bad_lexical_cast&)
+         {
+            std::cout << "Sorry, that was not valid input, try again [y/n]?";
+            std::getline(std::cin, line);
+            boost::algorithm::trim(line);
+            if(line == "y")
+               continue;
+            if(line == "n")
+               return false;
+            std::cout << "Sorry don't recognise that either, giving up...\n\n";
+            return false;
+         }
+#endif
+      }while(true);
+
+      do{
+         std::cout << "Data will be in the form a + r*2^b\n"
+            "for random value r,\n"
+            "enter value for the ending exponent b:";
+         std::getline(std::cin, line);
+         boost::algorithm::trim(line);
+#ifndef BOOST_NO_EXCEPTIONS
+         try{
+#endif
+            info.n2 = boost::lexical_cast<int>(line);
+            break;
+#ifndef BOOST_NO_EXCEPTIONS
+         }
+         catch(const boost::bad_lexical_cast&)
+         {
+            std::cout << "Sorry, that was not valid input, try again [y/n]?";
+            std::getline(std::cin, line);
+            boost::algorithm::trim(line);
+            if(line == "y")
+               continue;
+            if(line == "n")
+               return false;
+            std::cout << "Sorry don't recognise that either, giving up...\n\n";
+            return false;
+         }
+#endif
+      }while(true);
+
+      break;
+   default:
+      BOOST_ASSERT(0); // should never get here!!
+   }
+
+   return true;
+#ifdef BOOST_MSVC
+#  pragma warning(pop)
+#endif
+}
+
+template <class charT, class traits, class T>
+inline std::basic_ostream<charT, traits>& write_csv(std::basic_ostream<charT, traits>& os,
+                                             const test_data<T>& data)
+{
+   const charT defarg[] = { ',', ' ', '\0' };
+   return write_csv(os, data, defarg);
+}
+
+template <class charT, class traits, class T>
+std::basic_ostream<charT, traits>& write_csv(std::basic_ostream<charT, traits>& os,
+                                             const test_data<T>& data,
+                                             const charT* separator)
+{
+   typedef typename test_data<T>::const_iterator it_type;
+   typedef typename test_data<T>::value_type value_type;
+   typedef typename value_type::const_iterator value_type_iterator;
+   it_type a, b;
+   a = data.begin();
+   b = data.end();
+   while(a != b)
+   {
+      value_type_iterator x, y;
+      bool sep = false;
+      x = a->begin();
+      y = a->end();
+      while(x != y)
+      {
+         if(sep)
+            os << separator;
+         os << *x;
+         sep = true;
+         ++x;
+      }
+      os << std::endl;
+      ++a;
+   }
+   return os;
+}
+
+template <class T>
+std::ostream& write_code(std::ostream& os,
+                         const test_data<T>& data,
+                         const char* name)
+{
+   typedef typename test_data<T>::const_iterator it_type;
+   typedef typename test_data<T>::value_type value_type;
+   typedef typename value_type::const_iterator value_type_iterator;
+
+   BOOST_ASSERT(os.good());
+
+   it_type a, b;
+   a = data.begin();
+   b = data.end();
+   if(a == b)
+      return os;
+
+   os << "#ifndef SC_\n#  define SC_(x) static_cast<T>(BOOST_JOIN(x, L))\n#endif\n"
+   "   static const boost::array<boost::array<T, "
+   << a->size() << ">, " << data.size() << "> " << name << " = {{\n";
+
+   while(a != b)
+   {
+      if(a != data.begin())
+         os << ", \n";
+
+      value_type_iterator x, y;
+      x = a->begin();
+      y = a->end();
+      os << "      { ";
+      while(x != y)
+      {
+         if(x != a->begin())
+            os << ", ";
+         os << "SC_(" << *x << ")";
+         ++x;
+      }
+      os << " }";
+      ++a;
+   }
+   os << "\n   }};\n//#undef SC_\n\n";
+   return os;
+}
+
+} // namespace tools
+} // namespace math
+} // namespace boost
+
+#ifdef BOOST_MSVC
+#pragma warning(pop)
+#endif
+
+
+#endif // BOOST_MATH_TOOLS_TEST_DATA_HPP
+
+