在 C++ 中使用位掩码消除浮点错误
Eliminating floating point errors with bit masks in c++
假设我有一个问题,我想计算大量的双打,每次我计算一个新的双打,我想检查一下我以前是否见过这个双打。最好(最快)的方法是什么?
如果双打是准确的,我将创建一个 set
并检查成员资格,即 O(log n) 或 O(1) 的散列集。但它们并不准确,因此我需要遍历所有以前看到的双精度数并检查它们是否在新计算的双精度数的容差范围内,这很慢:O(n)。我的想法是只保留(比方说)double 的前 30 位,在集合中给出 2^-30 的精度并检查它们是否相等。这是个好主意还是有更好的方法?如何只保留 double 的 30 个最高有效位?
无论采用何种舍入模式,您总是会画一条线,落在线的一侧或另一侧的数字可以任意接近(在 Float 的情况下为 1 ulp)。
所以你需要至少使用两个舍入不同的集合,并检查一个舍入模式或另一个舍入模式的浮点数是否在集合中。
例如,对于公差 t
,它类似于:
is_close = ( round(x,t) is in set1) or ( round(x+t/2,t) is in set2 )
set1 add round(x,t)
set2 add round(x+t/2,t)
我建议使用 floor 或 ceiling 作为规律性的截断模式。
注1:n*t+t-eps
和n*t-t+eps
对于整数n
和小的0 < eps < t/100
,会被认为很近,距离接近2 t
,所以在上述公式中 t 应选择为半公差。
注2:这是绝对容差。如果目的是相对公差(截断有效数字的位 - 这会降低精度),那么处理二进制边界可能会涉及更多公式,但应该应用相同类型的算法(只是 t 会浮动...... .)
如果您需要确定您之前是否见过某个范围内的给定值(a
到 b
),那么您可以使用 std::set
所有以前看到的值。我们可以用upper_bound()
找到第一个大于a
的元素,然后测试是否小于b
:
// untested
bool have_seen(const std::set<double>& past, double a, double b)
{
auto it = past.upper_bound(a);
return it != past.end() && *it < b;
}
(如果要测试包含范围,请使用 lower_bound
和 <=
。)
当然,速度为 O(log n)。
您可以选择 std::nextafter
的多次迭代以获得适合您的值的 a
和 b
。
如果您确实需要降低 double 的精度(并且不清楚这是否是您特定问题的解决方案),那么如果您的系统使用通用标准二进制格式,这并不会太困难。
我们可以利用最低有效数字包含尾数这一事实,并将正确的数字设置为零:
#include <cstdint>
#include <limits>
double trunc(double d, int precision) noexcept
{
// This function requires IEEE-754 binary representations
static_assert(std::numeric_limits<double>::is_iec559);
static_assert(std::numeric_limits<double>::radix == 2);
using Integer = std::uint_fast32_t;
static const auto max_precision = std::numeric_limits<double>::digits;
if (precision < 1) {
// invalid
return 0;
}
if (precision >= max_precision) {
// no-op
return d;
}
auto& i = reinterpret_cast<Integer&>(d);
static_assert(sizeof i >= sizeof d);
auto mask = ~Integer{} << (max_precision - precision);
i &= mask;
return d;
}
示范:
#include <iostream>
int main()
{
for (double d: {100, 101, 102, 103, 104, 105, 106,
200, 202, 206, 207, 208, 209, 210}) {
std::cout << d << " -> " << trunc(d, 5) << '\n';
}
}
输出:
100 -> 100
101 -> 100
102 -> 100
103 -> 100
104 -> 104
105 -> 104
106 -> 104
200 -> 200
202 -> 200
206 -> 200
207 -> 200
208 -> 208
209 -> 208
210 -> 208
这是否真的是您需要的很难说。如果您想将结果收集到直方图桶中,这可能很有用:
#include <iostream>
#include <iomanip>
#include <map>
#include <random>
#include <string>
int main()
{
std::mt19937 gen(std::random_device{}());
std::normal_distribution<double> dist(0.5,0.1);
std::map<double,std::size_t> histogram;
for (int i = 0; i < 10000; ++i) {
auto d = trunc(dist(gen), 5);
++histogram[d];
}
for (auto const& [value, freq]: histogram) {
std::cout << std::fixed << std::setprecision(3) << std::setw(5)
<< value << ": " << std::string(freq/50, '*') << '\n';
}
}
0.203:
0.211:
0.219:
0.227:
0.234:
0.242:
0.250:
0.266:
0.281: *
0.297: *
0.312: **
0.328: ***
0.344: ***
0.359: *****
0.375: ******
0.391: *******
0.406: ********
0.422: *********
0.438: **********
0.453: ***********
0.469: ************
0.484: ************
0.500: *************************
0.531: **********************
0.562: ******************
0.594: **************
0.625: *********
0.656: *****
0.688: **
0.719: *
0.750:
0.781:
0.812:
0.844:
虽然这个正态分布 看起来 有偏差,但这只是因为当指数增加时桶大小在 0.5 处变化。
即使您负担不起将所有双打都保存在 RAM 中,也可以使用此方法。
- 收集双打名单。 O(N)
- 对它们进行排序。 O(N*logN)
- 浏览列表。对于每个项目,要么输出它,要么扔掉它,具体取决于它与最后一个输出值的接近程度。 O(n)
如果“'close' 的测试”成本很高,请注意该算法只执行 N 次测试。一些建议的算法是 O(N*M),其中 N 是输入计数,M 是输出计数。
进一步注意输出是不确定的。我的意思是,重新排列输入列表可能会改变 保留哪些 个双打,哪些被丢弃。此警告可能适用于任何和所有算法。
假设我有一个问题,我想计算大量的双打,每次我计算一个新的双打,我想检查一下我以前是否见过这个双打。最好(最快)的方法是什么?
如果双打是准确的,我将创建一个 set
并检查成员资格,即 O(log n) 或 O(1) 的散列集。但它们并不准确,因此我需要遍历所有以前看到的双精度数并检查它们是否在新计算的双精度数的容差范围内,这很慢:O(n)。我的想法是只保留(比方说)double 的前 30 位,在集合中给出 2^-30 的精度并检查它们是否相等。这是个好主意还是有更好的方法?如何只保留 double 的 30 个最高有效位?
无论采用何种舍入模式,您总是会画一条线,落在线的一侧或另一侧的数字可以任意接近(在 Float 的情况下为 1 ulp)。
所以你需要至少使用两个舍入不同的集合,并检查一个舍入模式或另一个舍入模式的浮点数是否在集合中。
例如,对于公差 t
,它类似于:
is_close = ( round(x,t) is in set1) or ( round(x+t/2,t) is in set2 )
set1 add round(x,t)
set2 add round(x+t/2,t)
我建议使用 floor 或 ceiling 作为规律性的截断模式。
注1:n*t+t-eps
和n*t-t+eps
对于整数n
和小的0 < eps < t/100
,会被认为很近,距离接近2 t
,所以在上述公式中 t 应选择为半公差。
注2:这是绝对容差。如果目的是相对公差(截断有效数字的位 - 这会降低精度),那么处理二进制边界可能会涉及更多公式,但应该应用相同类型的算法(只是 t 会浮动...... .)
如果您需要确定您之前是否见过某个范围内的给定值(a
到 b
),那么您可以使用 std::set
所有以前看到的值。我们可以用upper_bound()
找到第一个大于a
的元素,然后测试是否小于b
:
// untested
bool have_seen(const std::set<double>& past, double a, double b)
{
auto it = past.upper_bound(a);
return it != past.end() && *it < b;
}
(如果要测试包含范围,请使用 lower_bound
和 <=
。)
当然,速度为 O(log n)。
您可以选择 std::nextafter
的多次迭代以获得适合您的值的 a
和 b
。
如果您确实需要降低 double 的精度(并且不清楚这是否是您特定问题的解决方案),那么如果您的系统使用通用标准二进制格式,这并不会太困难。
我们可以利用最低有效数字包含尾数这一事实,并将正确的数字设置为零:
#include <cstdint>
#include <limits>
double trunc(double d, int precision) noexcept
{
// This function requires IEEE-754 binary representations
static_assert(std::numeric_limits<double>::is_iec559);
static_assert(std::numeric_limits<double>::radix == 2);
using Integer = std::uint_fast32_t;
static const auto max_precision = std::numeric_limits<double>::digits;
if (precision < 1) {
// invalid
return 0;
}
if (precision >= max_precision) {
// no-op
return d;
}
auto& i = reinterpret_cast<Integer&>(d);
static_assert(sizeof i >= sizeof d);
auto mask = ~Integer{} << (max_precision - precision);
i &= mask;
return d;
}
示范:
#include <iostream>
int main()
{
for (double d: {100, 101, 102, 103, 104, 105, 106,
200, 202, 206, 207, 208, 209, 210}) {
std::cout << d << " -> " << trunc(d, 5) << '\n';
}
}
输出:
100 -> 100
101 -> 100
102 -> 100
103 -> 100
104 -> 104
105 -> 104
106 -> 104
200 -> 200
202 -> 200
206 -> 200
207 -> 200
208 -> 208
209 -> 208
210 -> 208
这是否真的是您需要的很难说。如果您想将结果收集到直方图桶中,这可能很有用:
#include <iostream>
#include <iomanip>
#include <map>
#include <random>
#include <string>
int main()
{
std::mt19937 gen(std::random_device{}());
std::normal_distribution<double> dist(0.5,0.1);
std::map<double,std::size_t> histogram;
for (int i = 0; i < 10000; ++i) {
auto d = trunc(dist(gen), 5);
++histogram[d];
}
for (auto const& [value, freq]: histogram) {
std::cout << std::fixed << std::setprecision(3) << std::setw(5)
<< value << ": " << std::string(freq/50, '*') << '\n';
}
}
0.203:
0.211:
0.219:
0.227:
0.234:
0.242:
0.250:
0.266:
0.281: *
0.297: *
0.312: **
0.328: ***
0.344: ***
0.359: *****
0.375: ******
0.391: *******
0.406: ********
0.422: *********
0.438: **********
0.453: ***********
0.469: ************
0.484: ************
0.500: *************************
0.531: **********************
0.562: ******************
0.594: **************
0.625: *********
0.656: *****
0.688: **
0.719: *
0.750:
0.781:
0.812:
0.844:
虽然这个正态分布 看起来 有偏差,但这只是因为当指数增加时桶大小在 0.5 处变化。
即使您负担不起将所有双打都保存在 RAM 中,也可以使用此方法。
- 收集双打名单。 O(N)
- 对它们进行排序。 O(N*logN)
- 浏览列表。对于每个项目,要么输出它,要么扔掉它,具体取决于它与最后一个输出值的接近程度。 O(n)
如果“'close' 的测试”成本很高,请注意该算法只执行 N 次测试。一些建议的算法是 O(N*M),其中 N 是输入计数,M 是输出计数。
进一步注意输出是不确定的。我的意思是,重新排列输入列表可能会改变 保留哪些 个双打,哪些被丢弃。此警告可能适用于任何和所有算法。