Contest: fastest way to sort a big array of Gaussian-distributed data Following the interest in this question, I thought it would be interesting to make answers a bit more objective and quantitative by proposing a contest.

The idea is simple: I have generated a binary file containing 50 million gaussian-distributed doubles (average: 0, stdev 1). The goal is to make a program that will sort these in memory as fast as possible. A very simple reference implementation in python takes 1m4s to complete. How low can we go?

The rules are as follows: answer with a program that opens the file "gaussian.dat" and sorts the numbers in memory (no need to output them), and instructions for building and running the program. The program must be able to work on my Arch Linux machine (meaning you can use any programming language or library that is easily installable on this system).

The program must be reasonably readable, so that I can make sure it is safe to launch (no assembler-only solution, please!).

I will run the answers on my machine (quad core, 4 Gigabytes of RAM). The fastest solution will get the accepted answer and a 100 points bounty :)

The program used to generate the numbers:

#!/usr/bin/env python
import random
from array import array
from sys import argv
count=int(argv[1])
a=array('d',(random.gauss(0,1) for x in xrange(count)))
f=open("gaussian.dat","wb")
a.tofile(f)


The simple reference implementation:

#!/usr/bin/env python
from array import array
from sys import argv
count=int(argv[1])
a=array('d')
a.fromfile(open("gaussian.dat"),count)
print "sorting..."
b=sorted(a)


EDIT: only 4 GB of RAM, sorry

EDIT#2: Note that the point of the contest is to see if we can use prior information about the data. it is not supposed to be a pissing match between different programming language implementations!
Here is a solution in C++ which first partitions the numbers into buckets with same expected number of elements and then sorts each bucket separately. It precomputes a table of the cumulative distribution function based on some formulas from Wikipedia and then interpolates values from this table to get a fast approximation.

Several steps run in multiple threads to make use of the four cores.

#include &lt;cstdlib&gt;
#include &lt;math.h&gt;
#include &lt;stdio.h&gt;
#include &lt;algorithm&gt;

#include &lt;tbb/parallel_for.h&gt;

using namespace std;

typedef unsigned long long ull;

double signum(double x) {
    return (x&lt;0) ? -1 : (x&gt;0) ? 1 : 0;
}

const double fourOverPI = 4 / M_PI;

double erf(double x) {
    double a = 0.147;
    double x2 = x*x;
    double ax2 = a*x2;
    double f1 = -x2 * (fourOverPI + ax2) / (1 + ax2);
    double s1 = sqrt(1 - exp(f1));
    return signum(x) * s1;
}

const double sqrt2 = sqrt(2);

double cdf(double x) {
    return 0.5 + erf(x / sqrt2) / 2;
}

const int cdfTableSize = 200;
const double cdfTableLimit = 5;
double* computeCdfTable(int size) {
    double* res = new double[size];
    for (int i = 0; i &lt; size; ++i) {
        res[i] = cdf(cdfTableLimit * i / (size - 1));
    }
    return res;
}
const double* const cdfTable = computeCdfTable(cdfTableSize);

double cdfApprox(double x) {
    bool negative = (x &lt; 0);
    if (negative) x = -x;
    if (x &gt; cdfTableLimit) return negative ? cdf(-x) : cdf(x);
    double p = (cdfTableSize - 1) * x / cdfTableLimit;
    int below = (int) p;
    if (p == below) return negative ? -cdfTable[below] : cdfTable[below];
    int above = below + 1;
    double ret = cdfTable[below] +
            (cdfTable[above] - cdfTable[below])*(p - below);
    return negative ? 1 - ret : ret;
}

void print(const double* arr, int len) {
    for (int i = 0; i &lt; len; ++i) {
        printf("%e; ", arr[i]);
    }
    puts("");
}

void print(const int* arr, int len) {
    for (int i = 0; i &lt; len; ++i) {
        printf("%d; ", arr[i]);
    }
    puts("");
}

void fillBuckets(int N, int bucketCount,
        double* data, int* partitions,
        double* buckets, int* offsets) {
    for (int i = 0; i &lt; N; ++i) {
        ++offsets[partitions[i]];
    }

    int offset = 0;
    for (int i = 0; i &lt; bucketCount; ++i) {
        int t = offsets[i];
        offsets[i] = offset;
        offset += t;
    }
    offsets[bucketCount] = N;

    int next[bucketCount];
    memset(next, 0, sizeof(next));
    for (int i = 0; i &lt; N; ++i) {
        int p = partitions[i];
        int j = offsets[p] + next[p];
        ++next[p];
        buckets[j] = data[i];
    }
}

class Sorter {
public:
    Sorter(double* data, int* offsets) {
        this-&gt;data = data;
        this-&gt;offsets = offsets;
    }

