Computing the number of primitive sorting networks on n elements - seeking tiny optimizations

My goal here is to compute the 16th term of the seqeunce . Approximately 7 years ago I computed the 12th term and it took 250 GB of RAM and 24 hours. Now I've reduced computing the 12th term to 1 GB of RAM and 45 minutes with an iterative algorithm instead of recursive. The minimum amount of space I'll require to compute the 16th term is about 87 terabytes, which hopefully will become feasible sometime in the near future. However, with the fastest processors my most optimistic estimate for the time of the computation is about 4 years. My main reason for asking this question is trying to reduce that time to something manageable, like 6 months.

The main loop has $n!$ iterations and the inner loop then has as many as $2^{(n/2)}$ iterations. I reduced the time from 5.8 hours for $n=12$ to 45 minutes with tiny optimizations in the inner loop. I'm hoping someone else can find some additional tiny optimizations that I missed.

Here is the code.

#include <iostream>
#include <stdlib.h>
#include <sys/timeb.h>
#include <string.h>
#include <sstream>

using namespace std;

unsigned long *factorial,*fact2diff;
unsigned long long **pages;

int getMilliCount(){
timeb tb;
ftime(&tb);
return tb.millitm + (tb.time & 0xfffff) * 1000;
}

inline void push(unsigned long *stack, int &stackSize, unsigned long value)
{
stack[stackSize++] = value;
}

inline unsigned long pop(unsigned long *stack, int &stackSize)
{
return stack[--stackSize];
}

inline void push(int *stack, int &stackSize, int value)
{
stack[stackSize++] = value;
}

inline int pop(int *stack, int &stackSize)
{
return stack[--stackSize];
}

inline unsigned long long getAndDelete(unsigned long long **pages,unsigned long pageSize, unsigned long index)
{
unsigned long long value= pages[index/pageSize][index%pageSize];
if(index%pageSize==pageSize-1)
delete pages[index/pageSize];
return value;
}

inline void add(unsigned long long **pages, unsigned long pageSize, unsigned long index, unsigned long long value)
{
unsigned long page=index/pageSize;
unsigned long pageIndex = index%pageSize;

if(!pages[page])
{
pages[page] = new unsigned long long[pageSize];
memset(pages[page],0,sizeof(unsigned long long)*pageSize);
}
pages[page][pageIndex] += value;
}

inline unsigned long swapAscent(unsigned long cov,int i, int n)
{
unsigned long c2 = (i==n-2)?0:(cov%factorial[n-2-i])/factorial[n-3-i];
unsigned long c1 = (cov%factorial[n-1-i])/factorial[n-2-i];

if(c1<=c2)
{
return cov+factorial[n-2-i]+(c2-c1)*fact2diff[i];
}
}

int main(int argc,char **argv)
{
int n = atoi(argv[1]);
unsigned long pageSize = atol(argv[2]);
unsigned long fact=1;

factorial = new unsigned long[n];

factorial[0]=fact;
for(int i=1;i<n;i++)
{
fact*=i+1;
factorial[i]=fact;
}
fact2diff = new unsigned long[n-1];
for(int i=0;i<n-1;i++)
{
fact2diff[i] = factorial[n-2-i]-(i<n-2?factorial[n-3-i]:0);
}
unsigned long numPages = factorial[n-1]/pageSize+1;

pages = new unsigned long long*[numPages];
memset(pages,0,sizeof(unsigned long long*)*numPages);
cout<<"Pages "<<numPages<<endl;
cout<<"Pagesize "<<pageSize<<endl;
unsigned long *stack = new unsigned long[(1<<(n/2))];
int *spotStack = new int[(1<<(n/2))];

int stackSize = 0, spotStackSize=0;
time_t starttime=getMilliCount();
unsigned long long value=0;
unsigned long long v=0;

for(unsigned long cov=0;cov<factorial[n-1];cov++)
{
value = getAndDelete(pages,pageSize,cov);
for(int i=0;i<n-1;i++)
{
v=swapAscent(cov,i,n);
{
push(stack,stackSize,v);
push(spotStack,spotStackSize,i);
}
}
while(stackSize>0)
{
unsigned long f=pop(stack,stackSize);
int spot=pop(spotStack,spotStackSize);
int absspot = abs(spot);
if(spot<0)
{
for(int j=absspot+2;j<n-1;j++)
{
v=swapAscent(f,j,n);
{
push(stack,stackSize,v);
push(spotStack,spotStackSize,j);
}
}
}
else
{
for(int j=absspot+2;j<n-1;j++)
{
v=swapAscent(f,j,n);
{
push(stack,stackSize,v);
push(spotStack,spotStackSize,-j);
}
}
}
}
}
time_t endtime = getMilliCount();

cout<<(endtime-starttime)/1000.0<<endl;
cout<<value<<endl;
return 0;
}


