最准确的方法是在64位中进行组合乘法除法运算？

Question

对于在32位和64位程序(在Visual C++中)都能工作的64位整数,我能够进行乘法除法运算的最准确方法是什么？(如果溢出,我需要结果mod 2 ^64.)

(我正在寻找类似MulDiv64的东西,除了这个使用内联汇编,它只适用于32位程序.)

显然,可以投射到double后面,但是我想知道是否有更准确的方法并不太复杂.(即我不是在寻找任意精度的算术库!)

Answer 1

由于这是标记为Visual C++,我将提供一个滥用特定于MSVC的内在函数的解决方案.

这个例子相当复杂.它是GMP和java.math.BigInteger大分区使用的相同算法的高度简化版本.

虽然我有一个更简单的算法,它可能慢约30倍.

此解决方案具有以下约束/行为:

• 它需要x64.它不会在x86上编译.
• 商不为零.
• 如果它溢出64位,则商饱和.

请注意,这是针对无符号整数的情况.为此创建一个包装器以使其适用于已签名的案例是微不足道的.此示例还应生成正确的截断结果.

此代码未经过完全测试.但是,它已经通过了我抛出的所有测试用例.
(即使我故意构建以试图打破算法的情况.)

#include <intrin.h>

uint64_t muldiv2(uint64_t a, uint64_t b, uint64_t c){
    //  Normalize divisor
    unsigned long shift;
    _BitScanReverse64(&shift,c);
    shift = 63 - shift;

    c <<= shift;

    //  Multiply
    a = _umul128(a,b,&b);
    if (((b << shift) >> shift) != b){
        cout << "Overflow" << endl;
        return 0xffffffffffffffff;
    }
    b = __shiftleft128(a,b,shift);
    a <<= shift;


    uint32_t div;
    uint32_t q0,q1;
    uint64_t t0,t1;

    //  1st Reduction
    div = (uint32_t)(c >> 32);
    t0 = b / div;
    if (t0 > 0xffffffff)
        t0 = 0xffffffff;
    q1 = (uint32_t)t0;
    while (1){
        t0 = _umul128(c,(uint64_t)q1 << 32,&t1);
        if (t1 < b || (t1 == b && t0 <= a))
            break;
        q1--;
//        cout << "correction 0" << endl;
    }
    b -= t1;
    if (t0 > a) b--;
    a -= t0;

    if (b > 0xffffffff){
        cout << "Overflow" << endl;
        return 0xffffffffffffffff;
    }

    //  2nd reduction
    t0 = ((b << 32) | (a >> 32)) / div;
    if (t0 > 0xffffffff)
        t0 = 0xffffffff;
    q0 = (uint32_t)t0;

    while (1){
        t0 = _umul128(c,q0,&t1);
        if (t1 < b || (t1 == b && t0 <= a))
            break;
        q0--;
//        cout << "correction 1" << endl;
    }

//    //  (a - t0) gives the modulus.
//    a -= t0;

    return ((uint64_t)q1 << 32) | q0;
}

请注意,如果您不需要完美截断的结果,则可以完全删除最后一个循环.如果这样做,答案将不超过正确商数2.

测试用例:

cout << muldiv2(4984198405165151231,6132198419878046132,9156498145135109843) << endl;
cout << muldiv2(11540173641653250113, 10150593219136339683, 13592284235543989460) << endl;
cout << muldiv2(449033535071450778, 3155170653582908051, 4945421831474875872) << endl;
cout << muldiv2(303601908757, 829267376026, 659820219978) << endl;
cout << muldiv2(449033535071450778, 829267376026, 659820219978) << endl;
cout << muldiv2(1234568, 829267376026, 1) << endl;
cout << muldiv2(6991754535226557229, 7798003721120799096, 4923601287520449332) << endl;
cout << muldiv2(9223372036854775808, 2147483648, 18446744073709551615) << endl;
cout << muldiv2(9223372032559808512, 9223372036854775807, 9223372036854775807) << endl;
cout << muldiv2(9223372032559808512, 9223372036854775807, 12) << endl;
cout << muldiv2(18446744073709551615, 18446744073709551615, 9223372036854775808) << endl;

输出:

3337967539561099935
8618095846487663363
286482625873293138
381569328444
564348969767547451
1023786965885666768
11073546515850664288
1073741824
9223372032559808512
Overflow
18446744073709551615
Overflow
18446744073709551615

Answer 2

你只需要64位整数.有一些冗余操作但允许在调试器中使用10作为基础和步骤.

uint64_t const base = 1ULL<<32;
uint64_t const maxdiv = (base-1)*base + (base-1);

uint64_t multdiv(uint64_t a, uint64_t b, uint64_t c)
{
    // First get the easy thing
    uint64_t res = (a/c) * b + (a%c) * (b/c);
    a %= c;
    b %= c;
    // Are we done?
    if (a == 0 || b == 0)
        return res;
    // Is it easy to compute what remain to be added?
    if (c < base)
        return res + (a*b/c);
    // Now 0 < a < c, 0 < b < c, c >= 1ULL
    // Normalize
    uint64_t norm = maxdiv/c;
    c *= norm;
    a *= norm;
    // split into 2 digits
    uint64_t ah = a / base, al = a % base;
    uint64_t bh = b / base, bl = b % base;
    uint64_t ch = c / base, cl = c % base;
    // compute the product
    uint64_t p0 = al*bl;
    uint64_t p1 = p0 / base + al*bh;
    p0 %= base;
    uint64_t p2 = p1 / base + ah*bh;
    p1 = (p1 % base) + ah * bl;
    p2 += p1 / base;
    p1 %= base;
    // p2 holds 2 digits, p1 and p0 one

    // first digit is easy, not null only in case of overflow
    uint64_t q2 = p2 / c;
    p2 = p2 % c;

    // second digit, estimate
    uint64_t q1 = p2 / ch;
    // and now adjust
    uint64_t rhat = p2 % ch;
    // the loop can be unrolled, it will be executed at most twice for
    // even bases -- three times for odd one -- due to the normalisation above
    while (q1 >= base || (rhat < base && q1*cl > rhat*base+p1)) {
        q1--;
        rhat += ch;
    }
    // subtract 
    p1 = ((p2 % base) * base + p1) - q1 * cl;
    p2 = (p2 / base * base + p1 / base) - q1 * ch;
    p1 = p1 % base + (p2 % base) * base;

    // now p1 hold 2 digits, p0 one and p2 is to be ignored
    uint64_t q0 = p1 / ch;
    rhat = p1 % ch;
    while (q0 >= base || (rhat < base && q0*cl > rhat*base+p0)) {
        q0--;
        rhat += ch;
    }
    // we don't need to do the subtraction (needed only to get the remainder,
    // in which case we have to divide it by norm)
    return res + q0 + q1 * base; // + q2 *base*base
}