# Linear cryptanalysis of a substitution permutation network

I was translating a script I had that would automate linear cryptanalysis of SP networks from Python to C++. For those who aren't familiar with cryptography, linear cryptanalysis involves finding relationships between certain input and output bits of a $1 \rightarrow 1$ array.

What I'm doing is looping over all possible input bits and output bits, then seeing how many times there's a correlation between the bits. If this correlation is above a threshold value, I store it and move on.

My main problem is that running the substitute_rounds(rounds) function after rounds = setup_first_round() takes nearby 5 times as long to ran than the Python version did! Clearly, I wrote substitute_rounds() poorly!

Any suggestions would be very helpful, but I'm primarily trying to get a speed boost on substitute_rounds().

Two functions (Cartesian product and combinations) were just copy pasted from Stack Overflow, so I'm showing them after my code.

#include <vector>
#include <algorithm>
#include <iostream>
#include "math.h"

#define bits_in_sbox 8
#define subkey_bits 64
#define default_threshold 0.0625

using namespace std;

//return the i'th bit of s
inline bool ith_bit(int s, short i)
{
return (s >> i) & 1;
}

//I use shorts instead of ints here because of memory issues
//with each new round, the # of these I have to store goes up exponentially
struct linear_relationship
{
std::vector<short> inputs;
std::vector<short> outputs;
float bias;
};

//only used for debugging
void print_linear_relationship(linear_relationship & i)
{
for(auto &j:i.inputs)std::cout << j << " ";
std::cout << "| ";
for(auto &j:i.outputs)std::cout << j << " ";
std::cout << "| ";
std::cout << i.bias << std::endl;
}

//all linear relationships that use our input bits and a particular sbox that are above a certain threshold for the bias
std::vector<linear_relationship> bit(std::vector<short> & inputs, short * sbox)
{
std::vector<linear_relationship> ret;

//offset for shifting the output bits
int offset = inputs - (inputs % bits_in_sbox);

//make the inputs directly enterable as indexes of the sbox
//so bit 8 -> bit 0, bit 25 -> bit 1
std::vector<short> normalized_inputs;
for(short i: inputs)
normalized_inputs.push_back(i % bits_in_sbox);

std::vector<short> output_bits_index;
for(short i = 0; i < bits_in_sbox; ++i)
output_bits_index.push_back(i);

//loop through the number of output bits we'll be looking at
for(short number_of_output_bits = 1; number_of_output_bits <= bits_in_sbox; ++number_of_output_bits)
{

do
{
//the current combination of our output bit vector
std::vector<short> copy_of_output_bits( output_bits_index.begin(), output_bits_index.begin() + number_of_output_bits );

//loop through each possible input to the sbox
short number_of_times_linear_relationship_held = 0;

for(short i = 0; i < (1 << bits_in_sbox); ++i)
{
//linear sum (=xor) of our (parameter) input bits
bool xored_input_bits = 0;
for(short input_bit : normalized_inputs)
xored_input_bits ^= ith_bit(i, input_bit);

//linear sum of the output bits from the sbox
bool xored_output_bits = 0;
for(short output_bit : copy_of_output_bits)
xored_output_bits ^= ith_bit(sbox[i], output_bit);

//how well does the linear sum approximation hold for the entire range of inputs to the sbox?
if(xored_output_bits == xored_input_bits)
++number_of_times_linear_relationship_held;

}

//convert the number of times it held to the probabilities deviation from 0.5
//a perfectly random function would have bias = 0, ie it holds half the time
float bias = ((float)number_of_times_linear_relationship_held) / (1 << bits_in_sbox) - .5;

if( fabs(bias) >= default_threshold)
{
std::vector<short> copy_of_output_bits_plus_offset;
for(auto i: copy_of_output_bits) copy_of_output_bits_plus_offset.push_back(i + offset);

ret.push_back( (linear_relationship) {inputs, copy_of_output_bits_plus_offset, bias});
}

}

//loop through all possible combinations of our output bits (take k at a time)
//earlier, we tested each combination to find its bias
while(next_combination(output_bits_index.begin(),output_bits_index.begin() + number_of_output_bits,output_bits_index.end()));

}

return ret;
}

void permute_rounds(std::vector<linear_relationship> & rounds, short * pbox)
{
for(linear_relationship &line : rounds)
for(int i = 0; i < line.outputs.size(); ++i)
line.outputs[i] = pbox[line.outputs[i]];

}