The program takes two arguments. The first argument is n, and the second argument is the size of a page of memory (the unit is the size of unsigned long long, which I'm working with as 16 bytes). If anyone wants to run this code, my first test of any new optimization is to run

cccompute 10 1024


On my laptop this takes a little over 18 seconds at the moment, and on a faster computer it takes under 5 seconds.

A slightly more intensive test is

cccompute 11 131072


This takes 4 minutes on my laptop and about 70 seconds on a faster computer I've been using.

The code is based on my recurrence relation detailed in a paper that is linked on the sequence web page. I can explain it to anyone who wants an explanation, but I feel the kind of optimizations I'm looking for can likely be found without fully understanding the basis of the computation, which is why I didn't include an explanation here.

As an example, I just found an optimization that decreased the execution time by 25%. In calling swapAscent(cov,i+1, n) after swapAscent(cov,i, n) there's a redundancy in that c1 in the second call is the same as c2 in the first call. I expanded the function inline and eliminated this extra calculation.

• Comments are not for extended discussion; this conversation has been moved to chat. – Mathieu Guindon Dec 13 '17 at 22:04
• @snb index if it finds one. It was the most efficient way I found to enumerate the subsets. (Got a phone call as I was typing this.) – Matt Samuel Dec 13 '17 at 22:08

Names are key to maintainability

stack is a useful name for a datatype, but as a name for a variable it is opaque.

    unsigned long *stack = new unsigned long[(1<<(n/2))];


tells me that unsigned longs will be accessed in LIFO order, but doesn't tell me what each unsigned long means. Is it a count? An encoded permutation?

    int *spotStack = new int[(1<<(n/2))];


is a slight improvement, but without any hint as to what a spot is it still leaves me quite lost.

    ... unsigned long cov=0 ...
unsigned long long value=0;
unsigned long long v=0;


are all completely opaque.

But they also lead me into a subpoint: type names matter almost as much as variable names. You're using a language which has typedef: a typedef unsigned long long permutation_t; would go a long way to making the code easier to figure out. They might also make it slightly more obvious that the types are being mixed in ugly ways, which I suspect are unintentional and not noticed simply because your compiler reduces them to the same internal type:

   unsigned long *stack = new unsigned long[(1<<(n/2))];
...
unsigned long long v=0;
...
push(stack,stackSize,v);


Either stack should be unsigned long long* or v should be unsigned long. Given the signature of push I suspect it's the latter.

Finally on names,

(the unit is the size of unsigned long long, which I'm working with as 16 bytes)

Rather than using unsigned long and unsigned long long and making assumptions about the exact size of one of them, you're using a language which has <cstdint>. I would guess that your unsigned longs could be std::uint_fast32_t and your unsigned long longs could be std::uint64_t respectively, although it might be even more complicated than that. typedefing each logical type would not only provide clarity on the meaning of variables but also help you see where you can use fast types vs least types vs exact types.

Use unsigned only where appropriate

                add(pages,pageSize,f,-value);


The algorithm is doing an inclusion-exclusion, so that - really is a unary minus. Both value and pages should be using signed types, not unsigned.

Simplify edge cases

    unsigned long c2 = (i==n-2)?0:(cov%factorial[n-2-i])/factorial[n-3-i];

    factorial = new unsigned long[n];

factorial[0]=fact;
for(int i=1;i<n;i++)
{
fact*=i+1;
factorial[i]=fact;
}


Why not reindex factorial so that factorial[n] == n! (rather than, at present, factorial[n] == (n+1)!)? Then you can lose the special case in swapAscent. And as a bonus the variable name will be more transparent.

Optimise algorithms before loops

I've spent half an hour on the bus thinking more about the problem you're trying to solve, and I think that there are still algorithmic improvements which could be made before worrying about micro-optimisation.

            for(int j=$MINJ;j<n-1;j++) { v=swapAscent($PERM,j,n);
{
push(stack,stackSize,v);
push(spotStack,spotStackSize,j);
}
}


in essence occurs three times, although as I understand it that's really a recursive function which has been transformed to an iterative function with a stack. Conceptually viewing it as a single recursive function, I think it's a traversal through a directed acyclic graph. What is certainly the case is that the same recursive call is being made multiple times with different values of value (or, from the graph perspective, different values of value are being pushed down the same edge to be accumulated at the destination vertex).

Idea 1: if it's a DAG then a topological sort takes linear time (assuming that the structure of this particular graph doesn't allow parts of it to be done implicitly to get sublinear time), and then each edge of the graph only needs to be traversed once more. The memory cost here would be an additional $n! \;\textrm{sizeof}(n!)$.

Idea 2: whether it's a DAG or not, iterative deepening could be used to replace the inner loop (executed up to $2^{n/2}$ times) with an outer loop (executed up to, I think, $n/2$ times). The memory cost would be a doubling, to have a second paged array containing "accumulated value not yet pushed to my neighbours" per vertex.

Idea 3: if we were to label each vertex $v$ with a vector $d_v$ such that $d_v[i]$ is the number of paths from the root to $v$ of length $i$, I'm pretty sure we'd find isomorphic vertices (under a natural interpretation: same label and isomorphic neighbourhood). If there's an efficient way to identify isomorphic vertices or to canonicalise to a representative of the equivalence class then this would allow reducing the $n!$ of the outer loop to the number of equivalence classes (surely a saving of a least a factor of $n$), and in the best case would also save memory in the same proportion. [Edit: this turns out not to be as good as I hoped: the average equivalence class has only about 2.2 elements].

Idea 4: exploiting the particular structure of the graph, I think there's a kind of meet-in-the-middle which might be able to halve the memory usage. To be precise, I think that there's a symmetry which relates cov and factorial[n-1] - cov in such a way that if each of those pairs is swapped and each edge is reversed then the original graph is recovered. This is suggestive that it might be possible to work with only half of the graph by working forwards and then backwards.

• I need at least 48-bit integers to represent the permutations for n=16, and 128-bit integers to be used as values. – Matt Samuel Dec 13 '17 at 16:22
• Thanks for this - I'm also finding that the opaque names are a barrier to unravelling the algorithm. – Toby Speight Dec 13 '17 at 16:22
• @MattSamuel, of course, I should have calculated that 16 bytes is 128 bits. So you can't just use <cstdint>, but __uint128_t is still clearer than unsigned long long. – Peter Taylor Dec 13 '17 at 16:25
• I got it down to 178 seconds from 204 seconds for n=11 on my laptop with some optimizations that would probably make everyone here angry! :-P – Matt Samuel Dec 13 '17 at 17:06
• @MattSamuel That's not how it works, you should never have to sacrifice code readability for speed ups, and if you are doing "optimizations" that would make us angry that don't have to with code readability, in all likely hood you found the wrong place to optimize and there was probably another place that would have made that optimization obsolete. – opa Dec 13 '17 at 18:03

You can use whitespace to make the logic clearer; removing it won't make the program any faster.

#include <sys/timeb.h>


The man page for ftime(3) here says:

This function is obsolete. Don't use it.

Prefer <chrono> or <ctime> instead. A simple std::chrono timer looks like

#include <chrono>

    auto const start_time = std::chrono::steady_clock::now();

auto const result = count_crossings(n);

std::chrono::duration<double, std::chrono::seconds::period> time_taken
= end_time - start_time;

std::cout << time_taken.count() << " seconds" << std::endl;


These are specified for C compatibility, but deprecated:

#include <string.h>
#include <stdlib.h>


When writing new code, use the versions that declare names in the std namespace:

#include <cstring>
// No need for <cstdlib> - only used for ato¤(); use sto¤() instead
// No need for <sstream> - never used


Avoid using namespace

Never, ever write this:

using namespace std;


Prefer strongly-typed constants to macros:

constexpr auto BADPERM = ~0ul;


Avoid global variables

unsigned long *factorial,*fact2diff;
unsigned long long **pages;


The second of those is particularly problematic, as it is frequently shadowed, so the reader has to constantly check which is being used.

Don't reinvent wheels

inline void push(unsigned long *stack, int &stackSize, unsigned long value);
inline unsigned long pop(unsigned long *stack, int &stackSize);
inline void push(int *stack, int &stackSize, int value);
inline int pop(int *stack, int &stackSize);


That's why std::stack was invented. And the standard containers let you provide your own allocators if you need them; you may want create an allocator that uses mmap() rather than simply new[]; that will allow you to align to the hardware page size and may help you avoid fragmentation and other housekeeping if done well.

Pair new[] and delete[]

Valgrind picks up immediately that you need delete[] instead of delete here:

if (index%pageSize==pageSize-1)
delete[] pages[index/pageSize];


Better still, use a smart pointer or a collection rather than managing your own memory - doing so would solve the leaks of stack and spotstack.

Check the arguments

Don't just dereference argv[1] and argv[2] without checking:

if (argc < 3) {
std::cerr << "Usage: " << *argv << " <n> <pagesize>" << std::endl;
return 1;
}

int n = std::stoi(argv[1]);
size_t pageSize = std::stoul(argv[2]);
if (n <= 0 || pageSize == 0) {
std::cerr << "Usage: " << *argv << " <n> <pagesize>" << std::endl;
return 1;
}


Non-literal zero is not a null pointer

std::memset(pages,0,sizeof(unsigned long long*)*numPages);


is not a safe way to fill with null pointers. Prefer the correctly-typed algorithm:

std::fill_n(pages, numPages, nullptr);


Similarly, replace

    std::memset(pages[page],0,sizeof(unsigned long long)*pageSize);


with

    std::fill_n(pages[page], pageSize, 0ull);


You'll need to include <algorithm> instead of <cstring> to get std::fill_n().

• Considering your answer. One thing about the reinventing the wheel part - the thing is that my stack has at most exactly $2^{n/2}$ elements. I feel it's overkill to use an STL container when there's no reallocation required. Also, what do you recommend instead of global variables that would be faster? – Matt Samuel Dec 13 '17 at 11:34
• stack can quite obviously be in the scope of main(), since that's the only place it's used. factorial and fact2diff are used only by main() and swapAscent() (which is called only from main()), so they are easily passed by reference. Any good compiler will optimise the passing for you - don't micro-optimize without both profiling execution and inspecting the assembly code. – Toby Speight Dec 13 '17 at 12:01
• Nevermind, I googled it. – Matt Samuel Dec 13 '17 at 12:41
• I'm looking for examples on how to use mmap the way you mentioned and I'm only finding mapping files. Do you have a link handy on how to do this? – Matt Samuel Dec 13 '17 at 12:53
• Yes - look at MAP_ANONYMOUS in the mmap(2) man page. Consider also the Boost wrappers. – Toby Speight Dec 13 '17 at 13:11

I see this code in swapAscent():

unsigned long c2 = (i==n-2)?0:(cov%factorial[n-2-i])/factorial[n-3-i];
unsigned long c1 = (cov%factorial[n-1-i])/factorial[n-2-i];


Since n is constant during the run of the program, and i is just a loop index (either from i or j), you know exactly when the ternary's condition will be true: the last iteration of the for i/j loop. It will always return BADPERM. So you don't even need to call it then. Just change the termination test in each for i/j loop. Then get rid of the ternary in swapAscent().

I'm also wondering if it's possible to skip the divisions in swapAscent(), but doubtful as it's all integer division. In any event, I note the repeated pattern between the lines for c1 and c2 and wonder if any work can be saved.

• Yes, I just posted that I found an optimization in setting c1 to the old c2. I should probably reveal what these numbers are: they're encoded permutations of 1 through n. If you were to list the permutations in alphabetical order, cov is the index of the permutation on the list. swapAscent is checking if p(i)<p(i+1), and if it is then it swaps those two indices. – Matt Samuel Dec 13 '17 at 14:39
• I tried if (i==n-2) return BADPERM; there but got 0 as the program result - did I get that wrong? – Toby Speight Dec 13 '17 at 16:03
• @Toby If c1==0 then it shouldn't be BADPERM. – Matt Samuel Dec 13 '17 at 16:08