bgirgis/cpp-thirdparty
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

updated boost on windows

Browse Source
This commit is contained in:
Bassem Girgis
2019-08-13 21:48:48 -05:00
parent 7d77d485fd
commit b40a3bee82
5162 changed files with 473027 additions and 116452 deletions
Download Patch File Download Diff File Expand all files Collapse all files

93
winx64/include/boost/math/tools/atomic.hpp Normal file
View File

@@ -0,0 +1,93 @@
///////////////////////////////////////////////////////////////////////////////
// Copyright 2017 John Maddock
// 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)
#ifndef BOOST_MATH_ATOMIC_DETAIL_HPP
#define BOOST_MATH_ATOMIC_DETAIL_HPP
#include <boost/config.hpp>
#ifdef BOOST_HAS_THREADS
#ifndef BOOST_NO_CXX11_HDR_ATOMIC
# include <atomic>
# define BOOST_MATH_ATOMIC_NS std
namespace boost {
namespace math {
namespace detail {
#if ATOMIC_INT_LOCK_FREE == 2
typedef std::atomic<int> atomic_counter_type;
typedef std::atomic<unsigned> atomic_unsigned_type;
typedef int atomic_integer_type;
typedef unsigned atomic_unsigned_integer_type;
#elif ATOMIC_SHORT_LOCK_FREE == 2
typedef std::atomic<short> atomic_counter_type;
typedef std::atomic<unsigned short> atomic_unsigned_type;
typedef short atomic_integer_type;
typedef unsigned short atomic_unsigned_type;
#elif ATOMIC_LONG_LOCK_FREE == 2
typedef std::atomic<long> atomic_unsigned_integer_type;
typedef std::atomic<unsigned long> atomic_unsigned_type;
typedef unsigned long atomic_unsigned_type;
typedef long atomic_integer_type;
#elif ATOMIC_LLONG_LOCK_FREE == 2
typedef std::atomic<long long> atomic_unsigned_integer_type;
typedef std::atomic<unsigned long long> atomic_unsigned_type;
typedef long long atomic_integer_type;
typedef unsigned long long atomic_unsigned_integer_type;
#else
# define BOOST_MATH_NO_ATOMIC_INT
#endif
}
}}
#else // BOOST_NO_CXX11_HDR_ATOMIC
//
// We need Boost.Atomic, but on any platform that supports auto-linking we do
// not need to link against a separate library:
//
#define BOOST_ATOMIC_NO_LIB
#include <boost/atomic.hpp>
# define BOOST_MATH_ATOMIC_NS boost
namespace boost{ namespace math{ namespace detail{
//
// We need a type to use as an atomic counter:
//
#if BOOST_ATOMIC_INT_LOCK_FREE == 2
typedef boost::atomic<int> atomic_counter_type;
typedef boost::atomic<unsigned> atomic_unsigned_type;
typedef int atomic_integer_type;
typedef unsigned atomic_unsigned_integer_type;
#elif BOOST_ATOMIC_SHORT_LOCK_FREE == 2
typedef boost::atomic<short> atomic_counter_type;
typedef boost::atomic<unsigned short> atomic_unsigned_type;
typedef short atomic_integer_type;
typedef unsigned short atomic_unsigned_integer_type;
#elif BOOST_ATOMIC_LONG_LOCK_FREE == 2
typedef boost::atomic<long> atomic_counter_type;
typedef boost::atomic<unsigned long> atomic_unsigned_type;
typedef long atomic_integer_type;
typedef unsigned long atomic_unsigned_integer_type;
#elif BOOST_ATOMIC_LLONG_LOCK_FREE == 2
typedef boost::atomic<long long> atomic_counter_type;
typedef boost::atomic<unsigned long long> atomic_unsigned_type;
typedef long long atomic_integer_type;
typedef unsigned long long atomic_unsigned_integer_type;
#else
# define BOOST_MATH_NO_ATOMIC_INT
#endif
}}} // namespaces
#endif // BOOST_NO_CXX11_HDR_ATOMIC
#else // BOOST_HAS_THREADS
# define BOOST_MATH_NO_ATOMIC_INT
#endif // BOOST_HAS_THREADS
#endif // BOOST_MATH_ATOMIC_DETAIL_HPP

6
winx64/include/boost/math/tools/big_constant.hpp
View File

@@ -54,9 +54,9 @@ inline T make_big_value(largest_float, const char* s, mpl::false_ const&, mpl::f
}
#endif
template <class T>
inline BOOST_MATH_CONSTEXPR const char* make_big_value(largest_float, const char* s, mpl::false_ const&, mpl::true_ const&) BOOST_MATH_NOEXCEPT(T)
inline BOOST_MATH_CONSTEXPR T make_big_value(largest_float, const char* s, mpl::false_ const&, mpl::true_ const&) BOOST_MATH_NOEXCEPT(T)
{
return s;
return T(s);
}
//
@@ -78,7 +78,7 @@ inline BOOST_MATH_CONSTEXPR const char* make_big_value(largest_float, const char
#define BOOST_MATH_HUGE_CONSTANT(T, D, x)\
boost::math::tools::make_big_value<T>(0.0L, BOOST_STRINGIZE(x), \
mpl::bool_<is_floating_point<T>::value || (boost::math::tools::numeric_traits<T>::is_specialized && boost::math::tools::numeric_traits<T>::max_exponent <= boost::math::tools::numeric_traits<boost::math::tools::largest_float>::max_exponent && boost::math::tools::numeric_traits<T>::digits <= boost::math::tools::numeric_traits<boost::math::tools::largest_float>::digits)>(), \
boost::is_convertible<const char*, T>())
boost::is_constructible<const char*, T>())
}}} // namespaces

97
winx64/include/boost/math/tools/bivariate_statistics.hpp Normal file
View File

@@ -0,0 +1,97 @@
// (C) Copyright Nick Thompson 2018.
// 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_BIVARIATE_STATISTICS_HPP
#define BOOST_MATH_TOOLS_BIVARIATE_STATISTICS_HPP
#include <iterator>
#include <tuple>
#include <boost/assert.hpp>
#include <boost/multiprecision/detail/number_base.hpp>
namespace boost{ namespace math{ namespace tools {
template<class Container>
auto means_and_covariance(Container const & u, Container const & v)
{
using Real = typename Container::value_type;
using std::size;
BOOST_ASSERT_MSG(size(u) == size(v), "The size of each vector must be the same to compute covariance.");
BOOST_ASSERT_MSG(size(u) > 0, "Computing covariance requires at least one sample.");
// See Equation III.9 of "Numerically Stable, Single-Pass, Parallel Statistics Algorithms", Bennet et al.
Real cov = 0;
Real mu_u = u[0];
Real mu_v = v[0];
for(size_t i = 1; i < size(u); ++i)
{
Real u_tmp = (u[i] - mu_u)/(i+1);
Real v_tmp = v[i] - mu_v;
cov += i*u_tmp*v_tmp;
mu_u = mu_u + u_tmp;
mu_v = mu_v + v_tmp/(i+1);
}
return std::make_tuple(mu_u, mu_v, cov/size(u));
}
template<class Container>
auto covariance(Container const & u, Container const & v)
{
auto [mu_u, mu_v, cov] = boost::math::tools::means_and_covariance(u, v);
return cov;
}
template<class Container>
auto correlation_coefficient(Container const & u, Container const & v)
{
using Real = typename Container::value_type;
using std::size;
BOOST_ASSERT_MSG(size(u) == size(v), "The size of each vector must be the same to compute covariance.");
BOOST_ASSERT_MSG(size(u) > 0, "Computing covariance requires at least two samples.");
Real cov = 0;
Real mu_u = u[0];
Real mu_v = v[0];
Real Qu = 0;
Real Qv = 0;
for(size_t i = 1; i < size(u); ++i)
{
Real u_tmp = u[i] - mu_u;
Real v_tmp = v[i] - mu_v;
Qu = Qu + (i*u_tmp*u_tmp)/(i+1);
Qv = Qv + (i*v_tmp*v_tmp)/(i+1);
cov += i*u_tmp*v_tmp/(i+1);
mu_u = mu_u + u_tmp/(i+1);
mu_v = mu_v + v_tmp/(i+1);
}
// If both datasets are constant, then they are perfectly correlated.
if (Qu == 0 && Qv == 0)
{
return Real(1);
}
// If one dataset is constant and the other isn't, then they have no correlation:
if (Qu == 0 || Qv == 0)
{
return Real(0);
}
// Make sure rho in [-1, 1], even in the presence of numerical noise.
Real rho = cov/sqrt(Qu*Qv);
if (rho > 1) {
rho = 1;
}
if (rho < -1) {
rho = -1;
}
return rho;
}
}}}
#endif

57
winx64/include/boost/math/tools/complex.hpp Normal file
View File

@@ -0,0 +1,57 @@
// Copyright John Maddock 2018.
// 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)
//
// Tools for operator on complex as well as scalar types.
//
#include <boost/type_traits/is_complex.hpp>
namespace boost {
namespace math {
namespace tools {
//
// Speicalize this trait for user-defined complex types (ie Boost.Multiprecision):
//
template <class T>
struct is_complex_type : public boost::is_complex<T> {};
//
// Use this trait to typecast integer literals to something
// that will interoperate with T:
//
template <class T, bool = is_complex_type<T>::value>
struct integer_scalar_type
{
typedef int type;
};
template <class T>
struct integer_scalar_type<T, true>
{
typedef typename T::value_type type;
};
template <class T, bool = is_complex_type<T>::value>
struct unsigned_scalar_type
{
typedef unsigned type;
};
template <class T>
struct unsigned_scalar_type<T, true>
{
typedef typename T::value_type type;
};
template <class T, bool = is_complex_type<T>::value>
struct scalar_type
{
typedef T type;
};
template <class T>
struct scalar_type<T, true>
{
typedef typename T::value_type type;
};
} } }

139
winx64/include/boost/math/tools/condition_numbers.hpp Normal file
View File

@@ -0,0 +1,139 @@
// (C) Copyright Nick Thompson 2019.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_MATH_TOOLS_CONDITION_NUMBERS_HPP
#define BOOST_MATH_TOOLS_CONDITION_NUMBERS_HPP
#include <cmath>
#include <boost/math/differentiation/finite_difference.hpp>
namespace boost::math::tools {
template<class Real, bool kahan=true>
class summation_condition_number {
public:
summation_condition_number(Real const x = 0)
{
using std::abs;
m_l1 = abs(x);
m_sum = x;
m_c = 0;
}
void operator+=(Real const & x)
{
using std::abs;
// No need to Kahan the l1 calc; it's well conditioned:
m_l1 += abs(x);
if constexpr(kahan)
{
Real y = x - m_c;
Real t = m_sum + y;
m_c = (t-m_sum) -y;
m_sum = t;
}
else
{
m_sum += x;
}
}
inline void operator-=(Real const & x)
{
this->operator+=(-x);
}
// Is operator*= relevant? Presumably everything gets rescaled,
// (m_sum -> k*m_sum, m_l1->k*m_l1, m_c->k*m_c),
// but is this sensible? More important is it useful?
// In addition, it might change the condition number.
[[nodiscard]] Real operator()() const
{
using std::abs;
if (m_sum == Real(0) && m_l1 != Real(0))
{
return std::numeric_limits<Real>::infinity();
}
return m_l1/abs(m_sum);
}
[[nodiscard]] Real sum() const
{
// Higham, 1993, "The Accuracy of Floating Point Summation":
// "In [17] and [18], Kahan describes a variation of compensated summation in which the final sum is also corrected
// thus s=s+e is appended to the algorithm above)."
return m_sum + m_c;
}
[[nodiscard]] Real l1_norm() const
{
return m_l1;
}
private:
Real m_l1;
Real m_sum;
Real m_c;
};
template<class F, class Real>
auto evaluation_condition_number(F const & f, Real const & x)
{
using std::abs;
using std::isnan;
using std::sqrt;
using boost::math::differentiation::finite_difference_derivative;
Real fx = f(x);
if (isnan(fx))
{
return std::numeric_limits<Real>::quiet_NaN();
}
bool caught_exception = false;
Real fp;
try
{
fp = finite_difference_derivative(f, x);
}
catch(...)
{
caught_exception = true;
}
if (isnan(fp) || caught_exception)
{
// Check if the right derivative exists:
fp = finite_difference_derivative<decltype(f), Real, 1>(f, x);
if (isnan(fp))
{
// Check if a left derivative exists:
const Real eps = (std::numeric_limits<Real>::epsilon)();
Real h = - 2 * sqrt(eps);
h = boost::math::differentiation::detail::make_xph_representable(x, h);
Real yh = f(x + h);
Real y0 = f(x);
Real diff = yh - y0;
fp = diff / h;
if (isnan(fp))
{
return std::numeric_limits<Real>::quiet_NaN();
}
}
}
if (fx == 0)
{
if (x==0 || fp==0)
{
return std::numeric_limits<Real>::quiet_NaN();
}
return std::numeric_limits<Real>::infinity();
}
return abs(x*fp/fx);
}
}
#endif

