c++ - 计算许多数字的几何平均值的有效方法

Question

我需要计算大量数字的几何平均值，其值不受先验限制。天真的方法是

double geometric_mean(std::vector<double> const&data) // failure
{
  auto product = 1.0;
  for(auto x:data) product *= x;
  return std::pow(product,1.0/data.size());
}

但是，由于累积的 product 中的下溢或溢出，这很可能会失败(注意:long double 并不能真正避免这个问题)。因此，下一个选项是将对数相加:

double geometric_mean(std::vector<double> const&data)
{
  auto sumlog = 0.0;
  for(auto x:data) sum_log += std::log(x);
  return std::exp(sum_log/data.size());
}

这可行，但为每个元素调用 std::log()，这可能很慢。 我可以避免这种情况吗？例如，通过分别跟踪(相当于)累积 product 的指数和尾数？

Answer 1

“拆分指数和尾数”解决方案:

double geometric_mean(std::vector<double> const & data)
{
    double m = 1.0;
    long long ex = 0;
    double invN = 1.0 / data.size();

    for (double x : data)
    {
        int i;
        double f1 = std::frexp(x,&i);
        m*=f1;
        ex+=i;
    }

    return std::pow( std::numeric_limits<double>::radix,ex * invN) * std::pow(m,invN);
}

如果您担心 ex 可能会溢出，您可以将其定义为 double 而不是 long long，然后乘以 invN 的每一步，但使用这种方法可能会丢失很多精度。

编辑对于大型输入，我们可以将计算分成几个桶:

double geometric_mean(std::vector<double> const & data)
{
    long long ex = 0;
    auto do_bucket = [&data,&ex](int first,int last) -> double
    {
        double ans = 1.0;
        for ( ;first != last;++first)
        {
            int i;
            ans *= std::frexp(data[first],&i);
            ex+=i;
        }
        return ans;
    };

    const int bucket_size = -std::log2( std::numeric_limits<double>::min() );
    std::size_t buckets = data.size() / bucket_size;

    double invN = 1.0 / data.size();
    double m = 1.0;

    for (std::size_t i = 0;i < buckets;++i)
        m *= std::pow( do_bucket(i * bucket_size,(i+1) * bucket_size),invN );

    m*= std::pow( do_bucket( buckets * bucket_size, data.size() ),invN );

    return std::pow( std::numeric_limits<double>::radix,ex * invN ) * m;
}