11
\$\begingroup\$

My implementation of Karatsuba multiplication from Tim Roughgarden's Algorithms course. I've made a Integer class that holds an integer in string format. I've added operations to this class to add, subtract and multiply numbers. Karatsuba.cpp performs Gauss' trick (complex multiplication algorithm).

Integer.hpp

#pragma once
#include <string>
#include <functional>

class Integer{
    public:
        Integer();
        Integer(const std::string& input);
        ~Integer();
        Integer operator+(const Integer& other) const;
        Integer operator-(const Integer& other) const;
        Integer operator*(const Integer& other) const;
        size_t size() const{return fString.size();}
        void padRight(size_t num);
        void padLeft(size_t num);
        Integer substr(size_t pos, size_t length) const;
        void print() const;
        void stripLeadingZeros();
    private:
        std::string fString;
};
void equalizeLengthsPadLeft(Integer& first, 
                            Integer& second);
void equalizeLengthsPadRight(std::string& first, std::string& second);

Integer.cpp

#include "Integer.hpp"
#include <assert.h>
#include <string>
#include <algorithm>
#include <iostream>

namespace {
    int char2int(char c){
        return c - '0';
    }
}
Integer::Integer()
:
fString()
{
}

Integer::Integer(const std::string& input)
:
fString(input)
{
}

Integer::~Integer(){}

Integer Integer::substr(size_t pos, size_t len) const{
    return fString.substr(pos, len);
}

void equalizeLengthsPadLeft(Integer& first, 
                            Integer& second){
    if (first.size() < second.size()){
        first.padLeft(second.size()-first.size());
    }
    else if(first.size() > second.size()){
        second.padLeft(first.size() - second.size());
    }
}
void equalizeLengthsPadRight(std::string& first, std::string& second){
    if (first.size() < second.size()){
        first += std::string(second.size()-first.size(), '0');
    }
    else if(first.size() > second.size()){
        second += std::string(first.size() - second.size(), '0');
    }
}

Integer Integer::operator+(const Integer& other) const{
    // For the time being, just conver to integer and return
    std::string first = fString;

    std::reverse(first.begin(), first.end());
    std::string second = other.fString;
    std::reverse(second.begin(), second.end());
    equalizeLengthsPadRight(first,second);

    std::string::iterator first_it = first.begin();
    std::string::iterator second_it = second.begin();
    std::string resultStr;
    int carry = 0;
    while(first_it != first.end()){
        int sum = char2int(*first_it) + char2int(*second_it) + carry;
        carry = 0;
        if (sum >= 10){
            sum = sum%10;
            carry = 1;
        }
        resultStr += std::to_string(sum);        
        first_it++;second_it++;    
    }
    if (carry){
        resultStr += '1'; 
    }
    std::reverse(resultStr.begin(), resultStr.end());
    return resultStr;
}

Integer Integer::operator-(const Integer& other) const{

    std::string first = fString;
    std::reverse(first.begin(), first.end());
    std::string second = other.fString;
    std::reverse(second.begin(), second.end());

    // Equalize
    equalizeLengthsPadRight(first,second);

    std::string::iterator first_it = first.begin();
    std::string::iterator second_it = second.begin();
    std::string resultStr;
     while(first_it != first.end()){
         int up = char2int(*first_it);
         int down = char2int(*second_it);
         int localResult;
         if (up >= down){
            localResult = up-down;
         }
         else{
             // Keep searching forward until you get a non-zero value
            auto next_it = first_it+1;
            while(true){
                if (char2int(*next_it) > 0){
                    // Found the first non-zero number
                    break;
                }
                next_it++;
            }
            *next_it = std::to_string(char2int(*next_it)-1)[0];
            // Now chase back to first_it setting 9's 
            // on the way. Make sure everything was 0 
            // beforehand
            next_it--;
            while(next_it != first_it){
                assert(char2int(*next_it) == 0);
                *next_it = std::to_string(9)[0];
                next_it--;
            }

            localResult = up+10 -down;
         }
         assert(localResult>=0);
         resultStr += std::to_string(localResult);
        first_it++;
        second_it++;
    }
    std::reverse(resultStr.begin(), resultStr.end());
    return resultStr;

}

Integer Integer::operator*(const Integer& other) const {
    // Only allow multiplication when size is 1
    assert(size() == other.size() == 1);
    return std::to_string(std::stoi(fString)*std::stoi(other.fString));
}   

void Integer::padRight(size_t num){
    fString += std::string(num, '0');
}

void Integer::padLeft(size_t num){
    fString.insert(0,num,'0');
}

void Integer::print() const{
    std::cout << fString << std::endl;
}

void Integer::stripLeadingZeros(){
    // Don't strip if all are zeros - this will lead to an empty string
    if (std::all_of(fString.cbegin(), fString.cend(), [](char c){return ('0'== c); })){
        return;
    }

    fString.erase(0, fString.find_first_not_of('0'));
}

Karatsuba multiplication

#include <string>
#include <assert.h>
#include <cmath>
#include "Integer.hpp"

Integer multiply(const Integer& inp1, const Integer& inp2){
    Integer first = inp1;
    Integer second = inp2;

    equalizeLengthsPadLeft(first, second);

    assert(first.size()==second.size());
    size_t sz = first.size();

    // Base case 
    if (sz == 1){
        return first*second;
    }

    int n = sz/2;
    Integer A = first.substr(0,n);
    Integer B = first.substr(n, sz-n);
    Integer C = second.substr(0,n);
    Integer D = second.substr(n, sz-n);

    Integer AC = multiply(A, C);
    Integer BD = multiply(B, D);

    Integer A_plus_B = A+B;
    Integer C_plus_D = C+D;
    Integer sum =  multiply(A_plus_B, C_plus_D);
    Integer AD_plus_BC = sum - AC - BD;

    AC.padRight(2*(sz-n));
    AD_plus_BC.padRight(sz-n);

    Integer result = AC+ AD_plus_BC + BD;
    result.stripLeadingZeros();

    return result;

}

