Rewrite 12 consecutive ints based on bizarre logic

Question

Despite looking terribly ugly, it is pretty interesting due to some specific connections to cellular automata theory. I'm looking for ways of simplifying it further. This program takes 2.7 seconds to run on my machine. Is it possible, by merely changing the logic / control flow on magic, to speed it up? How fast can it be?

#include <stdio.h>

// Rewrites 12 consecutive ints based on bizarre logic.
void magic(int* mem){
    static int next = 1000000;
    int ak = mem[0+0], ax = mem[0+1], ay = mem[0+2], az = mem[0+3],
        bk = mem[4+0], bx = mem[4+1], by = mem[4+2], bz = mem[4+3],
        ck = mem[8+0], cx = mem[8+1], cy = mem[8+2], cz = mem[8+3], A, B, C, D;
    if (ak == 0 && bk == 0 && ck == 0) return;
    if (ak == 0 && bk == 0 && ck != 0) ak = ck, ax = cx, ay = cy, az = cz, ck =  0, cx =  0, cy =  0, cz =  0;
    if (ak != 0 && bk != 0 && ck == 0) ck = bk, cx = bx, cy = by, cz = bz, bk =  0, bx =  0, by =  0, bz =  0;
    if (ak == 1 && bk == 0 && ck != 0 && ax == -cx) cx = ay, ak = 0, ax = 0, ay = 0, az = 0;
    if (ak == 1 && bk == 0 && ck != 0 && ax == -cy) cy = ay, ak = 0, ax = 0, ay = 0, az = 0;
    if (ak == 1 && bk == 0 && ck != 0 && ax == -cz) cz = ay, ak = 0, ax = 0, ay = 0, az = 0;
    if (ck == 1 && bk == 0 && ck != 0 && cx == -ax) ax = cy, ck = 0, cx = 0, cy = 0, cz = 0;
    if (ck == 1 && bk == 0 && ck != 0 && cx == -ay) ay = cy, ck = 0, cx = 0, cy = 0, cz = 0;
    if (ck == 1 && bk == 0 && ck != 0 && cx == -az) az = cy, ck = 0, cx = 0, cy = 0, cz = 0;
    if (ay == -az && ay < 0) ay *= -1, az *= -1;
    if (by == -bz && by < 0) by *= -1, bz *= -1; 
    if (cy == -cz && cy < 0) cy *= -1, cz *= -1; 
    if (ak == 1 && ax == -ay) ak = 0, ax = 0, ay = 0, az = 0; 
    if (ck == 1 && cx == -cy) ck = 0, cx = 0, cy = 0, cz = 0; 
    if (ak > 0 && bk == 0 && ck > 0 && ((ax > 0 && cx < 0) || (ak == 1 && ck != 1 && ax > 0) || (ck == 1 && ak != 1 && cx < 0))){
        ak = ak+ck, ck = ak-ck, ak = ak-ck;
        ax = ax+cx, cx = ax-cx, ax = ax-cx;
        ay = ay+cy, cy = ay-cy, ay = ay-cy;
        az = az+cz, cz = az-cz, az = az-cz;
        if (ax == -cx) ax *= -1, cx *= -1;
        if (ax == -cy) ax *= -1, cy *= -1;
        if (ax == -cz) ax *= -1, cz *= -1;
        if (ay == -cx) ay *= -1, cx *= -1;
        if (ay == -cy) ay *= -1, cy *= -1;
        if (ay == -cz) ay *= -1, cz *= -1;
        if (az == -cx) az *= -1, cx *= -1;
        if (az == -cy) az *= -1, cy *= -1;
        if (az == -cz) az *= -1, cz *= -1;
    };
    if (ak < -1 && bk == 0)
        A = ak, B = ax, C = ay, D = az,
        ak = -A, ax = C, ay = B>0?B+0:-(-B+0), az = B>0?B+1:-(-B+2),
        bk = -A, bx = D, by = B>0?B+2:-(-B+1), bz = B>0?B+3:-(-B+3);
    if (ak > 1 && bk == 0 && ck > 1 && ax == -cx && ak == ck)
        ak = 1, ax = ay, ay = cy, ck = 1, cx = az, cy = cz, az = 0, cz = 0;
    if (ak > 1 && bk == 0 && ck > 1 && ax == -cx && ak != ck)
        A = (next+=6), ak = ak+ck, ck = ak-ck, ak = ak-ck, ak = -ak, ck = -ck, ax = A, cx = -A;
    mem[0+0] = ak, mem[0+1] = ax, mem[0+2] = ay, mem[0+3] = az,
    mem[4+0] = bk, mem[4+1] = bx, mem[4+2] = by, mem[4+3] = bz,
    mem[8+0] = ck, mem[8+1] = cx, mem[8+2] = cy, mem[8+3] = cz;
};

