Вычисление высоких 64 бит 64x64 int продукта в C

Question

Я бы хотел, чтобы моя функция C эффективно вычисляла 64-разрядные 64-разрядные данные из двух 64-битных подписных int. Я знаю, как это сделать в сборке x86-64, с imulq и вытягивая результат из% rdx. Но я не понимаю, как написать это в C вообще, не говоря уже о том, чтобы заставить компилятор сделать это эффективно.

У кого-нибудь есть предложения по написанию этого в C? Это чувствительно к производительности, поэтому «ручные методы» (например, русские крестьяне или библиотеки бигума) отсутствуют.

Этой функция тупой встроенного ассемблера я писал работу и примерно Codegen я после:

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

Общего ответ в том, что x * y может быть разбит на (a + b) * (c + d), где a и c являются высоким порядком частей.

Во-первых, расширить до ac + ad + bc + bd

Теперь вы умножать термины, как 32 битных чисел, хранящихся в long long (или еще лучше, uint64_t), и вы просто помните, что когда вы умножаются большее количество заказов, что вам нужно масштаб на 32 бита. Затем вы делаете добавления, не забывая обнаруживать перенос. Следите за знаками. Естественно, вам нужно делать добавления на куски.

Answer 2

Подождите, у вас есть совершенно хорошее, оптимизированное сборочное решение уже , работающее для этого, и вы хотите выполнить его резервное копирование и попытаться записать его в среде, которая не поддерживает 128-битную математику? Я не буду следовать.

Как вы, очевидно, знаете, эта операция представляет собой одну инструкцию на x86-64. Очевидно, что вы ничего не сделаете, это сделает работу лучше. Если вам действительно нужна портативная C, вам нужно сделать что-то вроде Код DigitalRoss выше и надеяться, что ваш оптимизатор определит, что вы делаете .

Если вам нужна архитектура портативность, но готовы ограничить себя для GCC платформ, есть __int128_t (и __uint128_t) типы в компилятора встроенных функций, которые будут делать то, что вы хотите.

Answer 3

Если вы используете относительно недавнюю GCC на x86_64:

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

В -O1 и выше, это компилирует к тому, что вы хотите:

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

Я считаю, что лязг и VC++ также имеют поддержку типа __int128_t, поэтому это также должно работать на этих платформах, с обычными оговорками о том, чтобы попробовать это самостоятельно.

Answer 4

Что касается вашего монтажного решения, не производите жесткие инструкции mov! Пусть компилятор сделает это за вас. Вот модифицированная версия коды:

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

Полезные ссылки: Machine Constraints

Answer 5

Так как вы сделали очень хорошую работу, решая свои собственные проблемы с машинным кодом, я понял, что вы заслужили помощь с портативной версией.Я бы оставил ifdef в том месте, где вы просто используете сборку, если в gnu на x86.

В любом случае, вот реализация ... Я уверен, что это правильно, но никаких гарантий, я просто ударил это вчера вечером ... вы, вероятно, должны избавиться от статики positive_result [] и result_negative, это только артефакты моего модульного теста ...

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