int main(){
    Integer first("3141592653589793238462643383279502884197169399375105820974944592");
    Integer second("2718281828459045235360287471352662497757247093699959574966967627");

    Integer ans = multiply(first, second);

    ans.print();
}
\$\endgroup\$
  • \$\begingroup\$ What is a Gauss trick, and how Karatsuba.cpp performs it? \$\endgroup\$ – vnp Feb 18 at 5:06
  • \$\begingroup\$ Karatsuba is div-conq multiplication. If \$x\$ and \$y\$ are \$n\$-digit strings in some base \$B\$, then for any positive integer \$m < n\$, you can write \$x = aB^m + b, y = cB^m + d\$ where \$b\$ and \$d\$ are less than \$B^m\$. The product \$x \cdot y = \ldots =acB^{2m} + (ad + bc)B^m + bd\$. Applying Gauss' complex multiplication algorithm, reuse the results of \$ac\$ and \$bd\$ to find \$(ad + bc)\$. Thus, \$(ad + bc) = (a + b)(c + d) - ac - bd\$, reducing the no. of recursive multiplications at the cost of extra additions. This is performed when AD_plus_BC is introduced in the code. \$\endgroup\$ – Snowhawk Feb 18 at 6:07
13
\$\begingroup\$

I could be wrong, but looking at Integer I think it is really NaturalNumber. It doesn't seem to support negative numbers.

There are a lot of reversals going on. It would almost certainly make more sense to use the reversed (little-endian) number internally, and only actually call std::reverse when converting between strings and Integers.

padLeft and stripLeadingZeros do not belong in the public API of the class (and I'm not even sure they should be required in the private API). An Integer should be an integer, not a string. padRight would be more appropriately named something like multiplyByPowerOfTen or shiftLeftDecimal.

If the Karatsuba multiplication is intended to be practical rather than purely an academic exercise, the base case should be when the multiplication can be done in a native type: probably when both numbers are less than 10 ** 9, so they fit in a (signed) 32-bit type and the result in a (signed) 64-bit type.

\$\endgroup\$
  • \$\begingroup\$ Absolutely right about renaming Integer to NaturalNumber. The assignment assumed positive numbers and that doesn't get implied anywhere in the code. Extending the class a little bit will open it to negative numbers too. I could store the number in a reverse manner but it makes it harder to reason about the code. When constructing an integer from a raw string, I could perhaps just reverse the input and store it in fString. On further thought, only 2 operations - operator+ and operator- need the string in reverse order, maybe I can store another reversed string as a member \$\endgroup\$ – cppprogrammer Feb 18 at 14:29
  • \$\begingroup\$ Only operator+ and operator- explicitly reverse, but the Karatsuba multiplication uses those operators... \$\endgroup\$ – Peter Taylor Feb 18 at 14:35
8
\$\begingroup\$

Your code is formatted inconsistently. Sometimes you placed spaces around operators, sometimes not, with no apparent patterns or rules. Just let your IDE format the code automatically.

Sometimes you use snake_case, sometimes camelCase. You should only use one of them, for consistency.

Instead of sum = sum%10 you can write sum -= 10 since the sum can never be greater than 19.

Having operator* restricted to a single digit is unexpected. That function should better be an internal implementation detail, and multiply should become the new operator*.

Instead of Integer::print, you should define the usual os << Integer operator. Don't use std::endl unless you really know you need it. A simple '\n' performs better.

\$\endgroup\$
8
\$\begingroup\$

I miss comments and use of a documentation tool like doxygen.

  • I did not expect separate digit handling loops in operator+() and operator-() (obviously, I expected handling of signed integers). Nor the asymmetry between the procedures.

In multiply(),
(I don't quite like sum for what it names (still undecided about sz):)

  • sz-n occurs just often enough to ponder introducing an abstraction, a name, a variable for it.
  • "Padding before adding" leads to two 2×sz additions: an alternative would be to exclude the least significant part (say, bd) of BD from addition, "appending" its most significant part to AC for a single sz+n addition and just append bd.

(The remaining remarks may lead to more code, code more difficult to maintain/read:
revisiting the test harness seems a good idea before coding any.)

  • The next step in avoiding complex digit operations would seem to be not padding in operator±() and using increment()/decrement() instead.)
  • With factors that differ widely in order of magnitude, there may be a lot of padding & stripping: consider an inner procedure doing just *three multiplications and five additions" (possibly with operand lengths passed in) and an outer one taking care of interface (like stripping leading zeroes).

(On second thought, I'm not sure where "tactics" should go: passing the shorter factor where it terminates recursion sooner, extra strategy ("long multiplication") when one factor reaches the base case, handling 0 and powers of "the base" separately (0th power = 1 being a prominent and trivial case).)

\$\endgroup\$
  • \$\begingroup\$ "Nor the asymmetry between the procedures." - How can this be improved? Aren't addition and subtraction fundamentally different? Inside the subtraction routine, trying to code the borrowing from the next power of 10 gave me a headache. Addition was a breeze \$\endgroup\$ – cppprogrammer Feb 18 at 14:41
  • \$\begingroup\$ A question looking so innocuous (just learned that I got the pronunciation wrong for decades) about something that tedious to implement. Faced with implementing both for integers in binary representation, someone came up with two's complement. Should work with any base greater than one. (Using a borrow in place of carry, I don't envision "decimal" subtraction any more difficult to implement on its own than addition?!) \$\endgroup\$ – greybeard Feb 19 at 8:01

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

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