//go through each combination of the input and output bits, and add a linear relation that holds 100% of the time
std::vector<linear_relationship> setup_first_round()
{
std::vector<linear_relationship> ret;

//each of the bits in an input to an sbox
std::vector<short> output_bits_index;
for(short i = 0; i < bits_in_sbox; ++i)
output_bits_index.push_back(i);

//each of the bits in an output from an sbox
std::vector<short> input_bits_index;
for(short i = 0; i < bits_in_sbox; ++i)
input_bits_index.push_back(i);

//go through each combination of the input and output bits, and add a linear relation that holds 100% of the time

//loop through the number of input bits we'll be choosing at a time
for(short number_of_input_bits = 1; number_of_input_bits <= bits_in_sbox; ++number_of_input_bits)
{
do
{
//the current combination of our output bit vector
std::vector<short> copy_of_input_bits( input_bits_index.begin(), input_bits_index.begin() + number_of_input_bits );

//loop through the number of output bits we'll be choosing at a time
for(short number_of_output_bits = 1; number_of_output_bits <= bits_in_sbox; ++number_of_output_bits)
{

do
{
//the current combination of our output bit vector
std::vector<short> copy_of_output_bits( output_bits_index.begin(), output_bits_index.begin() + number_of_output_bits );

for(int offset = 0; offset < subkey_bits; offset += bits_in_sbox)
{
std::vector<short> copy_of_output_bits_plus_offset;
for(auto i: copy_of_output_bits) copy_of_output_bits_plus_offset.push_back(i + offset);

std::vector<short> copy_of_input_bits_plus_offset;
for(auto i: copy_of_input_bits) copy_of_input_bits_plus_offset.push_back(i + offset);

ret.push_back( (linear_relationship) {copy_of_input_bits_plus_offset, copy_of_output_bits_plus_offset, .5});
}
}

//loop through all the actual output bits (combinatorics, pick k at a time)
while(next_combination(output_bits_index.begin(),output_bits_index.begin() + number_of_output_bits,output_bits_index.end()));

}
}

//loop through all the actual input bits (combinatorics, pick k at a time)
while(next_combination(input_bits_index.begin(),input_bits_index.begin() + number_of_input_bits,input_bits_index.end()));
}

return ret;
}

/*group the inputs by which sbox they are passed to
bits 0->7 go to sbox 0, bits 8->15 go to sbox 1, and so on...
{1,5, 8, 12, 20, 21, 22, 50} would become
{1, 5}, {8, 12}, {20, 21, 22}, {50}
*/
std::vector<std::vector<short>> group_inputs(linear_relationship & line)
{
std::vector<std::vector<short>> ret;
for(int index = 0; index < subkey_bits; index += bits_in_sbox)
{
std::vector<short> to_append;
for(int i: line.outputs)
{
if(index <= i && i < (index + bits_in_sbox))
to_append.push_back(i);
}

if(to_append.size() > 0)
{
std::sort(to_append.begin(), to_append.end());
ret.push_back(to_append);
}
}

return ret;
}

//to combine two linear_relationship's, you just multiply the biases, and then multiply that by two
//you combine all the output bits
linear_relationship merge_linear_relationships(Vi & output, linear_relationship base)
{
base.outputs.clear();
base.outputs.reserve(output.size());

for(auto i: output)
{
base.bias *= 2 * i.bias;
base.outputs.insert(base.outputs.end(), i.outputs.begin(), i.outputs.end());
}
return base;
}

std::vector<linear_relationship> substitute_rounds(std::vector<linear_relationship> & rounds, short ** sboxes)
{
std::vector<linear_relationship> ret;

for(auto line: rounds)
{
std::vector<std::vector<short>> grouped_inputs = group_inputs(line);
Vvi next_round_linearity;

for(auto i:grouped_inputs)
{
short * sbox = sboxes[ i/bits_in_sbox  ]; //picks the sbox to use based on the index of the grouped subsection of the input bits
next_round_linearity.push_back(bit(i, sbox));
}

Vvi output;
cart_product(output, next_round_linearity);

for(auto i: output)
ret.push_back(merge_linear_relationships(i, line));

}

return ret;
}