    static void radixSort(double* arr, int len) {
        ull* encoded = (ull*)arr;
        for (int i = 0; i &lt; len; ++i) {
            ull n = encoded[i];
            if (n &amp; signBit) {
                n ^= allBits;
            } else {
                n ^= signBit;
            }
            encoded[i] = n;
        }

        const int step = 11;
        const ull mask = (1ull &lt;&lt; step) - 1;
        int offsets[8][1ull &lt;&lt; step];
        memset(offsets, 0, sizeof(offsets));

        for (int i = 0; i &lt; len; ++i) {
            for (int b = 0, j = 0; b &lt; 64; b += step, ++j) {
                int p = (encoded[i] &gt;&gt; b) &amp; mask;
                ++offsets[j][p];
            }
        }

        int sum[8] = {0};
        for (int i = 0; i &lt;= mask; i++) {
            for (int b = 0, j = 0; b &lt; 64; b += step, ++j) {
                int t = sum[j] + offsets[j][i];
                offsets[j][i] = sum[j];
                sum[j] = t;
            }
        }

        ull* copy = new ull[len];
        ull* current = encoded;
        for (int b = 0, j = 0; b &lt; 64; b += step, ++j) {
            for (int i = 0; i &lt; len; ++i) {
                int p = (current[i] &gt;&gt; b) &amp; mask;
                copy[offsets[j][p]] = current[i];
                ++offsets[j][p];
            }

            ull* t = copy;
            copy = current;
            current = t;
        }

        if (current != encoded) {
            for (int i = 0; i &lt; len; ++i) {
                encoded[i] = current[i];
            }
        }

        for (int i = 0; i &lt; len; ++i) {
            ull n = encoded[i];
            if (n &amp; signBit) {
                n ^= signBit;
            } else {
                n ^= allBits;
            }
            encoded[i] = n;
        }
    }

    void operator() (tbb::blocked_range&lt;int&gt;&amp; range) const {
        for (int i = range.begin(); i &lt; range.end(); ++i) {
            double* begin = &amp;data[offsets[i]];
            double* end = &amp;data[offsets[i+1]];
            //std::sort(begin, end);
            radixSort(begin, end-begin);
        }
    }

private:
    double* data;
    int* offsets;
    static const ull signBit = 1ull &lt;&lt; 63;
    static const ull allBits = ~0ull;
};

void sortBuckets(int bucketCount, double* data, int* offsets) {
    Sorter sorter(data, offsets);
    tbb::blocked_range&lt;int&gt; range(0, bucketCount);
    tbb::parallel_for(range, sorter);
    //sorter(range);
}

class Partitioner {
public:
    Partitioner(int bucketCount, double* data, int* partitions) {
        this-&gt;data = data;
        this-&gt;partitions = partitions;
        this-&gt;bucketCount = bucketCount;
    }

    void operator() (tbb::blocked_range&lt;int&gt;&amp; range) const {
        for (int i = range.begin(); i &lt; range.end(); ++i) {
            double d = data[i];
            int p = (int) (cdfApprox(d) * bucketCount);
            partitions[i] = p;
        }
    }

private:
    double* data;
    int* partitions;
    int bucketCount;
};

const int bucketCount = 512;
int offsets[bucketCount + 1];

int main(int argc, char** argv) {
    if (argc != 2) {
        printf("Usage: %s N\n N = the size of the input\n", argv[0]);
        return 1;
    }

    puts("initializing...");
    int N = atoi(argv[1]);
    double* data = new double[N];
    double* buckets = new double[N];
    memset(offsets, 0, sizeof(offsets));
    int* partitions = new int[N];

    puts("loading data...");
    FILE* fp = fopen("gaussian.dat", "rb");
    if (fp == 0 || fread(data, sizeof(*data), N, fp) != N) {
        puts("Error reading data");
        return 1;
    }
    //print(data, N);

    puts("assigning partitions...");
    tbb::parallel_for(tbb::blocked_range&lt;int&gt;(0, N),
            Partitioner(bucketCount, data, partitions));

    puts("filling buckets...");
    fillBuckets(N, bucketCount, data, partitions, buckets, offsets);
    data = buckets;

    puts("sorting buckets...");
    sortBuckets(bucketCount, data, offsets);

    puts("done.");

    /*
    for (int i = 0; i &lt; N-1; ++i) {
        if (data[i] &gt; data[i+1]) {
            printf("error at %d: %e &gt; %e\n", i, data[i], data[i+1]);
        }
    }
    */

    //print(data, N);

    return 0;
}


To compile and run it, use this command:

g++ -O3 -ltbb -o gsort gsort.cpp &amp;&amp; time ./gsort 50000000


EDIT: All buckets are now placed into the same array to remove the need to copy the buckets back into the array. Also the size of the table with precomputed values was reduced, because the values are accurate enough. Still, if I change the number of buckets above 256, the program takes longer to run than with that number of buckets.

EDIT: The same algorithm, different programming language. I used C++ instead of Java and the running time reduced from ~3.2s to ~2.35s on my machine. The optimal number of buckets is still around 256 (again, on my computer).

By the way, tbb really is awesome.

EDIT: I was inspired by Alexandru's great solution and replaced the std::sort in the last phase by a modified version of his radix sort. I did use a different method to deal with the positive/negative numbers, even though it needs more passes through the array. I also decided to sort the array exactly and remove the insertion sort. I will later spend some time testing how do these changes influence performance and possibly revert them. However, by using radix sort, the time decreased from ~2.35s to ~1.63s.