In Visual Studio C++, I am trying to multiply two uint64_t numbers stored as little-endian in arrays. Unfortunately, there are two small carry propagation errors which I can't figure out where the problem is:

0x0306f7b285eead7d8cc88407a9f4c002 * 0xf9084d7a11528273377bf8560b3ffe =

0x2f1e00eebe083ca406e838b64ae89188195f931cca21eaaa66062c82cfffc correct
                              #               #
0x2f1e00eebe083ca406e838b64ae88188195f931cca21faaa66062c82cfffc my code

Here is my code:

#include <stdint.h>      // For uint64_t, uint8_t, ptrdiff_t
#include <stddef.h>      // For size_t, ptrdiff_t
#include <string.h>      // For memcpy, memset
#include <stdio.h>       // For printf
#include <intrin.h>      // For _umul128, _addcarry_u64

// Define uint128_t type using struct
typedef struct {
    uint64_t lo;
    uint64_t hi;
} uint128_t;

// Function to multiply two uint64_t numbers and get a uint128_t result
uint128_t multiply_uint64(uint64_t a, uint64_t b) {
    uint128_t result;
    result.lo = _umul128(a, b, &result.hi);
    return result;
}

size_t len = 2;
memset(result, 0, 2 * len * sizeof(uint64_t));
    for (size_t i = 0; i < len; ++i) {
        uint64_t carry = 0;
        for (size_t j = 0; j < len; ++j) {
            // Multiply x[i] and y[j]
            uint128_t product = multiply_uint64(x[i], y[j]);

            // Add product.lo to result[i + j] with carry
            uint8_t c1 = _addcarry_u64(0, product.lo, result[i + j], &result[i + j]);
            // Add carry from previous addition to result[i + j]
            uint8_t c2 = _addcarry_u64(c1, result[i + j], carry, &result[i + j]);

            // Update carry for next iteration
            uint64_t temp_carry = product.hi + c2;
            carry = temp_carry;
        }
        result[i + len] = carry;
    }

In Visual Studio C++, I am trying to multiply two uint64_t numbers stored as little-endian in arrays. Unfortunately, there are two small carry propagation errors which I can't figure out where the problem is:

0x0306f7b285eead7d8cc88407a9f4c002 * 0xf9084d7a11528273377bf8560b3ffe =

0x2f1e00eebe083ca406e838b64ae89188195f931cca21eaaa66062c82cfffc correct
                              #               #
0x2f1e00eebe083ca406e838b64ae88188195f931cca21faaa66062c82cfffc my code

Here is my code:

#include <stdint.h>      // For uint64_t, uint8_t, ptrdiff_t
#include <stddef.h>      // For size_t, ptrdiff_t
#include <string.h>      // For memcpy, memset
#include <stdio.h>       // For printf
#include <intrin.h>      // For _umul128, _addcarry_u64

// Define uint128_t type using struct
typedef struct {
    uint64_t lo;
    uint64_t hi;
} uint128_t;

// Function to multiply two uint64_t numbers and get a uint128_t result
uint128_t multiply_uint64(uint64_t a, uint64_t b) {
    uint128_t result;
    result.lo = _umul128(a, b, &result.hi);
    return result;
}

size_t len = 2;
memset(result, 0, 2 * len * sizeof(uint64_t));
    for (size_t i = 0; i < len; ++i) {
        uint64_t carry = 0;
        for (size_t j = 0; j < len; ++j) {
            // Multiply x[i] and y[j]
            uint128_t product = multiply_uint64(x[i], y[j]);

            // Add product.lo to result[i + j] with carry
            uint8_t c1 = _addcarry_u64(0, product.lo, result[i + j], &result[i + j]);
            // Add carry from previous addition to result[i + j]
            uint8_t c2 = _addcarry_u64(c1, result[i + j], carry, &result[i + j]);

            // Update carry for next iteration
            uint64_t temp_carry = product.hi + c2;
            carry = temp_carry;
        }
        result[i + len] = carry;
    }

• c++
• visual-studio
• bignum

• You may want to look into using a big integer library like GMP. – Jesper Juhl Commented Jan 17 at 19:47
• 3 I can't find out where is the problem -- Do the multiplication first by hand. Then compare those intermediate values with the values such as product, c1, by using a debugger. – PaulMcKenzie Commented Jan 17 at 19:56
• 1 On a side note: <stdint.h> <stddef.h> <string.h> <stdio.h> should be <cstdint> <cstddef> <cstring> <cstdio> in C++ – Remy Lebeau Commented Jan 17 at 20:01
• Imagine you have struct uint64_t { uint32_t lo; uint32_t hi; };, how would you multiply two of those together? Now think of it as struct uint128_t { uint32_t lo; uint32_t mlo; uint32_t mhi; uint32_t hi; };, how would you multiple two of those together? (And can you cheat, by using unsigned __int128?) – Eljay Commented Jan 17 at 20:05
• It's just a small piece of code which is Lucas-Lehmer Primality Testing. I am trying to prove 2^521-1 is a Marsenne Prime. At the end I will try to prove 2^136279841 - 1 for which you need just 17.5 MB to store it in uint64_t chunks ;-). – Caster Commented Jan 17 at 21:45

1 Answer 1

Fixed code by openai o1 preview free:

for (size_t i = 0; i < len; ++i) {
    uint64_t carry = 0;
    for (size_t j = 0; j < len; ++j) {
        // Multiply x[i] and y[j]
        uint128_t product = multiply_uint64(x[i], y[j]);

        // Add product.lo to result[i + j]
        uint64_t sum1;
        uint8_t c1 = _addcarry_u64(0, result[i + j], product.lo, &sum1);

        // Add carry to sum1
        uint64_t sum2;
        uint8_t c2 = _addcarry_u64(0, sum1, carry, &sum2);

        result[i + j] = sum2;

        // Compute new carry
        carry = product.hi + c1 + c2;
    }
    result[i + len] = carry;
}