15

I have following code where a sum is calculated, based on a very large series.

The series char *a is a char array, which contains digits only (0..9).

I wanted to ask if there is any possibility to make the code faster. It is currently a bottle neck in a distributed computing application.

A small reproduction code. Not the actual code, and more simplified.

int top = 999999999;

char *a;
a = (char*) calloc(top+1, sizeof(char));

// ... fill a with initial values ...

for (int i=0; i<10; ++i) {
    unsigned long long int sum = 0;

    for (m = 1, k = top; m < k; ++m, --k) {
        // Here is the bottle neck!!
        sum += a[m]*a[k];
    }

    printf("%d\n", sum);

    // ... Add something at the end of a, and increase top ...
}

I have already tried following:

  1. Optimizing the code with -O3 (gcc compiler). The compiler line is now:

    gcc -c -Wall -fopenmp -Wno-unused-function -O3 -std=c99 -g0 -march=native -pipe -D_FILE_OFFSET_BITS=64 -m64 -fwhole-program -fprefetch-loop-arrays -funsafe-loop-optimizations -Wunsafe-loop-optimizations -fselective-scheduling -fselective-scheduling2 -fsel-sched-pipelining -fsel-sched-pipelining-outer-loops -fgcse-sm -fgcse-lm -fgcse-las -fmodulo-sched -fgcse-after-reload -fsee -DLIBDIVIDE_USE_SSE2 -DLIBDIVIDE_USE_SSE4_1 xxx.c -o xxx.o
    
  2. Using of GNU openMP to split the for-loop to multiple cores

    unsigned long long int halfway = (top>>1) + 1; // = top/2 + 1
    // digits is defined as top+1
    
    #pragma omp parallel // firstprivate/*shared*/(a, digits, halfway)
    for (unsigned long long int m = 1; m < halfway; ++m) {
        sum += a[m] * a[digits-m];
    }
    

    Result: Much, much faster, but requires more cores, and I still would like to make it faster.

  3. Casting a[m] to unsigned long long int before multiplication

    sum += (unsigned long long int)a[m] * a[k];
    

    Result: A small performance boost.

  4. Using a multiplication lookup table, because an array-lookup is faster than the actual multiplication.

    sum += multiply_lookup[a[m]][a[k]]; // a[m]*a[k];
    

    Result: A small performance boost.

  5. I have tried to find a mathematical solution to reduce operations, but it seems like nothing can be optimized, mathematically seen.

I have following idea for optimization:

I have read that the multiplication of floats (asm fmul) is much faster than the multiplication of integers (asm mul). Just changing int to float doesn't help -- but I think the code might become much more performant if the work is done using MMX or SSE instruction sets, or if the work is done by the FPU. Although I have some assembler knowledge, I have no knowledge about these topics.

However, if you have additional ideas how to optimize it, I am glad to hear them.

Update Some additional information:

  • The series grows by 1 element after each loop.
  • While the series grows, top gets increased.
  • When top is reaching the array limit, a will get increased by 100000 bytes using realloc().
  • Platform: Debian Linux Jessie x64, on an Intel(R) Xeon(R) CPU X3440 @ 2.53GHz

Additional off-topic question: Do you know the mathematical name of this sum, where the pairs of elements of the series are multiplied from outside to inside?

