Karatsuba Implementation Issues

Question

Couple things:

• I am a university student
• This is not a homework problem (I'm curious about how bignum libraries work)
• All the code snippets here can be compiled with gcc -lgmp

Here is my implementation which keeps returning incorrect answers. I've tried going through with gdb & checking individual ops with GMP and marking which operations are definitely failing in the code, but I haven't managed to find specifically what's going wrong.

Beyond that, I can't identify where the issue is.

uint64_t *multiply(u256 c, u256 a, uint64_t alen, u256 b, uint64_t blen) {
    if (alen == 1 || blen == 1) {
        u2_set(c, 0);
        c[0] = a[0] * b[0];
        c[1] = c[0] >> 32LU;
        c[0] &= MASK;
        return c;
    } else {
        uint64_t al, bl;
        uint64_t M = min(alen, blen);
        al = ceil(alen, 2);
        bl = ceil(blen, 2);
        uint64_t m = M / 2;
        u256 low1, low2, high1, high2, z0, z1, z2;
        u2_set(low1, 0);
        u2_set(low2, 0);
        u2_set(high1, 0);
        u2_set(high2, 0);
        memcpy(low1, a, sizeof(uint64_t) * al);
        memcpy(low2, b, sizeof(uint64_t) * bl);
        memcpy(high1, &(a[al]), sizeof(uint64_t) * (alen - al));
        memcpy(high2, &(b[bl]), sizeof(uint64_t) * (blen - bl));

        /* Starts ok, goes bad immediately. Mostly ok. */
        multiply(z0, low1, al, low2, bl);

        u2_add(low1, low1, high1); /* CHECKED */
        u2_add(low2, low2, high2); /* CHECKED */

        multiply(z1, low1, al, low2, bl); /* Starts ok, goes bad immediately. */

        /* Bad but starts & usually is ok */
        multiply(z2, high1, (alen - al), high2, (blen - bl));

        /* ALL GOOD */
        u256 o;
        u2_sub(z1, z1, z2);
        u2_sub(z1, z1, z0);
        u2_shl(z1, z1, m);
        u2_shl(z2, z2, m * 2);
        u2_add(o, z2, u2_add(z1, z1, z0));

        u2_copy(c, o);
        return c;
    }
}

To be honest, the worst part is not having test vectors for the intermediate values. Beyond the ones in the wikipedia page example (which I seem to pass), I can't find anything.

u256 is an array of 8 uint64_t's. I'm only using the lower 32 bits of those integers to make masking and carries a bit easier. The idea is to do the Karatsuba multiplication base 2³².

The rest of the code can be found here.

Any help would be appreciated.

EDIT:

Example of what's wrong:

a = 89637787616193205010487395076171907549998432494691382753495063501471509582488
b = 62570744880332518477645686292054860709644349224438728618405612469614908709348

Multiplying the 2 should return

16439926136035242788212517705012640494860769275971090459163708006828967815008

Instead it returns

38696864898835518778416

That specific test can be run by compiling this code: here

with gcc and -lgmp

ceil is defined as so:

#define ceil(a, b) (((a) + (b) - 1LU) / (b))

using this test code, which checks stuff a bit more in-depth, I believe I narrowed the problem down some.

It occurs on the m=4 multiply of z0 = low1 * low2. The two values multiplied are 100000000000 and 100000000000, and the expected result is 10000000000000000000000 while it ends up being 1478053504202008644. Running the above test code will reproduce this. I assume this means my karatsuba implementation is simply broken, as the components aren't throwing any errors and it's specifically failing on the m=2 multiplications. I'm checking the m=1 multiplications and they seem fine.