c++ - 具有大值跨度的 double 组求和 : proper algorithm

Question

我有一个算法，我需要对 e-40 到 e+40 范围内的双数求和(很多时间)。

数组示例(从实际应用程序中随机转储):

-2.06991e-05 
7.58132e-06 
-3.91367e-06 
7.38921e-07 
-5.33143e-09
-4.13195e-11 
4.01724e-14 
6.03221e-17 
-4.4202e-20
6.58873 
-1.22257
-0.0606178 
0.00036508 
2.67599e-07 
0
-627.061
-59.048 
5.92985 
0.0885884
0.000276455 
-2.02579e-07

不言而喻，我知道这会导致舍入效应，我正在努力控制它:最终结果不应该在 double 的小数部分有任何缺失的信息，或者，如果不是可以避免的结果应该至少是 n 位准确(定义了 n)。最终结果需要 5 位数字加上指数。

经过一番深思熟虑，我最终得出了以下算法:

• 对数组进行排序，使绝对值最大的排在最前面，最接近零的排在最后。
• 在循环中添加所有内容

这个想法是，在这种情况下，任何大值(负值和正值)的取消都不会影响后面的较小值。 简而言之:

• (10e40 - 10e40) + 1 = 1 : 结果符合预期
• (1 + 10e-40) - 10e40 = 0 :不好

我最终使用 std::multiset(我的 PC 上的基准测试给出了 20% 的速度，与普通 double 相比，长 double - 我对 double 分辨率没问题)和使用 std:fabs 的自定义排序功能。

它仍然很慢(完成整个过程需要 5 秒)而且我仍然有这种“你在算法中遗漏了一些东西”的感觉。任何建议:

1. 用于速度优化。有没有更好的方法来分拣中间产品？对一组 40 个中间结果进行排序(通常)需要大约 70% 的总执行时间。
2. 针对遗漏的问题。是否仍有机会丢失关键数据(本应位于最终结果的小数部分的数据)？

---

从更大的角度来看，我正在实现纯虚变量(电阻抗:Z(jw))的实系数多项式类。 Z是代表用户定义系统的大多项式，系数指数范围很远。
“大”来自于将 Zc1 = 1/jC1w 添加到 Zc2 = 1/jC2w 中:
Zc1 + Zc2 = (C1C2(jw)^2 + 0(jw))/(C1+C2)(jw)<​​/p>

在这种情况下，C1 和 C2 在纳法拉 (10e-9) 中，C1C2 已经在 10e-18 中(而且它才刚刚开始......)

我的排序函数使用复数变量的曼哈顿距离(因为，我的是纯实数或纯虚数):

struct manhattan_complex_distance
{
        bool operator() (std::complex<long double> a, std::complex<long double> b)
        {
            return std::fabs(std::real(a) + std::imag(a)) > std::fabs(std::real(b) + std::imag(b));
        }
};

我的多组 Action :

std:complex<long double> get_value(std::vector<std::complex<long double>>& frequency_vector)
{
    //frequency_vector is precalculated once for all to have at index n the value (jw)^n. 
    std::multiset<std::complex<long double>, manhattan_distance> temp_list;   
    for (int i=0; i<m_coeficients.size(); ++i)
    {
        //   element of :       ℝ         *         ℂ
        temp_list.insert(m_coeficients[i] * frequency_vector[i]);
    }
    std::complex<long double> ret=0;
    for (auto i:temp_list)
    {
        // it is VERY important to start adding the big values before adding the small ones.
        // in informatics, 10^60 - 10^60 + 1 = 1; while 1 + 10^60 - 10^60 = 0. Of course you'd expected to get 1, not 0.
        ret += i;
    }
    return ret;
}

我的项目启用了c++11(主要是为了改进数学库和复数工具)

ps:我重构了代码以使其易于阅读，实际上所有复数和长双名称都是模板:我可以立即更改多项式类型或使用 ℝ 的正则多项式类

Answer 1

作为GuyGreer建议，可以用Kahan summation :

double sum = 0.0;
double c = 0.0;
for (double value : values) {
    double y = value - c;
    double t = sum + y;
    c = (t - sum) - y;
    sum = t;
}

编辑:您还应该考虑使用 Horner's method计算多项式。

double value = coeffs[degree];
for (auto i = degree; i-- > 0;) {
    value *= x;
    value += coeffs[i];
}