void compute(int *mem, int len){
    for (int j=0; j<3; ++j)
        for (int i=j; i<len-2; i+=3)
            magic(mem+i*4);
};

const int program[] = {3,2,3,-1,0,0,0,0,3,4,5,-2,0,0,0,0,3,-4,6,7,0,0,0,0,62,-6,8,9,0,0,0,0,3,-8,10,11,0,0,0,0,3,-7,-10,12,0,0,0,0,3,13,14,-12,0,0,0,0,4,15,16,-13,0,0,0,0,5,17,18,-15,0,0,0,0,6,19,20,-17,0,0,0,0,7,21,22,-19,0,0,0,0,8,23,24,-21,0,0,0,0,9,25,26,-23,0,0,0,0,10,27,28,-25,0,0,0,0,11,29,30,-27,0,0,0,0,12,31,32,-29,0,0,0,0,13,33,34,-31,0,0,0,0,14,35,36,-33,0,0,0,0,15,37,38,-35,0,0,0,0,16,39,40,-37,0,0,0,0,17,41,42,-39,0,0,0,0,18,43,44,-41,0,0,0,0,19,45,46,-43,0,0,0,0,20,47,48,-45,0,0,0,0,21,49,50,-47,0,0,0,0,22,51,52,-49,0,0,0,0,23,53,54,-51,0,0,0,0,24,55,56,-53,0,0,0,0,25,57,58,-55,0,0,0,0,26,59,60,-57,0,0,0,0,27,61,62,-59,0,0,0,0,28,63,64,-61,0,0,0,0,29,65,66,-63,0,0,0,0,30,67,68,-65,0,0,0,0,31,69,70,-67,0,0,0,0,32,71,72,-69,0,0,0,0,33,73,74,-71,0,0,0,0,34,75,76,-73,0,0,0,0,35,77,78,-75,0,0,0,0,36,79,80,-77,0,0,0,0,37,81,82,-79,0,0,0,0,38,83,84,-81,0,0,0,0,39,85,86,-83,0,0,0,0,40,87,88,-85,0,0,0,0,41,89,90,-87,0,0,0,0,42,91,92,-89,0,0,0,0,43,93,94,-91,0,0,0,0,44,95,96,-93,0,0,0,0,45,97,98,-95,0,0,0,0,46,99,100,-97,0,0,0,0,47,101,102,-99,0,0,0,0,48,103,104,-101,0,0,0,0,49,105,106,-103,0,0,0,0,50,107,108,-105,0,0,0,0,51,109,110,-107,0,0,0,0,52,111,112,-109,0,0,0,0,53,113,114,-111,0,0,0,0,54,115,116,-113,0,0,0,0,55,117,118,-115,0,0,0,0,56,119,120,-117,0,0,0,0,57,121,122,-119,0,0,0,0,58,123,124,-121,0,0,0,0,59,125,126,-123,0,0,0,0,60,127,128,-125,0,0,0,0,61,-9,129,-127,0,0,0,0,3,-129,-11,130,0,0,0,0,3,-128,-130,131,0,0,0,0,3,-126,-131,132,0,0,0,0,3,-124,-132,133,0,0,0,0,3,-122,-133,134,0,0,0,0,3,-120,-134,135,0,0,0,0,3,-118,-135,136,0,0,0,0,3,-116,-136,137,0,0,0,0,3,-114,-137,138,0,0,0,0,3,-112,-138,139,0,0,0,0,3,-110,-139,140,0,0,0,0,3,-108,-140,141,0,0,0,0,3,-106,-141,142,0,0,0,0,3,-104,-142,143,0,0,0,0,3,-102,-143,144,0,0,0,0,3,-100,-144,145,0,0,0,0,3,-98,-145,146,0,0,0,0,3,-96,-146,147,0,0,0,0,3,-94,-147,148,0,0,0,0,3,-92,-148,149,0,0,0,0,3,-90,-149,150,0,0,0,0,3,-88,-150,151,0,0,0,0,3,-86,-151,152,0,0,0,0,3,-84,-152,153,0,0,0,0,3,-82,-153,154,0,0,0,0,3,-80,-154,155,0,0,0,0,3,-78,-155,156,0,0,0,0,3,-76,-156,157,0,0,0,0,3,-74,-157,158,0,0,0,0,3,-72,-158,159,0,0,0,0,3,-70,-159,160,0,0,0,0,3,-68,-160,161,0,0,0,0,3,-66,-161,162,0,0,0,0,3,-64,-162,163,0,0,0,0,3,-62,-163,164,0,0,0,0,3,-60,-164,165,0,0,0,0,3,-58,-165,166,0,0,0,0,3,-56,-166,167,0,0,0,0,3,-54,-167,168,0,0,0,0,3,-52,-168,169,0,0,0,0,3,-50,-169,170,0,0,0,0,3,-48,-170,171,0,0,0,0,3,-46,-171,172,0,0,0,0,3,-44,-172,173,0,0,0,0,3,-42,-173,174,0,0,0,0,3,-40,-174,175,0,0,0,0,3,-38,-175,176,0,0,0,0,3,-36,-176,177,0,0,0,0,3,-34,-177,178,0,0,0,0,3,-32,-178,179,0,0,0,0,3,-30,-179,180,0,0,0,0,3,-28,-180,181,0,0,0,0,3,-26,-181,182,0,0,0,0,3,-24,-182,183,0,0,0,0,3,-22,-183,184,0,0,0,0,3,-20,-184,185,0,0,0,0,3,-18,-185,186,0,0,0,0,3,-16,-186,-14,0,0,0,0,3,-5,187,188,0,0,0,0,121,-187,189,190,0,0,0,0,3,-189,191,192,0,0,0,0,3,-188,-191,193,0,0,0,0,3,194,195,-193,0,0,0,0,63,196,197,-194,0,0,0,0,64,198,199,-196,0,0,0,0,65,200,201,-198,0,0,0,0,66,202,203,-200,0,0,0,0,67,204,205,-202,0,0,0,0,68,206,207,-204,0,0,0,0,69,208,209,-206,0,0,0,0,70,210,211,-208,0,0,0,0,71,212,213,-210,0,0,0,0,72,214,215,-212,0,0,0,0,73,216,217,-214,0,0,0,0,74,218,219,-216,0,0,0,0,75,220,221,-218,0,0,0,0,76,222,223,-220,0,0,0,0,77,224,225,-222,0,0,0,0,78,226,227,-224,0,0,0,0,79,228,229,-226,0,0,0,0,80,230,231,-228,0,0,0,0,81,232,233,-230,0,0,0,0,82,234,235,-232,0,0,0,0,83,236,237,-234,0,0,0,0,84,238,239,-236,0,0,0,0,85,240,241,-238,0,0,0,0,86,242,243,-240,0,0,0,0,87,244,245,-242,0,0,0,0,88,246,247,-244,0,0,0,0,89,248,249,-246,0,0,0,0,90,250,251,-248,0,0,0,0,91,252,253,-250,0,0,0,0,92,254,255,-252,0,0,0,0,93,256,257,-254,0,0,0,0,94,258,259,-256,0,0,0,0,95,260,261,-258,0,0,0,0,96,262,263,-260,0,0,0,0,97,264,265,-262,0,0,0,0,98,266,267,-264,0,0,0,0,99,268,269,-266,0,0,0,0,100,270,271,-268,0,0,0,0,101,272,273,-270,0,0,0,0,102,274,275,-272,0,0,0,0,103,276,277,-274,0,0,0,0,104,278,279,-276,0,0,0,0,105,280,281,-278,0,0,0,0,106,282,283,-280,0,0,0,0,107,284,285,-282,0,0,0,0,108,286,287,-284,0,0,0,0,109,288,289,-286,0,0,0,0,110,290,291,-288,0,0,0,0,111,292,293,-290,0,0,0,0,112,294,295,-292,0,0,0,0,113,296,297,-294,0,0,0,0,114,298,299,-296,0,0,0,0,115,300,301,-298,0,0,0,0,116,302,303,-300,0,0,0,0,117,304,305,-302,0,0,0,0,118,306,307,-304,0,0,0,0,119,308,309,-306,0,0,0,0,120,-190,310,-308,0,0,0,0,3,-310,-192,311,0,0,0,0,3,-309,-311,312,0,0,0,0,3,-307,-312,313,0,0,0,0,3,-305,-313,314,0,0,0,0,3,-303,-314,315,0,0,0,0,3,-301,-315,316,0,0,0,0,3,-299,-316,317,0,0,0,0,3,-297,-317,318,0,0,0,0,3,-295,-318,319,0,0,0,0,3,-293,-319,320,0,0,0,0,3,-291,-320,321,0,0,0,0,3,-289,-321,322,0,0,0,0,3,-287,-322,323,0,0,0,0,3,-285,-323,324,0,0,0,0,3,-283,-324,325,0,0,0,0,3,-281,-325,326,0,0,0,0,3,-279,-326,327,0,0,0,0,3,-277,-327,328,0,0,0,0,3,-275,-328,329,0,0,0,0,3,-273,-329,330,0,0,0,0,3,-271,-330,331,0,0,0,0,3,-269,-331,332,0,0,0,0,3,-267,-332,333,0,0,0,0,3,-265,-333,334,0,0,0,0,3,-263,-334,335,0,0,0,0,3,-261,-335,336,0,0,0,0,3,-259,-336,337,0,0,0,0,3,-257,-337,338,0,0,0,0,3,-255,-338,339,0,0,0,0,3,-253,-339,340,0,0,0,0,3,-251,-340,341,0,0,0,0,3,-249,-341,342,0,0,0,0,3,-247,-342,343,0,0,0,0,3,-245,-343,344,0,0,0,0,3,-243,-344,345,0,0,0,0,3,-241,-345,346,0,0,0,0,3,-239,-346,347,0,0,0,0,3,-237,-347,348,0,0,0,0,3,-235,-348,349,0,0,0,0,3,-233,-349,350,0,0,0,0,3,-231,-350,351,0,0,0,0,3,-229,-351,352,0,0,0,0,3,-227,-352,353,0,0,0,0,3,-225,-353,354,0,0,0,0,3,-223,-354,355,0,0,0,0,3,-221,-355,356,0,0,0,0,3,-219,-356,357,0,0,0,0,3,-217,-357,358,0,0,0,0,3,-215,-358,359,0,0,0,0,3,-213,-359,360,0,0,0,0,3,-211,-360,361,0,0,0,0,3,-209,-361,362,0,0,0,0,3,-207,-362,363,0,0,0,0,3,-205,-363,364,0,0,0,0,3,-203,-364,365,0,0,0,0,3,-201,-365,366,0,0,0,0,3,-199,-366,367,0,0,0,0,3,-197,-367,-195,0,0,0,0,3,-3,368,-368};
const int node_count = sizeof program / sizeof(int) / 4;
const int space = 4000;
int memory[space*4];
int main(){
    for (int i=0; i < space*4; ++i)
        memory[i] = i < node_count*4 ? program[i] : 0;

    // If this halts, `magic` is probably implemented correctly.
    while (memory[1]!=-1 || memory[2]!=-memory[3])
        compute(memory, space);

    printf("Pass!");
}