17
winx64/include/boost/math/tools/config.hpp
View File

@@ -11,7 +11,7 @@
#endif
#include <boost/config.hpp>
#include <boost/predef.h>
#include <boost/predef/architecture/x86.h>
#include <boost/cstdint.hpp> // for boost::uintmax_t
#include <boost/detail/workaround.hpp>
#include <boost/type_traits/is_integral.hpp>
@@ -31,7 +31,7 @@
#if (defined(__CYGWIN__) || defined(__FreeBSD__) || defined(__NetBSD__) \
|| (defined(__hppa) && !defined(__OpenBSD__)) || (defined(__NO_LONG_DOUBLE_MATH) && (DBL_MANT_DIG != LDBL_MANT_DIG))) \
&& !defined(BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS)
//# define BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
# define BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
#endif
#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
//
@@ -209,7 +209,7 @@
// constexpr support, early GCC implementations can't cope so disable
// constexpr for them:
//
#if !defined(__clang) && defined(__GNUC__)
#if !defined(__clang__) && defined(__GNUC__)
#if (__GNUC__ * 100 + __GNUC_MINOR__) < 490
# define BOOST_MATH_DISABLE_CONSTEXPR
#endif
@@ -451,6 +451,17 @@ namespace boost{ namespace math{
# define BOOST_MATH_THREAD_LOCAL
#endif
//
// Can we have constexpr tables?
//
#if (!defined(BOOST_NO_CXX11_HDR_ARRAY) && !defined(BOOST_NO_CXX14_CONSTEXPR)) || BOOST_WORKAROUND(BOOST_MSVC, >= 1910)
#define BOOST_MATH_HAVE_CONSTEXPR_TABLES
#define BOOST_MATH_CONSTEXPR_TABLE_FUNCTION constexpr
#else
#define BOOST_MATH_CONSTEXPR_TABLE_FUNCTION
#endif
#endif // BOOST_MATH_TOOLS_CONFIG_HPP

40
winx64/include/boost/math/tools/detail/is_const_iterable.hpp Normal file
View File

@@ -0,0 +1,40 @@
// (C) Copyright John Maddock 2018.
// 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_IS_CONST_ITERABLE_HPP
#define BOOST_MATH_TOOLS_IS_CONST_ITERABLE_HPP
#if !defined(BOOST_NO_CXX14_VARIABLE_TEMPLATES) && !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_SFINAE_EXPR)
#define BOOST_MATH_HAS_IS_CONST_ITERABLE
#include <boost/type_traits/is_detected.hpp>
#include <utility>
namespace boost {
namespace math {
namespace tools {
namespace detail {
template<class T>
using begin_t = decltype(std::declval<const T&>().begin());
template<class T>
using end_t = decltype(std::declval<const T&>().end());
template<class T>
using const_iterator_t = typename T::const_iterator;
template <class T>
struct is_const_iterable
: public boost::integral_constant<bool,
boost::is_detected<begin_t, T>::value
&& boost::is_detected<end_t, T>::value
&& boost::is_detected<const_iterator_t, T>::value
> {};
} } } }
#endif
#endif // BOOST_MATH_TOOLS_IS_CONST_ITERABLE_HPP

51
winx64/include/boost/math/tools/fraction.hpp
View File

@@ -15,6 +15,7 @@
#include <boost/type_traits/integral_constant.hpp>
#include <boost/mpl/if.hpp>
#include <boost/math/tools/precision.hpp>
#include <boost/math/tools/complex.hpp>
namespace boost{ namespace math{ namespace tools{
@@ -68,6 +69,22 @@ namespace detail
{
};
template <class T, bool = is_complex_type<T>::value>
struct tiny_value
{
static T get() {
return tools::min_value<T>();
}
};
template <class T>
struct tiny_value<T, true>
{
typedef typename T::value_type value_type;
static T get() {
return tools::min_value<value_type>();
}
};
} // namespace detail
//
@@ -93,14 +110,19 @@ inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(G
typedef detail::fraction_traits<Gen> traits;
typedef typename traits::result_type result_type;
typedef typename traits::value_type value_type;
typedef typename integer_scalar_type<result_type>::type integer_type;
typedef typename scalar_type<result_type>::type scalar_type;
result_type tiny = tools::min_value<result_type>();
integer_type const zero(0), one(1);
result_type tiny = detail::tiny_value<result_type>::get();
scalar_type terminator = abs(factor);
value_type v = g();
result_type f, C, D, delta;
f = traits::b(v);
if(f == 0)
if(f == zero)
f = tiny;
C = f;
D = 0;
@@ -110,15 +132,15 @@ inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(G
do{
v = g();
D = traits::b(v) + traits::a(v) * D;
if(D == 0)
if(D == result_type(0))
D = tiny;
C = traits::b(v) + traits::a(v) / C;
if(C == 0)
if(C == zero)
C = tiny;
D = 1/D;
D = one/D;
delta = C*D;
f = f * delta;
}while((fabs(delta - 1) > factor) && --counter);
}while((abs(delta - one) > terminator) && --counter);
max_terms = max_terms - counter;
@@ -183,15 +205,20 @@ inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(G
typedef detail::fraction_traits<Gen> traits;
typedef typename traits::result_type result_type;
typedef typename traits::value_type value_type;
typedef typename integer_scalar_type<result_type>::type integer_type;
typedef typename scalar_type<result_type>::type scalar_type;
result_type tiny = tools::min_value<result_type>();
integer_type const zero(0), one(1);
result_type tiny = detail::tiny_value<result_type>::get();
scalar_type terminator = abs(factor);
value_type v = g();
result_type f, C, D, delta, a0;
f = traits::b(v);
a0 = traits::a(v);
if(f == 0)
if(f == zero)
f = tiny;
C = f;
D = 0;
@@ -201,15 +228,15 @@ inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(G
do{
v = g();
D = traits::b(v) + traits::a(v) * D;
if(D == 0)
if(D == zero)
D = tiny;
C = traits::b(v) + traits::a(v) / C;
if(C == 0)
if(C == zero)
C = tiny;
D = 1/D;
D = one/D;
delta = C*D;
f = f * delta;
}while((fabs(delta - 1) > factor) && --counter);
}while((abs(delta - one) > terminator) && --counter);
max_terms = max_terms - counter;

640
winx64/include/boost/math/tools/norms.hpp Normal file
View File

@@ -0,0 +1,640 @@
// (C) Copyright Nick Thompson 2018.
// 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_NORMS_HPP
#define BOOST_MATH_TOOLS_NORMS_HPP
#include <algorithm>
#include <iterator>
#include <boost/type_traits/is_complex.hpp>
#include <boost/assert.hpp>
#include <boost/multiprecision/detail/number_base.hpp>
namespace boost::math::tools {
// Mallat, "A Wavelet Tour of Signal Processing", equation 2.60:
template<class ForwardIterator>
auto total_variation(ForwardIterator first, ForwardIterator last)
{
using T = typename std::iterator_traits<ForwardIterator>::value_type;
using std::abs;
BOOST_ASSERT_MSG(first != last && std::next(first) != last, "At least two samples are required to compute the total variation.");
auto it = first;
if constexpr (std::is_unsigned<T>::value)
{
T tmp = *it;
double tv = 0;
while (++it != last)
{
if (*it > tmp)
{
tv += *it - tmp;
}
else
{
tv += tmp - *it;
}
tmp = *it;
}
return tv;
}
else if constexpr (std::is_integral<T>::value)
{
double tv = 0;
double tmp = *it;
while(++it != last)
{
double tmp2 = *it;
tv += abs(tmp2 - tmp);
tmp = *it;
}
return tv;
}
else
{
T tmp = *it;
T tv = 0;
while (++it != last)
{
tv += abs(*it - tmp);
tmp = *it;
}
return tv;
}
}
template<class Container>
inline auto total_variation(Container const & v)
{
return total_variation(v.cbegin(), v.cend());
}
template<class ForwardIterator>
auto sup_norm(ForwardIterator first, ForwardIterator last)
{
BOOST_ASSERT_MSG(first != last, "At least one value is required to compute the sup norm.");
using T = typename std::iterator_traits<ForwardIterator>::value_type;
using std::abs;
if constexpr (boost::is_complex<T>::value ||
boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex)
{
auto it = std::max_element(first, last, [](T a, T b) { return abs(b) > abs(a); });
return abs(*it);
}
else if constexpr (std::is_unsigned<T>::value)
{
return *std::max_element(first, last);
}
else
{
auto pair = std::minmax_element(first, last);
if (abs(*pair.first) > abs(*pair.second))
{
return abs(*pair.first);
}
else
{
return abs(*pair.second);
}
}
}
template<class Container>
inline auto sup_norm(Container const & v)
{
return sup_norm(v.cbegin(), v.cend());
}
template<class ForwardIterator>
auto l1_norm(ForwardIterator first, ForwardIterator last)
{
using T = typename std::iterator_traits<ForwardIterator>::value_type;
using std::abs;
if constexpr (std::is_unsigned<T>::value)
{
double l1 = 0;
for (auto it = first; it != last; ++it)
{
l1 += *it;
}
return l1;
}
else if constexpr (std::is_integral<T>::value)
{
double l1 = 0;
for (auto it = first; it != last; ++it)
{
double tmp = *it;
l1 += abs(tmp);
}
return l1;
}
else
{
decltype(abs(*first)) l1 = 0;
for (auto it = first; it != last; ++it)
{
l1 += abs(*it);
}
return l1;
}
}
template<class Container>
inline auto l1_norm(Container const & v)
{
return l1_norm(v.cbegin(), v.cend());
}
template<class ForwardIterator>
auto l2_norm(ForwardIterator first, ForwardIterator last)
{
using T = typename std::iterator_traits<ForwardIterator>::value_type;
using std::abs;
using std::norm;
using std::sqrt;
using std::is_floating_point;
if constexpr (boost::is_complex<T>::value ||
boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex)
{
typedef typename T::value_type Real;
Real l2 = 0;
for (auto it = first; it != last; ++it)
{
l2 += norm(*it);
}
Real result = sqrt(l2);
if (!isfinite(result))
{
Real a = sup_norm(first, last);
l2 = 0;
for (auto it = first; it != last; ++it)
{
l2 += norm(*it/a);
}
return a*sqrt(l2);
}
return result;
}
else if constexpr (is_floating_point<T>::value ||
boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_floating_point)
{
T l2 = 0;
for (auto it = first; it != last; ++it)
{
l2 += (*it)*(*it);
}
T result = sqrt(l2);
// Higham, Accuracy and Stability of Numerical Algorithms,
// Problem 27.5 presents a different algorithm to deal with overflow.
// The algorithm used here takes 3 passes *if* there is overflow.
// Higham's algorithm is 1 pass, but more requires operations than the no oveflow case.
// I'm operating under the assumption that overflow is rare since the dynamic range of floating point numbers is huge.
if (!isfinite(result))
{
T a = sup_norm(first, last);
l2 = 0;
for (auto it = first; it != last; ++it)
{
T tmp = *it/a;
l2 += tmp*tmp;
}
return a*sqrt(l2);
}
return result;
}
else
{
double l2 = 0;
for (auto it = first; it != last; ++it)
{
double tmp = *it;
l2 += tmp*tmp;
}
return sqrt(l2);
}
}
template<class Container>
inline auto l2_norm(Container const & v)
{
return l2_norm(v.cbegin(), v.cend());
}
template<class ForwardIterator>
size_t l0_pseudo_norm(ForwardIterator first, ForwardIterator last)
{
using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type;
size_t count = 0;
for (auto it = first; it != last; ++it)
{
if (*it != RealOrComplex(0))
{
++count;
}
}
return count;
}
template<class Container>
inline size_t l0_pseudo_norm(Container const & v)
{
return l0_pseudo_norm(v.cbegin(), v.cend());
}
template<class ForwardIterator>
size_t hamming_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2)
{
size_t count = 0;
auto it1 = first1;
auto it2 = first2;
while (it1 != last1)
{
if (*it1++ != *it2++)
{
++count;
}
}
return count;
}
template<class Container>
inline size_t hamming_distance(Container const & v, Container const & w)
{
return hamming_distance(v.cbegin(), v.cend(), w.cbegin());
}
template<class ForwardIterator>
auto lp_norm(ForwardIterator first, ForwardIterator last, unsigned p)
{
using std::abs;
using std::pow;
using std::is_floating_point;
using std::isfinite;
using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type;
if constexpr (boost::is_complex<RealOrComplex>::value ||
boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_complex)
{
using std::norm;
using Real = typename RealOrComplex::value_type;
Real lp = 0;
for (auto it = first; it != last; ++it)
{
lp += pow(abs(*it), p);
}
auto result = pow(lp, Real(1)/Real(p));
if (!isfinite(result))
{
auto a = boost::math::tools::sup_norm(first, last);
Real lp = 0;
for (auto it = first; it != last; ++it)
{
lp += pow(abs(*it)/a, p);
}
result = a*pow(lp, Real(1)/Real(p));
}
return result;
}
else if constexpr (is_floating_point<RealOrComplex>::value ||
boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_floating_point)
{
BOOST_ASSERT_MSG(p >= 0, "For p < 0, the lp norm is not a norm");
RealOrComplex lp = 0;
for (auto it = first; it != last; ++it)
{
lp += pow(abs(*it), p);
}
RealOrComplex result = pow(lp, RealOrComplex(1)/RealOrComplex(p));
if (!isfinite(result))
{
RealOrComplex a = boost::math::tools::sup_norm(first, last);
lp = 0;
for (auto it = first; it != last; ++it)
{
lp += pow(abs(*it)/a, p);
}
result = a*pow(lp, RealOrComplex(1)/RealOrComplex(p));
}
return result;
}
else
{
double lp = 0;
for (auto it = first; it != last; ++it)
{
double tmp = *it;
lp += pow(abs(tmp), p);
}
double result = pow(lp, 1.0/double(p));
if (!isfinite(result))
{
double a = boost::math::tools::sup_norm(first, last);
lp = 0;
for (auto it = first; it != last; ++it)
{
double tmp = *it;
lp += pow(abs(tmp)/a, p);
}
result = a*pow(lp, double(1)/double(p));
}
return result;
}
}
template<class Container>
inline auto lp_norm(Container const & v, unsigned p)
{
return lp_norm(v.cbegin(), v.cend(), p);
}
template<class ForwardIterator>
auto lp_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2, unsigned p)
{
using std::pow;
using std::abs;
using std::is_floating_point;
using std::isfinite;
using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type;
auto it1 = first1;
auto it2 = first2;
if constexpr (boost::is_complex<RealOrComplex>::value ||
boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_complex)
{
using Real = typename RealOrComplex::value_type;
using std::norm;
Real dist = 0;
while(it1 != last1)
{
auto tmp = *it1++ - *it2++;
dist += pow(abs(tmp), p);
}
return pow(dist, Real(1)/Real(p));
}
else if constexpr (is_floating_point<RealOrComplex>::value ||
boost::multiprecision::number_category<RealOrComplex>::value == boost::multiprecision::number_kind_floating_point)
{
RealOrComplex dist = 0;
while(it1 != last1)
{
auto tmp = *it1++ - *it2++;
dist += pow(abs(tmp), p);
}
return pow(dist, RealOrComplex(1)/RealOrComplex(p));
}
else
{
double dist = 0;
while(it1 != last1)
{
double tmp1 = *it1++;
double tmp2 = *it2++;
// Naively you'd expect the integer subtraction to be faster,
// but this can overflow or wraparound:
//double tmp = *it1++ - *it2++;
dist += pow(abs(tmp1 - tmp2), p);
}
return pow(dist, 1.0/double(p));
}
}
template<class Container>
inline auto lp_distance(Container const & v, Container const & w, unsigned p)
{
return lp_distance(v.cbegin(), v.cend(), w.cbegin(), p);
}
template<class ForwardIterator>
auto l1_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2)
{
using std::abs;
using std::is_floating_point;
using std::isfinite;
using T = typename std::iterator_traits<ForwardIterator>::value_type;
auto it1 = first1;
auto it2 = first2;
if constexpr (boost::is_complex<T>::value ||
boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex)
{
using Real = typename T::value_type;
Real sum = 0;
while (it1 != last1) {
sum += abs(*it1++ - *it2++);
}
return sum;
}
else if constexpr (is_floating_point<T>::value ||
boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_floating_point)
{
T sum = 0;
while (it1 != last1)
{
sum += abs(*it1++ - *it2++);
}
return sum;
}
else if constexpr (std::is_unsigned<T>::value)
{
double sum = 0;
while(it1 != last1)
{
T x1 = *it1++;
T x2 = *it2++;
if (x1 > x2)
{
sum += (x1 - x2);
}
else
{
sum += (x2 - x1);
}
}
return sum;
}
else if constexpr (std::is_integral<T>::value)
{
double sum = 0;
while(it1 != last1)
{
double x1 = *it1++;
double x2 = *it2++;
sum += abs(x1-x2);
}
return sum;
}
else
{
BOOST_ASSERT_MSG(false, "Could not recognize type.");
}
}
template<class Container>
auto l1_distance(Container const & v, Container const & w)
{
using std::size;
BOOST_ASSERT_MSG(size(v) == size(w),
"L1 distance requires both containers to have the same number of elements");
return l1_distance(v.cbegin(), v.cend(), w.begin());
}
template<class ForwardIterator>
auto l2_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2)
{
using std::abs;
using std::norm;
using std::sqrt;
using std::is_floating_point;
using std::isfinite;
using T = typename std::iterator_traits<ForwardIterator>::value_type;
auto it1 = first1;
auto it2 = first2;
if constexpr (boost::is_complex<T>::value ||
boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex)
{
using Real = typename T::value_type;
Real sum = 0;
while (it1 != last1) {
sum += norm(*it1++ - *it2++);
}
return sqrt(sum);
}
else if constexpr (is_floating_point<T>::value ||
boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_floating_point)
{
T sum = 0;
while (it1 != last1)
{
T tmp = *it1++ - *it2++;
sum += tmp*tmp;
}
return sqrt(sum);
}
else if constexpr (std::is_unsigned<T>::value)
{
double sum = 0;
while(it1 != last1)
{
T x1 = *it1++;
T x2 = *it2++;
if (x1 > x2)
{
double tmp = x1-x2;
sum += tmp*tmp;
}
else
{
double tmp = x2 - x1;
sum += tmp*tmp;
}
}
return sqrt(sum);
}
else
{
double sum = 0;
while(it1 != last1)
{
double x1 = *it1++;
double x2 = *it2++;
double tmp = x1-x2;
sum += tmp*tmp;
}
return sqrt(sum);
}
}
template<class Container>
auto l2_distance(Container const & v, Container const & w)
{
using std::size;
BOOST_ASSERT_MSG(size(v) == size(w),
"L2 distance requires both containers to have the same number of elements");
return l2_distance(v.cbegin(), v.cend(), w.begin());
}
template<class ForwardIterator>
auto sup_distance(ForwardIterator first1, ForwardIterator last1, ForwardIterator first2)
{
using std::abs;
using std::norm;
using std::sqrt;
using std::is_floating_point;
using std::isfinite;
using T = typename std::iterator_traits<ForwardIterator>::value_type;
auto it1 = first1;
auto it2 = first2;
if constexpr (boost::is_complex<T>::value ||
boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex)
{
using Real = typename T::value_type;
Real sup_sq = 0;
while (it1 != last1) {
Real tmp = norm(*it1++ - *it2++);
if (tmp > sup_sq) {
sup_sq = tmp;
}
}
return sqrt(sup_sq);
}
else if constexpr (is_floating_point<T>::value ||
boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_floating_point)
{
T sup = 0;
while (it1 != last1)
{
T tmp = *it1++ - *it2++;
if (sup < abs(tmp))
{
sup = abs(tmp);
}
}
return sup;
}
else // integral values:
{
double sup = 0;
while(it1 != last1)
{
T x1 = *it1++;
T x2 = *it2++;
double tmp;
if (x1 > x2)
{
tmp = x1-x2;
}
else
{
tmp = x2 - x1;
}
if (sup < tmp) {
sup = tmp;
}
}
return sup;
}
}
template<class Container>
auto sup_distance(Container const & v, Container const & w)
{
using std::size;
BOOST_ASSERT_MSG(size(v) == size(w),
"sup distance requires both containers to have the same number of elements");
return sup_distance(v.cbegin(), v.cend(), w.begin());
}
}
#endif

12
winx64/include/boost/math/tools/numerical_differentiation.hpp Normal file
View File

@@ -0,0 +1,12 @@
// (C) Copyright Nick Thompson 2018.
// 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_NUMERICAL_DIFFERENTIATION_HPP
#define BOOST_MATH_TOOLS_NUMERICAL_DIFFERENTIATION_HPP
#include <boost/math/differentiation/finite_difference.hpp>
#include <boost/config/header_deprecated.hpp>
BOOST_HEADER_DEPRECATED("<boost/math/differentiation/finite_difference.hpp>");
#endif

365
winx64/include/boost/math/tools/polynomial.hpp
View File

@@ -15,14 +15,16 @@
#include <boost/assert.hpp>
#include <boost/config.hpp>
#include <boost/config/suffix.hpp>
#include <boost/function.hpp>
#ifdef BOOST_NO_CXX11_LAMBDAS
#include <boost/lambda/lambda.hpp>
#endif
#include <boost/math/tools/rational.hpp>
#include <boost/math/tools/real_cast.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/special_functions/binomial.hpp>
#include <boost/operators.hpp>
#include <boost/core/enable_if.hpp>
#include <boost/type_traits/is_convertible.hpp>
#include <boost/math/tools/detail/is_const_iterable.hpp>
#include <vector>
#include <ostream>
@@ -197,7 +199,7 @@ division(polynomial<T> u, const polynomial<T>& v)
BOOST_ASSERT(u);
typedef typename polynomial<T>::size_type N;
N const m = u.size() - 1, n = v.size() - 1;
N k = m - n;
polynomial<T> q;
@@ -213,13 +215,35 @@ division(polynomial<T> u, const polynomial<T>& v)
return std::make_pair(q, u);
}
template <class T>
struct identity
//
// These structures are the same as the void specializations of the functors of the same name
// in the std lib from C++14 onwards:
//
struct negate
{
T operator()(T const &x) const
{
return x;
}
template <class T>
T operator()(T const &x) const
{
return -x;
}
};
struct plus
{
template <class T, class U>
T operator()(T const &x, U const& y) const
{
return x + y;
}
};
struct minus
{
template <class T, class U>
T operator()(T const &x, U const& y) const
{
return x - y;
}
};
} // namespace detail
@@ -255,12 +279,7 @@ quotient_remainder(const polynomial<T>& dividend, const polynomial<T>& divisor)
template <class T>
class polynomial :
equality_comparable< polynomial<T>,
dividable< polynomial<T>,
dividable2< polynomial<T>, T,
modable< polynomial<T>,
modable2< polynomial<T>, T > > > > >
class polynomial
{
public:
// typedefs:
@@ -284,13 +303,26 @@ public:
normalize();
}
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
polynomial(std::vector<T>&& p) : m_data(std::move(p))
{
normalize();
}
#endif
template <class U>
explicit polynomial(const U& point)
explicit polynomial(const U& point, typename boost::enable_if<boost::is_convertible<U, T> >::type* = 0)
{
if (point != U(0))
m_data.push_back(point);
}
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
// move:
polynomial(polynomial&& p) BOOST_NOEXCEPT
: m_data(std::move(p.m_data)) { }
#endif
// copy:
polynomial(const polynomial& p)
: m_data(p.m_data) { }
@@ -298,17 +330,24 @@ public:
template <class U>
polynomial(const polynomial<U>& p)
{
m_data.resize(p.size());
for(unsigned i = 0; i < p.size(); ++i)
{
m_data.push_back(boost::math::tools::real_cast<T>(p[i]));
m_data[i] = boost::math::tools::real_cast<T>(p[i]);
}
}
#ifdef BOOST_MATH_HAS_IS_CONST_ITERABLE
template <class Range>
explicit polynomial(const Range& r, typename boost::enable_if<boost::math::tools::detail::is_const_iterable<Range> >::type* = 0)
: polynomial(r.begin(), r.end())
{
}
#endif
#if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST) && !BOOST_WORKAROUND(BOOST_GCC_VERSION, < 40500)
polynomial(std::initializer_list<T> l) : polynomial(std::begin(l), std::end(l))
{
}
polynomial&
operator=(std::initializer_list<T> l)
{
@@ -320,26 +359,32 @@ public:
// access:
size_type size()const { return m_data.size(); }
size_type degree()const
size_type size() const { return m_data.size(); }
size_type degree() const
{
if (size() == 0)
throw std::logic_error("degree() is undefined for the zero polynomial.");
return m_data.size() - 1;
}
}
value_type& operator[](size_type i)
{
return m_data[i];
}
const value_type& operator[](size_type i)const
const value_type& operator[](size_type i) const
{
return m_data[i];
}
T evaluate(T z)const
T evaluate(T z) const
{
return m_data.size() > 0 ? boost::math::tools::evaluate_polynomial(&m_data[0], z, m_data.size()) : 0;
return this->operator()(z);
}
std::vector<T> chebyshev()const
T operator()(T z) const
{
return m_data.size() > 0 ? boost::math::tools::evaluate_polynomial(&m_data[0], z, m_data.size()) : T(0);
}
std::vector<T> chebyshev() const
{
return polynomial_to_chebyshev(m_data);
}
@@ -354,9 +399,48 @@ public:
return m_data;
}
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
polynomial<T> prime() const
{
if (m_data.size() == 0)
{
return polynomial<T>({});
}
std::vector<T> p_data(m_data.size() - 1);
for (size_t i = 0; i < p_data.size(); ++i) {
p_data[i] = m_data[i+1]*static_cast<T>(i+1);
}
return polynomial<T>(std::move(p_data));
}
polynomial<T> integrate() const
{
std::vector<T> i_data(m_data.size() + 1);
// Choose integration constant such that P(0) = 0.
i_data[0] = T(0);
for (size_t i = 1; i < i_data.size(); ++i)
{
i_data[i] = m_data[i-1]/static_cast<T>(i);
}
return polynomial<T>(std::move(i_data));
}
// operators:
polynomial& operator =(polynomial&& p) BOOST_NOEXCEPT
{
m_data = std::move(p.m_data);
return *this;
}
#endif
polynomial& operator =(const polynomial& p)
{
m_data = p.m_data;
return *this;
}
template <class U>
polynomial& operator +=(const U& value)
typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator +=(const U& value)
{
addition(value);
normalize();
@@ -364,7 +448,7 @@ public:
}
template <class U>
polynomial& operator -=(const U& value)
typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator -=(const U& value)
{
subtraction(value);
normalize();
@@ -372,7 +456,7 @@ public:
}
template <class U>
polynomial& operator *=(const U& value)
typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator *=(const U& value)
{
multiplication(value);
normalize();
@@ -380,7 +464,7 @@ public:
}
template <class U>
polynomial& operator /=(const U& value)
typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator /=(const U& value)
{
division(value);
normalize();
@@ -388,7 +472,7 @@ public:
}
template <class U>
polynomial& operator %=(const U& /*value*/)
typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator %=(const U& /*value*/)
{
// We can always divide by a scalar, so there is no remainder:
this->set_zero();
@@ -410,20 +494,25 @@ public:
normalize();
return *this;
}
template <typename U, typename V>
void multiply(const polynomial<U>& a, const polynomial<V>& b) {
if (!a || !b)
{
this->set_zero();
return;
}
std::vector<T> prod(a.size() + b.size() - 1, T(0));
for (unsigned i = 0; i < a.size(); ++i)
for (unsigned j = 0; j < b.size(); ++j)
prod[i+j] += a.m_data[i] * b.m_data[j];
m_data.swap(prod);
}
template <class U>
polynomial& operator *=(const polynomial<U>& value)
{
// TODO: FIXME: use O(N log(N)) algorithm!!!
if (!value)
{
this->set_zero();
return *this;
}
std::vector<T> prod(size() + value.size() - 1, T(0));
for (size_type i = 0; i < value.size(); ++i)
for (size_type j = 0; j < size(); ++j)
prod[i+j] += m_data[j] * value[i];
m_data.swap(prod);
this->multiply(*this, value);
return *this;
}
@@ -456,13 +545,13 @@ public:
normalize();
return *this;
}
// Convenient and efficient query for zero.
bool is_zero() const
{
return m_data.empty();
}
// Conversion to bool.
#ifdef BOOST_NO_CXX11_EXPLICIT_CONVERSION_OPERATORS
typedef bool (polynomial::*unmentionable_type)() const;
@@ -483,74 +572,84 @@ public:
{
m_data.clear();
}
/** Remove zero coefficients 'from the top', that is for which there are no
* non-zero coefficients of higher degree. */
void normalize()
{
#ifndef BOOST_NO_CXX11_LAMBDAS
m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), [](const T& x)->bool { return x != T(0); }).base(), m_data.end());
#else
using namespace boost::lambda;
m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), _1 != T(0)).base(), m_data.end());
#endif
}
private:
template <class U, class R1, class R2>
polynomial& addition(const U& value, R1 sign, R2 op)
template <class U, class R>
polynomial& addition(const U& value, R op)
{
if(m_data.size() == 0)
m_data.push_back(sign(value));
else
m_data[0] = op(m_data[0], value);
m_data.resize(1, 0);
m_data[0] = op(m_data[0], value);
return *this;
}
template <class U>
polynomial& addition(const U& value)
{
return addition(value, detail::identity<U>(), std::plus<U>());
return addition(value, detail::plus());
}
template <class U>
polynomial& subtraction(const U& value)
{
return addition(value, std::negate<U>(), std::minus<U>());
return addition(value, detail::minus());
}
template <class U, class R1, class R2>
polynomial& addition(const polynomial<U>& value, R1 sign, R2 op)
template <class U, class R>
polynomial& addition(const polynomial<U>& value, R op)
{
size_type s1 = (std::min)(m_data.size(), value.size());
for(size_type i = 0; i < s1; ++i)
if (m_data.size() < value.size())
m_data.resize(value.size(), 0);
for(size_type i = 0; i < value.size(); ++i)
m_data[i] = op(m_data[i], value[i]);
for(size_type i = s1; i < value.size(); ++i)
m_data.push_back(sign(value[i]));
return *this;
}
template <class U>
polynomial& addition(const polynomial<U>& value)
{
return addition(value, detail::identity<U>(), std::plus<U>());
return addition(value, detail::plus());
}
template <class U>
polynomial& subtraction(const polynomial<U>& value)
{
return addition(value, std::negate<U>(), std::minus<U>());
return addition(value, detail::minus());
}
template <class U>
polynomial& multiplication(const U& value)
{
#ifndef BOOST_NO_CXX11_LAMBDAS
std::transform(m_data.begin(), m_data.end(), m_data.begin(), [&](const T& x)->T { return x * value; });
#else
using namespace boost::lambda;
std::transform(m_data.begin(), m_data.end(), m_data.begin(), ret<T>(_1 * value));
#endif
return *this;
}
template <class U>
polynomial& division(const U& value)
{
#ifndef BOOST_NO_CXX11_LAMBDAS
std::transform(m_data.begin(), m_data.end(), m_data.begin(), [&](const T& x)->T { return x / value; });
#else
using namespace boost::lambda;
std::transform(m_data.begin(), m_data.end(), m_data.begin(), ret<T>(_1 / value));
#endif
return *this;
}
@@ -565,6 +664,26 @@ inline polynomial<T> operator + (const polynomial<T>& a, const polynomial<T>& b)
result += b;
return result;
}
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
template <class T>
inline polynomial<T> operator + (polynomial<T>&& a, const polynomial<T>& b)
{
a += b;
return a;
}
template <class T>
inline polynomial<T> operator + (const polynomial<T>& a, polynomial<T>&& b)
{
b += a;
return b;
}
template <class T>
inline polynomial<T> operator + (polynomial<T>&& a, polynomial<T>&& b)
{
a += b;
return a;
}
#endif
template <class T>
inline polynomial<T> operator - (const polynomial<T>& a, const polynomial<T>& b)
@@ -573,61 +692,101 @@ inline polynomial<T> operator - (const polynomial<T>& a, const polynomial<T>& b)
result -= b;
return result;
}
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
template <class T>
inline polynomial<T> operator - (polynomial<T>&& a, const polynomial<T>& b)
{
a -= b;
return a;
}
template <class T>
inline polynomial<T> operator - (const polynomial<T>& a, polynomial<T>&& b)
{
b -= a;
return -b;
}
template <class T>
inline polynomial<T> operator - (polynomial<T>&& a, polynomial<T>&& b)
{
a -= b;
return a;
}
#endif
template <class T>
inline polynomial<T> operator * (const polynomial<T>& a, const polynomial<T>& b)
{
polynomial<T> result(a);
result *= b;
polynomial<T> result;
result.multiply(a, b);
return result;
}
template <class T>
inline polynomial<T> operator / (const polynomial<T>& a, const polynomial<T>& b)
{
return quotient_remainder(a, b).first;
}
template <class T>
inline polynomial<T> operator % (const polynomial<T>& a, const polynomial<T>& b)
{
return quotient_remainder(a, b).second;
}
template <class T, class U>
inline polynomial<T> operator + (const polynomial<T>& a, const U& b)
inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator + (polynomial<T> a, const U& b)
{
polynomial<T> result(a);
result += b;
return result;
a += b;
return a;
}
template <class T, class U>
inline polynomial<T> operator - (const polynomial<T>& a, const U& b)
inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator - (polynomial<T> a, const U& b)
{
polynomial<T> result(a);
result -= b;
return result;
a -= b;
return a;
}
template <class T, class U>
inline polynomial<T> operator * (const polynomial<T>& a, const U& b)
inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator * (polynomial<T> a, const U& b)
{
polynomial<T> result(a);
result *= b;
return result;
a *= b;
return a;
}
template <class T, class U>
inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator / (polynomial<T> a, const U& b)
{
a /= b;
return a;
}
template <class T, class U>
inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator % (const polynomial<T>&, const U&)
{
// Since we can always divide by a scalar, result is always an empty polynomial:
return polynomial<T>();
}
template <class U, class T>
inline polynomial<T> operator + (const U& a, const polynomial<T>& b)
inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator + (const U& a, polynomial<T> b)
{
polynomial<T> result(b);
result += a;
return result;
b += a;
return b;
}
template <class U, class T>
inline polynomial<T> operator - (const U& a, const polynomial<T>& b)
inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator - (const U& a, polynomial<T> b)
{
polynomial<T> result(a);
result -= b;
return result;
b -= a;
return -b;
}
template <class U, class T>
inline polynomial<T> operator * (const U& a, const polynomial<T>& b)
inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator * (const U& a, polynomial<T> b)
{
polynomial<T> result(b);
result *= a;
return result;
b *= a;
return b;
}
template <class T>
@@ -636,27 +795,31 @@ bool operator == (const polynomial<T> &a, const polynomial<T> &b)
return a.data() == b.data();
}
template <typename T, typename U>
polynomial<T> operator >> (const polynomial<T>& a, const U& b)
template <class T>
bool operator != (const polynomial<T> &a, const polynomial<T> &b)
{
polynomial<T> result(a);
result >>= b;
return result;
return a.data() != b.data();
}
template <typename T, typename U>
polynomial<T> operator << (const polynomial<T>& a, const U& b)
polynomial<T> operator >> (polynomial<T> a, const U& b)
{
polynomial<T> result(a);
result <<= b;
return result;
a >>= b;
return a;
}
template <typename T, typename U>
polynomial<T> operator << (polynomial<T> a, const U& b)
{
a <<= b;
return a;
}
// Unary minus (negate).
template <class T>
polynomial<T> operator - (polynomial<T> a)
{
std::transform(a.data().begin(), a.data().end(), a.data().begin(), std::negate<T>());
std::transform(a.data().begin(), a.data().end(), a.data().begin(), detail::negate());
return a;
}
@@ -713,7 +876,9 @@ inline std::basic_ostream<charT, traits>& operator << (std::basic_ostream<charT,
} // namespace math
} // namespace boost
//
// Polynomial specific overload of gcd algorithm:
//
#include <boost/math/tools/polynomial_gcd.hpp>
#endif // BOOST_MATH_TOOLS_POLYNOMIAL_HPP

