求和具有大值跨度的双精度数组:适当的算法

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

不言而喻,我知道这将导致舍入效应,我试图控制它:最终结果不应该在双倍的小数部分中有任何缺失的信息,或者如果不可避免的结果应该在最小n位精确(n定义).最终结果需要5位数加上指数.

经过一番体面的思考,我最终得到了以下算法:

• 对数组进行排序,使最大的绝对值先到最后,最接近零.
• 在循环中添加所有内容

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

• (10e40 - 10e40)+ 1 = 1:结果与预期一致
• (1 + 10e-40) - 10e40 = 0:不好

我最终使用std :: multiset(我的PC上的基准测试使用std:fabs,使用自定义排序功能,使用长双倍比普通双打更快20% - 我很好用双倍分辨率).

它仍然很慢(完成整个过程需要5秒钟)而且我仍然有这种感觉"你错过了你的算法".任何建议:

1. 用于速度优化.有没有更好的方法来分类中间产品？对一组40个中间结果(通常)进行排序大约占总执行时间的70%.
2. 对于错过的问题.是否还有机会丢失关键数据(应该是最终结果的小数部分)？

---

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

在这种情况下,对于纳米法(10e-9)的C1和C2,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));
        }
};

和我的多重作用:

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:我重构代码使得易于阅读,实际上所有复合体和长双精度都是模板:我可以立即更改多项式类型或使用类作为regular的正则多项式

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;
}

编辑：您还应该考虑使用霍纳的方法来评估多项式。

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