1

I am trying to implement the recursive Karatsuba algorithm for multiplication of two polynomials (of the same degree). My code still doesn't work for polynomials with degree greater than 1. Can anyone explain to me what I am doing wrong in my function? It should work for even and also for odd numbers of coefficients.

The polynomials are stored in an array of long. Every element represents one coefficient, so 1 2 3 4 means 1 + 2x + 3x^2 + 4x^3 and the argument size is the number of coefficients.

long* karatsuba(const long *A, const long *B, int size) {
    long *lowA, *highA, *lowB, *highB, *midA, *midB;

    if (size == 1)
        return naive(A, B, size, size);

    int half = size / 2;

    lowA = new long[half];
    lowB = new long[half];
    midA = new long[half];
    midB = new long[half];
    highA = new long[half];
    highB = new long[half];

    // init low coefficients
    for(int i=0; i<half; i++){
        lowA[i] = A[i];
        lowB[i] = B[i];
    }

    // init high // init low coefficients
    for(int i=half; i<size; i++){
        highA[i-half] = A[i];
        highB[i-half] = B[i];
    }

    // init mid // init low coefficients
    for(int i=0; i<half; i++){
        midA[i] = lowA[i] + highA[i];
        midB[i] = lowB[i] + highB[i];
    }

    long *z0 = karatsuba(lowA, lowB, half);
    long *z1 = karatsuba(midA, midB, half);
    long *z2 = karatsuba(highA, highB, half);

    // compute the result
    long *result = new long[2*size-1];
    for(int i=0; i<half; i++){
        result[i + 2*half] = z2[i];
        result[i + half] = z1[i] - z0[i] - z2[i];
        result[i] = z0[i];
    }
    return result;
}

My main problem is that this function doesn't compute the correct results.

  • 1
    Heed my advice: Use [`std::vector`](http://en.cppreference.com/w/cpp/container/vector). I don't even want to know how many memory leaks are in here... – Arnav Borborah Mar 14 '18 at 21:58
  • 1
    thank you for your advice but I have to use arrays like these, because I will have to work with vector registers and parallelization with OpenMP. –  Mar 14 '18 at 22:03
  • Welcome to StackOverflow. Please read and follow the posting guidelines in the help documentation. [Minimal, complete, verifiable example](http://stackoverflow.com/help/mcve) applies here. We cannot effectively help you until you post your MCVE code and accurately describe the problem. We should be able to paste your posted code into a text file and reproduce the problem you described. "Doesn't work" is not a problem specification. – Prune Mar 14 '18 at 22:13
  • In C++11, `vector::data()` will return a pointer to the underlying data. – David Eisenstat Mar 15 '18 at 01:26
  • how about debug this thing .... on a specific input (not too big so you do not step hundred times over) and compare with sub-result on paper check for first inconsistency.... – Spektre Mar 15 '18 at 08:12
  • `int half = size / 2` may be the problem. What if `size` is not even? Looks like you missed some digits or overflow somewhere (loops that fill arrays). Did you tried with original data of size 2^n? It may work well. – Jean-Baptiste Yunès Mar 15 '18 at 08:29

1 Answers1

0

I noticed one-and-a-half issues concerning correctness:

  1. the index in "the combination loop" is limited to the recursive calls parameter length instead of the result length (almost twice the parameter length)
  2. when the limit is correct, the intermediate results "overlap", necessitating accumulation instead of assignment.

I never used "modern C++" "in anger" - if this works, but leaves you dissatisfied with the code, consider posting your concerns on Code Review.

greybeard
  • 2,249
  • 8
  • 30
  • 66