209
winx64/include/boost/math/tools/polynomial_gcd.hpp Normal file
View File

@@ -0,0 +1,209 @@
// (C) Copyright Jeremy William Murphy 2016.
// 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_POLYNOMIAL_GCD_HPP
#define BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP
#ifdef _MSC_VER
#pragma once
#endif
#include <boost/math/tools/polynomial.hpp>
#include <boost/integer/common_factor_rt.hpp>
#include <boost/type_traits/is_pod.hpp>
namespace boost{
namespace integer {
namespace gcd_detail {
template <class T>
struct gcd_traits;
template <class T>
struct gcd_traits<boost::math::tools::polynomial<T> >
{
inline static const boost::math::tools::polynomial<T>& abs(const boost::math::tools::polynomial<T>& val) { return val; }
static const method_type method = method_euclid;
};
}
}
namespace math{ namespace tools{
/* From Knuth, 4.6.1:
*
* We may write any nonzero polynomial u(x) from R[x] where R is a UFD as
*
* u(x) = cont(u) . pp(u(x))
*
* where cont(u), the content of u, is an element of S, and pp(u(x)), the primitive
* part of u(x), is a primitive polynomial over S.
* When u(x) = 0, it is convenient to define cont(u) = pp(u(x)) = O.
*/
template <class T>
T content(polynomial<T> const &x)
{
return x ? boost::integer::gcd_range(x.data().begin(), x.data().end()).first : T(0);
}
// Knuth, 4.6.1
template <class T>
polynomial<T> primitive_part(polynomial<T> const &x, T const &cont)
{
return x ? x / cont : polynomial<T>();
}
template <class T>
polynomial<T> primitive_part(polynomial<T> const &x)
{
return primitive_part(x, content(x));
}
// Trivial but useful convenience function referred to simply as l() in Knuth.
template <class T>
T leading_coefficient(polynomial<T> const &x)
{
return x ? x.data().back() : T(0);
}
namespace detail
{
/* Reduce u and v to their primitive parts and return the gcd of their
* contents. Used in a couple of gcd algorithms.
*/
template <class T>
T reduce_to_primitive(polynomial<T> &u, polynomial<T> &v)
{
using boost::integer::gcd;
T const u_cont = content(u), v_cont = content(v);
u /= u_cont;
v /= v_cont;
return gcd(u_cont, v_cont);
}
}
/**
* Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
* Algorithm 4.6.1C: Greatest common divisor over a unique factorization domain.
*
* The subresultant algorithm by George E. Collins [JACM 14 (1967), 128-142],
* later improved by W. S. Brown and J. F. Traub [JACM 18 (1971), 505-514].
*
* Although step C3 keeps the coefficients to a "reasonable" size, they are
* still potentially several binary orders of magnitude larger than the inputs.
* Thus, this algorithm should only be used where T is a multi-precision type.
*
* @tparam T Polynomial coefficient type.
* @param u First polynomial.
* @param v Second polynomial.
* @return Greatest common divisor of polynomials u and v.
*/
template <class T>
typename enable_if_c< std::numeric_limits<T>::is_integer, polynomial<T> >::type
subresultant_gcd(polynomial<T> u, polynomial<T> v)
{
using std::swap;
BOOST_ASSERT(u || v);
if (!u)
return v;
if (!v)
return u;
typedef typename polynomial<T>::size_type N;
if (u.degree() < v.degree())
swap(u, v);
T const d = detail::reduce_to_primitive(u, v);
T g = 1, h = 1;
polynomial<T> r;
while (true)
{
BOOST_ASSERT(u.degree() >= v.degree());
// Pseudo-division.
r = u % v;
if (!r)
return d * primitive_part(v); // Attach the content.
if (r.degree() == 0)
return d * polynomial<T>(T(1)); // The content is the result.
N const delta = u.degree() - v.degree();
// Adjust remainder.
u = v;
v = r / (g * detail::integer_power(h, delta));
g = leading_coefficient(u);
T const tmp = detail::integer_power(g, delta);
if (delta <= N(1))
h = tmp * detail::integer_power(h, N(1) - delta);
else
h = tmp / detail::integer_power(h, delta - N(1));
}
}
/**
* @brief GCD for polynomials with unbounded multi-precision integral coefficients.
*
* The multi-precision constraint is enforced via numeric_limits.
*
* Note that intermediate terms in the evaluation can grow arbitrarily large, hence the need for
* unbounded integers, otherwise numeric loverflow would break the algorithm.
*
* @tparam T A multi-precision integral type.
*/
template <typename T>
typename enable_if_c<std::numeric_limits<T>::is_integer && !std::numeric_limits<T>::is_bounded, polynomial<T> >::type
gcd(polynomial<T> const &u, polynomial<T> const &v)
{
return subresultant_gcd(u, v);
}
// GCD over bounded integers is not currently allowed:
template <typename T>
typename enable_if_c<std::numeric_limits<T>::is_integer && std::numeric_limits<T>::is_bounded, polynomial<T> >::type
gcd(polynomial<T> const &u, polynomial<T> const &v)
{
BOOST_STATIC_ASSERT_MSG(sizeof(v) == 0, "GCD on polynomials of bounded integers is disallowed due to the excessive growth in the size of intermediate terms.");
return subresultant_gcd(u, v);
}
// GCD over polynomials of floats can go via the Euclid algorithm:
template <typename T>
typename enable_if_c<!std::numeric_limits<T>::is_integer && (std::numeric_limits<T>::min_exponent != std::numeric_limits<T>::max_exponent) && !std::numeric_limits<T>::is_exact, polynomial<T> >::type
gcd(polynomial<T> const &u, polynomial<T> const &v)
{
return boost::integer::gcd_detail::Euclid_gcd(u, v);
}
}
//
// Using declaration so we overload the default implementation in this namespace:
//
using boost::math::tools::gcd;
}
namespace integer
{
//
// Using declaration so we overload the default implementation in this namespace:
//
using boost::math::tools::gcd;
}
} // namespace boost::math::tools
#endif

