最近在思考一个很有趣的问题，float和double两种数据结构，一个32位、一个64位，占用内存差一半，精度会有明显差异吗？看到过一个结论，在大多数问题中，float的精度已经足够。而double精度则非常恐怖，例如用double来存储π计算太阳系的周长，和真实的周长相比，double计算出来的周长误差仅仅为3微米，而float计算出来的周长误差约为1公里。
为了验证float和double的精度问题，我编写了以下代码来模拟计算调和级数
#include <iostream>
#include <iomanip>
#include <cmath>
// 设置求和项数
const size_t N = 100000000; // 1e8 项，可根据需要调整
// float 精度求和
float harmonic_sum_float(size_t n) {
float sum = 0.0f;
for (size_t i = 1; i <= n; ++i) {
sum += 1.0f / i;
}
return sum;
}
// double 精度求和
double harmonic_sum_double(size_t n) {
double sum = 0.0;
for (size_t i = 1; i <= n; ++i) {
sum += 1.0 / i;
}
return sum;
}
// 渐近估算值：H_n ≈ ln(n) + γ（欧拉-马歇罗尼常数）
double harmonic_approx(size_t n) {
const double gamma = 0.57721566490153286060;
return std::log(n) + gamma;
}
int main() {
std::cout << std::fixed << std::setprecision(10);
std::cout << "Calculating harmonic sum H_" << N << "...\n";
float sum_f = harmonic_sum_float(N);
double sum_d = harmonic_sum_double(N);
double sum_ref = harmonic_approx(N);
std::cout << "\n--- float ---\n";
std::cout << "Sum : " << sum_f << "\n";
std::cout << "Abs Error : " << std::fabs(sum_f - sum_ref) << "\n";
std::cout << "\n--- double ---\n";
std::cout << "Sum : " << sum_d << "\n";
std::cout << "Abs Error : " << std::fabs(sum_d - sum_ref) << "\n";
std::cout << "\n--- Reference ---\n";
std::cout << "ln(n)+γ : " << sum_ref << "\n";
return 0;
}
它的输出如下：
Calculating harmonic sum H_100000000...
--- float ---
Sum : 15.4036827087
Abs Error : 3.5942137001
--- double ---
Sum : 18.9978964139
Abs Error : 0.0000000050
--- Reference ---
ln(n)+γ : 18.9978964089当继续加大N的数量级，我们会发现float收敛在15.4左右，而double加法持续增加，没有明显的误差。而我们知道调和级数是是不收敛的，当k趋向于无穷大，和也会趋于无穷大。
那么为什么使用float计算调和级数，会产生“假收敛”呢？因为发生了“大数吃小数”。举一个例子：一亿加一等于多少？二年级的小朋友都知道是一亿零一，我们看看计算机的答案：
#include <iomanip>
#include <iostream>
int main() {
constexpr float a = 1.0e8f;
constexpr double b = 1.0e8;
std::cout << std::setprecision(12) << a + 1.0f << std::endl;
std::cout << std::setprecision(12) << b + 1.0 << std::endl;
return 0;
}
// 输出：
100000000
100000001诶，我们发现使用float，对一亿加一还是一亿，这意味着，无论我重复多少次给一个大数加上一个小数，它都不会继续变大了。这就是为什么float调和级数到一定项数之后，它不再增加了。
那么，除了使用double之外，我们有什么办法可以解决这个问题呢？我首先想到的是归并求和，即小数之间互相两两求和，再和大数相加。这种方法是可以解决问题的，但是如果面对的不是调和级数，而是乱序的数，可能需要先排序再求和。而GPT教我了一种不需要排序，可以在O(n)时间复杂度内，经过一次遍历就能求和的算法：Kahan 求和。它的代码很短：
float kahanSum(vector<float> nums) {
float sum = 0.0f;
float c = 0.0f;
for (auto num : nums) {
float y = num - c;
float t = sum + y;
c = (t - sum) - y;
sum = t;
}
return sum;
}我尝试用这种方法计算float调和级数前一亿项，结果如下：
--- float kahan---
Sum : 18.9978961945
Abs Error : 0.0000002144取得了和double相当的精度！
简单解释其原理：c代表当前累计误差，将c从当前数字num上减除，得到修正后的当前项y。将修正后的当前项与sum相加，得到新的sum，记位t。接下来是最关键的一步，计算新的c。我们看到c = (t – sum) – y。又t = sum + y ，这里是比较容易疑惑的，难道 c = (sum + y – sum) – y，那不就是0吗？其实不然，原因是在t = sum + y这一步发生了浮点精度损失，导致t 和 真正的sum+y之间存在误差。举一个例子，假如这里sum = 100’000’000，y=1，根据我们刚刚的实验，t=100’000’000，1并没有被加进去。那么c = (t-sum)-y=(100’000’000-100’000’000) – y = -1。我们将-1保存了下来，存入了新的c当中。如果下一次加法，新的加数又是1，那我们会将这个-1从新的加数上减去，新的加数会变为2。这就是可汗求和的核心思想，将求和后的结果与原结果做差，算出sum的增量，将sum的增量与加数做差，得到增量的误差量，作为下一轮的修正量。如此一来我们就成功经过一次遍历就算出了误差极小的数组和。
发表回复