Подтвердить что ты не робот

Могу ли я использовать модули AVX FMA для выполнения битовых точных умножений на 52 бит?

AXV2 не имеет целочисленных умножений с источниками, превышающими 32 бит. Он предлагает 32 x 32 → 32, а также 32 x 32 → 64 умножает 1 но ничего с 64-разрядными источниками.

Скажем, мне нужно, чтобы без знака умножалось с входами более 32 бит, но меньше или равно 52 бит - могу ли я просто использовать плавучую точку DP умножить или FMA, и будет ли бит быть точным, когда целые входы и результаты могут быть представлены в 52 или менее битах (т.е. В диапазоне [0, 2 ^ 52-1])?

Как насчет более общего случая, когда я хочу все 104 бита продукта? Или случай, когда целочисленный продукт занимает более 52 бит (т.е. Продукт имеет ненулевые значения в битовых индексах > 52) - но я хочу только низкие 52 бита? В этом последнем случае MUL собирается дать мне более высокие биты и округлить некоторые из младших бит (возможно, с чем помогает IFMA?).

Изменить: на самом деле, возможно, он мог бы сделать что угодно до 2 ^ 53, на основе этого ответа - я забыл, что подразумеваемый ведущий 1, прежде чем мантисса эффективно даст вам другой бит.

1 Интересно, что операция 64-разрядного продукта PMULDQ имеет половину задержки и вдвое большую пропускную способность 32-разрядной версии PMULLD, так как Mystical объясняет в комментариях.

4b9b3361

ОТВЕТЫ

Ответ 1

Да, это возможно. Но с AVX2 он вряд ли будет лучше скалярных подходов с MULX/ADCX/ADOX.

Существует практически неограниченное количество вариантов этого подхода для разных областей ввода/вывода. Я расскажу только о трех из них, но их легко обобщить, как только вы знаете, как они работают.

Отказ от ответственности:

  • Все решения здесь предполагают, что режим округления округлен до четности.
  • Использование флагов быстрой математической оптимизации не рекомендуется, так как эти решения зависят от строгой IEEE.

Подписанные удвоения в диапазоне: [-2 51 2 51]

//  A*B = L + H*2^52
//  Input:  A and B are in the range [-2^51, 2^51]
//  Output: L and H are in the range [-2^51, 2^51]
void mul52_signed(__m256d& L, __m256d& H, __m256d A, __m256d B){
    const __m256d ROUND = _mm256_set1_pd(30423614405477505635920876929024.);    //  3 * 2^103
    const __m256d SCALE = _mm256_set1_pd(1. / 4503599627370496);                //  1 / 2^52

    //  Multiply and add normalization constant. This forces the multiply
    //  to be rounded to the correct number of bits.
    H = _mm256_fmadd_pd(A, B, ROUND);

    //  Undo the normalization.
    H = _mm256_sub_pd(H, ROUND);

    //  Recover the bottom half of the product.
    L = _mm256_fmsub_pd(A, B, H);

    //  Correct the scaling of H.
    H = _mm256_mul_pd(H, SCALE);
}

Это самый простой и единственный, который конкурирует со скалярными подходами. Окончательное масштабирование необязательно в зависимости от того, что вы хотите делать с выходами. Таким образом, это можно рассматривать только 3 инструкции. Но это также наименее полезно, поскольку оба входа и выхода являются значениями с плавающей запятой.

Абсолютно важно, чтобы оба FMA оставались сплавленными. И именно здесь быстрые математические оптимизации могут нарушать вещи. Если первый FMA разбит, то L больше не будет находиться в диапазоне [-2^51, 2^51]. Если вторая FMA будет разбита, L будет полностью неправильной.

Подписанные целые числа в диапазоне: [-2 51 2 51]

//  A*B = L + H*2^52
//  Input:  A and B are in the range [-2^51, 2^51]
//  Output: L and H are in the range [-2^51, 2^51]
void mul52_signed(__m256i& L, __m256i& H, __m256i A, __m256i B){
    const __m256d CONVERT_U = _mm256_set1_pd(6755399441055744);     //  3*2^51
    const __m256d CONVERT_D = _mm256_set1_pd(1.5);

    __m256d l, h, a, b;

    //  Convert to double
    A = _mm256_add_epi64(A, _mm256_castpd_si256(CONVERT_U));
    B = _mm256_add_epi64(B, _mm256_castpd_si256(CONVERT_D));
    a = _mm256_sub_pd(_mm256_castsi256_pd(A), CONVERT_U);
    b = _mm256_sub_pd(_mm256_castsi256_pd(B), CONVERT_D);

    //  Get top half. Convert H to int64.
    h = _mm256_fmadd_pd(a, b, CONVERT_U);
    H = _mm256_sub_epi64(_mm256_castpd_si256(h), _mm256_castpd_si256(CONVERT_U));

    //  Undo the normalization.
    h = _mm256_sub_pd(h, CONVERT_U);

    //  Recover bottom half.
    l = _mm256_fmsub_pd(a, b, h);

    //  Convert L to int64
    l = _mm256_add_pd(l, CONVERT_D);
    L = _mm256_sub_epi64(_mm256_castpd_si256(l), _mm256_castpd_si256(CONVERT_D));
}

