diff --git a/uint128_t.cpp b/uint128_t.cpp index ce5b74e..3dd2868 100644 --- a/uint128_t.cpp +++ b/uint128_t.cpp @@ -250,53 +250,96 @@ uint128_t & uint128_t::operator-=(const uint128_t & rhs){ return *this; } -uint128_t uint128_t::operator*(const uint128_t & rhs) const{ - // split values into 4 32-bit parts - uint64_t top[4] = {UPPER >> 32, UPPER & 0xffffffff, LOWER >> 32, LOWER & 0xffffffff}; - uint64_t bottom[4] = {rhs.UPPER >> 32, rhs.UPPER & 0xffffffff, rhs.LOWER >> 32, rhs.LOWER & 0xffffffff}; - uint64_t products[4][4]; - - // multiply each component of the values - for(int y = 3; y > -1; y--){ - for(int x = 3; x > -1; x--){ - products[3 - x][y] = top[x] * bottom[y]; - } - } - - // first row - uint64_t fourth32 = (products[0][3] & 0xffffffff); - uint64_t third32 = (products[0][2] & 0xffffffff) + (products[0][3] >> 32); - uint64_t second32 = (products[0][1] & 0xffffffff) + (products[0][2] >> 32); - uint64_t first32 = (products[0][0] & 0xffffffff) + (products[0][1] >> 32); - - // second row - third32 += (products[1][3] & 0xffffffff); - second32 += (products[1][2] & 0xffffffff) + (products[1][3] >> 32); - first32 += (products[1][1] & 0xffffffff) + (products[1][2] >> 32); - - // third row - second32 += (products[2][3] & 0xffffffff); - first32 += (products[2][2] & 0xffffffff) + (products[2][3] >> 32); - - // fourth row - first32 += (products[3][3] & 0xffffffff); - - // move carry to next digit - third32 += fourth32 >> 32; - second32 += third32 >> 32; - first32 += second32 >> 32; - - // remove carry from current digit - fourth32 &= 0xffffffff; - third32 &= 0xffffffff; - second32 &= 0xffffffff; - first32 &= 0xffffffff; - - // combine components - return uint128_t((first32 << 32) | second32, (third32 << 32) | fourth32); -} - -uint128_t & uint128_t::operator*=(const uint128_t & rhs){ +// Algorithm summary: +// +// First we do a 64-bit to 128-bit long multiply for the low bits, then we use a normal 64-bit +// multiply on the high bits. +// This allows us to take advantage of not only compiler intrinsics but native 64-bit arithmetic. + +// First we define the generic multlong64 methods. These will all do basically what _umul128 does. + +// MSVC _umul128 +#if _UINT128_T_MULT_TYPE == _UINT128_T_MULT_MSVC +#include +_UINT128_T_MULT_TARGET uint64_t uint128_t::multlong64(uint64_t lhs, uint64_t rhs, uint64_t *high){ + return _umul128(lhs, rhs, high); +} + +// GCC __uint128_t +#elif _UINT128_T_MULT_TYPE == _UINT128_T_MULT_GCC +_UINT128_T_MULT_TARGET uint64_t uint128_t::multlong64(uint64_t lhs, uint64_t rhs, uint64_t *high){ + __uint128_t product = static_cast<__uint128_t>(lhs) * static_cast<__uint128_t>(rhs); + *high = static_cast(product >> 64); + return static_cast(product & 0xFFFFFFFFFFFFFFFF); +} + +// Portable version +#else +// The double cast helps MSVC +_UINT128_T_MULT_TARGET static inline uint64_t lower32(uint64_t val){ + return static_cast(static_cast(val)); +} +_UINT128_T_MULT_TARGET static inline uint64_t upper32(uint64_t val){ + return static_cast(static_cast(val >> 32)); +} + +_UINT128_T_MULT_TARGET uint64_t uint128_t::multlong64(uint64_t lhs, uint64_t rhs, uint64_t *high){ + // This is a fast yet simple grade school 2x2 long multiply. + // The way we add the cross products avoids the need to track 64-bit carries due to the properties + // of multiplying by 11 (technically 0x100000001) capping the sums at 0xFFFFFFFFFFFFFFFF, and it + // tries to match the powerful ARMv6's UMAAL function which was explicitly designed for + // multiprecision multiplication: + // + // void umaal(uint32_t &RdLo, uint32_t &RdHi, const uint32_t Rn, const uint32_t Rm){ + // uint64_t product = static_cast(Rn) * static_cast(Rm); + // product += RdLo; + // product += RdHi; + // RdLo = static_cast(product & 0xFFFFFFFF); + // RdHi = static_cast(product >> 32); + // } + // + // This allows a 64-bit to 128-bit multiply to be calculated in 4 instructions, ~3 cycles each. + // + // It is still fast for other platforms, though. + // + // TODO: Use better variable names + + // Calculate the cross products... + uint64_t lo_lo = lower32(lhs) * lower32(rhs); + uint64_t hi_lo = upper32(lhs) * lower32(rhs); + uint64_t lo_hi = lower32(lhs) * upper32(rhs); + uint64_t hi_hi = upper32(lhs) * upper32(rhs); + + // then add them together. + uint64_t cross = upper32(lo_lo) + lower32(hi_lo) + lo_hi; + uint64_t top = upper32(hi_lo) + upper32(cross) + hi_hi; + + // Done + *high = top; + return (cross << 32) | (lo_lo & 0xFFFFFFFF); +} +#endif + +// Now we do the full 128-bit multiply. +// +// This is based on the 64-bit multiply idiom on ARM, only for 128-bit integers instead of 64-bit. +// +// @ {r0, r1} * {r2, r3} (little endian) +// umull r12, lr, r2, r0 @ {r12, lr} = static_cast(r2) * static_cast(r0) +// mla r4, r2, r1, lr @ r4 = r2 * r1 + lr; +// mla r1, r3, r0, r4 @ r1 = r3 * r0 + r4 +// @ result is in {r12, r1} + +_UINT128_T_MULT_TARGET uint128_t uint128_t::operator*(const uint128_t & rhs) const{ + uint64_t high; + uint64_t low = multlong64(LOWER, rhs.LOWER, &high); + uint128_t acc(high, low); + acc.UPPER += LOWER * rhs.UPPER; + acc.UPPER += UPPER * rhs.LOWER; + return acc; +} + +_UINT128_T_MULT_TARGET uint128_t & uint128_t::operator*=(const uint128_t & rhs){ *this = *this * rhs; return *this; } diff --git a/uint128_t.include b/uint128_t.include index 64066e4..d1d2816 100644 --- a/uint128_t.include +++ b/uint128_t.include @@ -264,18 +264,23 @@ class uint128_t{ *this = *this - rhs; return *this; } - - uint128_t operator*(const uint128_t & rhs) const; + // Note: _UINT128_T_MULTI_TARGET is for disabling SSE2 and switching to ARM mode from Thumb to greatly + // improve the performance. +private: + // XXX: make this public? + _UINT128_T_MULT_TARGET static uint64_t multlong64(uint64_t lhs, uint64_t rhs, uint64_t *high); +public: + _UINT128_T_MULT_TARGET uint128_t operator*(const uint128_t & rhs) const; template ::value, T>::type > - uint128_t operator*(const T & rhs) const{ + _UINT128_T_MULT_TARGET uint128_t operator*(const T & rhs) const{ return *this * uint128_t(rhs); } - uint128_t & operator*=(const uint128_t & rhs); + _UINT128_T_MULT_TARGET uint128_t & operator*=(const uint128_t & rhs); template ::value, T>::type > - uint128_t & operator*=(const T & rhs){ + _UINT128_T_MULT_TARGET uint128_t & operator*=(const T & rhs){ *this = *this * uint128_t(rhs); return *this; } diff --git a/uint128_t_config.include b/uint128_t_config.include index 66a6082..600c7bd 100644 --- a/uint128_t_config.include +++ b/uint128_t_config.include @@ -15,5 +15,54 @@ #define _UINT128_T_EXPORT __attribute__((visibility("default"))) #define _UINT128_T_IMPORT __attribute__((visibility("default"))) #endif + + // Multiply stuff. The algorithm is usually pretty efficient on its own, but we can do better. + // Notably this includes using target intrinsics and switching to ARM mode on Thumb-1 targets. + + // Portable grade school long multiply + #define _UINT128_T_MULT_PORTABLE 0 + // _umul128 + #define _UINT128_T_MULT_MSVC 1 + // __uint128_t + #define _UINT128_T_MULT_GCC 2 + + #ifndef _UINT128_T_MULT_TYPE + #if defined(_MSC_VER) && (defined(_M_IX64) || defined(_M_AMD64)) + #define _UINT128_T_MULT_TYPE _UINT128_T_MULT_MSVC + // Clang defines __uint128_t on WASM and asm.js even though it has to use builtins for multiplication. + // As a result, the algorithm is slower than it would be if it was done manually. + #elif defined(__GNUC__) && defined(__SIZEOF_INT128__) && !defined(__wasm__) && !defined(__asmjs__) + #define _UINT128_T_MULT_TYPE _UINT128_T_MULT_GCC + #else + #define _UINT128_T_MULT_TYPE _UINT128_T_MULT_PORTABLE + #endif + #endif + + // Some feature flags to optimize the manual algorithm. + #if defined(__GNUC__) && _UINT128_T_MULT_TYPE == _UINT128_T_MULT_PORTABLE + // GCC for x86 loves to use SSE2 in the multiply code, but it is inefficient because it uses shifts and shuffles which aren't 'free' + // like normal register swapping. Clang doesn't need this flag, as it never uses this. + + #if !defined(__clang__) && defined(__SSE2__) + #define _UINT128_T_MULT_TARGET __attribute__((__target__("no-sse2"))) + + // In Thumb-1, the multiply algorithm is heavily crippled because the powerful UMULL, UMLAL, and UMAAL are inaccesible. + // Even if reading 32-bit instructions is slower, it is almost always faster to switch to ARM mode here. + // If you are compiling for the M profile and this is falsely being triggered, define _UINT128_T_FORCE_THUMB or use the + // correct -march flag. + // + // Note: Clang sometimes emits warnings like this: + // '+soft-float-abi' is not a recognized feature for this target (ignoring feature) + // These can safely be ignored. + + #elif defined(__thumb__) && !defined(__thumb2__) && defined(__ARM_ARCH_ISA_ARM) && !defined(_UINT128_T_FORCE_THUMB) + #define _UINT128_T_MULT_TARGET __attribute__((__target__("arm"))) + #else + #define _UINT128_T_MULT_TARGET + #endif + #else + #define _UINT128_T_MULT_TARGET + #endif + #endif