Answer 1

You could probably reduce the number of branches in that code and get additional performance using SIMD, like SSE or AVX. All of the variables in magic() would fit in a few registers.

It is likely for the compiler to be using a lot of branches in that code. I can't tell from here but I would guess you have a lot of unpredictable branches. Unpredictable branches kill performance, because the processor is wasting its time trying to speculatively execute the wrong instructions 50% of the time, then flushing the pipeline and restarting.

You can immediately eliminate some ifs, by multiplying by 1 when you don't need to multiply by -1. A ternary cond?trueval:falseval is likely to be a conditional move, if the instruction set has it.

To eliminate branches, you can use vector comparison intrinsics. They do not work like typical comparisons, instead, they generate a result with all binary 1's for true, and all binary 0's for false. This way, you can write the code to compute both sides of the if, then bitwise-and (_mm_and_si128) the true path against the comparison, bitwise-and the else path against the ones-complement of the comparison (_mm_andnot_si128), then bitwise-or (_mm_or_si128) together the two execution paths. This will leave you with the correct execution path's result in the result.

The reason: it is often far faster to compute both paths and select the results at instruction retirement than it is to constantly speculate through branches in the wrong direction and flush the pipeline and start over at the correct target. Modern CPUs aggressively compute far ahead of where the program is really executing, and mispredictions throw that work away.