int main()
{
short SBox[] = { 0x63, 0x7C, 0x77, 0x7B, 0xF2, 0x6B, 0x6F, 0xC5, 0x30, 0x01, 0x67, 0x2B, 0xFE, 0xD7, 0xAB, 0x76, 0xCA, 0x82, 0xC9, 0x7D, 0xFA, 0x59, 0x47, 0xF0, 0xAD, 0xD4, 0xA2, 0xAF, 0x9C, 0xA4, 0x72, 0xC0, 0xB7, 0xFD, 0x93, 0x26, 0x36, 0x3F, 0xF7, 0xCC, 0x34, 0xA5, 0xE5, 0xF1, 0x71, 0xD8, 0x31, 0x15, 0x04, 0xC7, 0x23, 0xC3, 0x18, 0x96, 0x05, 0x9A, 0x07, 0x12, 0x80, 0xE2, 0xEB, 0x27, 0xB2, 0x75, 0x09, 0x83, 0x2C, 0x1A, 0x1B, 0x6E, 0x5A, 0xA0, 0x52, 0x3B, 0xD6, 0xB3, 0x29, 0xE3, 0x2F, 0x84, 0x53, 0xD1, 0x00, 0xED, 0x20, 0xFC, 0xB1, 0x5B, 0x6A, 0xCB, 0xBE, 0x39, 0x4A, 0x4C, 0x58, 0xCF, 0xD0, 0xEF, 0xAA, 0xFB, 0x43, 0x4D, 0x33, 0x85, 0x45, 0xF9, 0x02, 0x7F, 0x50, 0x3C, 0x9F, 0xA8, 0x51, 0xA3, 0x40, 0x8F, 0x92, 0x9D, 0x38, 0xF5, 0xBC, 0xB6, 0xDA, 0x21, 0x10, 0xFF, 0xF3, 0xD2, 0xCD, 0x0C, 0x13, 0xEC, 0x5F, 0x97, 0x44, 0x17, 0xC4, 0xA7, 0x7E, 0x3D, 0x64, 0x5D, 0x19, 0x73, 0x60, 0x81, 0x4F, 0xDC, 0x22, 0x2A, 0x90, 0x88, 0x46, 0xEE, 0xB8, 0x14, 0xDE, 0x5E, 0x0B, 0xDB, 0xE0, 0x32, 0x3A, 0x0A, 0x49, 0x06, 0x24, 0x5C, 0xC2, 0xD3, 0xAC, 0x62, 0x91, 0x95, 0xE4, 0x79, 0xE7, 0xC8, 0x37, 0x6D, 0x8D, 0xD5, 0x4E, 0xA9, 0x6C, 0x56, 0xF4, 0xEA, 0x65, 0x7A, 0xAE, 0x08, 0xBA, 0x78, 0x25, 0x2E, 0x1C, 0xA6, 0xB4, 0xC6, 0xE8, 0xDD, 0x74, 0x1F, 0x4B, 0xBD, 0x8B, 0x8A, 0x70, 0x3E, 0xB5, 0x66, 0x48, 0x03, 0xF6, 0x0E, 0x61, 0x35, 0x57, 0xB9, 0x86, 0xC1, 0x1D, 0x9E, 0xE1, 0xF8, 0x98, 0x11, 0x69, 0xD9, 0x8E, 0x94, 0x9B, 0x1E, 0x87, 0xE9, 0xCE, 0x55, 0x28, 0xDF, 0x8C, 0xA1, 0x89, 0x0D, 0xBF, 0xE6, 0x42, 0x68, 0x41, 0x99, 0x2D, 0x0F, 0xB0, 0x54, 0xBB, 0x16};

//in this case, all the sboxes are identical
//However, I wanted to support the case where they're different for different bits
short * SBoxes[subkey_bits/bits_in_sbox] = {SBox,SBox,SBox,SBox,SBox,SBox,SBox,SBox};

short PBox[] = {0, 51, 54, 12, 1, 55, 20, 6, 59, 56, 7, 17, 18, 39, 32, 21, 50, 25, 46, 5, 40, 34, 35, 43, 37, 26, 33, 28, 4, 14, 29, 42, 58, 41, 38, 24, 19, 23, 10, 13, 11, 22, 49, 3, 9, 48, 45, 8, 30, 47, 63, 16, 60, 15, 36, 27, 52, 44, 31, 53, 57, 2, 61, 62};

auto rounds = setup_first_round();
rounds = substitute_rounds(rounds, SBoxes);
}


Combinations of a vector. I didn't write this and don't need it reviewed. Just including it for completeness!

//all combinations of the input bits (order does not matter)
//take #number at a time
//stolen from stack overflow
template <typename Iterator>
inline bool next_combination(const Iterator first, Iterator k, const Iterator last)
{
/* Credits: Thomas Draper */
if ((first == last) || (first == k) || (last == k))
return false;
Iterator itr1 = first;
Iterator itr2 = last;
++itr1;
if (last == itr1)
return false;
itr1 = last;
--itr1;
itr1 = k;
--itr2;
while (first != itr1)
{
if (*--itr1 < *itr2)
{
Iterator j = k;
while (!(*itr1 < *j)) ++j;
std::iter_swap(itr1,j);
++itr1;
++j;
itr2 = k;
std::rotate(itr1,j,last);
while (last != j)
{
++j;
++itr2;
}
std::rotate(k,itr2,last);
return true;
}
}
std::rotate(first,k,last);
return false;
}


