Это можно сделать без деления на несколько слов
Предположим, мы хотим сделать ⌊2 64 × x ⁄ y ⌋тогда мы можем преобразовать выражение следующим образом:
Первый член тривиально выполняется как ((-y)/y + 1)*x согласно этому вопросу Как вычислить 2⁶⁴/ n в C?
Второй член эквивалентен (2 64 % y) / y * x и немного сложнее.Я пробовал разные способы, но всем нужно 128-битное умножение и 128/64 деление, если используются только целочисленные операции.Это можно сделать с помощью алгоритмов для вычисления MulDiv64(a, b, c) = a*b/c в приведенных ниже вопросах
Однако они могут быть медленными, и если у вас есть эти функции, вывычислить целое выражение проще, например MulDiv64(x, UINT64_MAX, y) + x/y + something, не путаясь с приведенным выше преобразованием
Использование long double представляется наиболее простым способом, если оно имеет 64-битную точность или более.Так что теперь это можно сделать с помощью (2 64 % y) / (long double) y * x
uint64_t divHi64(uint64_t x, uint64_t y) {
uint64_t mod_y = UINT64_MAX % y + 1;
uint64_t result = ((-y)/y + 1)*x;
if (mod_y != y)
result += (uint64_t)((mod_y/(long double)y)*x);
return result;
}
Проверка переполнения была опущена для упрощения.Небольшое изменение потребуется, если вам нужно подписанное деление
Если вы нацелены на 64-битную Windows , но вы используете MSVC, у которого нет __int128, тогда теперь он имеет 64-битное внутреннее деление , что значительно упрощает работу без 128-битного целочисленного типа.Тем не менее, вам все равно нужно обрабатывать переполнение, потому что инструкция div вызовет исключение в этом случае
uint64_t divHi64(uint64_t x, uint64_t y) {
uint64_t high, remainder;
uint64_t low = _umul128(UINT64_MAX, y, &high);
if (x <= high /* && 0 <= low */)
return _udiv128(x, 0, y, &remainder);
// overflow case
errno = EOVERFLOW;
return 0;
}
Приведенная выше проверка переполнения может быть упрощена до проверки того, x = y, то результат будет переполнен
См. также
Исчерпывающие тесты 16/16Битовое деление показывает, что мое решение работает правильно для всех случаев.Однако вам нужно double, даже если float имеет точность более 16 бит, в противном случае иногда будет возвращаться результат меньше единицы.Это можно исправить, добавив значение epsilon перед усечением: (uint64_t)((mod_y/(long double)y)*x + epsilon).Это означает, что вам потребуется __float128 (или опция -m128bit-long-double ) в gcc для точного 64/64-битного вывода, если вы не исправите результат с помощью epsilon .Однако этот тип доступен для 32-битных целей , в отличие от __int128, который поддерживается только для 64-битных целей, поэтому жизнь будет немного проще.Конечно, вы можете использовать функцию как есть, если нужен только очень близкий результат
Ниже приведен код, который я использовал для проверки
#include <thread>
#include <iostream>
#include <limits>
#include <climits>
#include <mutex>
std::mutex print_mutex;
#define MAX_THREAD 8
#define NUM_BITS 27
#define CHUNK_SIZE (1ULL << NUM_BITS)
// typedef uint32_t T;
// typedef uint64_t T2;
// typedef double D;
typedef uint64_t T;
typedef unsigned __int128 T2; // the type twice as wide as T
typedef long double D;
// typedef __float128 D;
const D epsilon = 1e-14;
T divHi(T x, T y) {
T mod_y = std::numeric_limits<T>::max() % y + 1;
T result = ((-y)/y + 1)*x;
if (mod_y != y)
result += (T)((mod_y/(D)y)*x + epsilon);
return result;
}
void testdiv(T midpoint)
{
T begin = midpoint - CHUNK_SIZE/2;
T end = midpoint + CHUNK_SIZE/2;
for (T i = begin; i != end; i++)
{
T x = i & ((1 << NUM_BITS/2) - 1);
T y = CHUNK_SIZE/2 - (i >> NUM_BITS/2);
// if (y == 0)
// continue;
auto q1 = divHi(x, y);
T2 q2 = ((T2)x << sizeof(T)*CHAR_BIT)/y;
if (q2 != (T)q2)
{
// std::lock_guard<std::mutex> guard(print_mutex);
// std::cout << "Overflowed: " << x << '&' << y << '\n';
continue;
}
else if (q1 != q2)
{
std::lock_guard<std::mutex> guard(print_mutex);
std::cout << x << '/' << y << ": " << q1 << " != " << (T)q2 << '\n';
}
}
std::lock_guard<std::mutex> guard(print_mutex);
std::cout << "Done testing [" << begin << ", " << end << "]\n";
}
uint16_t divHi16(uint32_t x, uint32_t y) {
uint32_t mod_y = std::numeric_limits<uint16_t>::max() % y + 1;
int result = ((((1U << 16) - y)/y) + 1)*x;
if (mod_y != y)
result += (mod_y/(double)y)*x;
return result;
}
void testdiv16(uint32_t begin, uint32_t end)
{
for (uint32_t i = begin; i != end; i++)
{
uint32_t y = i & 0xFFFF;
if (y == 0)
continue;
uint32_t x = i & 0xFFFF0000;
uint32_t q2 = x/y;
if (q2 > 0xFFFF) // overflowed
continue;
uint16_t q1 = divHi16(x >> 16, y);
if (q1 != q2)
{
std::lock_guard<std::mutex> guard(print_mutex);
std::cout << x << '/' << y << ": " << q1 << " != " << q2 << '\n';
}
}
}
int main()
{
std::thread t[MAX_THREAD];
for (int i = 0; i < MAX_THREAD; i++)
t[i] = std::thread(testdiv, std::numeric_limits<T>::max()/MAX_THREAD*i);
for (int i = 0; i < MAX_THREAD; i++)
t[i].join();
std::thread t2[MAX_THREAD];
constexpr uint32_t length = std::numeric_limits<uint32_t>::max()/MAX_THREAD;
uint32_t begin, end = length;
for (int i = 0; i < MAX_THREAD - 1; i++)
{
begin = end;
end += length;
t2[i] = std::thread(testdiv16, begin, end);
}
t2[MAX_THREAD - 1] = std::thread(testdiv, end, UINT32_MAX);
for (int i = 0; i < MAX_THREAD; i++)
t2[i].join();
std::cout << "Done\n";
}