Obliczanie wysokiej 64 bitów 64x64 int produktu w C

Question

Chciałbym mojej funkcji C, aby skutecznie obliczyć wysokie 64 bitów produktu dwóch 64-bitowych podpisanych int. Wiem, jak to zrobić w montażu x86-64, z imulq i wyciągając wynik z% rdx. Ale nie jestem w stanie napisać tego w C, nie mówiąc już o tym, żeby kompilator działał wydajnie.

Czy ktoś ma jakieś sugestie dotyczące pisania tego w C? Jest to wrażliwe na wydajność, więc "metody ręczne" (takie jak rosyjski chłop czy biblioteki bignum) są obecnie niedostępne.

Ta funkcja dorky montaż inline pisałem działa i jest grubsza Codegen Jestem po:

static long mull_hi(long inp1, long inp2) { 
    long output = -1; 
    __asm__("movq %[inp1], %%rax;" 
      "imulq %[inp2];" 
      "movq %%rdx, %[output];" 
      : [output] "=r" (output) 
      : [inp1] "r" (inp1), [inp2] "r" (inp2) 
      :"%rax", "%rdx"); 
    return output; 
} 

Answer 1

Ogólna odpowiedź jest taka, że ​​x * y można rozłożyć na (a + b) * (c + d), gdzie a i c są wysokie zamówienie Części.

pierwsze, rozszerzenie do ac + ad + bc + bd

Teraz pomnożyć warunki jak numery 32 bitowych przechowywanych jako long long (albo jeszcze lepiej, uint64_t), i po prostu pamiętać, że jeśli pomnożyć wyższy numer zamówienia, trzeba skaluj przez 32 bity. Potem robisz dodatki, pamiętając o wykryciu carry. Śledź znak. Oczywiście, musisz dodawać fragmenty.

Answer 2

Czekaj, masz doskonale dobre, zoptymalizowane rozwiązanie montażowe już działa , a chcesz je złożyć i spróbować zapisać w środowisku, które nie obsługuje matematyki 128-bitowej? Nie idę za tobą.

Jak widać, ta operacja jest pojedynczą instrukcją na temat x86-64. Oczywiście nic, co robisz, nie usprawni działania. Jeśli naprawdę potrzebujesz przenośnego C, musisz zrobić coś w stylu powyżej DigitalRoss i mieć nadzieję, że twój optymalizator zorientuje się, co robisz.

Jeśli potrzebujesz architektura przenoszenia, ale są skłonni ograniczać się do platform gcc, istnieje __int128_t (i __uint128_t) typy w kompilatora intrinsics że będą robić to, co chcesz.

Answer 3

Jeśli używasz stosunkowo niedawnej GCC na x86_64:

int64_t mulHi(int64_t x, int64_t y) { 
    return (int64_t)((__int128_t)x*y >> 64); 
} 

Na -O1 i wyższej, to kompiluje się, co chcesz:

_mulHi: 
0000000000000000 movq %rsi,%rax 
0000000000000003 imulq %rdi 
0000000000000006 movq %rdx,%rax 
0000000000000009 ret 

wierzę, że szczęk i VC++ również mają wsparcie dla typu __int128_t, więc powinno to działać również na tych platformach, z typowymi zastrzeżeniami dotyczącymi samodzielnego testowania.

Answer 4

W odniesieniu do swojego rozwiązania montażowego, nie koduj mocno instrukcji mov! Niech kompilator zrobi to za Ciebie. Oto zmodyfikowana wersja kodu:

static long mull_hi(long inp1, long inp2) { 
    long output; 
    __asm__("imulq %2" 
      : "=d" (output) 
      : "a" (inp1), "r" (inp2)); 
    return output; 
} 

pomocna referencyjny: Machine Constraints

Answer 5

Skoro zrobił kawał dobrej roboty rozwiązywania swój problem z kodu maszynowego, pomyślałem zasłużyłeś pomocy w wersji przenośnej.Zostawiłbym ifdef w miejscu, gdzie po prostu użyjesz zestawu, jeśli jesteś w gnu na x86.

W każdym razie, tutaj jest implementacja ... Jestem prawie pewna, że ​​to jest poprawne, ale nie ma gwarancji, właśnie uderzyłem to wczoraj wieczorem ... prawdopodobnie powinieneś pozbyć się statystyk positive_result [] i result_negative, to tylko artefakty z mojego testu jednostkowego ...

#include <stdlib.h> 
#include <stdio.h> 

// stdarg.h doesn't help much here because we need to call llabs() 

typedef unsigned long long uint64_t; 
typedef signed long long int64_t; 

#define B32 0xffffffffUL 

static uint64_t positive_result[2]; // used for testing 
static int result_negative;   // used for testing 

static void mixed(uint64_t *result, uint64_t innerTerm) 
{ 
    // the high part of innerTerm is actually the easy part 

    result[1] += innerTerm >> 32; 

    // the low order a*d might carry out of the low order result 

    uint64_t was = result[0]; 

    result[0] += (innerTerm & B32) << 32; 

    if (result[0] < was) // carry! 
     ++result[1]; 
} 


static uint64_t negate(uint64_t *result) 
{ 
    uint64_t t = result[0] = ~result[0]; 
    result[1] = ~result[1]; 
    if (++result[0] < t) 
    ++result[1]; 
    return result[1]; 
} 

uint64_t higherMul(int64_t sx, int64_t sy) 
{ 
    uint64_t x, y, result[2] = { 0 }, a, b, c, d; 

    x = (uint64_t)llabs(sx); 
    y = (uint64_t)llabs(sy); 

    a = x >> 32; 
    b = x & B32; 
    c = y >> 32; 
    d = y & B32; 

    // the highest and lowest order terms are easy 

    result[1] = a * c; 
    result[0] = b * d; 

    // now have the mixed terms ad + bc to worry about 

    mixed(result, a * d); 
    mixed(result, b * c); 

    // now deal with the sign 

    positive_result[0] = result[0]; 
    positive_result[1] = result[1]; 
    result_negative = sx < 0^sy < 0; 
    return result_negative ? negate(result) : result[1]; 
}