c - 将 64 位整数相除，就像被除数左移 64 位一样，没有 128 位类型

Question

为令人困惑的标题道歉。我不知道如何更好地描述我想要完成的事情。我实际上是在尝试 在 C 中为平台获得 64 位乘法的高半部分相反

int64_t divHi64(int64_t dividend, int64_t divisor) {
    return ((__int128)dividend << 64) / (__int128)divisor;
}

由于缺乏对__int128.

Answer 1

这可以在没有多词划分的情况下完成

假设我们想做 ⌊2 ⁶⁴ × ^x ⁄<sub>y⌋ 那么我们可以像这样转换表达式

根据((-y)/y + 1)*x这个问题如何在 C 中计算 2⁶⁴/n？

第二项等价于 (2 ⁶⁴ % y)/y*x 并且有点棘手。我尝试了各种方法，但如果只使用整数运算，都需要 128 位乘法和 128/64 除法。这可以使用MulDiv64(a, b, c) = a*b/c以下问题中的算法来计算

• 在 64 位中进行组合乘除运算的最准确方法是什么？
• 如何在 C++ 中将 64 位整数乘以分数，同时最小化错误？
• (a * b) / c MulDiv 和处理中间乘法的溢出
• 如何准确地乘除 64 位整数？

但是它们可能很慢，如果您有这些函数，则可以更轻松地计算整个表达式，MulDiv64(x, UINT64_MAX, y) + x/y + something而不会弄乱上述转换

如果它具有 64 位或更高的精度，则使用long double似乎是最简单的方法。所以现在可以通过 (2 ⁶⁴ % 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那么现在它有一个 128 位/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，因为如果 x >= y 那么结果会溢出

---

也可以看看

• x86（非 64 位）上两个 128 位整数的高效乘法/除法
• 通过快速浮点倒数高效计算 2**64 / 除数

---

对 16/16 位除法的详尽测试表明，我的解决方案适用于所有情况。double但是，即使精度超过 16 位，您也确实需要float，否则偶尔会返回小于一个的结果。它可以通过在截断之前添加一个epsilon(uint64_t)((mod_y/(long double)y)*x + epsilon)值来修复： 。这意味着如果您不使用epsilon__float128更正结果，您将需要（或-m128bit-long-doubleoption）在 gcc 中获得精确的 64/64 位输出。但是，该类型在 32 位目标上可用，与仅在 64 位目标上支持不同，因此生活会更轻松一些。当然，如果只需要非常接近的结果，您可以按原样使用该函数__int128

下面是我用来验证的代码

#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";
}