diff --git a/include/CppCore.Interface.C/cppcore.h b/include/CppCore.Interface.C/cppcore.h index 9b7e9adc..21c70e40 100644 --- a/include/CppCore.Interface.C/cppcore.h +++ b/include/CppCore.Interface.C/cppcore.h @@ -19,6 +19,7 @@ CPPCORE_EXPORT void name ## _add (void* a, void* b, void* r); \ CPPCORE_EXPORT void name ## _sub (void* a, void* b, void* r); \ CPPCORE_EXPORT void name ## _mul (void* a, void* b, void* r); \ + CPPCORE_EXPORT void name ## _square(void* a, void* r); \ CPPCORE_EXPORT void name ## _divmod(void* a, void* b, void* q, void* r); \ CPPCORE_EXPORT void name ## _mulmod(void* a, void* b, void* m, void* r); \ CPPCORE_EXPORT void name ## _powmod(void* a, void* b, void* m, void* r); \ diff --git a/include/CppCore.Test/Math/Util.h b/include/CppCore.Test/Math/Util.h index cc54929e..27263de3 100644 --- a/include/CppCore.Test/Math/Util.h +++ b/include/CppCore.Test/Math/Util.h @@ -473,6 +473,45 @@ namespace CppCore { namespace Test { namespace Math return true; } + template + INLINE static bool usquare_template() + { + CppCore::Random::Default64 prng; + uint64_t a[NA]; + uint64_t r1[NR]; + uint64_t r2[NR]; + for (size_t i = 0; i < I; i++) + { + prng.fill(a); + CppCore::umul(a, a, r1); + CppCore::usquare(a, r2); + if (!CppCore::equal(r1, r2)) + return false; + } + return true; + } + + INLINE static bool usquare() + { + uint64_t a[4] = { 0xFFFFFFFFFFFFFFFFULL,0xFFFFFFFFFFFFFFFFULL,0xFFFFFFFFFFFFFFFFULL,0xFFFFFFFFFFFFFFFFULL }; + uint64_t r1[9] = { 0x0000000000000001ULL,0x0000000000000000ULL,0x0000000000000000ULL,0x0000000000000000ULL, + 0xFFFFFFFFFFFFFFFEULL,0xFFFFFFFFFFFFFFFFULL,0xFFFFFFFFFFFFFFFFULL,0xFFFFFFFFFFFFFFFFULL, + 0x0000000000000000ULL }; + uint64_t r2[9]; + CppCore::usquare(a, r2); + if (!CppCore::equal(r1, r2)) + return false; + + return + usquare_template<4, 2>() && + usquare_template<4, 4>() && + usquare_template<4, 8>() && + usquare_template<4, 10>() && + usquare_template<3, 2>() && + usquare_template<3, 3>() && + usquare_template<3, 4>(); + } + INLINE static bool overflowadd16() { uint16_t ru; @@ -1581,6 +1620,7 @@ namespace CppCore { namespace Test { namespace VS { namespace Math { TEST_METHOD(SUBBORROW32) { Assert::AreEqual(true, CppCore::Test::Math::Util::subborrow32()); } TEST_METHOD(SUBBORROW64) { Assert::AreEqual(true, CppCore::Test::Math::Util::subborrow64()); } TEST_METHOD(MUL128) { Assert::AreEqual(true, CppCore::Test::Math::Util::mul128()); } + TEST_METHOD(USQUARE) { Assert::AreEqual(true, CppCore::Test::Math::Util::usquare()); } TEST_METHOD(OVERFLOWADD16) { Assert::AreEqual(true, CppCore::Test::Math::Util::overflowadd16()); } TEST_METHOD(OVERFLOWADD32) { Assert::AreEqual(true, CppCore::Test::Math::Util::overflowadd32()); } TEST_METHOD(OVERFLOWADD64) { Assert::AreEqual(true, CppCore::Test::Math::Util::overflowadd64()); } diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index e87123ad..56dc8eb2 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2717,6 +2717,7 @@ namespace CppCore template INLINE static void umul(const UINT1& a, const UINT2& b, UINT3& r) { + static_assert(sizeof(UINT1) != 0 && sizeof(UINT2) != 0 && sizeof(UINT3) != 0); if constexpr (sizeof(UINT1) < sizeof(size_t)) { CppCore::umul((size_t)a, b, r); } else if constexpr (sizeof(UINT2) < sizeof(size_t)) { CppCore::umul(a, (size_t)b, r); } else if constexpr (sizeof(UINT3) < sizeof(size_t)) @@ -2898,6 +2899,155 @@ namespace CppCore r = a * b; } + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // SQUARING + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + /// + /// Simple Squaring (a*a=r) in O((n*n)/2) automatically selecting 64-Bit or 32-Bit operations and chunks. + /// Calculates all bits of r, performing a full wide multiplication if sizeof(r) larger-equal sizeof(a)+sizeof(a). + /// This exploits the fact that a large number of pairwise products (triangles) are identical when squaring. + /// + template + INLINE static void usquare(const UINT1& a, UINT2& r) + { + static_assert(sizeof(UINT1) != 0 && sizeof(UINT2) != 0); + if constexpr (sizeof(UINT1) < sizeof(size_t)) { CppCore::umul((size_t)a, (size_t)a, r); } + else if constexpr (sizeof(UINT2) < sizeof(size_t)) + { + size_t t; + CppCore::usquare(a, t); + CppCore::clone(r, *(UINT2*)&t); + } + else if constexpr (sizeof(UINT1) % sizeof(size_t) != 0) + { + Padded t(a); + CppCore::usquare(t, r); + } + else if constexpr (sizeof(UINT2) % sizeof(size_t) != 0) + { + Padded t; + CppCore::usquare(a, t); + CppCore::clone(r, t.v); + } + #if defined(CPPCORE_CPU_64BIT) + else if constexpr (sizeof(UINT1) == 16 && sizeof(UINT2) == 16) + { + uint64_t* ap = (uint64_t*)&a; + uint64_t* rp = (uint64_t*)&r; + CppCore::umul128(ap[0], ap[1], ap[0], ap[1], rp[0], rp[1]); + } + else if constexpr (sizeof(UINT1) % 8 == 0 && sizeof(UINT2) % 8 == 0) + { + // 64-Bit CPU and Multiples of 64-Bit + constexpr size_t NA = sizeof(UINT1) / 8; + constexpr size_t NR = sizeof(UINT2) / 8; + assert((void*)&a != (void*)&r); + uint64_t* ap = (uint64_t*)&a; + uint64_t* rp = (uint64_t*)&r; + uint64_t tl, th, k; + uint8_t c; + for (size_t i = NA+NA; i < NR; i++) + rp[i] = 0ULL; + // calculate one triangle + rp[0] = 0ULL; + CppCore::umul128(ap[0], ap[1], rp[1], k); + for (size_t j = 2; j < MIN(NA,NR); j++) + { + CppCore::umul128(ap[0], ap[j], tl, th); + c = 0; + CppCore::addcarry64(tl, k, rp[j], c); + CppCore::addcarry64(th, 0ULL, k, c); + } + if constexpr (NA < NR) + rp[NA] = k; + for (size_t i = 1; i < NA && 2*i+1 < NR; i++) + { + k = 0ULL; + for (size_t j = i+1; j < NA && i+j < NR; j++) + { + CppCore::umul128(ap[i], ap[j], tl, th); + c = 0; + CppCore::addcarry64(tl, rp[i+j], tl, c); + CppCore::addcarry64(th, 0ULL, th, c); + c = 0; + CppCore::addcarry64(tl, k, rp[i+j], c); + CppCore::addcarry64(th, 0ULL, k, c); + } + if (i+NA < NR) + rp[i+NA] = k; + } + // double the triangular sum + c = 0; + for (size_t i = 0; i < MIN(NA+NA,NR); i++) + CppCore::addcarry64(rp[i], rp[i], rp[i], c); + // add diagonal terms + c = 0; + for (size_t i=0, j=0; i