Multiplying every element of one array by every element of another array and sort the new very large array

Question

Disclaimer This is an exercise of my course, not from an on going contest.

Problem description

The problem description is very straight forward:

You are given two arrays, A and B, containing n and m elements, correspondingly. The numbers which you need to sort are Ai*Bj , for 1 <= i <= n and 1 <= j <= m. In simple words, every element of the first array should be multipled by every element of the second array.

Let C be a result of this sorting, being a non-decreasing sequence of element. Print the sum of every tenth element of this sequence, that is, C1 + C11 + C21 + ... .

1 <= n,m <= 6000

1 <= Ai,Bj <= 40000

Memory limit: 512MB

Time limit: 2 seconds

My solution so far

First I use Java, using Arrays.sort, given the largest n,m. We will need to sort an array with size of 36000000. Then go through every tenth element in the array to get the sum. This passes 23 test cases, and the rest got TLE.

Then I switch to C++, also use the builtin sort method, and the result is just a little bit better, passes 29 test cases.

My observation

Given this input

4 4
7 1 4 9
2 7 8 11

If we sort two array A and B first then multiply them together, we got

2 8 14 18 7 28 49 63 8 32 56 72 11 44 77 99

which is an array with m sorted subarrays. But I couldn't think of any good solution to merge all of these sorted subarray in O(mn) or somewhere around that. Or we need to look at the problem from a different angle, is there any special properties involve with multiplying every elements of two array together?

Update 1: - using MinHeap - not fast enough. [TLE]

Update 2: - using k ways merge - still not fast enough. [TLE]

Update 3: - I forgot to mention about the range of elements in A and B so I've just updated it.

Update 4: - Radix sort base 256 [Accepted]

Conclusion

Through out this problem, I know more about sorting in general and some useful information of sorting with libraries in Java and C++.

• Builtin sort methods in C++ like std::sort is not stable because it is basically a quicksort but when the data format is not favorable for quicksort, then it switches to merge sort, but in general it is the fastest builtin sort of C++ (beside qsort, stable_sort).

• For Java, there are 3 types of sort, one with Arrays.sort(primitive[]), which uses merge sort under the hood, Arrays.sort(Object[]) which uses Timsort and Collections.sort which basically calls Arrays.sort to do its heavy processing stuff.

Big thanks to @rcgldr for his radix sort base 256 C++ code, it works like a champ with worse case of 6000*6000 elements, maximum running time is 1.187s.

• Interestingly, the std::sort of C++ only failed on last 3 biggest test cases, it works fine with input of size 6000*3000.

Answer 1

The clue to your answer lies in your observation...

If we sort two array A and B first then multiply them together, we got 2 8 14 18 7 28 49 63 8 32 56 72 11 44 77 99 which is an array with m sorted subarrays.

So there are n sequences of data which are sorted, and the problem is using these to generate the answer.

Hint 1: Can you solve this using a priority queue. The number of elements in the queue would be the same as the number of sorted lists which are generated.

With

#include <vector>
#include <algorithm>
#include <random>
#include <queue>

Given the following structures (C++)

// helper to catch every tenth element.
struct Counter {
    int mCount;
    double mSum;
    Counter() : mCount(0), mSum(0) {}
    void push_back(int val)
    {
        if (mCount++ % 10 == 0)
        {
            mSum += val;
        }
    }
    double sum() { return mSum; }
};

// Storage in the priority queue for each of the sorted results.
struct Generator {
    int i_lhs;
    int i_rhs;
    int product;
    Generator() : i_lhs(0), i_rhs(0), product(0) {}
    Generator(size_t lhs, size_t rhs, int p) : i_lhs(lhs), i_rhs(rhs), product(p)
    {
    }
 };

// comparitor to get lowest value product from a priority_queue
struct MinHeap
{
    bool operator()(const Generator & lhs, const Generator & rhs)
    {
        if (lhs.product > rhs.product) return true;
        return false;
    }
};

I measured ....

double Faster(std::vector<int> lhs, std::vector<int>  rhs)
{
    Counter result;
    if (lhs.size() == 0 || rhs.size() == 0) return 0;

    std::sort(lhs.begin(), lhs.end());
    std::sort(rhs.begin(), rhs.end());
    if (lhs.size() < rhs.size()) {
        std::swap(lhs, rhs);
    }
    size_t l = 0;
    size_t r = 0;
    size_t lhs_size = lhs.size();
    size_t rhs_size = rhs.size();
    std::priority_queue<Generator, std::vector< Generator >, MinHeap > queue;
    for (size_t i = 0; i < lhs_size; i++) {
        queue.push(Generator(i, 0, lhs[i] * rhs[0]));
    }
    Generator curr;
    while (queue.size()) {
        curr = queue.top();
        queue.pop();
        result.push_back(curr.product);
        curr.i_rhs++;
        if( curr.i_rhs < rhs_size ){
            queue.push(Generator(curr.i_lhs, curr.i_rhs, lhs[curr.i_lhs] * rhs[curr.i_rhs]));
        }
    }
    return result.sum();
 }

