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

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

Сначала разверните до ac + ad + bc + bd

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

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

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

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

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

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

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

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

Если вы '

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