我怎樣才能精確地乘以64位的整數？

Question

我有一個C函數：

int64_t fn(int64_t a, int32_t b, int32_t c, int32_t d) 
{ 
    /* should return (a * b * c)/d */ 
} 

這是可能的一個是INT64_MAX附近，但對最終結果沒有溢出，例如，如果B = 1，C = d = 40。然而，我我無法弄清楚如何計算這個值，所以我從不會丟失數據來舍入（通過先進行除法）或者使中間結果溢出。

如果我有權訪問足夠大的數據類型以適應a，b和c的整個產品，那麼我只需在該類型中進行數學計算，然後截斷，但是有沒有什麼辦法可以做到這一點整數？

Answer 1

寫a = q*d + r與|r| < |d|（我假設d != 0，否則計算是沒有意義的反正）。然後(a*b*c)/d = q*b*c + (r*b*c)/d。如果q*b*c溢出，則無論如何整個計算都會溢出，因此無論您是不在意，還是必須檢查溢出。 r*b*c仍可能溢出，所以我們再次用同樣的方法，以避免溢出，

int64_t q = a/d, r = a%d; 
int64_t part1 = q*b*c; 
int64_t q1 = (r*b)/d, r1 = (r*b)%d; 
return part1 + q1*c + (r1*c)/d; 

Answer 2

很容易發現某些輸入會產生無法用返回類型int64_t表示的輸出。例如，fn(INT64_MAX, 2, 1, 1)。但是，以下方法 應允許您爲任何確實適合int64_t範圍內的輸入組合返回正確答案。

int64_t fn(int64_t a, int32_t b, int32_t c, int32_t d) 
{ 
    /* find the integer and remainder portions of a/d */ 
    int64_t leftI = a/d; 
    int64_t leftR = a % d; 

    /* multiply the integer portion of the result by b and c */ 
    int64_t resultI = leftI * b * c; 

    /* multiply the remainder portion by b */ 
    int64_t resultR = leftR * b; 
    resultI = resultI + (resultR/d) * c; 

    /* multiply the remainder portion by c */ 
    resultR = (resultR % d) * c; 

    return resultI + (resultR/d); 
} 

Answer 3

如果你在x86_64的工作那麼ASM支持128個整數：

int64_t fn(uint64_t a, uint64_t b, uint64_t c, uint64_t d) { 

    asm (
     "mulq %1\n"   // a *= b 
     "movq %%rbx, %%rdx\n"// rbx = upper 64 bit of the multiplication 
     "mulq %2\n"   // multiply the lower 64 bits by c 
     "push %%rax\n"  // temporarily save the lowest 64 bits on the stack 
     "mov %%rcx, %%rdx\n" // rcx = upper 64 bits of the multiplication 
     "movq %%rax, %%rbx\n"// 
     "mulq %2\n"   // multiply the upper 64 bits by c 
     "addq %%rax, %%rcx\n"// combine the middle 64 bits 
     "addcq %%rdx, $0\n" // transfer carry tp the higest 64 bits if present 
     "divq %3\n"   // divide the upper 128 (of 192) bits by d 
     "mov %%rbx, %%rax\n" // rbx = result 
     "pop %%rax\n" 
     "divq %3\n"   // divide remainder:lower 64 bits by d 
     : "+a" (a)   // assigns a to rax register as in/out 
     , "+b" (b)   // assigns b to rbx register 
     : "g" (c)   // assigns c to random register 
     , "g" (d)   // assigns d to random register 
     : "edx", "rdx"  // tells the compiler that edx/rdx will be used internally, but does not need any input 
    ); 

    // b now holds the upper 64 bit if (a * b * c/d) > UINT64_MAX 
    return a; 
} 

請注意，所有的輸入整數必須是相同的長度。工作長度將是輸入的兩倍。僅與未簽名一起使用。

x86上的原生div和mul指令正確地使用雙精度來允許溢出。令人遺憾的是，我不瞭解使用它們的內在編譯器。

Answer 4

我會建議找d的最大公約數，每個A，B和C的 分割出去的共同因素，當您去：

common = gcd(a,d) // You can implement GCD using Euclid's algorithm 

a=a/common 
d=d/common 

common = gcd(b,d) 
b=b/common 
d=d/common 

common = gcd(c,d) 
c=c/common 
d=d/common 

然後刪除了所有的常見因素計算a*b*c/d 。歐幾里德的GCD算法 運行在對數時間，所以這應該是相當有效的。