From 1a98fb7da13731aa45cb35ad8c019eed61330b26 Mon Sep 17 00:00:00 2001 From: "easyaspi314 (Devin)" Date: Sun, 8 Sep 2019 10:06:20 -0400 Subject: [PATCH] Use a much better multiply algorithm which takes advantage of modern CPU instructions. This multiply algorithm is optimized for multiple platforms, 32-bit and 64-bit. Notably, instead of doing a full 128-bit product manually, it does a 64->128 bit long multiply on the low bits followed by two normal 64-bit multiplies which are added to the high bits. This is to take advantage of native 64-bit arithmetic, and if it is supported, compiler intrinsics/extensions which do the long multiply for us. This will use _umul128 or __uint128_t (except on wasm) if they are available. These will expand to the single operand MULQ instruction on x86_64, or MUL/UMULH on aarch64. Otherwise, the long multiply will do a simple yet fast grade-school multiply for other targets. It is optimized to suggest the powerful UMAAL instruction on ARM, but it doesn't require it to be fast (really you only need long multiply + adc). Because it makes a large difference in performance, I added a check which will switch off SSE2 on GCC for x86 and attempt to switch to ARM mode on Thumb-1 targets if it is available. Switching off SSE2 in GCC for x86 prevents a laggy partial vectorization which is detrimental to the performance (mainly because instead of switching registers we have to do PSHUFD, PUNPCKLDQ, and PSLLQ). Clang doesn't need this. Switching to ARM mode instead of Thumb-1 is crucial to make this algorithm usable. ARM is well known for its powerful multiplier, with the already-excellent UMULL and UMLAL instructions, and in ARMv6, the mighty UMAAL which allows a long multply to be done in 4 instructions. Thumb-1 does not have access to any of these, only having a 32-bit->32-bit multiply. This would result in calling compiler functions for something a simple mode switch can do in a quarter of the time or less. This will be turned off if the M profile is detected or if _UINT128_T_FORCE_THUMB is defined. Signed-off-by: easyaspi314 (Devin) --- uint128_t.cpp | 137 +++++++++++++++++++++++++-------------- uint128_t.include | 15 +++-- uint128_t_config.include | 49 ++++++++++++++ 3 files changed, 149 insertions(+), 52 deletions(-) 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