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

Question

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

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

Эта чертова встроенная функция сборки, которую я написал, работает и является примерно тем кодеком, который мне нужен:

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

Если вы используете относительно недавний 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

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

Answer 2

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

Сначала расширьте до ac + ad + bc + bd

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

Код, реализующий вышеизложенное, см. Мой другой ответ .

Answer 3

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

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

Полезные ссылки: Машинные ограничения

Answer 4

Поскольку вы проделали довольно хорошую работу по решению собственной проблемы с машинным кодом, я подумал, что вы заслужили некоторую помощь с переносной версией. Я бы оставил 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];
}

Answer 5

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

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

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