static int next = 1000000;

Function-level static is error prone.

Replace

int ak = mem[0+0], ax = mem[0+1], ay = mem[0+2], az = mem[0+3],
    bk = mem[4+0], bx = mem[4+1], by = mem[4+2], bz = mem[4+3],
    ck = mem[8+0], cx = mem[8+1], cy = mem[8+2], cz = mem[8+3],
    A, B, C, D;

with:

__m128i a = _mm_loadu_si128(mem);
__m128i b = _mm_loadu_si128(mem+4);
__m128i c = _mm_loadu_si128(mem+8);
__m128i A;

Some of the operations across the vectors, such as checks of ak bk ck, are the worst case scenario, because it makes a decision across one member of multiple vectors.

Determine what is zero in parallel

__m128i zero = _mm_setzero_si128();
__m128i a_zero = _mm_cmpeq_epi32(a, zero);
__m128i b_zero = _mm_cmpeq_epi32(b, zero);
__m128i c_zero = _mm_cmpeq_epi32(c, zero);

It is convenient to have testable ints for the early outs

int a_zero_msk = _mm_movemask_epi8(a_zero);
int b_zero_msk = _mm_movemask_epi8(b_zero);
int c_zero_msk = _mm_movemask_epi8(c_zero);

Replace

if (ak == 0 && bk == 0 && ck == 0) return;