2
winx64/include/boost/math/tools/precision.hpp
View File

@@ -377,7 +377,7 @@ inline T forth_root_epsilon_imp(const T*, const mpl::int_<0>&)
template <class T>
struct root_epsilon_traits
{
typedef mpl::int_< (::std::numeric_limits<T>::radix == 2) ? std::numeric_limits<T>::digits : 0> tag_type;
typedef mpl::int_< (::std::numeric_limits<T>::radix == 2) && (::std::numeric_limits<T>::digits != INT_MAX) ? std::numeric_limits<T>::digits : 0> tag_type;
BOOST_STATIC_CONSTANT(bool, has_noexcept = (tag_type::value == 113) || (tag_type::value == 64) || (tag_type::value == 53) || (tag_type::value == 24));
};

288
winx64/include/boost/math/tools/roots.hpp
View File

@@ -9,7 +9,10 @@
#ifdef _MSC_VER
#pragma once
#endif
#include <boost/multiprecision/detail/number_base.hpp> // test for multiprecision types.
#include <boost/type_traits/is_complex.hpp> // test for complex types
#include <iostream>
#include <utility>
#include <boost/config/no_tr1/cmath.hpp>
#include <stdexcept>
@@ -65,7 +68,7 @@ inline void unpack_0(const Tuple& t, T& val) BOOST_MATH_NOEXCEPT(T)
{
using dummy::get;
// Rely on ADL to find the correct overload of get:
val = get<0>(t);
val = get<0>(t);
}
template <class T, class U, class V>
@@ -267,8 +270,17 @@ T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_
#endif
if(fabs(delta * 2) > fabs(delta2))
{
// last two steps haven't converged, try bisection:
delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
// last two steps haven't converged.
T shift = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
if ((result != 0) && (fabs(shift) > fabs(result)))
{
delta = sign(delta) * result; // protect against huge jumps!
}
else
delta = shift;
// reset delta1/2 so we don't take this branch next time round:
delta1 = 3 * delta;
delta2 = 3 * delta;
}
guess = result;
result -= delta;
@@ -350,12 +362,11 @@ namespace detail{
T f0(0), f1, f2;
T result = guess;
T factor = static_cast<T>(ldexp(1.0, 1 - digits));
T factor = ldexp(static_cast<T>(1.0), 1 - digits);
T delta = (std::max)(T(10000000 * guess), T(10000000)); // arbitarily large delta
T last_f0 = 0;
T delta1 = delta;
T delta2 = delta;
bool out_of_bounds_sentry = false;
#ifdef BOOST_MATH_INSTRUMENT
@@ -415,13 +426,14 @@ namespace detail{
T convergence = fabs(delta / delta2);
if((convergence > 0.8) && (convergence < 2))
{
// last two steps haven't converged, try bisection:
// last two steps haven't converged.
delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
if(fabs(delta) > result)
if ((result != 0) && (fabs(delta) > result))
delta = sign(delta) * result; // protect against huge jumps!
// reset delta2 so that this branch will *not* be taken on the
// next iteration:
delta2 = delta * 3;
delta1 = delta * 3;
BOOST_MATH_INSTRUMENT_VARIABLE(delta);
}
guess = result;
@@ -431,7 +443,10 @@ namespace detail{
// check for out of bounds step:
if(result < min)
{
T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) ? T(1000) : T(result / min);
T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min)))
? T(1000)
: (fabs(min) < 1) && (fabs(tools::max_value<T>() * min) < fabs(result))
? ((min < 0) != (result < 0)) ? -tools::max_value<T>() : tools::max_value<T>() : T(result / min);
if(fabs(diff) < 1)
diff = 1 / diff;
if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
@@ -509,9 +524,10 @@ namespace detail{
template <class T>
static T step(const T& x, const T& f0, const T& f1, const T& f2) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T))
{
using std::fabs;
T ratio = f0 / f1;
T delta;
if(ratio / x < 0.1)
if((x != 0) && (fabs(ratio / x) < 0.1))
{
delta = ratio + (f2 / (2 * f1)) * ratio * ratio;
// check second derivative doesn't over compensate:
@@ -554,10 +570,262 @@ inline T schroeder_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEP
return schroder_iterate(f, guess, min, max, digits, m);
}
#ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS
/*
* Why do we set the default maximum number of iterations to the number of digits in the type?
* Because for double roots, the number of digits increases linearly with the number of iterations,
* so this default should recover full precision even in this somewhat pathological case.
* For isolated roots, the problem is so rapidly convergent that this doesn't matter at all.
*/
template<class Complex, class F>
Complex complex_newton(F g, Complex guess, int max_iterations=std::numeric_limits<typename Complex::value_type>::digits)
{
typedef typename Complex::value_type Real;
using std::norm;
using std::abs;
using std::max;
// z0, z1, and z2 cannot be the same, in case we immediately need to resort to Muller's Method:
Complex z0 = guess + Complex(1,0);
Complex z1 = guess + Complex(0,1);
Complex z2 = guess;
do {
auto pair = g(z2);
if (norm(pair.second) == 0)
{
// Muller's method. Notation follows Numerical Recipes, 9.5.2:
Complex q = (z2 - z1)/(z1 - z0);
auto P0 = g(z0);
auto P1 = g(z1);
Complex qp1 = static_cast<Complex>(1)+q;
Complex A = q*(pair.first - qp1*P1.first + q*P0.first);
Complex B = (static_cast<Complex>(2)*q+static_cast<Complex>(1))*pair.first - qp1*qp1*P1.first +q*q*P0.first;
Complex C = qp1*pair.first;
Complex rad = sqrt(B*B - static_cast<Complex>(4)*A*C);
Complex denom1 = B + rad;
Complex denom2 = B - rad;
Complex correction = (z1-z2)*static_cast<Complex>(2)*C;
if (norm(denom1) > norm(denom2))
{
correction /= denom1;
}
else
{
correction /= denom2;
}
z0 = z1;
z1 = z2;
z2 = z2 + correction;
}
else
{
z0 = z1;
z1 = z2;
z2 = z2 - (pair.first/pair.second);
}
// See: https://math.stackexchange.com/questions/3017766/constructing-newton-iteration-converging-to-non-root
// If f' is continuous, then convergence of x_n -> x* implies f(x*) = 0.
// This condition approximates this convergence condition by requiring three consecutive iterates to be clustered.
Real tol = max(abs(z2)*std::numeric_limits<Real>::epsilon(), std::numeric_limits<Real>::epsilon());
bool real_close = abs(z0.real() - z1.real()) < tol && abs(z0.real() - z2.real()) < tol && abs(z1.real() - z2.real()) < tol;
bool imag_close = abs(z0.imag() - z1.imag()) < tol && abs(z0.imag() - z2.imag()) < tol && abs(z1.imag() - z2.imag()) < tol;
if (real_close && imag_close)
{
return z2;
}
} while(max_iterations--);
// The idea is that if we can get abs(f) < eps, we should, but if we go through all these iterations
// and abs(f) < sqrt(eps), then roundoff error simply does not allow that we can evaluate f to < eps
// This is somewhat awkward as it isn't scale invariant, but using the Daubechies coefficient example code,
// I found this condition generates correct roots, whereas the scale invariant condition discussed here:
// https://scicomp.stackexchange.com/questions/30597/defining-a-condition-number-and-termination-criteria-for-newtons-method
// allows nonroots to be passed off as roots.
auto pair = g(z2);
if (abs(pair.first) < sqrt(std::numeric_limits<Real>::epsilon()))
{
return z2;
}
return {std::numeric_limits<Real>::quiet_NaN(),
std::numeric_limits<Real>::quiet_NaN()};
}
#endif
#if !defined(BOOST_NO_CXX17_IF_CONSTEXPR)
// https://stackoverflow.com/questions/48979861/numerically-stable-method-for-solving-quadratic-equations/50065711
namespace detail
{
template<class T>
inline T discriminant(T const & a, T const & b, T const & c)
{
T w = 4*a*c;
T e = std::fma(-c, 4*a, w);
T f = std::fma(b, b, -w);
return f + e;
}
}
template<class T>
auto quadratic_roots(T const& a, T const& b, T const& c)
{
using std::copysign;
using std::sqrt;
if constexpr (std::is_integral<T>::value)
{
// What I want is to write:
// return quadratic_roots(double(a), double(b), double(c));
// but that doesn't compile.
double nan = std::numeric_limits<double>::quiet_NaN();
if(a==0)
{
if (b==0 && c != 0)
{
return std::pair<double, double>(nan, nan);
}
else if (b==0 && c==0)
{
return std::pair<double, double>(0,0);
}
return std::pair<double, double>(-c/b, -c/b);
}
if (b==0)
{
double x0_sq = -double(c)/double(a);
if (x0_sq < 0) {
return std::pair<double, double>(nan, nan);
}
double x0 = sqrt(x0_sq);
return std::pair<double, double>(-x0,x0);
}
double discriminant = detail::discriminant(double(a), double(b), double(c));
if (discriminant < 0)
{
return std::pair<double, double>(nan, nan);
}
double q = -(b + copysign(sqrt(discriminant), double(b)))/T(2);
double x0 = q/a;
double x1 = c/q;
if (x0 < x1) {
return std::pair<double, double>(x0, x1);
}
return std::pair<double, double>(x1, x0);
}
else if constexpr (std::is_floating_point<T>::value)
{
T nan = std::numeric_limits<T>::quiet_NaN();
if(a==0)
{
if (b==0 && c != 0)
{
return std::pair<T, T>(nan, nan);
}
else if (b==0 && c==0)
{
return std::pair<T, T>(0,0);
}
return std::pair<T, T>(-c/b, -c/b);
}
if (b==0)
{
T x0_sq = -c/a;
if (x0_sq < 0) {
return std::pair<T, T>(nan, nan);
}
T x0 = sqrt(x0_sq);
return std::pair<T, T>(-x0,x0);
}
T discriminant = detail::discriminant(a, b, c);
// Is there a sane way to flush very small negative values to zero?
// If there is I don't know of it.
if (discriminant < 0)
{
return std::pair<T, T>(nan, nan);
}
T q = -(b + copysign(sqrt(discriminant), b))/T(2);
T x0 = q/a;
T x1 = c/q;
if (x0 < x1)
{
return std::pair<T, T>(x0, x1);
}
return std::pair<T, T>(x1, x0);
}
else if constexpr (boost::is_complex<T>::value || boost::multiprecision::number_category<T>::value == boost::multiprecision::number_kind_complex)
{
typename T::value_type nan = std::numeric_limits<typename T::value_type>::quiet_NaN();
if(a.real()==0 && a.imag() ==0)
{
using std::norm;
if (b.real()==0 && b.imag() && norm(c) != 0)
{
return std::pair<T, T>({nan, nan}, {nan, nan});
}
else if (b.real()==0 && b.imag() && c.real() ==0 && c.imag() == 0)
{
return std::pair<T, T>({0,0},{0,0});
}
return std::pair<T, T>(-c/b, -c/b);
}
if (b.real()==0 && b.imag() == 0)
{
T x0_sq = -c/a;
T x0 = sqrt(x0_sq);
return std::pair<T, T>(-x0, x0);
}
// There's no fma for complex types:
T discriminant = b*b - T(4)*a*c;
T q = -(b + sqrt(discriminant))/T(2);
return std::pair<T, T>(q/a, c/q);
}
else // Most likely the type is a boost.multiprecision.
{ //There is no fma for multiprecision, and in addition it doesn't seem to be useful, so revert to the naive computation.
T nan = std::numeric_limits<T>::quiet_NaN();
if(a==0)
{
if (b==0 && c != 0)
{
return std::pair<T, T>(nan, nan);
}
else if (b==0 && c==0)
{
return std::pair<T, T>(0,0);
}
return std::pair<T, T>(-c/b, -c/b);
}
if (b==0)
{
T x0_sq = -c/a;
if (x0_sq < 0) {
return std::pair<T, T>(nan, nan);
}
T x0 = sqrt(x0_sq);
return std::pair<T, T>(-x0,x0);
}
T discriminant = b*b - 4*a*c;
if (discriminant < 0)
{
return std::pair<T, T>(nan, nan);
}
T q = -(b + copysign(sqrt(discriminant), b))/T(2);
T x0 = q/a;
T x1 = c/q;
if (x0 < x1)
{
return std::pair<T, T>(x0, x1);
}
return std::pair<T, T>(x1, x0);
}
}
#endif
} // namespace tools
} // namespace math
} // namespace boost
#endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP

2
winx64/include/boost/math/tools/series.hpp
View File

@@ -35,7 +35,7 @@ inline typename Functor::result_type sum_series(Functor& func, const U& factor,
next_term = func();
result += next_term;
}
while((fabs(factor * result) < fabs(next_term)) && --counter);
while((abs(factor * result) < abs(next_term)) && --counter);
// set max_terms to the actual number of terms of the series evaluated:
max_terms = max_terms - counter;

346
winx64/include/boost/math/tools/signal_statistics.hpp Normal file
View File

@@ -0,0 +1,346 @@
// (C) Copyright Nick Thompson 2018.
// 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_SIGNAL_STATISTICS_HPP
#define BOOST_MATH_TOOLS_SIGNAL_STATISTICS_HPP
#include <algorithm>
#include <iterator>
#include <boost/type_traits/is_complex.hpp>
#include <boost/assert.hpp>
#include <boost/multiprecision/detail/number_base.hpp>
#include <boost/math/tools/roots.hpp>
#include <boost/math/tools/univariate_statistics.hpp>
namespace boost::math::tools {
template<class ForwardIterator>
auto absolute_gini_coefficient(ForwardIterator first, ForwardIterator last)
{
using std::abs;
using RealOrComplex = typename std::iterator_traits<ForwardIterator>::value_type;
BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Gini coefficient requires at least two samples.");
std::sort(first, last, [](RealOrComplex a, RealOrComplex b) { return abs(b) > abs(a); });
decltype(abs(*first)) i = 1;
decltype(abs(*first)) num = 0;
decltype(abs(*first)) denom = 0;
for (auto it = first; it != last; ++it)
{
decltype(abs(*first)) tmp = abs(*it);
num += tmp*i;
denom += tmp;
++i;
}
// If the l1 norm is zero, all elements are zero, so every element is the same.
if (denom == 0)
{
decltype(abs(*first)) zero = 0;
return zero;
}
return ((2*num)/denom - i)/(i-1);
}
template<class RandomAccessContainer>
inline auto absolute_gini_coefficient(RandomAccessContainer & v)
{
return boost::math::tools::absolute_gini_coefficient(v.begin(), v.end());
}
template<class ForwardIterator>
auto sample_absolute_gini_coefficient(ForwardIterator first, ForwardIterator last)
{
size_t n = std::distance(first, last);
return n*boost::math::tools::absolute_gini_coefficient(first, last)/(n-1);
}
template<class RandomAccessContainer>
inline auto sample_absolute_gini_coefficient(RandomAccessContainer & v)
{
return boost::math::tools::sample_absolute_gini_coefficient(v.begin(), v.end());
}
// The Hoyer sparsity measure is defined in:
// https://arxiv.org/pdf/0811.4706.pdf
template<class ForwardIterator>
auto hoyer_sparsity(const ForwardIterator first, const ForwardIterator last)
{
using T = typename std::iterator_traits<ForwardIterator>::value_type;
using std::abs;
using std::sqrt;
BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Hoyer sparsity requires at least two samples.");
if constexpr (std::is_unsigned<T>::value)
{
T l1 = 0;
T l2 = 0;
size_t n = 0;
for (auto it = first; it != last; ++it)
{
l1 += *it;
l2 += (*it)*(*it);
n += 1;
}
double rootn = sqrt(n);
return (rootn - l1/sqrt(l2) )/ (rootn - 1);
}
else {
decltype(abs(*first)) l1 = 0;
decltype(abs(*first)) l2 = 0;
// We wouldn't need to count the elements if it was a random access iterator,
// but our only constraint is that it's a forward iterator.
size_t n = 0;
for (auto it = first; it != last; ++it)
{
decltype(abs(*first)) tmp = abs(*it);
l1 += tmp;
l2 += tmp*tmp;
n += 1;
}
if constexpr (std::is_integral<T>::value)
{
double rootn = sqrt(n);
return (rootn - l1/sqrt(l2) )/ (rootn - 1);
}
else
{
decltype(abs(*first)) rootn = sqrt(static_cast<decltype(abs(*first))>(n));
return (rootn - l1/sqrt(l2) )/ (rootn - 1);
}
}
}
template<class Container>
inline auto hoyer_sparsity(Container const & v)
{
return boost::math::tools::hoyer_sparsity(v.cbegin(), v.cend());
}
template<class Container>
auto oracle_snr(Container const & signal, Container const & noisy_signal)
{
using Real = typename Container::value_type;
BOOST_ASSERT_MSG(signal.size() == noisy_signal.size(),
"Signal and noisy_signal must be have the same number of elements.");
if constexpr (std::is_integral<Real>::value)
{
double numerator = 0;
double denominator = 0;
for (size_t i = 0; i < signal.size(); ++i)
{
numerator += signal[i]*signal[i];
denominator += (noisy_signal[i] - signal[i])*(noisy_signal[i] - signal[i]);
}
if (numerator == 0 && denominator == 0)
{
return std::numeric_limits<double>::quiet_NaN();
}
if (denominator == 0)
{
return std::numeric_limits<double>::infinity();
}
return numerator/denominator;
}
else if constexpr (boost::is_complex<Real>::value ||
boost::multiprecision::number_category<Real>::value == boost::multiprecision::number_kind_complex)
{
using std::norm;
typename Real::value_type numerator = 0;
typename Real::value_type denominator = 0;
for (size_t i = 0; i < signal.size(); ++i)
{
numerator += norm(signal[i]);
denominator += norm(noisy_signal[i] - signal[i]);
}
if (numerator == 0 && denominator == 0)
{
return std::numeric_limits<typename Real::value_type>::quiet_NaN();
}
if (denominator == 0)
{
return std::numeric_limits<typename Real::value_type>::infinity();
}
return numerator/denominator;
}
else
{
Real numerator = 0;
Real denominator = 0;
for (size_t i = 0; i < signal.size(); ++i)
{
numerator += signal[i]*signal[i];
denominator += (signal[i] - noisy_signal[i])*(signal[i] - noisy_signal[i]);
}
if (numerator == 0 && denominator == 0)
{
return std::numeric_limits<Real>::quiet_NaN();
}
if (denominator == 0)
{
return std::numeric_limits<Real>::infinity();
}
return numerator/denominator;
}
}
template<class Container>
auto mean_invariant_oracle_snr(Container const & signal, Container const & noisy_signal)
{
using Real = typename Container::value_type;
BOOST_ASSERT_MSG(signal.size() == noisy_signal.size(), "Signal and noisy signal must be have the same number of elements.");
Real mu = boost::math::tools::mean(signal);
Real numerator = 0;
Real denominator = 0;
for (size_t i = 0; i < signal.size(); ++i)
{
Real tmp = signal[i] - mu;
numerator += tmp*tmp;
denominator += (signal[i] - noisy_signal[i])*(signal[i] - noisy_signal[i]);
}
if (numerator == 0 && denominator == 0)
{
return std::numeric_limits<Real>::quiet_NaN();
}
if (denominator == 0)
{
return std::numeric_limits<Real>::infinity();
}
return numerator/denominator;
}
template<class Container>
auto mean_invariant_oracle_snr_db(Container const & signal, Container const & noisy_signal)
{
using std::log10;
return 10*log10(boost::math::tools::mean_invariant_oracle_snr(signal, noisy_signal));
}
// Follows the definition of SNR given in Mallat, A Wavelet Tour of Signal Processing, equation 11.16.
template<class Container>
auto oracle_snr_db(Container const & signal, Container const & noisy_signal)
{
using std::log10;
return 10*log10(boost::math::tools::oracle_snr(signal, noisy_signal));
}
// A good reference on the M2M4 estimator:
// D. R. Pauluzzi and N. C. Beaulieu, "A comparison of SNR estimation techniques for the AWGN channel," IEEE Trans. Communications, Vol. 48, No. 10, pp. 1681-1691, 2000.
// A nice python implementation:
// https://github.com/gnuradio/gnuradio/blob/master/gr-digital/examples/snr_estimators.py
template<class ForwardIterator>
auto m2m4_snr_estimator(ForwardIterator first, ForwardIterator last, decltype(*first) estimated_signal_kurtosis=1, decltype(*first) estimated_noise_kurtosis=3)
{
BOOST_ASSERT_MSG(estimated_signal_kurtosis > 0, "The estimated signal kurtosis must be positive");
BOOST_ASSERT_MSG(estimated_noise_kurtosis > 0, "The estimated noise kurtosis must be positive.");
using Real = typename std::iterator_traits<ForwardIterator>::value_type;
using std::sqrt;
if constexpr (std::is_floating_point<Real>::value ||
boost::multiprecision::number_category<Real>::value == boost::multiprecision::number_kind_floating_point)
{
// If we first eliminate N, we obtain the quadratic equation:
// (ka+kw-6)S^2 + 2M2(3-kw)S + kw*M2^2 - M4 = 0 =: a*S^2 + bs*N + cs = 0
// If we first eliminate S, we obtain the quadratic equation:
// (ka+kw-6)N^2 + 2M2(3-ka)N + ka*M2^2 - M4 = 0 =: a*N^2 + bn*N + cn = 0
// I believe these equations are totally independent quadratics;
// if one has a complex solution it is not necessarily the case that the other must also.
// However, I can't prove that, so there is a chance that this does unnecessary work.
// Future improvements: There are algorithms which can solve quadratics much more effectively than the naive implementation found here.
// See: https://stackoverflow.com/questions/48979861/numerically-stable-method-for-solving-quadratic-equations/50065711#50065711
auto [M1, M2, M3, M4] = boost::math::tools::first_four_moments(first, last);
if (M4 == 0)
{
// The signal is constant. There is no noise:
return std::numeric_limits<Real>::infinity();
}
// Change to notation in Pauluzzi, equation 41:
auto kw = estimated_noise_kurtosis;
auto ka = estimated_signal_kurtosis;
// A common case, since it's the default:
Real a = (ka+kw-6);
Real bs = 2*M2*(3-kw);
Real cs = kw*M2*M2 - M4;
Real bn = 2*M2*(3-ka);
Real cn = ka*M2*M2 - M4;
auto [S0, S1] = boost::math::tools::quadratic_roots(a, bs, cs);
if (S1 > 0)
{
auto N = M2 - S1;
if (N > 0)
{
return S1/N;
}
if (S0 > 0)
{
N = M2 - S0;
if (N > 0)
{
return S0/N;
}
}
}
auto [N0, N1] = boost::math::tools::quadratic_roots(a, bn, cn);
if (N1 > 0)
{
auto S = M2 - N1;
if (S > 0)
{
return S/N1;
}
if (N0 > 0)
{
S = M2 - N0;
if (S > 0)
{
return S/N0;
}
}
}
// This happens distressingly often. It's a limitation of the method.
return std::numeric_limits<Real>::quiet_NaN();
}
else
{
BOOST_ASSERT_MSG(false, "The M2M4 estimator has not been implemented for this type.");
return std::numeric_limits<Real>::quiet_NaN();
}
}
template<class Container>
inline auto m2m4_snr_estimator(Container const & noisy_signal, typename Container::value_type estimated_signal_kurtosis=1, typename Container::value_type estimated_noise_kurtosis=3)
{
return m2m4_snr_estimator(noisy_signal.cbegin(), noisy_signal.cend(), estimated_signal_kurtosis, estimated_noise_kurtosis);
}
template<class ForwardIterator>
inline auto m2m4_snr_estimator_db(ForwardIterator first, ForwardIterator last, decltype(*first) estimated_signal_kurtosis=1, decltype(*first) estimated_noise_kurtosis=3)
{
using std::log10;
return 10*log10(m2m4_snr_estimator(first, last, estimated_signal_kurtosis, estimated_noise_kurtosis));
}
template<class Container>
inline auto m2m4_snr_estimator_db(Container const & noisy_signal, typename Container::value_type estimated_signal_kurtosis=1, typename Container::value_type estimated_noise_kurtosis=3)
{
using std::log10;
return 10*log10(m2m4_snr_estimator(noisy_signal, estimated_signal_kurtosis, estimated_noise_kurtosis));
}
}
#endif