Cartesian product function

//cartesion product of vector of vectors
//from stackoverflow
typedef std::vector<linear_relationship> Vi;
typedef std::vector<Vi> Vvi;
struct Digits {
Vi::const_iterator begin;
Vi::const_iterator end;
Vi::const_iterator me;
};
typedef std::vector<Digits> Vd;

void cart_product(Vvi & out,  Vvi & in)

{
Vd vd;

// Start all of the iterators at the beginning.
for(Vvi::const_iterator it = in.begin();
it != in.end();
++it) {
Digits d = {(*it).begin(), (*it).end(), (*it).begin()};
vd.push_back(d);
}

while(1) {

// Construct your first product vector by pulling
// out the element of each vector via the iterator.
Vi result;
for(Vd::const_iterator it = vd.begin();
it != vd.end();
it++) {
result.push_back(*(it->me));
}
out.push_back(result);

// Increment the rightmost one, and repeat.

// When you reach the end, reset that one to the beginning and
// increment the next-to-last one. You can get the "next-to-last"
// iterator by pulling it out of the neighboring element in your
// vector of iterators.
for(Vd::iterator it = vd.begin(); ; ) {
// okay, I started at the left instead. sue me
++(it->me);
if(it->me == it->end) {
if(it+1 == vd.end()) {
// I'm the last digit, and I'm about to roll
return;
} else {
it->me = it->begin;
++it;
}
} else {
// normal
break;
}
}
}
}


Try this version. I've made only two changes:

• Pass function parameters by const reference, rather than by (modifiable) reference, if you're promising not to modify them.

• In ranged for loops, always use for (auto&& elt : container) unless you know exactly what you're doing. Using bare auto means "make a copy of the element", which is the cause of a lot of your slowness: you're making copies of some pretty big vectors, if I understand correctly.

std::vector<linear_relationship> substitute_rounds(
const std::vector<linear_relationship>& rounds,
const short **sboxes)
{
std::vector<linear_relationship> ret;

for (auto&& line: rounds) {
std::vector<std::vector<short>> grouped_inputs = group_inputs(line);
Vvi next_round_linearity;

for (auto&& i : grouped_inputs) {
// picks the sbox to use based on the index of the grouped subsection of the input bits
const short *sbox = sboxes[i / bits_in_sbox];
next_round_linearity.push_back(bit(i, sbox));
}

Vvi output;
cart_product(output, next_round_linearity);

for (auto&& i : output) {
ret.push_back(merge_linear_relationships(i, line));
}
}

return ret;
}


For more on auto&& and what it means exactly (and why it's the proper default when you don't care about the details), see Stephan T. Lavavej's N3994 "Ranged For Loops: The Next Generation". Basically, it'll make a reference with the proper degree of constness and rvalueness for whatever the thing is that you're trying to iterate over.

auto& or const auto& would also work in this case, but not necessarily in all cases. The big thing in this case is to stay away from reference-less auto, because that makes copies.

• Thanks, this made it a lot faster. However, it's somehow still slower than pure Python! I changed all the parameters (except one) to const, and the ranged for loops to auto&&. Are there any improvements I can make to bit() or group_inputs() ? – robertkin Jun 27 '15 at 18:04

In response to this comment: Yes.

The first refinement you could make would be to avoid push_back (which might need to resize the result vector): just give the result vector the right size initially, and then index into it.

Also in this revision, I've used std::move to indicate that we want to transfer ownership of the vector to_append; we don't need to keep a copy for ourselves. This saves another big vector copy.

std::vector<std::vector<short>> group_inputs(const linear_relationship& line)
{
std::vector<std::vector<short>> result(subkey_bits);
for(int index = 0; index < subkey_bits; index += bits_in_sbox) {
std::vector<short> to_append;
for (int i : line.outputs) {
if (index <= i && i < (index + bits_in_sbox)) {
to_append.push_back(i);
}
}

std::sort(to_append.begin(), to_append.end());
result[index] = std::move(to_append);
}

return result;
}


However, we can do better by changing the algorithm from O(kn) to O(n):

std::vector<std::vector<short>> group_inputs(const linear_relationship& line)
{
std::vector<std::vector<short>> result(subkey_bits);
for (int i : line.outputs) {
result[i / bits_in_sbox].push_back(i);
}
for (auto&& v : result) {
std::sort(v.begin(), v.end());
}
return result;
}


If line.outputs.size() is expected to be much greater than bits_in_sbox * subkey_bits — let's say, greater than 5000 — then I would switch your algorithm again to get rid of the std::sort. Look up counting sort and think about how to apply it here.