快速获取接近2次幂的数字(浮点）
Fast way to get a close power-of-2 number (floating-point)
在数值计算中，通常需要将数字缩放到安全范围内。
例如，计算欧几里得距离：sqrt(a^2+b^2)。这里，如果a或b的大小太小/太大，则可能发生下溢/上溢。
解决这一问题的一种常见方法是将数字除以最大震级。然而，这个解决方案是：
- 慢(除法很慢)
- 导致一点额外的不准确性
所以我想，与其除以最大的星等数，不如乘以一个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(我不想使用低于标准的scaler，因为低于标准的数字可能很慢)。
什么是快速的方法？这可以用纯浮点运算完成吗？或者，我应该从value中提取指数，并使用简单的if来进行逻辑运算？是否有某种技巧可以快速与(-)1022进行比较(因为范围是对称的)？
注意：scaler不需要是最接近的2次方。如果某些逻辑需要它，scaler可以是离最接近值2的某个小幂。
函数s = get_scale(z)计算"2的闭合幂"。由于s的分数位为零，s的倒数只是一个(廉价的)整数减法：参见函数inv_of_scale。
在x86上，CCD_ 23和CCD_。编译器clang将三元运算符分别转换为CCD_ 25和CCD_，另请参阅Peter Cordes的评论。有了gcc将这些函数转换为x86内部函数代码(get_scale_x86和inv_of_scale_x86)，请参阅Godbolt。
注意C明确允许类型punning通过联合，而C++(C++11)没有这样的权限虽然gcc8.2和clang7.0没有抱怨联合，但您可以改进用CCD_ 29技巧代替工会把戏。对代码的这种修改应该是微不足道的。代码应该正确处理子规范。
#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 %14en", 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 %14en", 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和vminpd指令。
根据wim的回答，这里有另一个解决方案，它可以更快，因为它少了一条指令。输出有点不同，但仍然满足要求。
其想法是使用位运算来修复边界情况：将01放在指数的lsb中，无论其值如何。因此，指数：
- 0变为1(-1023变为-1022)
- 2046变为2045(1023变为1022)
- 其他指数也进行了修改，但只是略有修改：与wim的解决方案相比，这个数字可能会大两倍(当指数lsb从
00变为01时)，或者减半(当10->01时)或1/4(当11->01时
所以，这个修改后的例程可以工作(我认为只需2条快速asm指令就可以解决这个问题，这很酷)：
#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;
uint64_t and_i;
uint64_t or_i;
/* 0xFEDCBA9876543210 */
and_i = 0x7FD0000000000000ull;
or_i = 0x0010000000000000ull;
x.d = t;
x.i = (x.i & and_i)|or_i; /* Set fraction bits to zero, take absolute value */
return x.d;
}
double get_scale_x86(double t){
__m128d x = _mm_set_sd(t);
__m128d x_and = _mm_castsi128_pd(_mm_set1_epi64x(0x7FD0000000000000ull));
__m128d x_or = _mm_castsi128_pd(_mm_set1_epi64x(0x0010000000000000ull));
x = _mm_and_pd(x, x_and);
x = _mm_or_pd(x, x_or);
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 %14en", 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 %14en", z, s, u, s*u, z*u);
}
return 0;
}
您可以使用
double frexp (double x, int* exp);
返回值是x的小数部分，exp是指数(减去偏移量)。
或者，下面的代码获取double的指数部分。
int get_exp(double *d) {
long long *l = (long long *) d;
return ((*l & (0x7ffLL << 52) )>> 52)-1023 ;
}
相关文章：
- 从.txt文件中读取浮点型数字并在公式中使用它们
- 给定数字的浮点分辨率
- 用C++将浮点数字转换为本地化字符串
- C++中的 Json:将数字解析为字符串以避免浮点不准确
- printf如何从浮点数字中提取数字
- printf 的浮点格式标志 (%f) 仅适用于英文数字格式
- 浮点值未显示准确数字的数字C
- 浮点最低和最大的数字限制
- 不应该浮点 6 和双 15 可用数字吗？
- 快速获取接近2次幂的数字(浮点）
- php和c++中数字格式/浮点之间的区别
- 在Linux上的C++中，当在一行上打印5000个浮点数字时，换行
- 关于C中浮点算术的异常，从本身中减去数字
- 编译器为什么要将浮点数字的位数固定为6
- C++中是否有一个浮点文字后缀来使数字的精度翻倍
- 完全递增/递减浮点型数字
- 有没有一种方法可以设置浮点变量=某个数字*字符串中的变量
- 在C++中将字符串转换为带完整有效数字的浮点值[修复了在开发人员C++中工作的问题]
- 使用数组查找最接近的数字C++，浮点
- 浮点运算是如何在一个大数字上加一的