to be faster than the following naive implementation

double Naive(std::vector<int> lhs, std::vector<int>  rhs)
{
    std::vector<int> result;
    result.reserve(lhs.size() * rhs.size());
    for (size_t i = 0; i < lhs.size(); i++) {
        for (size_t j = 0; j < rhs.size(); j++) {
            result.push_back(lhs[i] * rhs[j]);
        }
    }
    std::sort(result.begin(), result.end());
    Counter aCount;
    for (size_t i = 0; i < result.size(); i++) {
        aCount.push_back(result[i]);
    }
    return aCount.sum();
}

Sorting the input vectors is much faster than the output vector. For each row we create a generator, which will iterate over all the columns. The current product is added as the priority value to the queue, and once we have all generators made, we read them out of the queue.

Then if there is another column left for each generator, we add it back to the queue. This is from the observation that there were m subarrays of size n in the output of a pre-sorted input. The queue holds all m current minimum value for each sub-array, and the smallest of that set, is the smallest remaining of the whole list. When a Generator is removed and re-added, it ensures the top value is the next smallest item of the results.

The loop is still O(nm), as each generator is created once, reading the smallest value is O(1), and inserting into the Queue is O( log n). Which we do once for each row, so O( nm * log n + nm) which simplifies to O( nm log n ).

The Naive solution is O(nm log nm).

The performance bottle neck I found from the above solution, was the cost of inserting into the queue, and I had a performance speed up for that, but I don't think it is algorithmically" much faster.

Answer 2

merge all of these sorted subarray in O(mn)

The products are < 2^31, so 32 bit integers are enough, and a radix sort base 256, will work. The sum of every 10th item could need 64 bits.

Update - you didn't mention a memory limit of 256MB in your comments, I just noticed this. The input array size is 6000*6000*4 = 137.33MB. Allocate a working array half the size of the original array (rounded up: work_size = (1+original_size)/2 ), worst case, 3000*6000 elements (< 210MB total space needed). Treat the original (product) array as two halves and use radix sort to sort the two halves of the original array. Move the lower sorted half into the working array, then merge the working array with the upper half of the original array back into the original array. On my system (Intel 3770K 3.5 ghz, Win 7 Pro 64 bit), the 2 radix sorts will take less than 0.4 seconds (~0.185 seconds each) and the one time merge of 3000*6000 integers will take about 0.16 seconds, less than 0.6 seconds for the sort part. With this approach, there's no need to sort A or B before doing the multiply.

Are you allowed to use SIMD / xmm registers to do the outer product multiply of A and B (A o.x B) ?

Example C++ code for base 256 radix sort:

//  a is input array, b is working array
uint32_t * RadixSort(uint32_t * a, uint32_t *b, size_t count)
{
size_t mIndex[4][256] = {0};            // count / index matrix
size_t i,j,m,n;
uint32_t u;
    for(i = 0; i < count; i++){         // generate histograms
        u = a[i];
        for(j = 0; j < 4; j++){
            mIndex[j][(size_t)(u & 0xff)]++;
            u >>= 8;
        }       
    }
    for(j = 0; j < 4; j++){             // convert to indices
        m = 0;
        for(i = 0; i < 256; i++){
            n = mIndex[j][i];
            mIndex[j][i] = m;
            m += n;
        }       
    }
    for(j = 0; j < 4; j++){             // radix sort
        for(i = 0; i < count; i++){     //  sort by current lsb
            u = a[i];
            m = (size_t)(u>>(j<<3))&0xff;
            b[mIndex[j][m]++] = u;
        }
        std::swap(a, b);                //  swap ptrs
    }
    return(a);
}

---

Merge sort could be used, but it's slower. Assuming m >= n, then a conventional 2 way merge sort will take O(mn ⌈log2(n)⌉) to sort n sorted runs, each of of size m. On my system, sorting 6000 runs of 6000 integers takes about 1.7 seconds, and I don't know how long the matrix multiply would take.

Using a heap or other form of priority queue will just add to the overhead. A conventional 2-way merge sort would be faster that k-way merge sort with heap.

On a system with 16 registers, 8 of which are used as working and ending indexes or pointers to runs, a 4 way merge sort (without heap) will probably be a bit faster (around 15%), it's the same total number of operations, 1.5 x number of compares, but 0.5 x number of moves, which is a bit more cache friendly.