18
  • 2
    I'm voting to close this question as off-topic because SO is no code review site. Commented Dec 13, 2015 at 18:00
  • 4
    It is not a code review. I want to know how a fast 8-bit multiplication in assembler is possible. Commented Dec 13, 2015 at 18:01
  • You have running code and ask for optimisation. This requires very well a code-review. And the multiplications will not be performed in 8 bits anyway. This is C. Commented Dec 13, 2015 at 18:01
  • 8
    @Olaf: Optimization questions are generally on topic for SO. There isn't enough code here for a proper post to CodeReview.SE. This is a perfectly acceptable question on SO.
    – Cornstalks
    Commented Dec 13, 2015 at 18:04
  • 5
    @Olaf: Look at the tour. This question definitely fits in the "on-topic" box on that page (and definitely doesn't fit in the "off-topic" box).
    – Cornstalks
    Commented Dec 13, 2015 at 18:08

3 Answers 3

7

You can use the little-known PMADDUBSW (Multiply and Add Packed Signed and Unsigned Bytes) for this. The signed/unsigned business doesn't matter here, everything is in the interval [0 .. 9] anyway. The add is saturating, but that doesn't matter here because 9*9 is only 81. With intrinsics that's _mm_maddubs_epi16. Because the k index goes down, you have to byte-reverse it, which you can do with PSHUFB (_mm_shuffle_epi8). An annoying thing happens when the indexes "meet" in the middle, you can do that part one by one..

Here's a try, only slightly tested:

__m128i sum = _mm_setzero_si128();
int m, k;
for (m = 1, k = top - 15; m + 15 < k; m += 16, k -= 16) {
   __m128i am = _mm_loadu_si128((__m128i*)(a + m));
   __m128i ak = _mm_loadu_si128((__m128i*)(a + k));
   ak = _mm_shuffle_epi8(ak, _mm_set_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 ,15));
   sum = _mm_add_epi16(sum, _mm_maddubs_epi16(am, ak));
}
// could use phaddw, but I do this the long way to avoid overflow slightly longer
sum = _mm_add_epi32(_mm_unpacklo_epi16(sum, _mm_setzero_si128()),
                    _mm_unpackhi_epi16(sum, _mm_setzero_si128()));
sum = _mm_hadd_epi32(sum, sum);
sum = _mm_hadd_epi32(sum, sum);
int s = _mm_cvtsi128_si32(sum);
// this is for the "tail"
k += 15;
for (; m < k; ++m, --k)
    s += a[m] * a[k];

Also I ignore overflow. You can do this for (216-1)/(2*81) = 404 iterations and still definitely have no overflow. If you need more, periodically add this to a 32bit result.

In a quick benchmark, this is about 7 times as fast as the simple way (tested with 2KB of random data on a 4770K, taking the best out of a hundred runs for each).

Using pointers as suggested by an other answer improves it further, to about 9 times as fast as the simple way. With indices there was some weird sign-extension going on.

int foobar(char* a, int top)
{
    __m128i sum = _mm_setzero_si128();

    char *m, *k;
    for (m = a + 1, k = a + top - 15; m + 15 < k; m += 16, k -= 16) {
       __m128i am = _mm_loadu_si128((__m128i*)(m));
       __m128i ak = _mm_loadu_si128((__m128i*)(k));
       ak = _mm_shuffle_epi8(ak, _mm_set_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15));
       sum = _mm_add_epi16(sum, _mm_maddubs_epi16(am, ak));
    }

    sum = _mm_add_epi32(_mm_unpacklo_epi16(sum, _mm_setzero_si128()),
                        _mm_unpackhi_epi16(sum, _mm_setzero_si128()));
    sum = _mm_hadd_epi32(sum, sum);
    sum = _mm_hadd_epi32(sum, sum);
    int s = _mm_cvtsi128_si32(sum);

    k += 15;
    for (; m < k; ++m, --k)
        s += *m * *k;

    return s;
}

Split up in parts, still about 9 times as fast as the original despite the extra logic:

int foobar(char* a, int top)
{
    int s = 0;
    char *m, *k;
    for (m = a + 1, k = a + top - 15; m + 15 < k;) {
        __m128i sum = _mm_setzero_si128();
        for (int i = 0; i < 404 && m + 15 < k; m += 16, k -= 16, ++i) {
           __m128i am = _mm_loadu_si128((__m128i*)(m));
           __m128i ak = _mm_loadu_si128((__m128i*)(k));
           ak = _mm_shuffle_epi8(ak, _mm_set_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 ,15));
           sum = _mm_add_epi16(sum, _mm_maddubs_epi16(am, ak));
        }
        sum = _mm_add_epi32(_mm_unpacklo_epi16(sum, _mm_setzero_si128()),
                            _mm_unpackhi_epi16(sum, _mm_setzero_si128()));
        sum = _mm_hadd_epi32(sum, sum);
        sum = _mm_hadd_epi32(sum, sum);
        s += _mm_cvtsi128_si32(sum);
    }

    k += 15;
    for (; m < k; ++m, --k)
        s += *m * *k;

    return s;
}
18
  • Wow, you seem to have a lot of knowledge with that! But it doesn't work at the moment. I think something relating to SSE is not included or not available in gcc? error: unknown type name '__m128i' . Commented Dec 13, 2015 at 18:35
  • The compiler line is gcc -c -Wall -fopenmp -Wno-unused-function -O3 -std=c99 -g0 -march=native -pipe -D_FILE_OFFSET_BITS=64 -m64 -fwhole-program -fprefetch-loop-arrays -funsafe-loop-optimizations -Wunsafe-loop-optimizations -fselective-scheduling -fselective-scheduling2 -fsel-sched-pipelining -fsel-sched-pipelining-outer-loops -fgcse-sm -fgcse-lm -fgcse-las -fmodulo-sched -fgcse-after-reload -fsee -DLIBDIVIDE_USE_SSE2 -DLIBDIVIDE_USE_SSE4_1 -DWITH_SSE2=ON -S -masm=intel -fverbose-asm xx.c -o xx.asm . Commented Dec 13, 2015 at 18:35
  • You need something like #include <x86intrin.h> for this
    – user555045
    Commented Dec 13, 2015 at 18:35
  • @DanielMarschall More than half of the flags you pass are implied by -O3.
    – fuz
    Commented Dec 13, 2015 at 18:36
  • @harold Thanks. It compiles, but I think something is still wrong. Here are some compiler warnings and my code pastebin.com/d5w2A4Fh . Commented Dec 13, 2015 at 18:50
3

This loop,

for (m = 1, k = top; m < k; ++m, --k) {
    // Here is the bottle neck!!
    sum += a[m]*a[k];
}

could gain benefit from changing to:

char *b = a + top;
a++;
for (; a < b; ) 
{
    sum += ( *a++ ) * ( *b--);
}

By removing the [] you are saving alot of arithmetic for every access into the array. This halves the theoretical number of address computations from: 4 with ++m --k and a[m] a[k] to 2 with *a++ *b--

Simple pointer increment is cheaper and usually faster overall, as array access is not always optimized perfectly using the [].

Hope this helps

10
  • 1
    Address computations are generally free on modern architectures. x86 processors actually have dedicated ALUs just for address computation. If anything, this is going to confuse the compiler even more.
    – fuz
    Commented Dec 13, 2015 at 18:21
  • 1
    @FUZxxl: I thought that too, but I just tested his code and this answer gave a ~5.2x speedup.
    – Cornstalks
    Commented Dec 13, 2015 at 18:22
  • You may be right, but I doubt it is perfectly 'free'. Perhaps the SO can try it and see.
    – Grantly
    Commented Dec 13, 2015 at 18:22
  • @Cornstalks Can you show me the code you benchmarked with? This is interesting.
    – fuz
    Commented Dec 13, 2015 at 18:22
  • @FUZxxl: this is the code I profiled with. Dead simple (too simple, really). I think the subscript access really threw off the optimizer. With subscript accessing, things take longer if I change i<10 to i<20. However, without subscript accessing (i.e. this answer), it's always constant time. I suspect the optimizer has an easier time optimizing the outer for-loop (and just changing it to a constant multiply) away when incremental pointers are used.
    – Cornstalks
    Commented Dec 13, 2015 at 18:31
2

The operation you want to perform is called a discrete convolution, it appears when multiplying large numbers. The naïve algorithm you use has a complexity of O(n2), but an O(n log n) solution obtains using a discrete Fourier transform.

