boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Vectorize mean calculation.

Browse Source
This commit is contained in:
Nick
2019-04-13 17:22:24 -04:00
parent 748172f279
commit 165cf9e480
2 changed files with 67 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

43
include/boost/math/tools/univariate_statistics.hpp
View File

@@ -30,12 +30,47 @@ auto mean(ForwardIterator first, ForwardIterator last)
}
return mu;
}
else if constexpr (std::is_same_v<typename std::iterator_traits<ForwardIterator>::iterator_category, std::random_access_iterator_tag>)
{
size_t elements = std::distance(first, last);
Real mu0 = 0;
Real mu1 = 0;
Real mu2 = 0;
Real mu3 = 0;
Real i = 1;
auto end = last - (elements % 4);
for(auto it = first; it != end; it += 4) {
Real inv = Real(1)/i;
Real tmp0 = (*it - mu0);
Real tmp1 = (*(it+1) - mu1);
Real tmp2 = (*(it+2) - mu2);
Real tmp3 = (*(it+3) - mu3);
// please generate a vectorized fma here
mu0 += tmp0*inv;
mu1 += tmp1*inv;
mu2 += tmp2*inv;
mu3 += tmp3*inv;
i += 1;
}
Real num1 = Real(elements - (elements %4))/Real(4);
Real num2 = num1 + Real(elements % 4);
for (auto it = end; it != last; ++it)
{
mu3 += (*it-mu3)/i;
i += 1;
}
return (num1*(mu0+mu1+mu2) + num2*mu3)/Real(elements);
}
else
{
Real mu = 0;
Real i = 1;
for(auto it = first; it != last; ++it) {
mu = mu + (*it - mu)/i;
auto it = first;
Real mu = *it;
Real i = 2;
while(++it != last)
{
mu += (*it - mu)/i;
i += 1;
}
return mu;

28
test/univariate_statistics_test.cpp
View File

@@ -124,6 +124,17 @@ void test_integer_mean()
BOOST_TEST(abs(m1 - m2) < tol*abs(m1));
}
template<class RandomAccessContainer>
auto naive_mean(RandomAccessContainer const & v)
{
typename RandomAccessContainer::value_type sum = 0;
for (auto & x : v)
{
sum += x;
}
return sum/v.size();
}
template<class Real>
void test_mean()
{
@@ -172,6 +183,23 @@ void test_mean()
}
Real m2 = boost::math::tools::mean(v);
BOOST_TEST(abs(m1 - m2) < tol*abs(m1));
// Stress test:
for (size_t i = 1; i < 30; ++i)
{
v = generate_random_vector<Real>(i, 12803);
auto naive_ = naive_mean(v);
auto higham_ = boost::math::tools::mean(v);
if (abs(higham_ - naive_) >= 100*tol*abs(naive_))
{
std::cout << std::hexfloat;
std::cout << "Terms = " << v.size() << "\n";
std::cout << "higham = " << higham_ << "\n";
std::cout << "naive_ = " << naive_ << "\n";
}
BOOST_TEST(abs(higham_ - naive_) < 100*tol*abs(naive_));
}
}
template<class Complex>