Вычисление отношения (64-разрядное целое число без знака) * (64-разрядное целое число без знака), деленное на 2^64

Question

Я хочу вычислить частное от "умножения двух 64-битных целых чисел без знака", деленное на 2^64.

Некоторые люди говорят, что могут использовать 128-битные целые числа, но в некоторых языках программирования или на некоторых платформах (например, в Visual Studio C++) 128-битные целые числа не поддерживаются как встроенные.

Но мы не хотим использовать разделение, потому что разделение занимает слишком много времени. Я думаю, что это можно сделать с помощью сложения / вычитания, умножения и побитовой операции (например, сдвиг битов).

Answer 1

Разбейте свои числа на две части (используя битовые сдвиги и битовые маски) и примените некоторую алгебру.

• Первый номер: A*2^32 + C, где A а также C каждый меньше 2^32.
• Второй номер: B*2^32 + D, где B а также D каждый меньше 2^32.
• (A*2^32 + C) * (B*2^32 + D) = (A*B)*2^64 + (A*D)*2^32 + (B*C)*2^32 + (C*D)
• Разделите на 2^64: (A*B) + (A*D)/2^32 + (B*C)/2^32 + (C*D)/2^64

Так что ответ почти (A*B) + (A*D)>>32 + (B*C)>>32, но это может допустить ошибку округления. В чем ошибка? Вычтите этот почти-ответ из реального (с плавающей запятой) частного:

• (A*D)&0xFFFFFFFF/2^32 + (B*C)&0xFFFFFFFF/2^32 + (C*D)/2^64 (пожалуйста, рассмотрите деления как "реальные" или с плавающей точкой).
• = [(A*D)&0xFFFFFFFF + (B*C)&0xFFFFFFFF + (C*D)/2^32] / 2^32 (опять реальное разделение)
• = [(A*D)&0xFFFFFFFF + (B*C)&0xFFFFFFFF + (C*D)>>32] >> 32 плюс что-то меньше 1.

Таким образом, вы можете получить желаемый номер с(A*B) + (A*D)>>32 + (B*C)>>32 + [(A*D)&0xFFFFFFFF + (B*C)&0xFFFFFFFF + (C*D)>>32] >> 32

Answer 2

Разложите ваши числа a и b на 32-битные части:

a = a1 * 2**32 + a0
b = b1 * 2**32 + b0
a * b = (a1 * b1) * 2**64 + (a1 * b0 + a0 * b1) * 2**32 + a0 * b0

Чтобы получить 64 старших бита результата, часть a1 * b1 очевидна, другая часть должна быть разложена: сначала добавьте 32 старших бита a1 * b0 к 32 старшим битам a0 * b1; во-вторых, добавьте 32 младших бита (a1 * b0 + a0 * b1) к 32 старшим битам a0 * b0 и сохраните 32 старших бита этого промежуточного результата (фактически 1 значащий бит), чтобы учесть переполнение от младшего биты, прежде чем выбросить их.

Ниже приведен код с базовой проверкой результата.

#include <cstdint>
#include <iostream>

using namespace std;

inline uint64_t high32(uint64_t x) {
    return x >> 32;
}

inline uint64_t low32(uint64_t x) {
    return static_cast<uint32_t>(x);
}

uint64_t mul64(uint64_t a, uint64_t b)
{
    uint64_t a1 = high32(a);
    uint64_t a0 = low32(a);
    uint64_t b1 = high32(b);
    uint64_t b0 = low32(b);

    uint64_t a1_b0 = a1 * b0;
    uint64_t a0_b1 = a0 * b1;

    return a1 * b1 + high32(a1_b0) + high32(a0_b1)
         + high32(low32(a1_b0 + a0_b1) + high32(a0 * b0));
}

int main()
{
    cout << mul64(UINT64_MAX, UINT64_MAX) << endl;
    cout << UINT64_MAX - 1 << endl;

    cout << endl;
    cout << mul64(UINT64_MAX - 1, UINT64_MAX) << endl;
    cout << UINT64_MAX - 2 << endl;

    cout << endl;
    cout << mul64(UINT64_MAX - 2, UINT64_MAX) << endl;
    cout << UINT64_MAX - 3 << endl;
}