Как я могу точно умножить и разделить 64-битные числа?

Question

У меня есть функция C:

int64_t fn(int64_t a, int32_t b, int32_t c, int32_t d)
{
    /* should return (a * b * c)/d */   
}

Возможно, что значение a будет близко к INT64_MAX, но конечный результат не будет переполнен, например, если b = 1, c = d = 40. Однако у меня возникают проблемы с выяснением, как вычислить это, чтобы я никогда не терял данные к округлению (выполнив сначала деление) или переполнению промежуточного результата.

Если бы у меня был доступ к достаточно большому типу данных, чтобы он соответствовал всем продуктам a, b и c, я просто делал бы математические операции в этом типе и затем усекал их, но есть ли способ, которым я могу сделать это без больших целых чисел?

Answer 1

Написать a = q*d + r с |r| < |d| (Я предполагаю, что d != 0иначе все равно вычисления не имеют смысла). затем (a*b*c)/d = q*b*c + (r*b*c)/d, Если q*b*c переполнение, все вычисления будут переполнены в любом случае, так что либо вас это не волнует, либо вы должны проверить переполнение. r*b*c может все еще переполниться, поэтому мы снова используем тот же метод, чтобы избежать переполнения,

int64_t q = a/d, r = a%d;
int64_t part1 = q*b*c;
int64_t q1 = (r*b)/d, r1 = (r*b)%d;
return part1 + q1*c + (r1*c)/d;

Answer 2

Легко видеть, что некоторые входные данные будут производить выходные данные, которые не могут быть представлены возвращаемым типом int64_t, Например, fn(INT64_MAX, 2, 1, 1), Тем не менее, следующий подход должен позволить вам вернуть правильный ответ для любой комбинации входных данных, которая в действительности соответствует диапазону int64_t,

int64_t fn(int64_t a, int32_t b, int32_t c, int32_t d)
{
    /* find the integer and remainder portions of a/d */
    int64_t leftI = a / d;
    int64_t leftR = a % d;

    /* multiply the integer portion of the result by b and c */
    int64_t resultI = leftI * b * c;

    /* multiply the remainder portion by b */
    int64_t resultR = leftR * b;
    resultI = resultI + (resultR / d) * c;

    /* multiply the remainder portion by c */
    resultR = (resultR % d) * c;

    return resultI + (resultR / d);
}

Answer 3

Я бы предложил найти наибольший общий делитель d для каждого из a, b и c, разделив общие факторы по мере продвижения:

common = gcd(a,d) // You can implement GCD using Euclid's algorithm

a=a/common
d=d/common

common = gcd(b,d)
b=b/common
d=d/common

common = gcd(c,d)
c=c/common
d=d/common

Затем рассчитайте a*b*c/d со всеми общими факторами удалены. Алгоритм Евклида GCD выполняется в логарифмическом времени, поэтому он должен быть достаточно эффективным.

Answer 4

Если вы работаете с x86_64, то asm поддерживает 128-битные целые числа:

int64_t fn(uint64_t a, uint64_t b, uint64_t c, uint64_t d) {

    asm (
        "mulq %1\n"          // a *= b
        "movq %%rbx, %%rdx\n"// rbx = upper 64 bit of the multiplication
        "mulq %2\n"          // multiply the lower 64 bits by c
        "push %%rax\n"       // temporarily save the lowest 64 bits on the stack
        "mov %%rcx, %%rdx\n" // rcx = upper 64 bits of the multiplication
        "movq %%rax, %%rbx\n"// 
        "mulq %2\n"          // multiply the upper 64 bits by c
        "addq %%rax, %%rcx\n"// combine the middle 64 bits
        "addcq %%rdx, $0\n"  // transfer carry tp the higest 64 bits if present
        "divq %3\n"          // divide the upper 128 (of 192) bits by d
        "mov %%rbx, %%rax\n" // rbx = result
        "pop %%rax\n"
        "divq %3\n"          // divide remainder:lower 64 bits by d
        : "+a" (a)           // assigns a to rax register as in/out
        , "+b" (b)           // assigns b to rbx register
        : "g" (c)            // assigns c to random register
        , "g" (d)            // assigns d to random register
        : "edx", "rdx"       // tells the compiler that edx/rdx will be used internally, but does not need any input
    );

    // b now holds the upper 64 bit if (a * b * c / d) > UINT64_MAX
    return a;
}

Обратите внимание, что все входные целые должны быть одинаковой длины. Рабочая длина будет вдвое больше. Работает только с неподписанным.

Родной div а также mul инструкции на x86 работают на двойной длине, чтобы учесть переполнение. К сожалению, я не знаю, какой компилятор присущ им.