usquare

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); \

include/CppCore.Test/Math/Util.h

@@ -473,6 +473,45 @@ namespace CppCore { namespace Test { namespace Math
          return true;
       }
 
+      template<size_t NA, size_t NR, size_t I=1000>
+      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()); }

include/CppCore/Math/Util.h

@@ -2717,6 +2717,7 @@ namespace CppCore
    template<typename UINT1, typename UINT2, typename UINT3>
    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
+   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+   /// <summary>
+   /// 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.
+   /// </summary>
+   template<typename UINT1, typename UINT2>
+   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<UINT1> t(a);
+         CppCore::usquare(t, r);
+      }
+      else if constexpr (sizeof(UINT2) % sizeof(size_t) != 0)
+      {
+         Padded<UINT2> 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<NA && j<NR; i++, j+=2)
+         {
+            CppCore::umul128(ap[i], ap[i], tl, th);
+            CppCore::addcarry64(rp[j], tl, rp[j], c);
+            if (j+1 < NR)
+               CppCore::addcarry64(rp[j+1], th, rp[j+1], c);
+         }
+      }
+   #endif
+      else if constexpr (sizeof(UINT1) % 4 == 0 && sizeof(UINT2) % 4 == 0)
+      {
+         // 32-Bit CPU or only Multiples of 32-Bit
+         constexpr size_t NA = sizeof(UINT1) / 4;
+         constexpr size_t NR = sizeof(UINT2) / 4;
+         assert((void*)&a != (void*)&r);
+         uint32_t* ap = (uint32_t*)&a;
+         uint32_t* rp = (uint32_t*)&r;
+         uint32_t  tl, th, k;
+         uint8_t   c;
+         for (size_t i = NA+NA; i < NR; i++)
+            rp[i] = 0U;
+         // calculate one triangle
+         rp[0] = 0U;
+         CppCore::umul64(ap[0], ap[1], rp[1], k);
+         for (size_t j = 2; j < MIN(NA,NR); j++)
+         {
+            CppCore::umul64(ap[0], ap[j], tl, th);
+            c = 0;
+            CppCore::addcarry32(tl, k, rp[j], c);
+            CppCore::addcarry32(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::umul64(ap[i], ap[j], tl, th);
+               c = 0;
+               CppCore::addcarry32(tl, rp[i+j], tl, c);
+               CppCore::addcarry32(th, 0ULL, th, c);
+               c = 0;
+               CppCore::addcarry32(tl, k, rp[i+j], c);
+               CppCore::addcarry32(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::addcarry32(rp[i], rp[i], rp[i], c);
+         // add diagonal terms
+         c = 0;
+         for (size_t i=0, j=0; i<NA && j<NR; i++, j+=2)
+         {
+            CppCore::umul64(ap[i], ap[i], tl, th);
+            CppCore::addcarry32(rp[j], tl, rp[j], c);
+            if (j+1 < NR)
+               CppCore::addcarry32(rp[j+1], th, rp[j+1], c);
+         }
+      }
+   }
+
    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    // UNSIGNED DIVISION+MODULO BY CONSTANTS
    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

src/CppCore.Interface.C/cppcore.cpp

@@ -63,6 +63,7 @@ void  cppcore_free (void* ptr)         { CPPCORE_ALIGNED_FREE(ptr); }
   void name ## _add   (void* a, void* b, void* r)         { CppCore::uadd(*(cname*)a, *(cname*)b, *(cname*)r);}            \
   void name ## _sub   (void* a, void* b, void* r)         { CppCore::usub(*(cname*)a, *(cname*)b, *(cname*)r);}            \
   void name ## _mul   (void* a, void* b, void* r)         { CppCore::umul(*(cname*)a, *(cname*)b, *(cname*)r);}            \
+  void name ## _square(void* a, void* r)                  { CppCore::usquare(*(cname*)a, *(cname*)r);}                     \
   void name ## _divmod(void* a, void* b, void* q, void* r){ CppCore::udivmod(*(cname*)q,*(cname*)r,*(cname*)a,*(cname*)b);}\
   void name ## _mulmod(void* a, void* b, void* m, void* r){ CppCore::umulmod(*(cname*)a,*(cname*)b,*(cname*)m,*(cname*)r);}\
   void name ## _powmod(void* a, void* b, void* m, void* r){ CppCore::upowmod(*(cname*)a,*(cname*)b,*(cname*)m,*(cname*)r);}\

src/CppCore.Test/Test.cpp

@@ -365,6 +365,7 @@ int main()
    TEST(CppCore::Test::Math::Util::subborrow32,     "subborrow32:      ", std::endl);
    TEST(CppCore::Test::Math::Util::subborrow64,     "subborrow64:      ", std::endl);
    TEST(CppCore::Test::Math::Util::mul128,          "mul128:           ", std::endl);
+   TEST(CppCore::Test::Math::Util::usquare,         "usquare:          ", std::endl);
    TEST(CppCore::Test::Math::Util::overflowadd16,   "overflowadd16:    ", std::endl);
    TEST(CppCore::Test::Math::Util::overflowadd32,   "overflowadd32:    ", std::endl);
    TEST(CppCore::Test::Math::Util::overflowadd64,   "overflowadd64:    ", std::endl);