I am trying to multiply two uint64_t
s and store the result to uint64_t
. I found an existing answer on Stackoverflow which splits the inputs in to their four uint32_t
s and joins the result later:
https://stackoverflow.com/a/28904636/1107474
I have created a full example using the code and pasted it below.
However, for 37 x 5
I am getting the result 0
instead of 185
?
#include <iostream>
int main()
{
uint64_t a = 37; // Input 1
uint64_t b = 5; // Input 2
uint64_t a_lo = (uint32_t)a;
uint64_t a_hi = a >> 32;
uint64_t b_lo = (uint32_t)b;
uint64_t b_hi = b >> 32;
uint64_t a_x_b_hi = a_hi * b_hi;
uint64_t a_x_b_mid = a_hi * b_lo;
uint64_t b_x_a_mid = b_hi * a_lo;
uint64_t a_x_b_lo = a_lo * b_lo;
uint64_t carry_bit = ((uint64_t)(uint32_t)a_x_b_mid
(uint64_t)(uint32_t)b_x_a_mid
(a_x_b_lo >> 32) ) >> 32;
uint64_t multhi = a_x_b_hi
(a_x_b_mid >> 32) (b_x_a_mid >> 32)
carry_bit;
std::cout << multhi << std::endl; // Outputs 0 instead of 185?
}
CodePudding user response:
I'm merging your code with another answer in the original link.
#include <iostream>
int main()
{
uint64_t a = 37; // Input 1
uint64_t b = 5; // Input 2
uint64_t a_lo = (uint32_t)a;
uint64_t a_hi = a >> 32;
uint64_t b_lo = (uint32_t)b;
uint64_t b_hi = b >> 32;
uint64_t a_x_b_hi = a_hi * b_hi;
uint64_t a_x_b_mid = a_hi * b_lo;
uint64_t b_x_a_mid = b_hi * a_lo;
uint64_t a_x_b_lo = a_lo * b_lo;
/*
This is implementing schoolbook multiplication:
x1 x0
X y1 y0
-------------
00 LOW PART
-------------
00
10 10 MIDDLE PART
01
-------------
01
11 11 HIGH PART
-------------
*/
// 64-bit product two 32-bit values
uint64_t middle = a_x_b_mid (a_x_b_lo >> 32) uint32_t(b_x_a_mid);
// 64-bit product two 32-bit values
uint64_t carry = a_x_b_hi (middle >> 32) (b_x_a_mid >> 32);
// Add LOW PART and lower half of MIDDLE PART
uint64_t result = (middle << 32) | uint32_t(a_x_b_lo);
std::cout << result << std::endl;
std::cout << carry << std::endl;
}
This results in
Program stdout
185
0
Godbolt link: https://godbolt.org/z/97xhMvY53
Or you could use __uint128_t which is non-standard but widely available.
static inline void mul64(uint64_t a, uint64_t b, uint64_t& result, uint64_t& carry) {
__uint128_t va(a);
__uint128_t vb(b);
__uint128_t vr = va * vb;
result = uint64_t(vr);
carry = uint64_t(vr >> 64);
}
CodePudding user response:
In the title of this question, you said you wanted to multiply two integers. But the code you found on that other Q&A (Getting the high part of 64 bit integer multiplication) isn't trying to do that, it's only trying to get the high half of the full product. For a 64x64 => 128-bit product, the high half is product >> 64
.
37 x 5 = 185
185 >> 64 = 0
It's correctly emulating multihi = (37 * (unsigned __int128)5) >> 64
, and you're forgetting about the >>64
part.
__int128
is a GNU C extension; it's much more efficient than emulating it manually with pure ISO C, but only supported on 64-bit targets by current compilers. See my answer on the same question. (ISO C23 is expected to have _BitInt(128)
or whatever width you specify.)
In comments you were talking about floating-point mantissas. In an FP multiply, you have two n-bit mantissas (usually with their leading bits set), so the high half of the 2n-bit product will have n significant bits (more or less; maybe actually one place to the right IIRC).
Something like 37 x 5 would only happen with tiny subnormal floats, where the product would indeed underflow to zero. But in that case, it would be because you only get subnormals at the limits of the exponent range, and (37 * 2^-1022) * (5 * 2^-1022)
would be 186 * 2^-2044
, an exponent way too small to be represented in an FP format like IEEE binary64 aka double
where -1022
was the minimum exponent.
You're using integers where a>>63
isn't 1
, in fact they're both less than 2^32 so there are no significant bits outside the low 64 bits of the full 128-bit product.