118
winx64/include/boost/math/tools/test_value.hpp Normal file
View File

@@ -0,0 +1,118 @@
// Copyright Paul A. Bristow 2017.
// Copyright John Maddock 2017.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
// test_value.hpp
#ifndef TEST_VALUE_HPP
#define TEST_VALUE_HPP
// BOOST_MATH_TEST_VALUE is used to create a test value of suitable type from a decimal digit string.
// Two parameters, both a floating-point literal double like 1.23 (not long double so no suffix L)
// and a decimal digit string const char* like "1.23" must be provided.
// The decimal value represented must be the same of course, with at least enough precision for long double.
// Note there are two gotchas to this approach:
// * You need all values to be real floating-point values
// * and *MUST* include a decimal point (to avoid confusion with an integer literal).
// * It's slow to compile compared to a simple literal.
// Speed is not an issue for a few test values,
// but it's not generally usable in large tables
// where you really need everything to be statically initialized.
// Macro BOOST_MATH_INSTRUMENT_CREATE_TEST_VALUE provides a global diagnostic value for create_type.
#include <boost/cstdfloat.hpp> // For float_64_t, float128_t. Must be first include!
#include <boost/lexical_cast.hpp>
#include <boost/type_traits/is_constructible.hpp>
#include <boost/type_traits/is_convertible.hpp>
#ifdef BOOST_MATH_INSTRUMENT_CREATE_TEST_VALUE
// global int create_type(0); must be defined before including this file.
#endif
#ifdef BOOST_HAS_FLOAT128
typedef __float128 largest_float;
#define BOOST_MATH_TEST_LARGEST_FLOAT_SUFFIX(x) x##Q
#define BOOST_MATH_TEST_LARGEST_FLOAT_DIGITS 113
#else
typedef long double largest_float;
#define BOOST_MATH_TEST_LARGEST_FLOAT_SUFFIX(x) x##L
#define BOOST_MATH_TEST_LARGEST_FLOAT_DIGITS std::numeric_limits<long double>::digits
#endif
template <class T, class T2>
inline T create_test_value(largest_float val, const char*, const boost::mpl::true_&, const T2&)
{ // Construct from long double or quad parameter val (ignoring string/const char* str).
// (This is case for MPL parameters = true_ and T2 == false_,
// and MPL parameters = true_ and T2 == true_ cpp_bin_float)
// All built-in/fundamental floating-point types,
// and other User-Defined Types that can be constructed without loss of precision
// from long double suffix L (or quad suffix Q),
//
// Choose this method, even if can be constructed from a string,
// because it will be faster, and more likely to be the closest representation.
// (This is case for MPL parameters = mpl::true_ and T2 == mpl::true_).
#ifdef BOOST_MATH_INSTRUMENT_CREATE_TEST_VALUE
create_type = 1;
#endif
return static_cast<T>(val);
}
template <class T>
inline T create_test_value(largest_float, const char* str, const boost::mpl::false_&, const boost::mpl::true_&)
{ // Construct from decimal digit string const char* @c str (ignoring long double parameter).
// For example, extended precision or other User-Defined types which ARE constructible from a string
// (but not from double, or long double without loss of precision).
// (This is case for MPL parameters = mpl::false_ and T2 == mpl::true_).
#ifdef BOOST_MATH_INSTRUMENT_CREATE_TEST_VALUE
create_type = 2;
#endif
return T(str);
}
template <class T>
inline T create_test_value(largest_float, const char* str, const boost::mpl::false_&, const boost::mpl::false_&)
{ // Create test value using from lexical cast of decimal digit string const char* str.
// For example, extended precision or other User-Defined types which are NOT constructible from a string
// (NOR constructible from a long double).
// (This is case T1 = mpl::false and T2 == mpl::false).
#ifdef BOOST_MATH_INSTRUMENT_CREATE_TEST_VALUE
create_type = 3;
#endif
return boost::lexical_cast<T>(str);
}
// T real type, x a decimal digits representation of a floating-point, for example: 12.34.
// It must include a decimal point (or it would be interpreted as an integer).
// x is converted to a long double by appending the letter L (to suit long double fundamental type), 12.34L.
// x is also passed as a const char* or string representation "12.34"
// (to suit most other types that cannot be constructed from long double without possible loss).
// BOOST_MATH_TEST_LARGEST_FLOAT_SUFFIX(x) makes a long double or quad version, with
// suffix a letter L (or Q) to suit long double (or quad) fundamental type, 12.34L or 12.34Q.
// #x makes a decimal digit string version to suit multiprecision and fixed_point constructors, "12.34".
// (Constructing from double or long double (or quad) could lose precision for multiprecision or fixed-point).
// The matching create_test_value function above is chosen depending on the T1 and T2 mpl bool truths.
// The string version from #x is used if the precision of T is greater than long double.
// Example: long double test_value = BOOST_MATH_TEST_VALUE(double, 1.23456789);
#define BOOST_MATH_TEST_VALUE(T, x) create_test_value<T>(\
BOOST_MATH_TEST_LARGEST_FLOAT_SUFFIX(x),\
#x,\
boost::mpl::bool_<\
std::numeric_limits<T>::is_specialized &&\
(std::numeric_limits<T>::radix == 2)\
&& (std::numeric_limits<T>::digits <= BOOST_MATH_TEST_LARGEST_FLOAT_DIGITS)\
&& boost::is_convertible<largest_float, T>::value>(),\
boost::mpl::bool_<\
boost::is_constructible<T, const char*>::value>()\
)
#endif // TEST_VALUE_HPP