with

if ((a_zero_msk & b_zero_msk & c_zero_msk) & 0xFF) return;

You want to avoid short-circuit evaluation, using bitwise operations disables any attempt to jump over any "unnecessary" comparisons.

The rest of the code consists of decisions about whether to multiply several things by either 1 or -1. You should use branchless merging of the conditions to select a 1 or -1 multiplier. Then unconditionally multiply.

You don't have to completely eliminate branching, just reduce it to a level where the CPU has a nice chunk of work to do between mispredictions. It will have sufficient time to get ahead and hide latencies.

---

For an example of converting some code to a branchless alternative, consider the operator< for a std::map or std::set, it has completely random branching, and benefits from using a branchless implementation.

Branchy code

    // Changing shaders is most expensive
    if (shaderState < r.shaderState)
        return true;
    if (shaderState > r.shaderState)
        return false;

    // Changing textures is slightly expensive
    if (textureState < r.textureState)
        return true;
    if (textureState > r.textureState)
        return false;

    // Changing geometry is least expensive
    if (geometryState < r.geometryState)
        return true;
    if (geometryState > r.geometryState)
        return false;

    // Sort by drawing mode
    if (drawCall.mode < r.drawCall.mode)
        return true;
    if (drawCall.mode > r.drawCall.mode)
        return false;

    // Sort biggest calls first
    if (drawCall.count < r.drawCall.count)
        return false;
    if (drawCall.count > r.drawCall.count)
        return true;

    // Then by original order
    if (drawCall.indices < r.drawCall.indices)
        return true;
    if (drawCall.indices > r.drawCall.indices)
        return false;

    return false;
}

