Compiler Explorer

Source code

#include <stdint.h>

#ifdef _MSC_VER
#  include <intrin.h>
#else
// MSVC doesn't optimize 32x32 => 64b multiplication without its intrinsic
// But good compilers can just use this to get a single mul instruction
static inline
uint64_t __emulu(uint32_t x, uint32_t y) {
     return x * (uint64_t)y;
}
#endif

uint64_t __umul128(uint64_t top, uint64_t bot) {

uint64_t BD = (top & 0xFFFFFFFF) * (bot & 0xFFFFFFFF);
    uint64_t AC = (top >> 32)        * (bot >> 32);

uint64_t BC = (top & 0xFFFFFFFF) * (bot >> 32);
    uint64_t AD = (top >> 32)        * (bot & 0xFFFFFFFF);

uint64_t sum1 = (BC & 0xFFFFFFFF);
    uint64_t sum2 = (AD & 0xFFFFFFFF);

sum1 += (BD >> 32);

sum2 += sum1;

sum1 = sum2 >> 32;
    sum2 <<= 32;

sum1 += (BD & 0xFFFFFFFF);
    sum2 += AC;

sum1 += (BC >> 32);
    sum2 += (AD >> 32);

sum2 += sum1;

return sum2;
}

uint64_t big128_mul(uint64_t x, uint64_t y)
{
	/* x * y = (z2 << 64) + (z1 << 32) + z0
	 * where z2 = x1 * y1
	 *       z1 = x0 * y1 + x1 * y0
	 *       z0 = x0 * y0
	 */
	
	uint64_t x0 =	x & 0xFFFFFFFF, x1 = x >> 32, y0 = y  & 0xFFFFFFFF, y1 = y >> 32;
	uint64_t z0 = x0 * y0;
	uint64_t z1a = x1 * y0;
	uint64_t z1b = x0 * y1;
	uint64_t z2 = x1 * y1;

uint32_t z0l = z0 & UINT32_MAX;
	uint32_t z0h = z0 >> 32u;
	
	uint64_t z1al = z1a & UINT32_MAX;
	uint64_t z1bl = z1b & UINT32_MAX;
	uint64_t z1l = z1al + z1bl + z0h;
	
	uint64_t z1h = (z1a >> 32u) + (z1b >> 32u) + (z1l >> 32u);
	z2 += z1h;

return ((z0l & 0xFFFFFFFF) | (z1l & 0xFFFFFFFF)) + ((z2 & 0xFFFFFFFF) | (z2 >> 32u));
}

// This is still pretty ugly with MSVC, branching on the carry
//  and using XMM store / integer reload to zero a register!
// But at least it inlines 4 mul instructions
//  instead of calling a generic 64x64 => 64b multiply helper function
uint64_t mul128_msvc(uint64_t multiplier, uint64_t multiplicand) 
{
    // multiplier   = ab = a * 2^32 + b
    // multiplicand = cd = c * 2^32 + d
    // ab * cd = a * c * 2^64 + (a * d + b * c) * 2^32 + b * d
    uint64_t a = multiplier >> 32;
    uint64_t b = (uint32_t)multiplier; // & 0xFFFFFFFF;
    uint64_t c = multiplicand >> 32;
    uint64_t d = (uint32_t)multiplicand; // & 0xFFFFFFFF;

//uint64_t ac = __emulu(a, c);
    uint64_t ad = __emulu(a, d);
    //uint64_t bc = __emulu(b, c);
    uint64_t bd = __emulu(b, d);

uint64_t adbc = ad + __emulu(b , c);
    uint64_t adbc_carry = (adbc < ad); // ? 1 : 0;
    // MSVC gets confused by the ternary and makes worse code than using a boolean in an integer context for 1 : 0

// multiplier * multiplicand = product_hi * 2^64 + product_lo
    uint64_t product_lo = bd + (adbc << 32);
    uint64_t product_lo_carry = (product_lo < bd); // ? 1 : 0;
    return product_lo + __emulu(a , c) + (adbc >> 32) + (adbc_carry << 32) + product_lo_carry;

}

// The above compiles badly in 64-bit mode
// This compiles to a single mul or mulx
#if defined(__x86_64__) && defined(__GNUC__)
uint64_t umul128_gcc(uint64_t multiplier, uint64_t multiplicand, 
    uint64_t *product_hi) {
    unsigned __int128 a = multiplier, b = multiplicand;
    unsigned __int128 prod = a*b;
    *product_hi = prod >> 64;
    return prod;
}
#endif