一种快速的方法来获得一个接近的2次方数（浮点）

Question

例如，计算欧氏距离: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和Scaler互惠，两者都是2的幂数，其中Scaler接近于value，而Scaler互惠是Scaler的倒数。

scaler/scalereciprocal允许的最大指数为-1022..1022（我不想使用次常规的scaler，因为次常规的数字可能很慢）。

什么是一个快速的方法来做到这一点？这能用纯浮点运算完成吗？或者我应该从value中提取指数，并使用简单的if来执行逻辑？有什么方法可以快速与（-）1022进行比较吗（因为范围是对称的）？

注意:scaler不需要尽可能接近2的幂。如果某些逻辑需要它，scaler可以离最接近的值有一个很小的2次方。

Answer 1

函数s=get_scale(z)计算“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)的效率略高一些，请参阅Godbolt。

请注意，C明确允许通过联合使用类型双关，而C++(C++11)没有这种权限尽管GCC8.2和CLANG7.0没有抱怨联合，但您可以通过使用memcpy技巧而不是联合技巧来方便地改进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应使用支持自动向量化的编译器向量化。下面的代码很好地利用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]);
    }
}