# Rewrite 12 consecutive ints based on bizarre logic

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!");
}

• It's really interesting ... to you... Because you know about what it does. We don't. Could you explain? Because that might make us more interested :) May 5 '16 at 21:33
• No, this is not interesting. This is like the evil, upgraded twin brother of an obfuscated duff's device. Please tell us what it does and whether it does so like expected. Comments through the code wouldn't hurt either in this case.
– Mast
May 5 '16 at 21:38
• The problem is that it would take a whole book to explain what it does. There is no way to fit it into a single post. The whole point of the function is that it does many things with as few instructions as possible - it is turing complete, for one, and is capable of executing functional programs optimally. But there are years of study behind its logic. I don't know how to start. My hope was that by just analyzing the code and with some C knowledge one could find some ways to make it more efficient. But I guess this isn't the kind of question for SO. May 5 '16 at 21:51
• Please post a version with magic() indexing mem instead of using ak to cz. Compare the machine code an "optimizing" compiler produces for both. May 6 '16 at 5:22

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 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);


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
return true;
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(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
{
}

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

private:
T total;
};

bool InstancingManager::Instance::operator <(Instance const& r) const
{
BranchlessCompare cmp;
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.

• I can't express my gratitude for providing such an insightful answer even when the question itself was regarded as bad by the community. That's exactly what I was expecting and now you gave me a lot of paths to explore that I wouldn't otherwise. Thank you. May 6 '16 at 9:49

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.