Multiply two uint64_ts and store result to uint64_t doesnt seem to work?

Question

I am trying to multiply two uint64_ts and store the result to uint64_t. I found an existing answer on Stackoverflow which splits the inputs in to their four uint32_ts 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?
}

Answer 1

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);
}

Answer 2

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.