diff options
| -rw-r--r-- | lib/math/div64.c | 108 |
1 files changed, 65 insertions, 43 deletions
diff --git a/lib/math/div64.c b/lib/math/div64.c index 191761b1b623..b7fc75246399 100644 --- a/lib/math/div64.c +++ b/lib/math/div64.c @@ -186,55 +186,77 @@ EXPORT_SYMBOL(iter_div_u64_rem); #ifndef mul_u64_u64_div_u64 u64 mul_u64_u64_div_u64(u64 a, u64 b, u64 c) { - u64 res = 0, div, rem; - int shift; + if (ilog2(a) + ilog2(b) <= 62) + return div64_u64(a * b, c); - /* can a * b overflow ? */ - if (ilog2(a) + ilog2(b) > 62) { - /* - * Note that the algorithm after the if block below might lose - * some precision and the result is more exact for b > a. So - * exchange a and b if a is bigger than b. - * - * For example with a = 43980465100800, b = 100000000, c = 1000000000 - * the below calculation doesn't modify b at all because div == 0 - * and then shift becomes 45 + 26 - 62 = 9 and so the result - * becomes 4398035251080. However with a and b swapped the exact - * result is calculated (i.e. 4398046510080). - */ - if (a > b) - swap(a, b); +#if defined(__SIZEOF_INT128__) + + /* native 64x64=128 bits multiplication */ + u128 prod = (u128)a * b; + u64 n_lo = prod, n_hi = prod >> 64; + +#else + + /* perform a 64x64=128 bits multiplication manually */ + u32 a_lo = a, a_hi = a >> 32, b_lo = b, b_hi = b >> 32; + u64 x, y, z; + + x = (u64)a_lo * b_lo; + y = (u64)a_lo * b_hi + (u32)(x >> 32); + z = (u64)a_hi * b_hi + (u32)(y >> 32); + y = (u64)a_hi * b_lo + (u32)y; + z += (u32)(y >> 32); + x = (y << 32) + (u32)x; + + u64 n_lo = x, n_hi = z; + +#endif + + int shift = __builtin_ctzll(c); + /* try reducing the fraction in case the dividend becomes <= 64 bits */ + if ((n_hi >> shift) == 0) { + u64 n = (n_lo >> shift) | (n_hi << (64 - shift)); + + return div64_u64(n, c >> shift); /* - * (b * a) / c is equal to - * - * (b / c) * a + - * (b % c) * a / c - * - * if nothing overflows. Can the 1st multiplication - * overflow? Yes, but we do not care: this can only - * happen if the end result can't fit in u64 anyway. - * - * So the code below does - * - * res = (b / c) * a; - * b = b % c; + * The remainder value if needed would be: + * res = div64_u64_rem(n, c >> shift, &rem); + * rem = (rem << shift) + (n_lo - (n << shift)); */ - div = div64_u64_rem(b, c, &rem); - res = div * a; - b = rem; - - shift = ilog2(a) + ilog2(b) - 62; - if (shift > 0) { - /* drop precision */ - b >>= shift; - c >>= shift; - if (!c) - return res; - } } - return res + div64_u64(a * b, c); + if (n_hi >= c) { + /* overflow: result is unrepresentable in a u64 */ + return -1; + } + + /* Do the full 128 by 64 bits division */ + + shift = __builtin_clzll(c); + c <<= shift; + + int p = 64 + shift; + u64 res = 0; + bool carry; + + do { + carry = n_hi >> 63; + shift = carry ? 1 : __builtin_clzll(n_hi); + if (p < shift) + break; + p -= shift; + n_hi <<= shift; + n_hi |= n_lo >> (64 - shift); + n_lo <<= shift; + if (carry || (n_hi >= c)) { + n_hi -= c; + res |= 1ULL << p; + } + } while (n_hi); + /* The remainder value if needed would be n_hi << p */ + + return res; } EXPORT_SYMBOL(mul_u64_u64_div_u64); #endif |