From dd8228581a4cc74d4bff0e0b116372891ee8f82f Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sat, 25 Jul 2026 18:28:23 +0200 Subject: [PATCH 01/22] usquare wip --- include/CppCore.Test/Math/Util.h | 29 +++++++ include/CppCore/Math/Util.h | 131 +++++++++++++++++++++++++++++++ src/CppCore.Test/Test.cpp | 1 + 3 files changed, 161 insertions(+) diff --git a/include/CppCore.Test/Math/Util.h b/include/CppCore.Test/Math/Util.h index cc54929e..9e5aa457 100644 --- a/include/CppCore.Test/Math/Util.h +++ b/include/CppCore.Test/Math/Util.h @@ -473,6 +473,34 @@ namespace CppCore { namespace Test { namespace Math return true; } + INLINE static bool usquare() + { + CppCore::Random::Default64 prng; + + uint64_t a[4]; + uint64_t r1_1[4]; + uint64_t r1_2[4]; + uint64_t r2_1[8]; + uint64_t r2_2[8]; + + for (size_t i = 0; i < 100; i++) + { + prng.fill(a); + + CppCore::umul(a, a, r1_1); + CppCore::usquare(a, r1_2); + if (!CppCore::equal(r1_1, r1_2)) + return false; + + CppCore::umul(a, a, r2_1); + CppCore::usquare(a, r2_2); + if (!CppCore::equal(r2_1, r2_2)) + return false; + + } + return true; + } + INLINE static bool overflowadd16() { uint16_t ru; @@ -1581,6 +1609,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..0497f51f 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2898,6 +2898,137 @@ namespace CppCore r = a * b; } + + + /// + /// + /// + template + INLINE static void usquare(const UINT1& a, UINT2& r) + { + 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; + + // limbs assumed little-endian in memory (word 0 = least significant) + const uint64_t* A = reinterpret_cast(&a); + uint64_t* R = reinterpret_cast(&r); + + for (size_t k = 0; k < NR; ++k) + R[k] = 0; + for (size_t i = 0; i < NA; ++i) + { + // smallest column this row can touch is i+(i+1) = 2i+1; once that's + // out of range, this and every later row contribute nothing + if (2 * i + 1 >= NR) + break; + + uint64_t carry64 = 0; // full 64-bit inter-term carry for this row + size_t j = i + 1; + + for (; j < NA; ++j) + { + size_t idx = i + j; + if (idx >= NR) + break; + + uint64_t lo, hi; + umul128(A[i], A[j], lo, hi); + + // fold in the carry from the previous term (its "hi") first + uint64_t combined; + uint8_t c1 = 0; + uint8_t c2 = 0; + addcarry64(lo, carry64, combined, c1); + + // then accumulate into the destination limb + addcarry64(R[idx], combined, R[idx], c2); + + // hi <= 2^64-2 always, so this plain add never overflows 64 bits + carry64 = hi + static_cast(c1) + static_cast(c2); + } + + // loop ran to completion (never broke early) => the pending carry64 + // (this row's last term's "hi", effectively) lands at i+NA + if (j == NA && i + NA < NR) + { + uint8_t carry = 0; + addcarry64(R[i + NA], carry64, R[i + NA], carry); + for (size_t k = i + NA + 1; k < NR; ++k) + addcarry64(R[k], 0, R[k], carry); + } + } + + // ---- pass 2: double the triangular sum in one sweep (single carry chain) ---- + { + uint8_t carry = 0; + for (size_t k = 0; k < NR; ++k) + addcarry64(R[k], R[k], R[k], carry); + } + + // ---- pass 3: diagonal terms a[i]*a[i]; positions (2i,2i+1,2i+2,...) + // never overlap between consecutive i (unlike pass 1's rows), so a + // simple 1-bit carry chain is exact here ---- + { + uint8_t carry = 0; + size_t i = 0; + + for (; i < NA; ++i) + { + size_t idx = 2 * i; + if (idx >= NR) + break; + + uint64_t lo, hi; + umul128(A[i], A[i], lo, hi); + + addcarry64(R[idx], lo, R[idx], carry); + + if (idx + 1 < NR) + addcarry64(R[idx + 1], hi, R[idx + 1], carry); + } + + if (NR > 2 * NA) + addcarry64(R[2 * NA], 0, R[2 * NA], carry); + } + } + #endif + else if constexpr (sizeof(UINT1) % 4 == 0 && sizeof(UINT2) % 4 == 0) + { + assert(false); + } + } + + + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // UNSIGNED DIVISION+MODULO BY CONSTANTS //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/src/CppCore.Test/Test.cpp b/src/CppCore.Test/Test.cpp index b2f23f39..04ea9a00 100644 --- a/src/CppCore.Test/Test.cpp +++ b/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); From ae096b0b2bbe64013320ba7f782baed4f6d98b97 Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sat, 25 Jul 2026 19:16:37 +0200 Subject: [PATCH 02/22] wip --- include/CppCore.Interface.C/cppcore.h | 1 + include/CppCore/Math/Util.h | 19 +++++++++---------- src/CppCore.Interface.C/cppcore.cpp | 1 + 3 files changed, 11 insertions(+), 10 deletions(-) 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/Math/Util.h b/include/CppCore/Math/Util.h index 0497f51f..fff0f49a 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2950,7 +2950,7 @@ namespace CppCore if (2 * i + 1 >= NR) break; - uint64_t carry64 = 0; // full 64-bit inter-term carry for this row + uint64_t k = 0; // full 64-bit inter-term carry for this row size_t j = i + 1; for (; j < NA; ++j) @@ -2959,20 +2959,19 @@ namespace CppCore if (idx >= NR) break; - uint64_t lo, hi; - umul128(A[i], A[j], lo, hi); + uint64_t tl, th; + umul128(A[i], A[j], tl, th); // fold in the carry from the previous term (its "hi") first - uint64_t combined; + //uint64_t combined; uint8_t c1 = 0; uint8_t c2 = 0; - addcarry64(lo, carry64, combined, c1); + addcarry64(tl, k, tl, c1); + addcarry64(th, 0, k, c1); // then accumulate into the destination limb - addcarry64(R[idx], combined, R[idx], c2); - - // hi <= 2^64-2 always, so this plain add never overflows 64 bits - carry64 = hi + static_cast(c1) + static_cast(c2); + addcarry64(tl, R[idx], R[idx], c2); + addcarry64(k, 0, k, c2); } // loop ran to completion (never broke early) => the pending carry64 @@ -2980,7 +2979,7 @@ namespace CppCore if (j == NA && i + NA < NR) { uint8_t carry = 0; - addcarry64(R[i + NA], carry64, R[i + NA], carry); + addcarry64(R[i + NA], k, R[i + NA], carry); for (size_t k = i + NA + 1; k < NR; ++k) addcarry64(R[k], 0, R[k], carry); } diff --git a/src/CppCore.Interface.C/cppcore.cpp b/src/CppCore.Interface.C/cppcore.cpp index 451f36f3..85fabab0 100644 --- a/src/CppCore.Interface.C/cppcore.cpp +++ b/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);}\ From bfe58cf70937d3901db22e83258f3a84a57fcc2b Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sat, 25 Jul 2026 19:22:56 +0200 Subject: [PATCH 03/22] wip --- include/CppCore.Test/Math/Util.h | 2 +- include/CppCore/Math/Util.h | 14 +++++--------- 2 files changed, 6 insertions(+), 10 deletions(-) diff --git a/include/CppCore.Test/Math/Util.h b/include/CppCore.Test/Math/Util.h index 9e5aa457..858499ba 100644 --- a/include/CppCore.Test/Math/Util.h +++ b/include/CppCore.Test/Math/Util.h @@ -483,7 +483,7 @@ namespace CppCore { namespace Test { namespace Math uint64_t r2_1[8]; uint64_t r2_2[8]; - for (size_t i = 0; i < 100; i++) + for (size_t i = 0; i < 10000000; i++) { prng.fill(a); diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index fff0f49a..cd1a1b67 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2960,18 +2960,14 @@ namespace CppCore break; uint64_t tl, th; - umul128(A[i], A[j], tl, th); - // fold in the carry from the previous term (its "hi") first - //uint64_t combined; + CppCore::umul128(A[i], A[j], tl, th); uint8_t c1 = 0; + CppCore::addcarry64(tl, k, tl, c1); + CppCore::addcarry64(th, 0, th, c1); uint8_t c2 = 0; - addcarry64(tl, k, tl, c1); - addcarry64(th, 0, k, c1); - - // then accumulate into the destination limb - addcarry64(tl, R[idx], R[idx], c2); - addcarry64(k, 0, k, c2); + CppCore::addcarry64(tl, R[idx], R[idx], c2); + CppCore::addcarry64(th, 0, k, c2); } // loop ran to completion (never broke early) => the pending carry64 From fd7cd85a5f74aaf8cba13a703c17d33f92be08ae Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sat, 25 Jul 2026 19:31:55 +0200 Subject: [PATCH 04/22] wip --- include/CppCore/Math/Util.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index cd1a1b67..166a8968 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2963,11 +2963,11 @@ namespace CppCore CppCore::umul128(A[i], A[j], tl, th); uint8_t c1 = 0; - CppCore::addcarry64(tl, k, tl, c1); - CppCore::addcarry64(th, 0, th, c1); + CppCore::addcarry64(tl, R[idx], tl, c1); + CppCore::addcarry64(th, 0ULL, th, c1); uint8_t c2 = 0; - CppCore::addcarry64(tl, R[idx], R[idx], c2); - CppCore::addcarry64(th, 0, k, c2); + CppCore::addcarry64(tl, k, R[idx], c2); + CppCore::addcarry64(th, 0ULL, k, c2); } // loop ran to completion (never broke early) => the pending carry64 From 373b23504cd3e21ee7999090d619d0cf8f39097a Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sat, 25 Jul 2026 19:57:13 +0200 Subject: [PATCH 05/22] wip --- include/CppCore/Math/Util.h | 49 ++++++++++++++++++------------------- 1 file changed, 24 insertions(+), 25 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index 166a8968..2264cd5f 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2936,13 +2936,14 @@ namespace CppCore // 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); - // limbs assumed little-endian in memory (word 0 = least significant) const uint64_t* A = reinterpret_cast(&a); uint64_t* R = reinterpret_cast(&r); for (size_t k = 0; k < NR; ++k) R[k] = 0; + for (size_t i = 0; i < NA; ++i) { // smallest column this row can touch is i+(i+1) = 2i+1; once that's @@ -2972,47 +2973,45 @@ namespace CppCore // loop ran to completion (never broke early) => the pending carry64 // (this row's last term's "hi", effectively) lands at i+NA - if (j == NA && i + NA < NR) + if (i+NA < NR) + R[i+NA] = k; + + /*if (j == NA && i + NA < NR) { uint8_t carry = 0; addcarry64(R[i + NA], k, R[i + NA], carry); for (size_t k = i + NA + 1; k < NR; ++k) addcarry64(R[k], 0, R[k], carry); - } + }*/ } - // ---- pass 2: double the triangular sum in one sweep (single carry chain) ---- - { - uint8_t carry = 0; - for (size_t k = 0; k < NR; ++k) - addcarry64(R[k], R[k], R[k], carry); - } + // double the triangular sum in one sweep (single carry chain) + uint8_t carry = 0; + for (size_t k = 0; k < NR; ++k) + CppCore::addcarry64(R[k], R[k], R[k], carry); // ---- pass 3: diagonal terms a[i]*a[i]; positions (2i,2i+1,2i+2,...) // never overlap between consecutive i (unlike pass 1's rows), so a // simple 1-bit carry chain is exact here ---- - { - uint8_t carry = 0; - size_t i = 0; - for (; i < NA; ++i) - { - size_t idx = 2 * i; - if (idx >= NR) - break; + carry = 0; + for (size_t i = 0; i < NA; ++i) + { + size_t idx = 2 * i; + if (idx >= NR) + break; - uint64_t lo, hi; - umul128(A[i], A[i], lo, hi); + uint64_t lo, hi; - addcarry64(R[idx], lo, R[idx], carry); + CppCore::umul128(A[i], A[i], lo, hi); - if (idx + 1 < NR) - addcarry64(R[idx + 1], hi, R[idx + 1], carry); - } + CppCore::addcarry64(R[idx], lo, R[idx], carry); - if (NR > 2 * NA) - addcarry64(R[2 * NA], 0, R[2 * NA], carry); + if (idx + 1 < NR) + CppCore::addcarry64(R[idx + 1], hi, R[idx + 1], carry); } + if (NR > 2 * NA) + CppCore::addcarry64(R[2 * NA], 0, R[2 * NA], carry); } #endif else if constexpr (sizeof(UINT1) % 4 == 0 && sizeof(UINT2) % 4 == 0) From 4346d9566bb91cba8588740ffeea7852981b2c2c Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sat, 25 Jul 2026 20:10:17 +0200 Subject: [PATCH 06/22] wip --- include/CppCore/Math/Util.h | 28 +++++++++------------------- 1 file changed, 9 insertions(+), 19 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index 2264cd5f..c8496a29 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2944,6 +2944,8 @@ namespace CppCore for (size_t k = 0; k < NR; ++k) R[k] = 0; + uint64_t tl, th; + for (size_t i = 0; i < NA; ++i) { // smallest column this row can touch is i+(i+1) = 2i+1; once that's @@ -2960,8 +2962,6 @@ namespace CppCore if (idx >= NR) break; - uint64_t tl, th; - CppCore::umul128(A[i], A[j], tl, th); uint8_t c1 = 0; CppCore::addcarry64(tl, R[idx], tl, c1); @@ -2985,30 +2985,20 @@ namespace CppCore }*/ } - // double the triangular sum in one sweep (single carry chain) + // double the triangular sum uint8_t carry = 0; for (size_t k = 0; k < NR; ++k) CppCore::addcarry64(R[k], R[k], R[k], carry); - // ---- pass 3: diagonal terms a[i]*a[i]; positions (2i,2i+1,2i+2,...) - // never overlap between consecutive i (unlike pass 1's rows), so a - // simple 1-bit carry chain is exact here ---- - + // add diagonal terms carry = 0; - for (size_t i = 0; i < NA; ++i) + //size_t idx = 0; + for (size_t i=0, idx=0; i= NR) - break; - - uint64_t lo, hi; - - CppCore::umul128(A[i], A[i], lo, hi); - - CppCore::addcarry64(R[idx], lo, R[idx], carry); - + CppCore::umul128(A[i], A[i], tl, th); + CppCore::addcarry64(R[idx], tl, R[idx], carry); if (idx + 1 < NR) - CppCore::addcarry64(R[idx + 1], hi, R[idx + 1], carry); + CppCore::addcarry64(R[idx + 1], th, R[idx + 1], carry); } if (NR > 2 * NA) CppCore::addcarry64(R[2 * NA], 0, R[2 * NA], carry); From d75b5716fba0098092c6faeed997639f2e7a63af Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sat, 25 Jul 2026 20:17:15 +0200 Subject: [PATCH 07/22] wip --- include/CppCore.Test/Math/Util.h | 21 +++++++++++++++++---- include/CppCore/Math/Util.h | 9 ++++----- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/include/CppCore.Test/Math/Util.h b/include/CppCore.Test/Math/Util.h index 858499ba..07d56329 100644 --- a/include/CppCore.Test/Math/Util.h +++ b/include/CppCore.Test/Math/Util.h @@ -478,25 +478,38 @@ namespace CppCore { namespace Test { namespace Math CppCore::Random::Default64 prng; uint64_t a[4]; + uint64_t r0_1[2]; + uint64_t r0_2[2]; uint64_t r1_1[4]; uint64_t r1_2[4]; uint64_t r2_1[8]; uint64_t r2_2[8]; + uint64_t r3_1[12]; + uint64_t r3_2[12]; - for (size_t i = 0; i < 10000000; i++) + for (size_t i = 0; i < 1000000; i++) { prng.fill(a); - + // smaller + CppCore::umul(a, a, r0_1); + CppCore::usquare(a, r0_2); + if (!CppCore::equal(r0_1, r0_2)) + return false; + // same size CppCore::umul(a, a, r1_1); CppCore::usquare(a, r1_2); if (!CppCore::equal(r1_1, r1_2)) return false; - + // full wide mul CppCore::umul(a, a, r2_1); CppCore::usquare(a, r2_2); if (!CppCore::equal(r2_1, r2_2)) return false; - + // zero padded + CppCore::umul(a, a, r3_1); + CppCore::usquare(a, r3_2); + if (!CppCore::equal(r3_1, r3_2)) + return false; } return true; } diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index c8496a29..e06acd7d 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2992,13 +2992,12 @@ namespace CppCore // add diagonal terms carry = 0; - //size_t idx = 0; - for (size_t i=0, idx=0; i 2 * NA) CppCore::addcarry64(R[2 * NA], 0, R[2 * NA], carry); From e0da4da4b444bdf2dced4fcf4fc2ab58fca13d18 Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sat, 25 Jul 2026 20:24:04 +0200 Subject: [PATCH 08/22] wip --- include/CppCore/Math/Util.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index e06acd7d..e0368a3d 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2999,8 +2999,8 @@ namespace CppCore if (j+1 < NR) CppCore::addcarry64(R[j+1], th, R[j+1], carry); } - if (NR > 2 * NA) - CppCore::addcarry64(R[2 * NA], 0, R[2 * NA], carry); + //if (NR > 2 * NA) + // CppCore::addcarry64(R[2 * NA], 0, R[2 * NA], carry); } #endif else if constexpr (sizeof(UINT1) % 4 == 0 && sizeof(UINT2) % 4 == 0) From 4698f2f6f6157195bce8c1105e8dde9ebf528be8 Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sat, 25 Jul 2026 20:35:06 +0200 Subject: [PATCH 09/22] wip --- include/CppCore.Test/Math/Util.h | 2 +- include/CppCore/Math/Util.h | 7 +++---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/include/CppCore.Test/Math/Util.h b/include/CppCore.Test/Math/Util.h index 07d56329..0cf57c1b 100644 --- a/include/CppCore.Test/Math/Util.h +++ b/include/CppCore.Test/Math/Util.h @@ -487,7 +487,7 @@ namespace CppCore { namespace Test { namespace Math uint64_t r3_1[12]; uint64_t r3_2[12]; - for (size_t i = 0; i < 1000000; i++) + for (size_t i = 0; i < 100000; i++) { prng.fill(a); // smaller diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index e0368a3d..3de4f51c 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2954,9 +2954,8 @@ namespace CppCore break; uint64_t k = 0; // full 64-bit inter-term carry for this row - size_t j = i + 1; - for (; j < NA; ++j) + for (size_t j = i+1; j < NA; j++) { size_t idx = i + j; if (idx >= NR) @@ -2987,8 +2986,8 @@ namespace CppCore // double the triangular sum uint8_t carry = 0; - for (size_t k = 0; k < NR; ++k) - CppCore::addcarry64(R[k], R[k], R[k], carry); + for (size_t i = 0; i < NR; i++) + CppCore::addcarry64(R[i], R[i], R[i], carry); // add diagonal terms carry = 0; From 6339a47478e001837fc0aac5e9332e6686834f30 Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sat, 25 Jul 2026 21:53:14 +0200 Subject: [PATCH 10/22] wip --- include/CppCore/Math/Util.h | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index 3de4f51c..d010b7be 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2946,15 +2946,9 @@ namespace CppCore uint64_t tl, th; - for (size_t i = 0; i < NA; ++i) + for (size_t i = 0; i < NA && 2*i+1 < NR; i++) { - // smallest column this row can touch is i+(i+1) = 2i+1; once that's - // out of range, this and every later row contribute nothing - if (2 * i + 1 >= NR) - break; - - uint64_t k = 0; // full 64-bit inter-term carry for this row - + uint64_t k = 0; for (size_t j = i+1; j < NA; j++) { size_t idx = i + j; From a4a2745d16fc4e37b555014246df7fd470aff2dd Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sat, 25 Jul 2026 21:57:41 +0200 Subject: [PATCH 11/22] wip --- include/CppCore/Math/Util.h | 16 +--------------- 1 file changed, 1 insertion(+), 15 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index d010b7be..ae467bf6 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2946,6 +2946,7 @@ namespace CppCore uint64_t tl, th; + // calculate one triangle for (size_t i = 0; i < NA && 2*i+1 < NR; i++) { uint64_t k = 0; @@ -2963,26 +2964,13 @@ namespace CppCore CppCore::addcarry64(tl, k, R[idx], c2); CppCore::addcarry64(th, 0ULL, k, c2); } - - // loop ran to completion (never broke early) => the pending carry64 - // (this row's last term's "hi", effectively) lands at i+NA if (i+NA < NR) R[i+NA] = k; - - /*if (j == NA && i + NA < NR) - { - uint8_t carry = 0; - addcarry64(R[i + NA], k, R[i + NA], carry); - for (size_t k = i + NA + 1; k < NR; ++k) - addcarry64(R[k], 0, R[k], carry); - }*/ } - // double the triangular sum uint8_t carry = 0; for (size_t i = 0; i < NR; i++) CppCore::addcarry64(R[i], R[i], R[i], carry); - // add diagonal terms carry = 0; for (size_t i=0, j=0; i 2 * NA) - // CppCore::addcarry64(R[2 * NA], 0, R[2 * NA], carry); } #endif else if constexpr (sizeof(UINT1) % 4 == 0 && sizeof(UINT2) % 4 == 0) From 1e554f0548913fe4d1314e53d2e5d52e58ed7e99 Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sat, 25 Jul 2026 22:08:15 +0200 Subject: [PATCH 12/22] wip --- include/CppCore/Math/Util.h | 34 +++++++++++++++------------------- 1 file changed, 15 insertions(+), 19 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index ae467bf6..ef814b92 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2944,41 +2944,37 @@ namespace CppCore for (size_t k = 0; k < NR; ++k) R[k] = 0; - uint64_t tl, th; - + uint64_t tl, th, k; + uint8_t c; // calculate one triangle for (size_t i = 0; i < NA && 2*i+1 < NR; i++) { - uint64_t k = 0; - for (size_t j = i+1; j < NA; j++) + k = 0ULL; + for (size_t j = i+1; j < NA && i+j < NR; j++) { - size_t idx = i + j; - if (idx >= NR) - break; - CppCore::umul128(A[i], A[j], tl, th); - uint8_t c1 = 0; - CppCore::addcarry64(tl, R[idx], tl, c1); - CppCore::addcarry64(th, 0ULL, th, c1); - uint8_t c2 = 0; - CppCore::addcarry64(tl, k, R[idx], c2); - CppCore::addcarry64(th, 0ULL, k, c2); + c = 0; + CppCore::addcarry64(tl, R[i+j], tl, c); + CppCore::addcarry64(th, 0ULL, th, c); + c = 0; + CppCore::addcarry64(tl, k, R[i+j], c); + CppCore::addcarry64(th, 0ULL, k, c); } if (i+NA < NR) R[i+NA] = k; } // double the triangular sum - uint8_t carry = 0; + c = 0; for (size_t i = 0; i < NR; i++) - CppCore::addcarry64(R[i], R[i], R[i], carry); + CppCore::addcarry64(R[i], R[i], R[i], c); // add diagonal terms - carry = 0; + c = 0; for (size_t i=0, j=0; i Date: Sun, 26 Jul 2026 11:11:35 +0200 Subject: [PATCH 13/22] wip --- include/CppCore/Math/Util.h | 32 +++++++++++++++----------------- 1 file changed, 15 insertions(+), 17 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index ef814b92..a16cbc3b 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2937,44 +2937,42 @@ namespace CppCore constexpr size_t NA = sizeof(UINT1) / 8; constexpr size_t NR = sizeof(UINT2) / 8; assert((void*)&a != (void*)&r); - - const uint64_t* A = reinterpret_cast(&a); - uint64_t* R = reinterpret_cast(&r); - + uint64_t* ap = (uint64_t*)&a; + uint64_t* rp = (uint64_t*)&r; + uint64_t tl, th, k; + uint8_t c; + // TODO: only clear starting at 2NA, write the others before reading for (size_t k = 0; k < NR; ++k) - R[k] = 0; - - uint64_t tl, th, k; - uint8_t c; + rp[k] = 0; // calculate one triangle for (size_t i = 0; i < NA && 2*i+1 < NR; i++) { k = 0ULL; for (size_t j = i+1; j < NA && i+j < NR; j++) { - CppCore::umul128(A[i], A[j], tl, th); + CppCore::umul128(ap[i], ap[j], tl, th); c = 0; - CppCore::addcarry64(tl, R[i+j], tl, c); + CppCore::addcarry64(tl, rp[i+j], tl, c); CppCore::addcarry64(th, 0ULL, th, c); c = 0; - CppCore::addcarry64(tl, k, R[i+j], c); + CppCore::addcarry64(tl, k, rp[i+j], c); CppCore::addcarry64(th, 0ULL, k, c); } if (i+NA < NR) - R[i+NA] = k; + rp[i+NA] = k; } // double the triangular sum - c = 0; + c = 0; //TODO: This doesn't need to go beyond 2NA? for (size_t i = 0; i < NR; i++) - CppCore::addcarry64(R[i], R[i], R[i], c); + CppCore::addcarry64(rp[i], rp[i], rp[i], c); // add diagonal terms c = 0; for (size_t i=0, j=0; i Date: Sun, 26 Jul 2026 12:00:29 +0200 Subject: [PATCH 14/22] wip --- include/CppCore/Math/Util.h | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index a16cbc3b..f94b5a30 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2942,10 +2942,20 @@ namespace CppCore uint64_t tl, th, k; uint8_t c; // TODO: only clear starting at 2NA, write the others before reading - for (size_t k = 0; k < NR; ++k) - rp[k] = 0; + for (size_t i = 0; i < NR; i++) + rp[i] = 0; // calculate one triangle - for (size_t i = 0; i < NA && 2*i+1 < NR; i++) + k = 0ULL; + for (size_t j = 1; j < NA && j < 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++) From 393a382aea3c3b7a2d76e37634523346816df44e Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sun, 26 Jul 2026 12:14:45 +0200 Subject: [PATCH 15/22] wip --- include/CppCore/Math/Util.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index f94b5a30..1e14944c 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2945,7 +2945,7 @@ namespace CppCore for (size_t i = 0; i < NR; i++) rp[i] = 0; // calculate one triangle - k = 0ULL; + /*k = 0ULL; for (size_t j = 1; j < NA && j < NR; j++) { CppCore::umul128(ap[0], ap[j], tl, th); @@ -2954,8 +2954,8 @@ namespace CppCore 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++) + rp[NA] = k;*/ + for (size_t i = 0; i < NA && 2*i+1 < NR; i++) { k = 0ULL; for (size_t j = i+1; j < NA && i+j < NR; j++) From e7fd0fb5c3ee5a09eaf6e5a2efbb67447c5c0b9f Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sun, 26 Jul 2026 12:23:48 +0200 Subject: [PATCH 16/22] wip --- include/CppCore.Test/Math/Util.h | 53 +++++++++++++------------------- 1 file changed, 21 insertions(+), 32 deletions(-) diff --git a/include/CppCore.Test/Math/Util.h b/include/CppCore.Test/Math/Util.h index 0cf57c1b..0b285c43 100644 --- a/include/CppCore.Test/Math/Util.h +++ b/include/CppCore.Test/Math/Util.h @@ -473,47 +473,36 @@ namespace CppCore { namespace Test { namespace Math return true; } - INLINE static bool usquare() + template + INLINE static bool usquare_template() { CppCore::Random::Default64 prng; - - uint64_t a[4]; - uint64_t r0_1[2]; - uint64_t r0_2[2]; - uint64_t r1_1[4]; - uint64_t r1_2[4]; - uint64_t r2_1[8]; - uint64_t r2_2[8]; - uint64_t r3_1[12]; - uint64_t r3_2[12]; - - for (size_t i = 0; i < 100000; i++) + uint64_t a[NA]; + uint64_t r1[NR]; + uint64_t r2[NR]; + for (size_t i = 0; i < I; i++) { prng.fill(a); - // smaller - CppCore::umul(a, a, r0_1); - CppCore::usquare(a, r0_2); - if (!CppCore::equal(r0_1, r0_2)) - return false; - // same size - CppCore::umul(a, a, r1_1); - CppCore::usquare(a, r1_2); - if (!CppCore::equal(r1_1, r1_2)) - return false; - // full wide mul - CppCore::umul(a, a, r2_1); - CppCore::usquare(a, r2_2); - if (!CppCore::equal(r2_1, r2_2)) - return false; - // zero padded - CppCore::umul(a, a, r3_1); - CppCore::usquare(a, r3_2); - if (!CppCore::equal(r3_1, r3_2)) + CppCore::umul(a, a, r1); + CppCore::usquare(a, r2); + if (!CppCore::equal(r1, r2)) return false; } return true; } + INLINE static bool usquare() + { + 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; From cc02490b45d611f1e789fdbfa5a605df2ee5ffdd Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sun, 26 Jul 2026 12:35:27 +0200 Subject: [PATCH 17/22] wip --- include/CppCore.Test/Math/Util.h | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/include/CppCore.Test/Math/Util.h b/include/CppCore.Test/Math/Util.h index 0b285c43..27263de3 100644 --- a/include/CppCore.Test/Math/Util.h +++ b/include/CppCore.Test/Math/Util.h @@ -493,6 +493,15 @@ namespace CppCore { namespace Test { namespace Math 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>() && From f92df7f46aab79a013ed8142bb2200d6b2116fa8 Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sun, 26 Jul 2026 12:42:42 +0200 Subject: [PATCH 18/22] wip --- include/CppCore/Math/Util.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index 1e14944c..ca1ddb62 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2941,11 +2941,11 @@ namespace CppCore uint64_t* rp = (uint64_t*)&r; uint64_t tl, th, k; uint8_t c; - // TODO: only clear starting at 2NA, write the others before reading - for (size_t i = 0; i < NR; i++) - rp[i] = 0; + rp[0] = 0ULL; + for (size_t i = NA+NA; i < NR; i++) + rp[i] = 0ULL; // calculate one triangle - /*k = 0ULL; + k = 0ULL; for (size_t j = 1; j < NA && j < NR; j++) { CppCore::umul128(ap[0], ap[j], tl, th); @@ -2954,8 +2954,8 @@ namespace CppCore CppCore::addcarry64(th, 0ULL, k, c); } if constexpr (NA < NR) - rp[NA] = k;*/ - for (size_t i = 0; i < NA && 2*i+1 < NR; i++) + 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++) From 514551ef3e4f59774a3a2c981406c60fe5c0d6de Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sun, 26 Jul 2026 12:52:07 +0200 Subject: [PATCH 19/22] wip --- include/CppCore/Math/Util.h | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index ca1ddb62..cc9e44a0 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,14 +2899,19 @@ 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)) { From 553bb03e6622549d88ddfb653153847dc7712bdc Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sun, 26 Jul 2026 13:00:18 +0200 Subject: [PATCH 20/22] wip --- include/CppCore/Math/Util.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index cc9e44a0..68096af1 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2952,7 +2952,7 @@ namespace CppCore rp[i] = 0ULL; // calculate one triangle k = 0ULL; - for (size_t j = 1; j < NA && j < NR; j++) + for (size_t j = 1; j < MIN(NA,NR); j++) { CppCore::umul128(ap[0], ap[j], tl, th); c = 0; @@ -2978,8 +2978,8 @@ namespace CppCore rp[i+NA] = k; } // double the triangular sum - c = 0; //TODO: This doesn't need to go beyond 2NA? - for (size_t i = 0; i < NR; i++) + 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; From 30278f6073d4690eb957638617de5ff3e890427b Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sun, 26 Jul 2026 13:14:33 +0200 Subject: [PATCH 21/22] wip --- include/CppCore/Math/Util.h | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index 68096af1..3f730cfd 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2947,12 +2947,12 @@ namespace CppCore uint64_t* rp = (uint64_t*)&r; uint64_t tl, th, k; uint8_t c; - rp[0] = 0ULL; for (size_t i = NA+NA; i < NR; i++) rp[i] = 0ULL; // calculate one triangle - k = 0ULL; - for (size_t j = 1; j < MIN(NA,NR); j++) + 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; @@ -2998,9 +2998,6 @@ namespace CppCore } } - - - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // UNSIGNED DIVISION+MODULO BY CONSTANTS //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// From 355fa98bdaa6f3666841c780e497c1e389e0466d Mon Sep 17 00:00:00 2001 From: Clint Banzhaf Date: Sun, 26 Jul 2026 13:30:53 +0200 Subject: [PATCH 22/22] wip --- include/CppCore/Math/Util.h | 52 ++++++++++++++++++++++++++++++++++++- 1 file changed, 51 insertions(+), 1 deletion(-) diff --git a/include/CppCore/Math/Util.h b/include/CppCore/Math/Util.h index 3f730cfd..56dc8eb2 100644 --- a/include/CppCore/Math/Util.h +++ b/include/CppCore/Math/Util.h @@ -2994,7 +2994,57 @@ namespace CppCore #endif else if constexpr (sizeof(UINT1) % 4 == 0 && sizeof(UINT2) % 4 == 0) { - assert(false); + // 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