В первом примере мы объединим его с обобщенной версией тэга fast double <-> int64.

Это более полезно, так как вы работаете с целыми числами. Но даже с быстрым конверсионным трюком большая часть времени будет потрачена на конверсии. К счастью, вы можете устранить некоторые из входных преобразований, если вы умножаете один и тот же операнд несколько раз.

Целые числа без знака в диапазоне: [0, 2 52)

//  A*B = L + H*2^52
//  Input:  A and B are in the range [0, 2^52)
//  Output: L and H are in the range [0, 2^52)
void mul52_unsigned(__m256i& L, __m256i& H, __m256i A, __m256i B){
    const __m256d CONVERT_U = _mm256_set1_pd(4503599627370496);     //  2^52
    const __m256d CONVERT_D = _mm256_set1_pd(1);
    const __m256d CONVERT_S = _mm256_set1_pd(1.5);

    __m256d l, h, a, b;

    //  Convert to double
    A = _mm256_or_si256(A, _mm256_castpd_si256(CONVERT_U));
    B = _mm256_or_si256(B, _mm256_castpd_si256(CONVERT_D));
    a = _mm256_sub_pd(_mm256_castsi256_pd(A), CONVERT_U);
    b = _mm256_sub_pd(_mm256_castsi256_pd(B), CONVERT_D);

    //  Get top half. Convert H to int64.
    h = _mm256_fmadd_pd(a, b, CONVERT_U);
    H = _mm256_xor_si256(_mm256_castpd_si256(h), _mm256_castpd_si256(CONVERT_U));

    //  Undo the normalization.
    h = _mm256_sub_pd(h, CONVERT_U);

    //  Recover bottom half.
    l = _mm256_fmsub_pd(a, b, h);

    //  Convert L to int64
    l = _mm256_add_pd(l, CONVERT_S);
    L = _mm256_sub_epi64(_mm256_castpd_si256(l), _mm256_castpd_si256(CONVERT_S));

    //  Make Correction
    H = _mm256_sub_epi64(H, _mm256_srli_epi64(L, 63));
    L = _mm256_and_si256(L, _mm256_set1_epi64x(0x000fffffffffffff));
}

Наконец, мы получим ответ на исходный вопрос. Это завершает работу целочисленного решения со знаком, регулируя конверсии и добавляя шаг исправления.

Но на данный момент мы находимся в 13 инструкциях, половина из которых содержит инструкции с высокой задержкой, не считая многочисленных задержек байпаса FP <-> int. Поэтому маловероятно, что это приведет к выигрышу в любых тестах. Для сравнения, умножение 64 x 64 -> 128-bit SIMD может быть выполнено в 16 инструкциях (14, если вы предварительно обрабатываете входы.)

Шаг коррекции может быть опущен, если режим округления округлен или округлен до нуля. Единственная инструкция, где это имеет значение h = _mm256_fmadd_pd(a, b, CONVERT_U);. Таким образом, на AVX512 вы можете переопределить округление для этой команды и оставить режим округления самостоятельно.

Заключительные мысли:

Стоит отметить, что диапазон операций 2 52 может быть уменьшен путем корректировки магических констант. Это может быть полезно для первого решения (с плавающей точкой), поскольку оно дает вам дополнительную мантиссу для накопления в плавающей точке. Это позволяет обойти необходимость постоянно конвертировать назад и вперед между int64 и double, как в последних двух решениях.

В то время как три примера здесь вряд ли будут лучше скалярных методов, AVX512 почти наверняка опрокинет баланс. Рыцари Посадка, в частности, имеет низкую пропускную способность для ADCX и ADOX.

И, конечно, все это спорно, когда AVX512-IFMA выходит. Это уменьшает полный 52 x 52 -> 104-bit продукт до 2 инструкций и дает накопление бесплатно.

Ответ 2

Один из способов сделать многословную целочисленную арифметику - double-double_arithmetic. Начнем с кода с двойным двойным умножением