8
winx64/include/boost/math/tools/toms748_solve.hpp
View File

@@ -302,6 +302,12 @@ std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax, const
static const char* function = "boost::math::tools::toms748_solve<%1%>";
//
// Sanity check - are we allowed to iterate at all?
//
if (max_iter == 0)
return std::make_pair(ax, bx);
boost::uintmax_t count = max_iter;
T a, b, fa, fb, c, u, fu, a0, b0, d, fd, e, fe;
static const T mu = 0.5f;
@@ -477,6 +483,8 @@ inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax
template <class F, class T, class Tol, class Policy>
inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
{
if (max_iter <= 2)
return std::make_pair(ax, bx);
max_iter -= 2;
std::pair<T, T> r = toms748_solve(f, ax, bx, f(ax), f(bx), tol, max_iter, pol);
max_iter += 2;

393
winx64/include/boost/math/tools/univariate_statistics.hpp Normal file
View File

@@ -0,0 +1,393 @@
// (C) Copyright Nick Thompson 2018.
// 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_UNIVARIATE_STATISTICS_HPP
#define BOOST_MATH_TOOLS_UNIVARIATE_STATISTICS_HPP
#include <algorithm>
#include <iterator>
#include <boost/type_traits/is_complex.hpp>
#include <boost/assert.hpp>
#include <boost/multiprecision/detail/number_base.hpp>
namespace boost::math::tools {
template<class ForwardIterator>
auto mean(ForwardIterator first, ForwardIterator last)
{
using Real = typename std::iterator_traits<ForwardIterator>::value_type;
BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute the mean.");
if constexpr (std::is_integral<Real>::value)
{
double mu = 0;
double i = 1;
for(auto it = first; it != last; ++it) {
mu = mu + (*it - mu)/i;
i += 1;
}
return mu;
}
else
{
Real mu = 0;
Real i = 1;
for(auto it = first; it != last; ++it) {
mu = mu + (*it - mu)/i;
i += 1;
}
return mu;
}
}
template<class Container>
inline auto mean(Container const & v)
{
return mean(v.cbegin(), v.cend());
}
template<class ForwardIterator>
auto variance(ForwardIterator first, ForwardIterator last)
{
using Real = typename std::iterator_traits<ForwardIterator>::value_type;
BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute mean and variance.");
// Higham, Accuracy and Stability, equation 1.6a and 1.6b:
if constexpr (std::is_integral<Real>::value)
{
double M = *first;
double Q = 0;
double k = 2;
for (auto it = std::next(first); it != last; ++it)
{
double tmp = *it - M;
Q = Q + ((k-1)*tmp*tmp)/k;
M = M + tmp/k;
k += 1;
}
return Q/(k-1);
}
else
{
Real M = *first;
Real Q = 0;
Real k = 2;
for (auto it = std::next(first); it != last; ++it)
{
Real tmp = (*it - M)/k;
Q += k*(k-1)*tmp*tmp;
M += tmp;
k += 1;
}
return Q/(k-1);
}
}
template<class Container>
inline auto variance(Container const & v)
{
return variance(v.cbegin(), v.cend());
}
template<class ForwardIterator>
auto sample_variance(ForwardIterator first, ForwardIterator last)
{
size_t n = std::distance(first, last);
BOOST_ASSERT_MSG(n > 1, "At least two samples are required to compute the sample variance.");
return n*variance(first, last)/(n-1);
}
template<class Container>
inline auto sample_variance(Container const & v)
{
return sample_variance(v.cbegin(), v.cend());
}
// Follows equation 1.5 of:
// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
template<class ForwardIterator>
auto skewness(ForwardIterator first, ForwardIterator last)
{
using Real = typename std::iterator_traits<ForwardIterator>::value_type;
BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute skewness.");
if constexpr (std::is_integral<Real>::value)
{
double M1 = *first;
double M2 = 0;
double M3 = 0;
double n = 2;
for (auto it = std::next(first); it != last; ++it)
{
double delta21 = *it - M1;
double tmp = delta21/n;
M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
M2 = M2 + tmp*(n-1)*delta21;
M1 = M1 + tmp;
n += 1;
}
double var = M2/(n-1);
if (var == 0)
{
// The limit is technically undefined, but the interpretation here is clear:
// A constant dataset has no skewness.
return double(0);
}
double skew = M3/(M2*sqrt(var));
return skew;
}
else
{
Real M1 = *first;
Real M2 = 0;
Real M3 = 0;
Real n = 2;
for (auto it = std::next(first); it != last; ++it)
{
Real delta21 = *it - M1;
Real tmp = delta21/n;
M3 += tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
M2 += tmp*(n-1)*delta21;
M1 += tmp;
n += 1;
}
Real var = M2/(n-1);
if (var == 0)
{
// The limit is technically undefined, but the interpretation here is clear:
// A constant dataset has no skewness.
return Real(0);
}
Real skew = M3/(M2*sqrt(var));
return skew;
}
}
template<class Container>
inline auto skewness(Container const & v)
{
return skewness(v.cbegin(), v.cend());
}
// Follows equation 1.5/1.6 of:
// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
template<class ForwardIterator>
auto first_four_moments(ForwardIterator first, ForwardIterator last)
{
using Real = typename std::iterator_traits<ForwardIterator>::value_type;
BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute the first four moments.");
if constexpr (std::is_integral<Real>::value)
{
double M1 = *first;
double M2 = 0;
double M3 = 0;
double M4 = 0;
double n = 2;
for (auto it = std::next(first); it != last; ++it)
{
double delta21 = *it - M1;
double tmp = delta21/n;
M4 = M4 + tmp*(tmp*tmp*delta21*((n-1)*(n*n-3*n+3)) + 6*tmp*M2 - 4*M3);
M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
M2 = M2 + tmp*(n-1)*delta21;
M1 = M1 + tmp;
n += 1;
}
return std::make_tuple(M1, M2/(n-1), M3/(n-1), M4/(n-1));
}
else
{
Real M1 = *first;
Real M2 = 0;
Real M3 = 0;
Real M4 = 0;
Real n = 2;
for (auto it = std::next(first); it != last; ++it)
{
Real delta21 = *it - M1;
Real tmp = delta21/n;
M4 = M4 + tmp*(tmp*tmp*delta21*((n-1)*(n*n-3*n+3)) + 6*tmp*M2 - 4*M3);
M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
M2 = M2 + tmp*(n-1)*delta21;
M1 = M1 + tmp;
n += 1;
}
return std::make_tuple(M1, M2/(n-1), M3/(n-1), M4/(n-1));
}
}
template<class Container>
inline auto first_four_moments(Container const & v)
{
return first_four_moments(v.cbegin(), v.cend());
}
// Follows equation 1.6 of:
// https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
template<class ForwardIterator>
auto kurtosis(ForwardIterator first, ForwardIterator last)
{
auto [M1, M2, M3, M4] = first_four_moments(first, last);
if (M2 == 0)
{
return M2;
}
return M4/(M2*M2);
}
template<class Container>
inline auto kurtosis(Container const & v)
{
return kurtosis(v.cbegin(), v.cend());
}
template<class ForwardIterator>
auto excess_kurtosis(ForwardIterator first, ForwardIterator last)
{
return kurtosis(first, last) - 3;
}
template<class Container>
inline auto excess_kurtosis(Container const & v)
{
return excess_kurtosis(v.cbegin(), v.cend());
}
template<class RandomAccessIterator>
auto median(RandomAccessIterator first, RandomAccessIterator last)
{
size_t num_elems = std::distance(first, last);
BOOST_ASSERT_MSG(num_elems > 0, "The median of a zero length vector is undefined.");
if (num_elems & 1)
{
auto middle = first + (num_elems - 1)/2;
std::nth_element(first, middle, last);
return *middle;
}
else
{
auto middle = first + num_elems/2 - 1;
std::nth_element(first, middle, last);
std::nth_element(middle, middle+1, last);
return (*middle + *(middle+1))/2;
}
}
template<class RandomAccessContainer>
inline auto median(RandomAccessContainer & v)
{
return median(v.begin(), v.end());
}
template<class RandomAccessIterator>
auto gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)
{
using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;
BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Gini coefficient requires at least two samples.");
std::sort(first, last);
if constexpr (std::is_integral<Real>::value)
{
double i = 1;
double num = 0;
double denom = 0;
for (auto it = first; it != last; ++it)
{
num += *it*i;
denom += *it;
++i;
}
// If the l1 norm is zero, all elements are zero, so every element is the same.
if (denom == 0)
{
return double(0);
}
return ((2*num)/denom - i)/(i-1);
}
else
{
Real i = 1;
Real num = 0;
Real denom = 0;
for (auto it = first; it != last; ++it)
{
num += *it*i;
denom += *it;
++i;
}
// If the l1 norm is zero, all elements are zero, so every element is the same.
if (denom == 0)
{
return Real(0);
}
return ((2*num)/denom - i)/(i-1);
}
}
template<class RandomAccessContainer>
inline auto gini_coefficient(RandomAccessContainer & v)
{
return gini_coefficient(v.begin(), v.end());
}
template<class RandomAccessIterator>
inline auto sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)
{
size_t n = std::distance(first, last);
return n*gini_coefficient(first, last)/(n-1);
}
template<class RandomAccessContainer>
inline auto sample_gini_coefficient(RandomAccessContainer & v)
{
return sample_gini_coefficient(v.begin(), v.end());
}
template<class RandomAccessIterator>
auto median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator last, typename std::iterator_traits<RandomAccessIterator>::value_type center=std::numeric_limits<typename std::iterator_traits<RandomAccessIterator>::value_type>::quiet_NaN())
{
using std::abs;
using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;
using std::isnan;
if (isnan(center))
{
center = boost::math::tools::median(first, last);
}
size_t num_elems = std::distance(first, last);
BOOST_ASSERT_MSG(num_elems > 0, "The median of a zero-length vector is undefined.");
auto comparator = [&center](Real a, Real b) { return abs(a-center) < abs(b-center);};
if (num_elems & 1)
{
auto middle = first + (num_elems - 1)/2;
std::nth_element(first, middle, last, comparator);
return abs(*middle);
}
else
{
auto middle = first + num_elems/2 - 1;
std::nth_element(first, middle, last, comparator);
std::nth_element(middle, middle+1, last, comparator);
return (abs(*middle) + abs(*(middle+1)))/abs(static_cast<Real>(2));
}
}
template<class RandomAccessContainer>
inline auto median_absolute_deviation(RandomAccessContainer & v, typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN())
{
return median_absolute_deviation(v.begin(), v.end(), center);
}
}
#endif