You can replace that with a helper that uses branchless code to carry forward the true or false with bit logic.

Branchless code

class BranchlessCompare
{
public:
    BranchlessCompare()
    {
        reset();
    }

    void reset() noexcept
    {
        force_false = 0;
        force_true = 0;
    }

    template<typename A>
    void compare(A& lhs, A& rhs) noexcept
    {
        // Force T if we are not forcing F and we would return T
        // Force F if we are not forcing T and we would return F
        force_true |= (force_false^1) & int(lhs < rhs);
        force_false |= (force_true^1) & int(rhs < lhs);
    }

    bool get() const noexcept
    {
        return force_true & (force_false^1);
    }

    operator bool() const noexcept
    {
        return get();
    }

private:
    int force_true;
    int force_false;
};

bool InstancingManager::Instance::operator <(Instance const& r) const
{
    BranchlessCompare cmp;
    cmp.compare(shaderState, r.shaderState);
    cmp.compare(textureState, r.textureState);
    cmp.compare(geometryState, r.geometryState);
    cmp.compare(drawCall.mode, r.drawCall.mode);
    cmp.compare(drawCall.count, r.drawCall.count);
    cmp.compare(drawCall.count, r.drawCall.count);
    cmp.compare(drawCall.indices, r.drawCall.indices);
    return cmp;
}

You can go another step and make the comparer use only one value to track the relative value during the comparison. This is limited to 31 comparisons for a 32-bit integer.

class BranchlessCompare
{
public:
    using T = std::int32_t;

    BranchlessCompare()
    {
        reset();
    }

    void reset() noexcept
    {
        total = 0;
    }

    template<typename A>
    void compare(A& lhs, A& rhs) noexcept
    {
        // Make earlier comparison more significant
        // And make total more positive if lhs < rhs
        // And make total more negative if rhs < lhs
        total = (total << 1) + T((lhs < rhs) - (lhs < rhs));
    }

    bool get() const noexcept
    {
        return total > 0;
    }

    operator bool() const noexcept
    {
        return get();
    }

private:
    T total;
};

bool InstancingManager::Instance::operator <(Instance const& r) const
{
    BranchlessCompare cmp;
    cmp.compare(shaderState, r.shaderState);
    cmp.compare(textureState, r.textureState);
    cmp.compare(geometryState, r.geometryState);
    cmp.compare(drawCall.mode, r.drawCall.mode);
    cmp.compare(drawCall.count, r.drawCall.count);
    cmp.compare(drawCall.count, r.drawCall.count);
    cmp.compare(drawCall.indices, r.drawCall.indices);
    return cmp;
}

Hopefully that gives you more ideas how you can use masking to make expressions that update the state incrementally, without branching around.

Answer 2

It looks as though you only operate on a single chunk of the array at a time, and once done you don't revisit an index. If that is the case, this is a good candidate for doing these operations in parallel with multiple threads.

Also, if your consecutive IFs are mutually exclusive ( I can't tell ), then you could possibly convert them to a switch statement and/or use Boolean operators rather than logical operators. At a minimum you could use else-ifs.

Finally, you could always re-write this in assembly, or perhaps fortran, if you are super concerned with performance.