#include <math.h>
typedef struct {
  double hi;
  double lo;
} doubledouble;

static doubledouble quick_two_sum(double a, double b) {
  double s = a + b;
  double e = b - (s - a);
  return (doubledouble){s, e};
}

static doubledouble two_prod(double a, double b) {
  double p = a*b;
  double e = fma(a, b, -p);
  return (doubledouble){p, e};
}

doubledouble df64_mul(doubledouble a, doubledouble b) {
  doubledouble p = two_prod(a.hi, b.hi);
  p.lo += a.hi*b.lo;
  p.lo += a.lo*b.hi;
  return quick_two_sum(p.hi, p.lo);
}

Функция two_prod может выполнять целое число 53bx53b → 106b в двух инструкциях. Функция df64_mul может выполнять целое число 106bx106b → 106b.

Сравним это с целым числом 128bx128b → 128b с целым аппаратным обеспечением.

__int128 mul128(__int128 a, __int128 b) {
  return a*b;
}

Сборка для mul128

imul    rsi, rdx
mov     rax, rdi
imul    rcx, rdi
mul     rdx
add     rcx, rsi
add     rdx, rcx

Сборка для df64_mul (скомпилирована с gcc -O3 -S i128.c -masm=intel -mfma -ffp-contract=off)

vmulsd      xmm4, xmm0, xmm2
vmulsd      xmm3, xmm0, xmm3
vmulsd      xmm1, xmm2, xmm1
vfmsub132sd xmm0, xmm4, xmm2
vaddsd      xmm3, xmm3, xmm0
vaddsd      xmm1, xmm3, xmm1
vaddsd      xmm0, xmm1, xmm4
vsubsd      xmm4, xmm0, xmm4
vsubsd      xmm1, xmm1, xmm4

mul128 выполняет три скалярных умножения и два скалярных дополнения/вычитания, тогда как df64_mul выполняет 3 SIMD-умножения, 1 SIMD FMA и 5 SIMD-добавлений/вычитаний. Я не профилировал эти методы, но мне не кажется необоснованным, что df64_mul может превзойти mul128 с использованием 4-парного числа на регистр AVX (изменить sd на pd и xmm на ymm).

Загадочно сказать, что проблема заключается в возвращении к целочисленному домену. Но почему это необходимо? Вы можете делать все в домене с плавающей точкой. Давайте посмотрим на некоторые примеры. Мне легче unit test с float, чем с double.

doublefloat two_prod(float a, float b) {
  float p = a*b;
  float e = fma(a, b, -p);
  return (doublefloat){p, e};
}

//3202129*4807935=15395628093615
x = two_prod(3202129,4807935)
int64_t hi = p, lo = e, s = hi+lo
//p = 1.53956280e+13, e = 1.02575000e+05  
//hi = 15395627991040, lo = 102575, s = 15395628093615

//1450779*1501672=2178594202488
y = two_prod(1450779, 1501672)
int64_t hi = p, lo = e, s = hi+lo 
//p = 2.17859424e+12, e = -4.00720000e+04
//hi = 2178594242560 lo = -40072, s = 2178594202488

Итак, мы заканчиваем разными диапазонами, а во втором случае ошибка (e) является даже отрицательной, но сумма по-прежнему верна. Мы могли бы даже добавить два значения doublefloat x и y вместе (как только мы знаем, как делать двойное двойное добавление - см. Код в конце) и получите 15395628093615+2178594202488. Нет необходимости нормализовать результаты.

Но добавление вызывает основную проблему с двойной двойной арифметикой. А именно, сложение/вычитание происходит медленно, например. 128b + 128b → 128b требуется не менее 11 с добавлением с плавающей запятой, тогда как с целыми числами требуется только два (add и adc).

Итак, если алгоритм тяжелый при умножении, но свет на добавлении, то выполнение многословных целых операций с двойным двойным может выиграть.

В качестве дополнительной заметки язык C является достаточно гибким, чтобы обеспечить реализацию, где целые числа реализованы полностью через аппаратные средства с плавающей точкой. int может быть 24-битным (с одной плавающей точкой), long может быть 54 бит. (из двойной с плавающей запятой), а long long может быть 106-битным (из double-double). C даже не требует двух комплиментов, и поэтому целые числа могут использовать знаковое значение для отрицательных чисел, как обычно, с плавающей точкой.

Здесь работает код C с двойным удвоением и добавлением (я не реализовал деление или другие операции, такие как sqrt, но есть документы, показывающие, как это сделать), если кто-то захочет поиграть с ним. Было бы интересно узнать, может ли это быть оптимизировано для целых чисел.

