Multiplication of 2 64-digit numbers in C++

Question

I have karatsuba multiplication algorithm implemented. I want to improve it in this way that I can multiply 2 64-digit numbers but I don't know how to do this. I was given a hint that both numbers contain number of digits that is a power of two but it suggests me nothing. Could you give any other hints? Either math hint or algorithm improvement hint.

#include <iostream>
#include <math.h>
using namespace std;

int getLength(long long value);
long long multiply(long long x, long long y);

int getLength(long long value)
{
    int counter = 0;
    while (value != 0)
    {
        counter++;
        value /= 10;
    }
    return counter;
}

long long multiply(long long x, long long y)
{
    int xLength = getLength(x);
    int yLength = getLength(y);

    // the bigger of the two lengths
    int N = (int)(fmax(xLength, yLength));

    // if the max length is small it's faster to just flat out multiply the two nums
    if (N < 10)
        return x * y;

    //max length divided and rounded up
    N = (N / 2) + (N % 2);

    long long multiplier = pow(10, N);

    long long b = x / multiplier;
    long long a = x - (b * multiplier);
    long long d = y / multiplier;
    long long c = y - (d * N);

    long long z0 = multiply(a, c);
    long long z1 = multiply(a + b, c + d);
    long long z2 = multiply(b, d);


    return z0 + ((z1 - z0 - z2) * multiplier) + (z2 * (long long)(pow(10, 2 * N)));

}

int main()
{
    long long a;
    long long b;
    cin >> a;
    cout << '\n';
    cin >> b;
    cout << '\n' << multiply(a, b) << endl;
    return 0;
}

Answer 1

Here's a hint:

(A + kB) * (C + kD) = AC + k(BC + AD) + k^2(BD)

This helps if k is a power of the base you are keeping your numbers in. For example, if k is 1'000'000'000 and your numbers are base-10, then the multiplications by k are done by simply shifting the numbers around (adding zeroes.)

Anyways, consider breaking your 64-digit numbers in two parts, each 32 digits and do the math like the above. To calculate AC, BC, AD, and BD you are multiplying a pair of 32-digit numbers which can be done similarly.

Since your number of digits is a power of two, you can keep breaking your numbers in half until you get to number sizes that are manageable (e.g. 1-digit numbers.)

BTW, it's not clear from your question whether you're talking about 64 bits or 64 decimal digits. If all you're looking for is multiplying 64-bit numbers, just do this:

// I haven't actually run this code, so...

typedef unsigned long long u64;

u64 high32 (u64 x) {return x >> 32;}
u64 low32  (u64 x) {return x & 0xFFFFFFFF;}

u64 add_with_carry (u64 a, u64 b, u64 * carry)
{
    u64 result = a + b;
    *carry = (result < a ? 1 : 0);
    return result;
}

void mul (u64 a, u64 b, u64 * result_low, u64 * result_high)
{
    u64 a0 = low32(a), a1 = high32(a);
    u64 b0 = low32(b), b1 = high32(b);

    u64 a0b0 = a0 * b0;
    u64 a0b1 = a0 * b1;
    u64 a1b0 = a1 * b0;
    u64 a1b1 = a1 * b1;

    u64 c0 = 0, c1 = 0;
    u64 mid_part = add_with_carry(a0b1, a1b0, &c1);

    *result_low  = add_with_carry(a0b0, (low32(mid_part) << 32, &c0);
    *result_high = high32(mid_part) + a1b1 + (c1 << 32) + c0; // this won't overflow
}

This implementation is the same idea outlined above. Since in standard C/C++, the largest number of meaningful bits we can have in the result of a multiplication is 64, then we can only multiply two 32-bit numbers at a time. Which is what we are doing here.

The final result will be 128 bits, which we return in two unsigned 64-bit numbers. We are doing a 64-bit by 64-bit multiply, by doing 4 32-bit multiplies and a few adds.

As a side note, this is one of those few case where writing this is assembly is usually waaaaay easier than C. For example, in x64 assembly, this is literally a single mul instruction that multiplies two 64-bit unsigned integers and return the 128-bit result in two 64-bit registers.

Even if you don't have a 64-bit to 128-bit multiplication instruction, still writing this in assembly is easier (because of adc or similar instructions, e.g the whole code above is just two movs, four muls, four adds, and four adcs in plain x86 assembly.) Even if you don't want to write in assembly, you should check your compiler's "intrinsics". It might have one for a large multiplication (but you'll be platform-dependent.)

Answer 2

To apply Karatsuba or any other multiplication on lower bit width arithmetics you need to divide your number into smaller "digits". Before anything else you need to access these "digits" so here is how to do it:

you got number 1234 and want to divide it to 10^1 digits so

1234 = 1*1000 + 2*100 + 3*10 + 4

you can obtain the digits like this:

x=1234;
a0=x%10; x/=10; // 4
a1=x%10; x/=10; // 3
a2=x%10; x/=10; // 2
a3=x%10; x/=10; // 1

if you want 10^2 digits then:

x=1234;
a0=x%100; x/=100; // 34
a1=x%100; x/=100; // 12

Now the problem is that to do this you need to have division on the full number which you do not. If you got the number as string then it is easily done but let assume you do not. Computers are based on binary computations so it is a good idea to use power of 2 as base for "digits" so:

x = 1234 = 0100 1101 0010 bin

Now if we want to have for example 2^4=16 base digits then:

a0=x%16; x/=16; // 0010
a1=x%16; x/=16; // 1101
a2=x%16; x/=16; // 0100

Now if you realize that dividing by power of 2 is just bit shift right and remainder can be expressed as AND then:

a0=x&15; x>>=4; // 0010
a1=x&15; x>>=4; // 1101
a2=x&15; x>>=4; // 0100

The bit shift can be stacked to any bitwidth numbers so now you got all you need. But that is not all if you select as "digit" for example 2^8 which is BYTE then you can use pointers instead for example:

DWORD x=0x12345678; // 32 bit number
BYTE *db=(BYTE*)(&x); // 8bit pointer that points to x

a0=db[0]; // 0x78
a1=db[1]; // 0x56
a2=db[2]; // 0x34
a3=db[3]; // 0x12

so you can directly access the digits or reconstruct the x from digits:

DWORD x; // 32 bit number
BYTE *db=(BYTE*)(&x); // 8bit pointer that points to x
db[0]=0x78;
db[1]=0x56;
db[2]=0x34;
db[3]=0x12;
// here x should be 0x12345678

Beware that the order depends on platform MSB or LSB first order. Now you can apply multiplication. For example 32*32=64 bit done on 16 bit multiplication is done like this with naive O(n^2) approach:

x(a0+a1<<16) * y(b0+b1<<16) = a0*b0 + a0*b1<<16 + a1*b0<<16 + a1*b1<<32

Where a0,a1,b0,b1 are the digits of operands. Note that result of each ai*bj multiplication is 2 digit wide so you need to split it into digits and store to the result digit addressed by the bit shift. Beware that adding can cause overflow to higher digit. To handle that you need either to use at least twice the arithmetics width for addition (16*16 bit mul -> 32bit add) or use carry flag. Sadly apart using assembly in C++ you have no access to the Carry flag. luckily it can be simulated see:

• Cant make value propagate through carry

And now you can construct your Karatsuba or more advanced multiplications for more info see:

• Fast bignum square computation