Skip to content

Use a much better multiply algorithm which takes advantage of modern CPUs #9

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
easyaspi314 wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Use a much better multiply algorithm which takes advantage of modern CPUs #9

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
137 changes: 90 additions & 47 deletions uint128_t.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <intrin.h>
_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<uint64_t>(product >> 64);
return static_cast<uint64_t>(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<uint64_t>(static_cast<uint32_t>(val));
}
_UINT128_T_MULT_TARGET static inline uint64_t upper32(uint64_t val){
return static_cast<uint64_t>(static_cast<uint32_t>(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<uint64_t>(Rn) * static_cast<uint64_t>(Rm);
// product += RdLo;
// product += RdHi;
// RdLo = static_cast<uint32_t>(product & 0xFFFFFFFF);
// RdHi = static_cast<uint32_t>(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<uint64_t>(r2) * static_cast<uint64_t>(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;
}
Expand Down
15 changes: 10 additions & 5 deletions uint128_t.include
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <typename T, typename = typename std::enable_if<std::is_integral<T>::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 <typename T, typename = typename std::enable_if<std::is_integral<T>::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;
}
Expand Down
49 changes: 49 additions & 0 deletions uint128_t_config.include
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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