Могу ли я использовать 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 версия, как объясняет мистик в комментариях.
3 ответа
Да, это возможно. Но что касается AVX2, он вряд ли будет лучше скалярных подходов с MULX/ADCX/ADOX.
Существует практически неограниченное количество вариантов этого подхода для разных доменов ввода / вывода. Я расскажу только о 3 из них, но их легко обобщить, если вы знаете, как они работают.
Отказ от ответственности:
- Все решения здесь предполагают, что режим округления округлен до четного.
- Использование флагов быстрой математической оптимизации не рекомендуется, поскольку эти решения основаны на строгом стандарте IEEE.
Подписано в парном разряде: [-251, 251]
// 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 будет совершенно не так.
Целые числа со знаком в диапазоне: [-251, 251]
// 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));
}
Основываясь на первом примере, мы объединяем его с обобщенной версией double <-> int64 конверсионный трюк.
Этот более полезен, так как вы работаете с целыми числами. Но даже с помощью быстрого конверсионного трюка, большая часть времени будет потрачена на конверсию. К счастью, вы можете исключить некоторые входные преобразования, если вы умножаете на один и тот же операнд несколько раз.
Целые числа без знака в диапазоне: [0, 252)
// 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 вы можете переопределить округление для этой инструкции и оставить режим округления в покое.
Последние мысли:
Стоит отметить, что диапазон действия 252 можно уменьшить, отрегулировав магические константы. Это может быть полезно для первого решения (с плавающей запятой), поскольку оно дает вам дополнительную мантиссу для использования для накопления с плавающей запятой. Это позволяет вам обходиться без необходимости постоянно конвертировать туда и обратно между int64 и double, как в последних 2 решениях.
Хотя 3 примера здесь вряд ли будут лучше скалярных методов, AVX512 почти наверняка перевесит баланс. В частности, Knights Landing имеет плохую пропускную способность для ADCX и ADOX.
И, конечно, все это спорно, когда AVX512-IFMA выходит. Это уменьшает полный 52 x 52 -> 104-bit Продукт по 2 инструкции и дает накопление бесплатно.
Один из способов сделать целочисленную арифметику из нескольких слов - использовать двойную арифметику. Давайте начнем с некоторого двойного-двойного кода умножения
#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).
Соблазнительно сказать, что проблема заключается в переключении обратно на целочисленную область. Но зачем это нужно? Вы можете делать все в домене с плавающей запятой. Давайте посмотрим на некоторые примеры. Я считаю, что с юнит-тестом легче 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 бит (из двойного-двойного). 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);
}
Что ж, вы, безусловно, можете выполнять операции FP-lane с целыми числами. И они всегда будут точными: хотя существуют инструкции SSE, которые не гарантируют надлежащую точность и округление IEEE-754, все они без исключения не имеют целочисленного диапазона, так что вы все равно не смотрите. Итог: сложение / вычитание / умножение всегда будет точным в целочисленной области, даже если вы делаете их на упакованных числах с плавающей запятой.
Что касается плавающих чисел с квадратичной точностью (>52-битная мантисса), нет, они не поддерживаются и, скорее всего, не появятся в обозримом будущем. Просто не так много их зовут. Они появились в нескольких архитектурах рабочих станций эпохи SPARC, но, честно говоря, они были просто поводом для неполного понимания разработчиками того, как писать численно устойчивые алгоритмы, и со временем они исчезли.
Операции с широкими целыми числами оказываются неподходящими для SSE. Я действительно пытался использовать это недавно, когда внедрял большую целочисленную библиотеку, и, честно говоря, это не помогло мне. x86 был разработан для многословной арифметики; Вы можете увидеть это в таких операциях, как ADC (который производит и потребляет бит переноса) и IDIV (который позволяет делителю быть в два раза шире, чем дивиденд, если частное не шире, чем дивиденд, ограничение, которое делает его бесполезно ни для чего, кроме деления нескольких слов). Но многословная арифметика по своей природе последовательна, а SSE по своей природе параллельна. Если вам повезло, что в ваших числах достаточно битов, чтобы поместиться в мантиссу FP, поздравляю. Но если у вас большие целые числа, SSE, вероятно, не будет вашим другом.