最准确的方式在64位做一个组合乘法和除法操作？ [英] Most accurate way to do a combined multiply-and-divide operation in 64-bit?

问题描述

什么是最准确的方式，我可以做，在32位和64位程序的工作原理（在Visual C ++）的64位整数的乘法和除法操作？ （在溢出的情况下，我需要的结果MOD ​​2 ⁶⁴ ）

What is the most accurate way I can do a multiply-and-divide operation for 64-bit integers that works in both 32-bit and 64-bit programs (in Visual C++)? (In case of overflow, I need the result mod 2⁶⁴.)

（我正在寻找像 MulDiv64 ，只是此人使用内嵌装配，它只能在32位程序。）

(I'm looking for something like MulDiv64, except that this one uses inline assembly, which only works in 32-bit programs.)

显然，铸造双击和背部是可能的，但我想知道如果有，是不是太复杂，一个更准确的方法。 （即我不是在寻找arbitrary- precision运算库在这里！）

Obviously, casting to double and back is possible, but I'm wondering if there's a more accurate way that isn't too complicated. (i.e. I'm not looking for arbitrary-precision arithmetic libraries here!)

推荐答案

由于这个被标记的Visual C ++我会给一个解决方案，滥用特定MSVC-内部函数。

这个例子是相当复杂的。这是相同的算法的一个高度简化的版本所使用的GMP和 java.math.BigInteger中的为大分工。

This example is fairly complicated. It's a highly simplified version of the same algorithm that is used by GMP and java.math.BigInteger for large division.

虽然我心里有一个简单的算法，它可能是30倍左右慢。

Although I have a simpler algorithm in mind, it's probably about 30x slower.

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

• 它需要64位。它不会编译在x86上。


• 的商不为零。


• 的商如果饱和溢出64位。

请注意，这是为无符号整数的情况下。这是微不足道的围绕这一个包装，使之成为签约的情况下正常工作。这个例子也应该产生正确截断结果。

Note that this is for the unsigned integer case. It's trivial to build a wrapper around this to make it work for signed cases as well. This example should also produce correctly truncated results.

这code未全面测试。然而，它通过了所有，我已经在它抛出的测试用例。结果（即使是我故意构造的情况下试图打破的算法。）

This code is not fully tested. However, it has passed all the tests cases that I've thrown at it.
(Even cases that I've intentionally constructed to try to break the algorithm.)

#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。

Note that if you don't need a perfectly truncated result, you can remove the last loop completely. If you do this, the answer will be no more than 2 larger than the correct quotient.

测试用例：

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

这篇关于最准确的方式在64位做一个组合乘法和除法操作？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！