Discrete Convolutions

The discrete convolution c = ab of two sequences a = a0, a1, …, an − 1 and b = b0, b1, …, bn − 1 with n elements each is a sequence of 2n − 1 elements defined for each k as the sum:

ci = ∑max(0, in + 1) ≤ j < min(n, i + 1) ajbij

If we assume that ai = bi = 0 for i ∉ {0, …, n − 1}, then we can simplify this and sum over all integer i:

ci = ∑j ajbij

Notice how this is the operation you want to perform: a = b = a[] and the i-th iteration of your loop, sum is just out ci.

The discrete convolution is well-studied and appears in a number of mathematical problems related to signal-processing, combinatorics and statistics. Fortunately, it can be computed in superlinear time O(n log n) instead of the naïve O(n2).

Cyclic Discrete Convolutions

We can expand b to bN such that

bk = bk mod n, or
bN = b0, b1, …, bn − 1, b0, …, bn − 1, b0, …

This is called the cyclic expansion of b to bN. The discrete convolution abN is called the discrete cyclic convolution of a and b.

Note that the non-cyclic discrete convolution of a and b can be computed from a cyclic convolution by appending sufficiently many zeroes to a and b before convolving such that the cyclic expansion of b doesn't alter the result. See this article for more details.

Discrete Fourier Transform

The discrete Fourier transform (DFT) transforms a sequence of samples a into a complex frequency spectrum ℱ(a) of the sampled signal. The Fourier transform is invertible and computable in superlinear time O(n log n) using various fast Fourier transform (FFT) algorithms. This transform has many applications in signal-processing. It has the useful property that

abN = ℱ−1(ℱ(a) × ℱ(b))

That is, the cyclic convolution of two sequences a and b is element-wise multiplication under Fourier transform. Since element-wise multiplication can be done in linear time, this makes it possible to compute the result you need in total time O(n log n).

Implementation

The FFTW project provides highly optimized implementations of the fast Fourier transform.

Problems

Fourier transforms work on complex numbers and requires a floating-point unit to be fast. The result is likely slightly imprecise. As far as I know, precise methods are possible with some modulo magic (this is what's used in the Schönhage-Strassen algorithm for large-integer multiplication) but I'm not sure about the details.

15
  • So it looks like i am searching Sum = (a*a) [top]? But the definition shows that m goes from -inf to +inf ; so the case that the sum stops once both sides meet each other, isn't there? I do not quite understand why you are referring that algorithm in re multiplying large numbers ; isn't the actual problem doing a lot of multiplications of digits into a huge sum? How are they related to a bigint multiplication? Commented Dec 18, 2015 at 19:49
  • @DanielMarschall The sum goes from −∞ to ∞ but your series is 0 on all but the entries 0 to n − 1, so all but at most n terms of the convolution are zero. Ever done long multiplication (schriftliches multiplizieren)? You basically do the kind of operation and then sum up all the sums shifted by a power of 10 to get a result. This page might be more useful to you.
    – fuz
    Commented Dec 18, 2015 at 19:54
  • If I understood correctly, this is the solution pastebin.com/C1WLbRXj . I did long multiplication long ago, but I don't understand how this solves the problem. Don't the carry values destroy the calculation? If I did multiply 123*654 using Base 81, then the resulting number {6, 10, 12} could be used to form the sum = 6+10+12 = 28. Commented Dec 18, 2015 at 20:20
  • Update: A better solution than the one I wrote in pastebin: a={1,2,3}; b={4,5,6}; (a*b)[2] = 28 It seems that this is a 1D convolution. (Actually, I worked with (2D) convolution before when I implemented a quick'n'dirty Photoshop filter. But it was slow, and used only simple multiplication). Commented Dec 18, 2015 at 20:56
  • 1
    The precise variant you allude to is the number-theoretic transform, it's really just an FFT over a finite field
    – user555045
    Commented Dec 20, 2015 at 0: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.