6
\$\begingroup\$

This is the Karatsuba algorithm I wrote:

#include <stdio.h>
#include <math.h>

long long karatsuba(long long, long long);
long long findn(long long);

int main() {
    long long x = 1234, y = 1234;
    long long result;

    result = karatsuba(x,y);

    printf("%lld * %lld = %lld\n", x, y, result);
}

long long karatsuba(long long x, long long y) {

    int n = fmax(findn(x), findn(y));
    int n2 = floor(n/2);

    if ( n == 1 )  return x * y;

    long long a = x / pow(10, n2);
    long long b = x % (long long) pow(10, n2);
    long long c = y / pow(10, n2);
    long long d = y % (long long) pow(10, n2);
    long long p = a + b;
    long long q = c + d;

    long long ac = karatsuba(a, c);
    long long bd = karatsuba(b, d);
    long long pq = karatsuba(p, q);

    long long adbc = pq - ac - bd;

    long long result = (long long) pow(10, n2*2) * ac + (long long) pow(10, n2) * adbc + bd;
    return result;

}

long long findn(long long x) {
    long long n = 1;

    while ( x/10 >= 1 ) {
        ++n;
        x /= 10;
    }
    return n;
}

Inspired by some comments on my post about the same algorithm on Stack Overflow in which I asked help for a modified version I wanted to code, I'd like to submit the original working algorithm.

I understand it would be better to ditch the pow() function, or working with binary numbers, but I'm not sure.

I'd really appreciate some guidance on how to improve performance.

\$\endgroup\$
1
  • 1
    \$\begingroup\$ int main() is as fine as it is in C23. So no comment on it. Though, you could define main() at the end of the function and avoid having to specify prototypes for other functions twice. \$\endgroup\$
    – Harith
    Commented May 18 at 15:28

1 Answer 1

6
\$\begingroup\$

Avoid mixing integer and floating point mathematics

You are mixing operations with different types of values. This can lead to subtle problems. C doesn't define exactly what precision long long and double have, this depends on the platform. Most likely, long long has a higher precision than double, so conversions between them could cause rounding errors. Apart from that, those conversions are also not free, and will reduce the performance of your code.

Luckily, you don't need any of the floating point operations in your code. fmax() is replaced by a simple if-statement, floor() wasn't necessary at all, and an integer pow() can be implemented quite easily.

Choice of base

You are implementing the Karatsuba algorithm in base 10. Unfortunately, computers don't work in base 10, they work in powers of 2. The simplest is to use base 2 (also known as binary). Your code would then look like:

long long a = x >> n2;
long long b = x - (a << n2);
long long c = y >> n2;
long long d = y - (c << n2);

And findn() can be replaced by ffs() or a similar function that counts the size of a number in bits. This then completely avoids using division and modulo operations, which are typically quite expensive for a CPU to perform.

However, you don't need to restrict yourself to binary, you can also work in base 16 (hexadecimal), base 256 (so each byte is a digit), or in the extreme, base \$2^{64}\$, which means each digit is a long long (on most 64-bit architectures). Of course, that means that if your inputs are just two long long numbers, your algorithm will just immediately return x * y. This will tell you that the Karatsuba algorithm is not very useful for any numbers that fit inside long longs. The algorithm does start to make a lot of sense though when you have much larger numbers, which have to be stored in multiple long longs. In that case, you make use of the fast long long operations the CPU can perform itself, and the algorithm will ensure you can expand this to numbers of arbitrary size.

\$\endgroup\$
10
  • \$\begingroup\$ Which function from stdc_first_leading_one*()/stdc_first_trailing_one*() is equivalent to this ffs()? \$\endgroup\$
    – Harith
    Commented May 18 at 15:12
  • 2
    \$\begingroup\$ Ah, if you can use C23, you'd use stdc_bit_width*(). \$\endgroup\$
    – G. Sliepen
    Commented May 18 at 15:33
  • \$\begingroup\$ @G. Sliepen, your answer seems really helpful, i'll be sure to reading it carefully. I guess that if I wanted to implement it in base 2 i would also need a function to perform operations on binary numbers ? Suppose x = 1010 (10 in decimal), y = 1000 (8 in decimal); eventually the recursion would calculate 2^4 * 100 + 2^2 * 100 + 0 which equals to 80 but the program would calculate2000, since it would count 100 as decimal and not binary. Is this correct? \$\endgroup\$
    – Talete
    Commented May 18 at 16:54
  • 1
    \$\begingroup\$ Operations are done on values, and while values can either be represented in decimal, binary or any other base, the resulting value of an operation will not be affected by the representation. \$\endgroup\$
    – G. Sliepen
    Commented May 18 at 17:48
  • 1
    \$\begingroup\$ Aside: "Computers don't work in base 10, they work in powers of 2." --> C2x offers optional decimal floating point types. \$\endgroup\$ Commented May 19 at 22:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.