//if compiling with -mfma you must also use -ffp-contract=off
//float-float is easier to debug. If you want double-double replace
//all float words with double and fmaf with fma 
#include <stdio.h>
#include <math.h>
#include <inttypes.h>
#include <x86intrin.h>
#include <stdlib.h>

//#include <float.h>

typedef struct {
  float hi;
  float lo;
} doublefloat;

typedef union {
  float f;
  int i;
  struct {
    unsigned mantisa : 23;
    unsigned exponent: 8;
    unsigned sign: 1;
  };
} float_cast;

void print_float(float_cast a) {
  printf("%.8e, 0x%x, mantisa 0x%x, exponent 0x%x, expondent-127 %d, sign %u\n", a.f, a.i, a.mantisa, a.exponent, a.exponent-127, a.sign);
}

void print_doublefloat(doublefloat a) {
  float_cast hi = {a.hi};
  float_cast lo = {a.lo};
  printf("hi: "); print_float(hi);
  printf("lo: "); print_float(lo);
}

doublefloat quick_two_sum(float a, float b) {
  float s = a + b;
  float e = b - (s - a);
  return (doublefloat){s, e};
  // 3 add
}

doublefloat two_sum(float a, float b) {
  float s = a + b;
  float v = s - a;
  float e = (a - (s - v)) + (b - v);
  return (doublefloat){s, e};
  // 6 add 
}

doublefloat df64_add(doublefloat a, doublefloat b) {
  doublefloat s, t;
  s = two_sum(a.hi, b.hi);
  t = two_sum(a.lo, b.lo);
  s.lo += t.hi;
  s = quick_two_sum(s.hi, s.lo);
  s.lo += t.lo;
  s = quick_two_sum(s.hi, s.lo);
  return s;
  // 2*two_sum, 2 add, 2*quick_two_sum = 2*6 + 2 + 2*3 = 20 add
}

doublefloat split(float a) {
  //#define SPLITTER (1<<27) + 1
#define SPLITTER (1<<12) + 1
  float t = (SPLITTER)*a;
  float hi = t - (t - a);
  float lo = a - hi;
  return (doublefloat){hi, lo};
  // 1 mul, 3 add
}

doublefloat split_sse(float a) {
  __m128 k = _mm_set1_ps(4097.0f);
  __m128 a4 = _mm_set1_ps(a);
  __m128 t = _mm_mul_ps(k,a4);
  __m128 hi4 = _mm_sub_ps(t,_mm_sub_ps(t, a4));
  __m128 lo4 = _mm_sub_ps(a4, hi4);
  float tmp[4];
  _mm_storeu_ps(tmp, hi4);
  float hi = tmp[0];
  _mm_storeu_ps(tmp, lo4);
  float lo = tmp[0];
  return (doublefloat){hi,lo};

}

float mult_sub(float a, float b, float c) {
  doublefloat as = split(a), bs = split(b);
  //print_doublefloat(as);
  //print_doublefloat(bs);
  return ((as.hi*bs.hi - c) + as.hi*bs.lo + as.lo*bs.hi) + as.lo*bs.lo;
  // 4 mul, 4 add, 2 split = 6 mul, 10 add
}

doublefloat two_prod(float a, float b) {
  float p = a*b;
  float e = mult_sub(a, b, p);
  return (doublefloat){p, e};
  // 1 mul, one mult_sub
  // 7 mul, 10 add
}

float mult_sub2(float a, float b, float c) {
  doublefloat as = split(a);
  return ((as.hi*as.hi -c ) + 2*as.hi*as.lo) + as.lo*as.lo;
}

doublefloat two_sqr(float a) {
  float p = a*a;
  float e = mult_sub2(a, a, p);
  return (doublefloat){p, e};
}

doublefloat df64_mul(doublefloat a, doublefloat b) {
  doublefloat p = two_prod(a.hi, b.hi);
  p.lo += a.hi*b.lo;
  p.lo += a.lo*b.hi;
  return quick_two_sum(p.hi, p.lo);
  //two_prod, 2 add, 2mul, 1 quick_two_sum = 9 mul, 15 add 
  //or 1 mul, 1 fma, 2add 2mul, 1 quick_two_sum = 3 mul, 1 fma, 5 add
}

doublefloat df64_sqr(doublefloat a) {
  doublefloat p = two_sqr(a.hi);
  p.lo += 2*a.hi*a.lo;
  return quick_two_sum(p.hi, p.lo);
}

