快速获得接近 2 的幂数（浮点数）的方法
Fast way to get a close power-of-2 number (floating-point)
在数值计算中，经常需要将数字缩放到安全范围内。
例如计算欧氏距离：sqrt(a^2+b^2)。在这里，如果 a 或 b 的大小太 small/large，则可能会发生 underflow/overflow。
解决这个问题的一种常见方法是将数字除以最大幅度的数字。然而，这个解决方案是：
- 慢（除法慢）
- 导致一些额外的不准确
所以我认为与其除以最大幅度的数，不如将它乘以一个接近的 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); }
}
这个函数应该 return 归一化 scaler 和 scalerReciprocal，两者都是 2 的幂数，其中 scaler 接近 value , scalerReciprocal 是 scaler.
的倒数
scaler/scaleReciprocal 的最大允许指数是 -1022..1022（我不想使用次正规 scaler，因为次正规数可能很慢） .
快速的方法是什么？这可以用纯浮点运算来完成吗？或者我应该从 value 中提取指数并使用简单的 ifs 来执行逻辑？是否有某种技巧可以快速与 (-)1022 进行比较（因为范围是对称的）？
注意：scaler 不需要是最接近的 2 的幂。如果某些逻辑需要它，scaler 可以与最接近的值相差一些小的 2 次幂。
您可以使用
double frexp (double x, int* exp);
返回值是x的小数部分，exp是指数（减去偏移量）。
或者，以下代码获取双精度数的指数部分。
int get_exp(double *d) {
long long *l = (long long *) d;
return ((*l & (0x7ffLL << 52) )>> 52)-1023 ;
}
函数 s = get_scale(z) 计算 "close power of 2"。由于 s 的小数位
为零，s 的倒数只是一个（廉价的）整数减法：参见函数 inv_of_scale。
在 x86 上 get_scale 和 inv_of_scale 使用 clang 编译为非常高效的汇编。
编译器 clang 将三元运算符转换为 minsd 和 maxsd，
另见 Peter Cordes 的 。
使用 gcc，效率稍微高一点
将这些函数转换为 x86 内在函数
代码（get_scale_x86 和 inv_of_scale_x86），see Godbolt.
请注意 C 明确允许类型双关
通过联合，而 C++ (c++11) 没有这样的权限
虽然 gcc 8.2 和 clang 7.0 并没有抱怨联合，但你可以改进
通过 using the memcpy trick 而不是可移植的 C++
工会把戏。这样的代码修改应该是微不足道的。
代码应该正确处理次正规。
#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 应使用支持自动矢量化的编译器进行矢量化。下面这一段
代码 vectorizes very well with 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 fast asm instructions 就可以解决这个问题非常酷）：
#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 %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;
}
在数值计算中，经常需要将数字缩放到安全范围内。
例如计算欧氏距离：sqrt(a^2+b^2)。在这里，如果 a 或 b 的大小太 small/large，则可能会发生 underflow/overflow。
解决这个问题的一种常见方法是将数字除以最大幅度的数字。然而，这个解决方案是：
- 慢（除法慢）
- 导致一些额外的不准确
所以我认为与其除以最大幅度的数，不如将它乘以一个接近的 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); }
}
这个函数应该 return 归一化 scaler 和 scalerReciprocal，两者都是 2 的幂数，其中 scaler 接近 value , scalerReciprocal 是 scaler.
scaler/scaleReciprocal 的最大允许指数是 -1022..1022（我不想使用次正规 scaler，因为次正规数可能很慢） .
快速的方法是什么？这可以用纯浮点运算来完成吗？或者我应该从 value 中提取指数并使用简单的 ifs 来执行逻辑？是否有某种技巧可以快速与 (-)1022 进行比较（因为范围是对称的）？
注意：scaler 不需要是最接近的 2 的幂。如果某些逻辑需要它，scaler 可以与最接近的值相差一些小的 2 次幂。
您可以使用
double frexp (double x, int* exp);
返回值是x的小数部分，exp是指数（减去偏移量）。
或者，以下代码获取双精度数的指数部分。
int get_exp(double *d) {
long long *l = (long long *) d;
return ((*l & (0x7ffLL << 52) )>> 52)-1023 ;
}
函数 s = get_scale(z) 计算 "close power of 2"。由于 s 的小数位
为零，s 的倒数只是一个（廉价的）整数减法：参见函数 inv_of_scale。
在 x86 上 get_scale 和 inv_of_scale 使用 clang 编译为非常高效的汇编。
编译器 clang 将三元运算符转换为 minsd 和 maxsd，
另见 Peter Cordes 的 get_scale_x86 和 inv_of_scale_x86），see Godbolt.
请注意 C 明确允许类型双关
通过联合，而 C++ (c++11) 没有这样的权限
虽然 gcc 8.2 和 clang 7.0 并没有抱怨联合，但你可以改进
通过 using the memcpy trick 而不是可移植的 C++
工会把戏。这样的代码修改应该是微不足道的。
代码应该正确处理次正规。
#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 应使用支持自动矢量化的编译器进行矢量化。下面这一段
代码 vectorizes very well with 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 fast asm instructions 就可以解决这个问题非常酷）：
#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 %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;
}