我有一个算法,我需要对 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 秒)而且我仍然有这种“你在算法中遗漏了一些东西”的感觉。任何建议:
- 用于速度优化。有没有更好的方法来分拣中间产品?对一组 40 个中间结果进行排序(通常)需要大约 70% 的总执行时间。
- 针对遗漏的问题。是否仍有机会丢失关键数据(本应位于最终结果的小数部分的数据)?
从更大的角度来看,我正在实现纯虚变量(电阻抗: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:我重构了代码以使其易于阅读,实际上所有复数和长双名称都是模板:我可以立即更改多项式类型或使用 ℝ 的正则多项式类
最佳答案
作为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];
}
关于c++ - 具有大值跨度的 double 组求和 : proper algorithm,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/37950122/