c++ - 获得接近 2 次幂数的快速方法(浮点数)
问题描述
在数值计算中,通常需要将数字缩放到安全范围内。
例如,计算欧几里得距离:sqrt(a^2+b^2)
。在这里,如果a
orb
的大小太小/太大,则可能发生下溢/溢出。
解决此问题的常用方法是将数字除以最大幅度数。但是,此解决方案是:
- 慢(除法很慢)
- 导致一些额外的不准确
所以我认为,与其除以最大幅度数,不如将它乘以一个接近的 2 次幂倒数。这似乎是一个更好的解决方案,因为:
- 乘法比除法快得多
- 更好的准确性,因为乘以 2 的幂数是准确的
所以,我想创建一个小的实用函数,它具有这样的逻辑(^
我的意思是求幂):
void getScaler(double value, double &scaler, double &scalerReciprocal) {
int e = <exponent of value>;
if (e<-1022) { scaler=2^-1022; scalerReciprocal = 2^1022; }
} else if (e>1022) { scaler=2^1022; scalerReciprocal = 2^-1022; }
} else { scaler=2^e; scalerReciprocal = 2^(2046-e); }
}
这个函数应该返回一个归一化的scaler
& scalerReciprocal
,两者都是 2 的幂数,其中scaler
接近value
,并且scalerReciprocal
是 的倒数scaler
。
scaler
/的最大允许指数scaleReciprocal
是-1022..1022
(我不想使用 subnormal scaler
,因为 subnormal 数字可能很慢)。
有什么快速的方法可以做到这一点?这可以通过纯浮点运算来完成吗?或者我应该从中提取指数value
,并使用简单if
的 s 来执行逻辑?是否有某种技巧可以快速与 (-)1022 进行比较(因为范围是对称的)?
注意:scaler
不需要是最接近的 2 次幂。如果某些逻辑需要它,scaler
可以与最接近的值相差一些小的 2 次方。
解决方案
函数s = get_scale(z)
计算“2 的关闭幂”。由于 的小数位s
为零,因此 的倒数s
只是一个(便宜的)整数减法:参见 function inv_of_scale
。
在 x86 上get_scale
并inv_of_scale
使用 clang 编译为非常有效的程序集。编译器 clang 将三元运算符转换为minsd
and maxsd
,另请参见 Peter Cordes 的评论。get_scale_x86
使用 gcc,将这些函数转换为 x86 内部代码 (和inv_of_scale_x86
)会更有效,请参阅 Godbolt。
请注意,C 明确允许通过联合进行类型双关,而 C++ (c++11) 没有这样的权限
虽然 gcc 8.2 和 clang 7.0 不会抱怨联合,但您可以通过使用memcpy
技巧而不是工会把戏。对代码的这种修改应该是微不足道的。代码应正确处理次规范。
#include<stdio.h>
#include<stdint.h>
#include<immintrin.h>
/* gcc -Wall -m64 -O3 -march=sandybridge dbl_scale.c */
union dbl_int64{
double d;
uint64_t i;
};
double get_scale(double t){
union dbl_int64 x;
union dbl_int64 x_min;
union dbl_int64 x_max;
uint64_t mask_i;
/* 0xFEDCBA9876543210 */
x_min.i = 0x0010000000000000ull;
x_max.i = 0x7FD0000000000000ull;
mask_i = 0x7FF0000000000000ull;
x.d = t;
x.i = x.i & mask_i; /* Set fraction bits to zero, take absolute value */
x.d = (x.d < x_min.d) ? x_min.d : x.d; /* If subnormal: set exponent to 1 */
x.d = (x.d > x_max.d) ? x_max.d : x.d; /* If exponent is very large: set exponent to 7FD, otherwise the inverse is a subnormal */
return x.d;
}
double get_scale_x86(double t){
__m128d x = _mm_set_sd(t);
__m128d x_min = _mm_castsi128_pd(_mm_set1_epi64x(0x0010000000000000ull));
__m128d x_max = _mm_castsi128_pd(_mm_set1_epi64x(0x7FD0000000000000ull));
__m128d mask = _mm_castsi128_pd(_mm_set1_epi64x(0x7FF0000000000000ull));
x = _mm_and_pd(x, mask);
x = _mm_max_sd(x, x_min);
x = _mm_min_sd(x, x_max);
return _mm_cvtsd_f64(x);
}
/* Compute the inverse 1/t of a double t with all zero fraction bits */
/* and exponent between the limits of function get_scale */
/* A single integer subtraction is much less expensive than a */
/* floating point division. */
double inv_of_scale(double t){
union dbl_int64 x;
/* 0xFEDCBA9876543210 */
uint64_t inv_mask = 0x7FE0000000000000ull;
x.d = t;
x.i = inv_mask - x.i;
return x.d;
}
double inv_of_scale_x86(double t){
__m128i inv_mask = _mm_set1_epi64x(0x7FE0000000000000ull);
__m128d x = _mm_set_sd(t);
__m128i x_i = _mm_sub_epi64(inv_mask, _mm_castpd_si128(x));
return _mm_cvtsd_f64(_mm_castsi128_pd(x_i));
}
int main(){
int n = 14;
int i;
/* Several example values, 4.94e-324 is the smallest subnormal */
double y[14] = { 4.94e-324, 1.1e-320, 1.1e-300, 1.1e-5, 0.7, 1.7, 123.1, 1.1e300,
1.79e308, -1.1e-320, -0.7, -1.7, -123.1, -1.1e307};
double z, s, u;
printf("Portable code:\n");
printf(" x pow_of_2 inverse pow2*inv x*inverse \n");
for (i = 0; i < n; i++){
z = y[i];
s = get_scale(z);
u = inv_of_scale(s);
printf("%14e %14e %14e %14e %14e\n", z, s, u, s*u, z*u);
}
printf("\nx86 specific SSE code:\n");
printf(" x pow_of_2 inverse pow2*inv x*inverse \n");
for (i = 0; i < n; i++){
z = y[i];
s = get_scale_x86(z);
u = inv_of_scale_x86(s);
printf("%14e %14e %14e %14e %14e\n", z, s, u, s*u, z*u);
}
return 0;
}
输出看起来不错:
Portable code:
x pow_of_2 inverse pow2*inv x*inverse
4.940656e-324 2.225074e-308 4.494233e+307 1.000000e+00 2.220446e-16
1.099790e-320 2.225074e-308 4.494233e+307 1.000000e+00 4.942713e-13
1.100000e-300 7.466109e-301 1.339386e+300 1.000000e+00 1.473324e+00
1.100000e-05 7.629395e-06 1.310720e+05 1.000000e+00 1.441792e+00
7.000000e-01 5.000000e-01 2.000000e+00 1.000000e+00 1.400000e+00
1.700000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.700000e+00
1.231000e+02 6.400000e+01 1.562500e-02 1.000000e+00 1.923437e+00
1.100000e+300 6.696929e+299 1.493222e-300 1.000000e+00 1.642544e+00
1.790000e+308 4.494233e+307 2.225074e-308 1.000000e+00 3.982882e+00
-1.099790e-320 2.225074e-308 4.494233e+307 1.000000e+00 -4.942713e-13
-7.000000e-01 5.000000e-01 2.000000e+00 1.000000e+00 -1.400000e+00
-1.700000e+00 1.000000e+00 1.000000e+00 1.000000e+00 -1.700000e+00
-1.231000e+02 6.400000e+01 1.562500e-02 1.000000e+00 -1.923437e+00
-1.100000e+307 5.617791e+306 1.780059e-307 1.000000e+00 -1.958065e+00
x86 specific SSE code:
x pow_of_2 inverse pow2*inv x*inverse
4.940656e-324 2.225074e-308 4.494233e+307 1.000000e+00 2.220446e-16
1.099790e-320 2.225074e-308 4.494233e+307 1.000000e+00 4.942713e-13
1.100000e-300 7.466109e-301 1.339386e+300 1.000000e+00 1.473324e+00
1.100000e-05 7.629395e-06 1.310720e+05 1.000000e+00 1.441792e+00
7.000000e-01 5.000000e-01 2.000000e+00 1.000000e+00 1.400000e+00
1.700000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.700000e+00
1.231000e+02 6.400000e+01 1.562500e-02 1.000000e+00 1.923437e+00
1.100000e+300 6.696929e+299 1.493222e-300 1.000000e+00 1.642544e+00
1.790000e+308 4.494233e+307 2.225074e-308 1.000000e+00 3.982882e+00
-1.099790e-320 2.225074e-308 4.494233e+307 1.000000e+00 -4.942713e-13
-7.000000e-01 5.000000e-01 2.000000e+00 1.000000e+00 -1.400000e+00
-1.700000e+00 1.000000e+00 1.000000e+00 1.000000e+00 -1.700000e+00
-1.231000e+02 6.400000e+01 1.562500e-02 1.000000e+00 -1.923437e+00
-1.100000e+307 5.617791e+306 1.780059e-307 1.000000e+00 -1.958065e+00
矢量化
函数get_scale
应该使用支持自动矢量化的编译器进行矢量化。以下代码可以很好地使用 clang 进行矢量化(无需编写 SSE/AVX 内部代码)。
/* Test how well get_scale vectorizes: */
void get_scale_vec(double * __restrict__ t, double * __restrict__ x){
int n = 1024;
int i;
for (i = 0; i < n; i++){
x[i] = get_scale(t[i]);
}
}
不幸的是 gcc 没有找到vmaxpd
andvminpd
说明。
推荐阅读
- python - Python请求Json数据,想从excel复制回csv文件
- apache-camel - Camel Bindy 在 Talend ESB Karaf 上给出未解决的依赖错误
- macos - Apple 权限地狱:osascript 和 Accessibility(辅助访问)
- node.js - NodeJs列出具有绝对路径的文件中的文件名
- r - 正则表达式删除不相关的键值 - R
- primefaces - 在 JSF 中使用 PrimeFaces 的 FontAwesome
- excel - VBA 目标变量复制问题
- r - 如何列出每个类别的所有唯一值
- javascript - JOI - 验证复杂对象
- kubernetes - 404 使用 Istio 从 K8s pod 中点击 `metadata.google.internal`