Быстрый способ получить близкое число степени 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 & scalerReciprocal, оба являются числами степени 2, где scaler близко к value, а scalerReciprocal равно обратной величине scaler.

Максимально допустимые показатели для scaler / scaleReciprocal равны -1022..1022 (я не хочу работать с субнормальными scaler, так как субнормальные числа могут быть медленными).

Каким будет быстрый способ сделать это?Можно ли это сделать с помощью операций с плавающей запятой?Или я должен извлечь показатель степени из value и использовать простые if s для выполнения логики?Есть ли какая-нибудь хитрость для быстрого сравнения с (-) 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, см. также комментарий Питера Кордеса . С gcc это немного эффективнее перевести эти функции в х86 код (get_scale_x86 и inv_of_scale_x86), см. Годболт .

Обратите внимание, что C явно разрешает наказание типа через объединение, тогда как C ++ (c ++ 11) не имеет такого разрешения Хотя gcc 8.2 и clang 7.0 не жалуются на объединение, вы можете улучшить C ++ переносится с использованием 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 и vminpd.

Answer 2

Основываясь на ответе wim, вот еще одно решение, которое может быть быстрее, так как содержит на одну инструкцию меньше. Вывод немного отличается, но все еще отвечает требованиям.

Идея состоит в том, чтобы использовать битовые операции для исправления граничных случаев: поместите 01 в lsb показателя степени, независимо от его значения. Итак, показатель степени:

• 0 становится 1 (-1023 становится -1022)
• 2046 становится 2045 (1023 становится 1022)
• другие показатели также изменены, но незначительно: число может стать в два раза больше по сравнению с решением Вима (когда показатель степени lsb изменяется с 00 на 01), или вдвое (при 10-> 01) или 1 / 4 (при 11-> 01)

Итак, эта модифицированная подпрограмма работает (и я думаю, что это здорово, что проблему можно решить только с помощью 2 инструкций fast 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 %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;
}

Answer 3

Вы можете использовать

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 ;
}