int float2int(float a) {
  int M = 0xc00000; //1100 0000 0000 0000 0000 0000
  a += M;
  float_cast x;
  x.f = a;
  return x.i - 0x4b400000;
}

doublefloat add22(doublefloat a, doublefloat b) {
  float r = a.hi + b.hi;
  float s = fabsf(a.hi) > fabsf(b.hi) ?
    (((a.hi - r) + b.hi) + b.lo ) + a.lo :
    (((b.hi - r) + a.hi) + a.lo ) + b.lo;
  return two_sum(r, s);  
  //11 add 
}

int main(void) {
  //print_float((float_cast){1.0f});
  //print_float((float_cast){-2.0f});
  //print_float((float_cast){0.0f});
  //print_float((float_cast){3.14159f});
  //print_float((float_cast){1.5f});
  //print_float((float_cast){3.0f});
  //print_float((float_cast){7.0f});
  //print_float((float_cast){15.0f});
  //print_float((float_cast){31.0f});

  //uint64_t t = 0xffffff;
  //print_float((float_cast){1.0f*t});
  //printf("%" PRId64 " %" PRIx64 "\n", t*t,t*t);

  /*
    float_cast t1;
    t1.mantisa = 0x7fffff;
    t1.exponent = 0xfe;
    t1.sign = 0;
    print_float(t1);
  */
  //doublefloat z = two_prod(1.0f*t, 1.0f*t);
  //print_doublefloat(z);
  //double z2 = (double)z.hi + (double)z.lo;
  //printf("%.16e\n", z2);
  doublefloat s = {0};
  int64_t si = 0;
  for(int i=0; i<100000; i++) {
    int ai = rand()%0x800, bi = rand()%0x800000;
    float a = ai, b = bi;
    doublefloat z = two_prod(a,b);
    int64_t zi = (int64_t)ai*bi;
    //print_doublefloat(z);
    //s = df64_add(s,z);
    s = add22(s,z);
    si += zi;
    print_doublefloat(z);
    printf("%d %d ", ai,bi);
    int64_t h = z.hi;
    int64_t l = z.lo;
    int64_t t = h+l;
    //if(t != zi) printf("%" PRId64 " %" PRId64 "\n", h, l);

    printf("%" PRId64 " %" PRId64 " %" PRId64 " %" PRId64 "\n", zi, h, l, h+l);

    h = s.hi;
    l = s.lo;
    t = h + l;
    //if(si != t) printf("%" PRId64 " %" PRId64 "\n", h, l);

    if(si > (1LL<<48)) {
      printf("overflow after %d iterations\n", i); break;
    }
  }

  print_doublefloat(s);
  printf("%" PRId64 "\n", si);
  int64_t x = s.hi;
  int64_t y = s.lo;
  int64_t z = x+y;
  //int hi = float2int(s.hi);
  printf("%" PRId64 " %" PRId64 " %" PRId64 "\n", z,x,y);
}

Ответ 3

Ну, вы, конечно же, можете выполнять операции FP-lane над целыми числами. И они всегда будут точными: хотя существуют инструкции SSE, которые не гарантируют правильную точность и округление IEEE-754, без исключения они являются теми, у которых нет целочисленного диапазона, поэтому не те, которые вы смотрите в любом случае. Итог: добавление/вычитание/умножение всегда будут точными в целочисленном домене, даже если вы делаете их на упакованных поплавках.

Что касается поплавков с четырьмя точками ( > 52 бит мантиссы), нет, они не поддерживаются и, вероятно, не будут в обозримом будущем. Просто не так много призыв к ним. Они появились в нескольких архитектурах рабочих станций SPARC, но, честно говоря, они были всего лишь повязкой над неполным пониманием разработчиками того, как писать численно устойчивые алгоритмы, и со временем они исчезли.

Широкоцелевые операции оказываются очень плохими для SSE. Я действительно пытался использовать его в последнее время, когда я реализовал библиотеку большого числа, и, честно говоря, это не помогло мне. x86 был разработан для многословной арифметики; вы можете увидеть его в таких операциях, как АЦП (который производит и потребляет бит переноса) и IDIV (что позволяет делителю быть вдвое шире, чем дивиденд, если коэффициент не шире дивиденда, ограничение, которое делает это бесполезно ни для чего, кроме разделения многословных). Но многословная арифметика по своей природе последовательна, а SSE по своей природе параллельна. Если вам повезло, что у ваших номеров достаточно бит, чтобы вписаться в мантису FP, поздравляю. Но если у вас большие целые числа, SSE, вероятно, не будет вашим другом.