text:
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)

Only 4 GB of RAM is available.
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!

I'm a bit lost concerning your implementation. I saved the first snippet as generator.py and did python generator.py 1000000. However, as I saved the second snippet as sort.py and entered in the console python sort.py 1000000, I got a EOFError. I know nearly nothing in Python, so if you could help... Edit : The gaussian.dat file is indeed at the right place, but the serialization made it look like garbage in my editor. Normal behaviour with python's tofile()?

@Fezvez: it's supposed to look as garbage since it's stored as raw binary data. However what you did is supposed to work... What version of python did you use? I used python 2.4, which is rather old. Also, how big is the resulting file?

I used python 2.7, and the file size was 8M. If anyone cared to check with this version of python, it'll help. Anyway, that's not the point, I'll generate the data myself, and sort it myself too. It's just that I wondered why I couldn't make the baseline work for me.

@Fezvez: I'll try tonight on a newer machine. In the meantime, can you report the exact size in bytes? If it's exactly 8,000,000 bytes, then the file is probably good.

Hehe : 7,62 MB (8 000 000 bytes)! In the meantime, I'll work on the sorting algorithm =)

@Fezvez: I just tried it again with Python 2.7, and it works fine...

I believe that using the Gaussian distribution needlessly complicates things, if you want to do smart calculations. If we have the uniform distribution, the point (using the knowledge of the distribution) stays the same, and the arithmetics gets much simpler...

@xofon: I picked the Gaussian distribution because it's what is most widely encountered in practice.

@xofon: I would be surprised if the direct indexing method described by Mike didn't hands down win such a method, while in theory with this method a strategy to avoid the cost of the phi method could be created.

Take each value and move it directly to its "expected" position, repeat for the displaced value. Not sure how to resolve a couple issue with that. When done, bubble sort until complete (a couple passes should do).

I will post a bucket sort solution tomorrow evening if this hasn't been closed by then :)

@static_rtti - as a heavy CG user, this is exactly the sort of thing "we" like to hack on at CG.SE. To any reading mods, move this to CG, don't close it.

Welcome to CodeGolf.SE! I've cleared a lot of the commentary from the SO original concerning where this does or does not belong, and re-tagged to be closer to the CodeGolf.SE mainstream.

The one tricky issue here is that we look for objective winning criteria, and "fastest" introduces platform dependencies...does a O(n^{1.2}) algorithm implemented on the cpython virtual machine beat a O(n^{1.3}) algorithm with a similar constant implemented in c? I generally suggest some discussion on the performance characteristics of each solution, as this may help people to judge what is going on.

Thanks for the welcome! The problem you're referring to is why I run the code on my computer, so we have a single reference. I agree different languages are going to have different performance results, that is why I think discussing the algorithms is as important as looking at the raw numbers.

Hopefully the winning solution (at least!) will be tested for correctness in addition to speed?  I can write sort; exit very, very fast... as far as you know!

@Ben Jackson: of course! This is becomming increasingly needed as solutions become more and more cryptic!

Note that I'll add a bounty, so the contest should go on for about a week.

@static_rtti: Would you like to split the contest into two sections: single-threaded (wall-time) vs multi-threaded (real-time)? I have access to a 48 core machine (gcc 4.3.4) which I can use it to time multi-threaded solutions.

Why don't you give the timings on that machine? I'm not sure it would make sense to start a second question, but the results would definitely be interesting. Maybe increase the problem size as well?

Here are some timings on a 16 core machine: link. The 48 core machine is currently busy, but none of the tested solutions can take advantage of more than 4 cpus. The testing methodology is in the git repository I linked to. Anybody is welcome to add more tests.

@Alexandru, the fork/join solution I provided should definitely be able to take advantage of more than 4 cpus. Did you test it?

@arjan: Unfortunately, I don't have Java 1.7 on that machine, but 1.6. I expected the TBB solution to use all CPUs, but it wasn't the case.

I tried to evaluate all new solutions, I hope I didn't miss any! Note that I didn't have much time, so I didn't look at the solutions in depth, but I will try to do so (and also check a few ones for correctness) later. Thanks to everyone for your very nice solutions!

@Alexandru - and you don't have download rights/capabilities on that 48 core machine? You don't have to install Java 7, but can simply put it in your home dir. I'll see if I can try it tomorrow on a simple (compared to your's ;)) 12 core dev machine.

@arjan: Done. Your solution is a second close.

@Alexandru: Thanks! :) What are specs of the machine btw? I also executed my solution on our 12 core dev machine, which is ~3.46Ghz dual Intel Nehalem hex core with some 96GB RAM and I think a dual 8x SSD. Hyperthreading was disabled. To get some idea of the scaling factor, I run the code with 1 core till 12 cores and another run with 10 million doubles. See my updated answer for a table with the results of that.

OK, if nobody objects, I will retest the fastest solutions tonight, check them for correctness, and award the bounty.

Nice contest, you should probably point out that a real implementation would have to estimate the mean and standard deviation first or somehow on the fly - might change the algorithm somewhat.

@Andrew, I agree, but the idea was to focus on the sorting. Estimation is another problem, and is already well understood.

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 <cstdlib>
#include <math.h>
#include <stdio.h>
#include <algorithm>

#include <tbb/parallel_for.h>

using namespace std;

typedef unsigned long long ull;

double signum(double x) {
    return (x<0) ? -1 : (x>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 < size; ++i) {
        res[i] = cdf(cdfTableLimit * i / (size - 1));
    }
    return res;
}
const double* const cdfTable = computeCdfTable(cdfTableSize);

double cdfApprox(double x) {
    bool negative = (x < 0);
    if (negative) x = -x;
    if (x > 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 < len; ++i) {
        printf("%e; ", arr[i]);
    }
    puts("");
}

void print(const int* arr, int len) {
    for (int i = 0; i < 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 < N; ++i) {
        ++offsets[partitions[i]];
    }

    int offset = 0;
    for (int i = 0; i < 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 < 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->data = data;
        this->offsets = offsets;
    }

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

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

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

        int sum[8] = {0};
        for (int i = 0; i <= mask; i++) {
            for (int b = 0, j = 0; b < 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 < 64; b += step, ++j) {
            for (int i = 0; i < len; ++i) {
                int p = (current[i] >> b) & 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 < len; ++i) {
                encoded[i] = current[i];
            }
        }

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

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

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

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

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

    void operator() (tbb::blocked_range<int>& range) const {
        for (int i = range.begin(); i < 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<int>(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 < N-1; ++i) {
        if (data[i] > data[i+1]) {
            printf("error at %d: %e > %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 && 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.

Nice.  I got 3.055 on mine.  Lowest I was able to get mine to was 6.3.  I'm picking through yours to get the stats better.  Why did you choose 256 as the number of buckets?  I tried 128 and 512, yet 256 worked the best.

Why did I choose 256 as the number of buckets? I tried 128 and 512, yet 256 worked the best. :) I found it empirically and I am not sure why increasing the number of buckets slows down the algorithm - memory allocation should not take that long. Maybe something related to cache size?

2.725s on my machine. Pretty nice for a java solution, taking into account the loading time of the JVM.

I switched your code over to use the nio packages, per my and Arjan's solution(used his syntax, since it was cleaner than mine) and was able to get it .3 seconds faster.  I've got an ssd, I wonder what the implications might be if not.  It also gets rid of some of your bit twiddling. Modded sections are here.

@Scott: I agree that my way of loading data is suboptimal. I have not used the java.nio.* classes before, so I have had no idea how to do the loading better. My approach was possibly also slower because I tried to load the file in parts and start partitioning the parts that were already loaded before all the data was read. On my machine (no ssd) the time went from ~3.5s to ~3.2s. Thanks for your suggestion, hope you do not mind it if I use it in my program?

Have at.  You might be able to merge the two behaviors(use nio and start a thread when you have the data), if you use a different call to do the reading.

This is the fastest parallel solution on my tests (16core cpu). 1.22s far from 1.94s second place.

1.188s on my machine! I think we have a winner! (although Alexandru deserves big kudos for his extremely fast sequential version)

Thanks. I agree that Alexandru is the one who came with the idea of using radix sort on doubles. My solution uses his approach but makes it faster by: a) parallelization, b) using the cache better by decreasing the size of the arrays on which the radix sort is run.

@static_rtti https://duck.ac/problem/1001/board They sort 1e8 int32s in ~750ms no parallel, was there so much hardware improve for such kind of question in these years?

Without getting smart, just to provide a much faster naïve sorter, here's one in C which should be pretty much equivalent to your Python one:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
int cmp(const void* av, const void* bv) {
    double a = *(const double*)av;
    double b = *(const double*)bv;
    return a < b ? -1 : a > b ? 1 : 0;
}
int main(int argc, char** argv) {
    if (argc <= 1)
        return puts("No argument!");
    unsigned count = atoi(argv[1]);

    double *a = malloc(count * sizeof *a);

    FILE *f = fopen("gaussian.dat", "rb");
    if (fread(a, sizeof *a, count, f) != count)
        return puts("fread failed!");
    fclose(f);

    puts("sorting...");
    double *b = malloc(count * sizeof *b);
    memcpy(b, a, count * sizeof *b);
    qsort(b, count, sizeof *b, cmp);
    return 0;
}

Compiled with gcc -O3, on my machine this takes more than a minute less than the Python: about 11 s compared to 87 s.

Took 10.086s on my machine, which makes you the current leader! But I'm pretty sure we can do better :)

Could you try to remove the second ternary operator and simply return 1 for that case because random doubles are like not equal to each other in these amount of data.

@Codism: I would add that we don't care about swapping locations of equivalent data so therefore even if we could get equivalent values it would be an appropriate simplification.

I partitioned into segments based on the standard deviation that should best break it down into 4ths.  Edit: Rewritten to partition based on x value in http://en.wikipedia.org/wiki/Error_function#Table_of_values
http://www.wolframalpha.com/input/?i=percentages+by++normal+distribution
I tried using smaller buckets, but it seemed to have little effect once  2 * beyond the number of cores available. Without any parallel collections, it would take 37 seconds on my box and 24 with the parallel collections.  If partitioning via distribution, you can't just use an array, so there is some more overhead.  I'm not clear on when a value would be boxed/unboxed in scala.
I'm using scala 2.9, for the parallel collection.  You can just download the tar.gz distribution of it.
To compile: scalac SortFile.scala (I just copied it directly in the scala/bin folder.
To run:     JAVA_OPTS="-Xmx4096M" ./scala SortFile
(I ran it with 2 gigs of ram and got about the same time)
Edit:  Removed allocateDirect, slower than just allocate.  Removed priming of initial size for array buffers.  Actually made it read the whole 50000000 values.  Rewrote to hopefully avoid autoboxing issues (still slower than naive c)
import java.io.FileInputStream;
import java.nio.ByteBuffer
import java.nio.ByteOrder
import scala.collection.mutable.ArrayBuilder

object SortFile {

//used partition numbers from Damascus' solution
val partList = List(0, 0.15731, 0.31864, 0.48878, 0.67449, 0.88715, 1.1503, 1.5341)

val listSize = partList.size * 2;
val posZero = partList.size;
val neg = partList.map( _ * -1).reverse.zipWithIndex
val pos = partList.map( _ * 1).zipWithIndex.reverse

def partition(dbl:Double): Int = {

//for each partition, i am running through the vals in order
//could make this a binary search to be more performant... but our list size is 4 (per side)

  if(dbl < 0) { return neg.find( dbl < _._1).get._2  }
  if(dbl > 0) { return posZero  + pos.find( dbl > _._1).get._2  }
      return posZero;

}

  def main(args: Array[String])
    {

    var l = 0
    val dbls = new Array[Double](50000000)
    val partList = new Array[Int](50000000)
    val pa = Array.fill(listSize){Array.newBuilder[Double]}
    val channel = new FileInputStream("../../gaussian.dat").getChannel()
    val bb = ByteBuffer.allocate(50000000 * 8)
    bb.order(ByteOrder.LITTLE_ENDIAN)
    channel.read(bb)
    bb.rewind
    println("Loaded" + System.currentTimeMillis())
    var dbl = 0.0
    while(bb.hasRemaining)
    {
      dbl = bb.getDouble
      dbls.update(l,dbl)

      l+=1
    }
    println("Beyond first load" + System.currentTimeMillis());

    for( i <- (0 to 49999999).par) { partList.update(i, partition(dbls(i)))}

    println("Partition computed" + System.currentTimeMillis() )
    for(i <- (0 to 49999999)) { pa(partList(i)) += dbls(i) }
    println("Partition completed " + System.currentTimeMillis())
    val toSort = for( i <- pa) yield i.result()
    println("Arrays Built" + System.currentTimeMillis());
    toSort.par.foreach{i:Array[Double] =>scala.util.Sorting.quickSort(i)};

    println("Read\t" + System.currentTimeMillis());

  }
}

8.185s! Nice for a scala solution, I guess... Also, bravo for providing the first solution that actually uses the Gaussian distribution in some way!

I was only aiming to compete with the c# solution.  Didn't figure I'd beat c/c++.  Also.. it's behaving much differently for you than for me.  I'm using openJDK on my end and it's a lot slower.  I wonder if adding more partitions would help in your env.

Just put this in a cs file and compile it with csc in theory: (Requires mono)
using System;
using System.IO;
using System.Threading;

namespace Sort
{
    class Program
    {
        const int count = 50000000;
        static double[][] doubles;
        static WaitHandle[] waiting = new WaitHandle[4];
        static AutoResetEvent[] events = new AutoResetEvent[4];

        static double[] Merge(double[] left, double[] right)
        {
            double[] result = new double[left.Length + right.Length];
            int l = 0, r = 0, spot = 0;
            while (l < left.Length && r < right.Length)
            {
                if (right[r] < left[l])
                    result[spot++] = right[r++];
                else
                    result[spot++] = left[l++];
            }
            while (l < left.Length) result[spot++] = left[l++];
            while (r < right.Length) result[spot++] = right[r++];
            return result;
        }

        static void ThreadStart(object data)
        {
            int index = (int)data;
            Array.Sort(doubles[index]);
            events[index].Set();
        }

        static void Main(string[] args)
        {
            System.Diagnostics.Stopwatch watch = new System.Diagnostics.Stopwatch();
            watch.Start();
            byte[] bytes = File.ReadAllBytes(@"..\..\..\SortGuassian\Data.dat");
            doubles = new double[][] { new double[count / 4], new double[count / 4], new double[count / 4], new double[count / 4] };
            for (int i = 0; i < 4; i++)
            {
                for (int j = 0; j < count / 4; j++)
                {
                    doubles[i][j] = BitConverter.ToDouble(bytes, i * count/4 + j * 8);
                }
            }
            Thread[] threads = new Thread[4];
            for (int i = 0; i < 4; i++)
            {
                threads[i] = new Thread(ThreadStart);
                waiting[i] = events[i] = new AutoResetEvent(false);
                threads[i].Start(i);
            }
            WaitHandle.WaitAll(waiting);
            double[] left = Merge(doubles[0], doubles[1]);
            double[] right = Merge(doubles[2], doubles[3]);
            double[] result = Merge(left, right);
            watch.Stop();
            Console.WriteLine(watch.Elapsed.ToString());
            Console.ReadKey();
        }
    }
}

Can I run your solutions with Mono? How should I do it?

Haven't used Mono, didn't think of that, you should be able to compile the F# and then run it.

Updated to use four threads to improve performance. Now gives me 6 sec. Note that this could be significantly improved (5 sec likely) if you use only one spare array and avoid initializing a ton of memory to zero, which is done by the CLR, since everything is being written to at least once.

9.598s on my machine! You're the current leader :)

(Note that I patched your code to remove the stopwatch and the console stuff: I measure time with time)

@static_rtti: That's fine I just like excluding bootup of CLR etc, but I can see where you would want an overall time for this kind of thing.

@Guvante: it's mostly because it's simpler to do it that way rather than to request that everyone implements a stopwatch mechanism and ensure that the comparison is fair.

My mother told me to stay away from guys with Mono!

Since you know what the distribution is, you can use a direct indexing O(N) sort.
(If you're wondering what that is, suppose you have a deck of 52 cards and you want to sort it. Just have 52 bins and toss each card into it's own bin.)
You have 5e7 doubles. Allocate a result array R of 5e7 doubles. Take each number x and get i = phi(x) * 5e7. Basically do R[i] = x. Have a way to handle collisions, such as moving the number it may be colliding with (as in simple hash coding). Alternatively, you could make R a few times larger, filled with a unique empty value.
At the end, you just sweep up the elements of R.
phi is just the gaussian cumulative distribution function.
It converts a gaussian distributed number between +/- infinity into a uniform distributed number between 0 and 1.
A simple way to calculate it is with table lookup and interpolation.

Be careful: you know the approximate distribution, not the exact distribution. You know the data was generated using a Gaussian law, but since it is finite, it doesn't exactly follow a Gaussian.

@static_rtti: In this case the necessary approximation of phi would create a bigger hassle than any irregularities in the data set IMO.

@static_rtti: it doesn't have to be exact. It only has to spread out the data so it is approximately uniform, so it doesn't bunch up too much in some places.

Suppose you have 5e7 doubles.  Why not just make each entry in R a vector of, say, 5e6 vectors of double.  Then, push_back each double in its appropriate vector.  Sort the vectors and you're done.  This should take linear time in the size of the input.

Actually, I see that mdkess already came up with that solution.

Here is another sequential solution:
#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
#include <ctime>

typedef unsigned long long ull;

int size;
double *dbuf, *copy;
int cnt[8][1 << 16];

void sort()
{
  const int step = 10;
  const int start = 24;
  ull mask = (1ULL << step) - 1;

  ull *ibuf = (ull *) dbuf;
  for (int i = 0; i < size; i++) {
    for (int w = start, v = 0; w < 64; w += step, v++) {
      int p = (~ibuf[i] >> w) & mask;
      cnt[v][p]++;
    }
  }

  int sum[8] = { 0 };
  for (int i = 0; i <= mask; i++) {
    for (int w = start, v = 0; w < 64; w += step, v++) {
      int tmp = sum[v] + cnt[v][i];
      cnt[v][i] = sum[v];
      sum[v] = tmp;
    }
  }

  for (int w = start, v = 0; w < 64; w += step, v++) {
    ull *ibuf = (ull *) dbuf;
    for (int i = 0; i < size; i++) {
      int p = (~ibuf[i] >> w) & mask;
      copy[cnt[v][p]++] = dbuf[i];
    }

    double *tmp = copy;
    copy = dbuf;
    dbuf = tmp;
  }

  for (int p = 0; p < size; p++)
    if (dbuf[p] >= 0.) {
      std::reverse(dbuf + p, dbuf + size);
      break;
    }

  // Insertion sort
  for (int i = 1; i < size; i++) {
    double value = dbuf[i];
    if (value < dbuf[i - 1]) {
      dbuf[i] = dbuf[i - 1];
      int p = i - 1;
      for (; p > 0 && value < dbuf[p - 1]; p--)
        dbuf[p] = dbuf[p - 1];
      dbuf[p] = value;
    }
  }
}

int main(int argc, char **argv) {
  size = atoi(argv[1]);
  dbuf = new double[size];
  copy = new double[size];

  FILE *f = fopen("gaussian.dat", "r");
  fread(dbuf, size, sizeof(double), f);
  fclose(f);

  clock_t c0 = clock();
  sort();
  printf("Finished after %.3f\n", (double) ((clock() - c0)) / CLOCKS_PER_SEC);
  return 0;
}

I doubt it beats the multi-threaded solution, but the timings on my i7 laptop are (stdsort is the C++ solution provided in another answer):
$ g++ -O3 mysort.cpp -o mysort && ./mysort 50000000
Finished after 2.10
$ g++ -O3 stdsort.cpp -o stdsort && ./stdsort
Finished after 7.12

Note that this solution has linear time complexity (because it uses the special representation of doubles).
EDIT: Fixed the order of elements to be increasing.
EDIT: Improved speed by almost half a second.
EDIT: Improved speed by another 0.7 seconds. Made the algorithm more cache friendly.
EDIT: Improved speed by another 1 second. Since there only 50.000.000 elements I can partially sort the mantissa and use insert sort (which is cache friendly) to fix out-of-place elements. This idea removes around two iterations from the last radix sorting loop.
EDIT: 0.16 less seconds. First std::reverse can be eliminated if the sorting order is reversed.

Now that is getting interesting! What kind of sort algorithm is it?

Least significant digit radix sort. You can sort the mantissa, then the exponent, then the sign. The algorithm presented here takes this idea one step further. It can be parallelized using a partitioning idea provided in a different answer.

Pretty fast for a single threaded solution: 2.552s! Do you think you could change your solution to make use of the fact that the data is normally distributed? You could probably do better than the current best multi-threaded solutions.

@static_rtti: I see that Damascus Steel already posted a multithreaded version of this implementation. I improved caching behavior of this algorithm, so you should get better timings now. Please test this new version.

1.525s!! Amazing, you're the leader now. It's nice to see that a sequential solution actually beats the best parallel solution.

1.459s in my latests tests. While this solution is not the winner per my rules, it really deserves big kudos. Congratulations!

Taking Christian Ammer's solution and parallelizing it with Intel's Threaded Building Blocks
#include <iostream>
#include <fstream>
#include <algorithm>
#include <vector>
#include <ctime>
#include <tbb/parallel_sort.h>

int main(void)
{
    std::ifstream ifs("gaussian.dat", std::ios::binary | std::ios::in);
    std::vector<double> values;
    values.reserve(50000000);
    double d;
    while (ifs.read(reinterpret_cast<char*>(&d), sizeof(double)))
    values.push_back(d);
    clock_t c0 = clock();
    tbb::parallel_sort(values.begin(), values.end());
    std::cout << "Finished after "
              << static_cast<double>((clock() - c0)) / CLOCKS_PER_SEC
              << std::endl;
}

If you have access to Intel's Performance Primitives (IPP) library, you can use its radix sort.  Just replace
#include <tbb/parallel_sort.h>

with
#include "ipps.h"

and
tbb::parallel_sort(values.begin(), values.end());

with
std::vector<double> copy(values.size());
ippsSortRadixAscend_64f_I(&values[0], &copy[0], values.size());

On my dual core laptop, the timings are
C               16.4 s
C#              20 s
C++ std::sort   7.2 s
C++ tbb         5 s
C++ ipp         4.5 s
python          too long

2.958s! TBB seems pretty cool and easy to use!

TBB is absurdly awesome. It's exactly the right level of abstraction for algorithmic work.

How about an implementation of parallel quicksort that chooses its pivot values based on the statistics of the distribution, thereby ensuring equal sized partitions?  The first pivot would be at the mean (zero in this case), the next pair would be at the 25th and 75th percentiles (+/- -0.67449 standard deviations), and so on, with each partition halving the remaining data set more or less perfectly.

That's effectively what I did on mine.. of course you got this post up before I could finish my writeup.

Very ugly (why use arrays when I can use variables ending with numbers), but fast code (my first try to std::threads), whole time (time real) on my system 1,8 s (comparing to std::sort() 4,8 s), compile with g++ -std=c++0x -O3 -march=native -pthread
Just pass data through stdin (works only for 50M).
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <thread>
using namespace std;
const size_t size=50000000;

void pivot(double* start,double * end, double middle,size_t& koniec){
    double * beg=start;
    end--;
    while (start!=end){
        if (*start>middle) swap (*start,*end--);
        else start++;
    }
    if (*end<middle) start+=1;
    koniec= start-beg;
}
void s(double * a, double* b){
    sort(a,b);
}
int main(){
    double *data=new double[size];
    FILE *f = fopen("gaussian.dat", "rb");
    fread(data,8,size,f);
    size_t end1,end2,end3,temp;
    pivot(data, data+size,0,end2);
    pivot(data, data+end2,-0.6745,end1);
    pivot(data+end2,data+size,0.6745,end3);
    end3+=end2;
    thread ts1(s,data,data+end1);
    thread ts2(s,data+end1,data+end2);
    thread ts3(s,data+end2,data+end3);
    thread ts4(s,data+end3,data+size);
    ts1.join(),ts2.join(),ts3.join(),ts4.join();
    //for (int i=0; i<size-1; i++){
    //  if (data[i]>data[i+1]) cerr<<"BLAD\n";
    //}
    fclose(f);
    //fwrite(data,8,size,stdout);
}

//Edit changed to read gaussian.dat file.

Could you change it to read gaussian.dat, as the above C++ solutions do?

I'll try later when I come home.

Very nice solution, you're the current leader (1.949s)! And nice use of the gaussian distribution :)

A C++ solution using std::sort (eventually faster than qsort, regarding Performance of qsort vs std::sort)
#include <iostream>
#include <fstream>
#include <algorithm>
#include <vector>
#include <ctime>

int main(void)
{
    std::ifstream ifs("C:\\Temp\\gaussian.dat", std::ios::binary | std::ios::in);
    std::vector<double> values;
    values.reserve(50000000);
    double d;
    while (ifs.read(reinterpret_cast<char*>(&d), sizeof(double)))
        values.push_back(d);
    clock_t c0 = clock();
    std::sort(values.begin(), values.end());
    std::cout << "Finished after "
              << static_cast<double>((clock() - c0)) / CLOCKS_PER_SEC
              << std::endl;
}

I can't reliable say how long it takes because I only have 1GB on my machine and with the given Python code I could only make a gaussian.dat file with only 25mio doubles (without getting a Memory error). But I'm very interested how long the std::sort algorithm runs.

6.425s! As expected, C++ shines :)

@static_rtti: I've tried swensons Timsort algorithm (as suggested from Matthieu M. in your first question). I had to make some changes to the sort.h file to compile it with C++. It was about twice as slow than std::sort. Don't know why, maybe because of compiler optimizations?

Here is a mix of Alexandru's radix sort with Zjarek's threaded smart pivoting.  Compile it with
g++ -std=c++0x -pthread -O3 -march=native sorter_gaussian_radix.cxx -o sorter_gaussian_radix

You can change the radix size by defining STEP (e.g. add -DSTEP=11).  I found the best for my laptop is 8 (the default).
By default, it splits the problem into 4 pieces and runs that on multiple threads.  You can change that by passing a depth parameter to the command line.  So if you have two cores, run it as
sorter_gaussian_radix 50000000 1

and if you have 16 cores
sorter_gaussian_radix 50000000 4

The max depth right now is 6 (64 threads).  If you put too many levels, you will just slow the code down.
One thing I also tried was the radix sort from the Intel Performance Primitives (IPP) library.  Alexandru's implementation soundly trounces IPP, with IPP being about 30% slower.  That variation is also included here (commented out).
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <ctime>
#include <iostream>
#include <thread>
#include <vector>
#include <boost/cstdint.hpp>
// #include "ipps.h"

#ifndef STEP
#define STEP 8
#endif

const int step = STEP;
const int start_step=24;
const int num_steps=(64-start_step+step-1)/step;
int size;
double *dbuf, *copy;

clock_t c1, c2, c3, c4, c5;

const double distrib[]={-2.15387,
                        -1.86273,
                        -1.67594,
                        -1.53412,
                        -1.4178,
                        -1.31801,
                        -1.22986,
                        -1.15035,
                        -1.07752,
                        -1.00999,
                        -0.946782,
                        -0.887147,
                        -0.830511,
                        -0.776422,
                        -0.724514,
                        -0.67449,
                        -0.626099,
                        -0.579132,
                        -0.53341,
                        -0.488776,
                        -0.445096,
                        -0.40225,
                        -0.36013,
                        -0.318639,
                        -0.27769,
                        -0.237202,
                        -0.197099,
                        -0.157311,
                        -0.11777,
                        -0.0784124,
                        -0.0391761,
                        0,
                        0.0391761,
                        0.0784124,
                        0.11777,
                        0.157311,
                        0.197099,
                        0.237202,
                        0.27769,
                        0.318639,
                        0.36013,
                        0.40225,
                        0.445097,
                        0.488776,
                        0.53341,
                        0.579132,
                        0.626099,
                        0.67449,
                        0.724514,
                        0.776422,
                        0.830511,
                        0.887147,
                        0.946782,
                        1.00999,
                        1.07752,
                        1.15035,
                        1.22986,
                        1.31801,
                        1.4178,
                        1.53412,
                        1.67594,
                        1.86273,
                        2.15387};

class Distrib
{
  const int value;
public:
  Distrib(const double &v): value(v) {}

  bool operator()(double a)
  {
    return a<value;
  }
};

void recursive_sort(const int start, const int end,
                    const int index, const int offset,
                    const int depth, const int max_depth)
{
  if(depth<max_depth)
    {
      Distrib dist(distrib[index]);
      const int middle=std::partition(dbuf+start,dbuf+end,dist) - dbuf;

      // const int middle=
      //   std::partition(dbuf+start,dbuf+end,[&](double a)
      //                  {return a<distrib[index];})
      //   - dbuf;

      std::thread lower(recursive_sort,start,middle,index-offset,offset/2,
                        depth+1,max_depth);
      std::thread upper(recursive_sort,middle,end,index+offset,offset/2,
                        depth+1,max_depth);
      lower.join(), upper.join();
    }
  else
    {
  // ippsSortRadixAscend_64f_I(dbuf+start,copy+start,end-start);

      c1=clock();

      double *dbuf_local(dbuf), *copy_local(copy);
      boost::uint64_t mask = (1 << step) - 1;
      int cnt[num_steps][mask+1];

      boost::uint64_t *ibuf = reinterpret_cast<boost::uint64_t *> (dbuf_local);

      for(int i=0;i<num_steps;++i)
        for(uint j=0;j<mask+1;++j)
          cnt[i][j]=0;

      for (int i = start; i < end; i++)
        {
          for (int w = start_step, v = 0; w < 64; w += step, v++)
            {
              int p = (~ibuf[i] >> w) & mask;
              (cnt[v][p])++;
            }
        }

      c2=clock();

      std::vector<int> sum(num_steps,0);
      for (uint i = 0; i <= mask; i++)
        {
          for (int w = start_step, v = 0; w < 64; w += step, v++)
            {
              int tmp = sum[v] + cnt[v][i];
              cnt[v][i] = sum[v];
              sum[v] = tmp;
            }
        }

      c3=clock();

      for (int w = start_step, v = 0; w < 64; w += step, v++)
        {
          ibuf = reinterpret_cast<boost::uint64_t *>(dbuf_local);

          for (int i = start; i < end; i++)
            {
              int p = (~ibuf[i] >> w) & mask;
              copy_local[start+((cnt[v][p])++)] = dbuf_local[i];
            }
          std::swap(copy_local,dbuf_local);
        }

      // Do the last set of reversals
      for (int p = start; p < end; p++)
        if (dbuf_local[p] >= 0.)
          {
            std::reverse(dbuf_local+p, dbuf_local + end);
            break;
          }

      c4=clock();

      // Insertion sort
      for (int i = start+1; i < end; i++) {
        double value = dbuf_local[i];
        if (value < dbuf_local[i - 1]) {
          dbuf_local[i] = dbuf_local[i - 1];
          int p = i - 1;
          for (; p > 0 && value < dbuf_local[p - 1]; p--)
            dbuf_local[p] = dbuf_local[p - 1];
          dbuf_local[p] = value;
        }
      }
      c5=clock();

    }
}

int main(int argc, char **argv) {
  size = atoi(argv[1]);
  copy = new double[size];

  dbuf = new double[size];
  FILE *f = fopen("gaussian.dat", "r");
  fread(dbuf, size, sizeof(double), f);
  fclose(f);

  clock_t c0 = clock();

  const int max_depth= (argc > 2) ? atoi(argv[2]) : 2;

  // ippsSortRadixAscend_64f_I(dbuf,copy,size);

  recursive_sort(0,size,31,16,0,max_depth);

  if(num_steps%2==1)
    std::swap(dbuf,copy);

  // for (int i=0; i<size-1; i++){
  //   if (dbuf[i]>dbuf[i+1])
  //     std::cout << "BAD "
  //               << i << " "
  //               << dbuf[i] << " "
  //               << dbuf[i+1] << " "
  //               << "\n";
  // }

  std::cout << "Finished after "
            << (double) (c1 - c0) / CLOCKS_PER_SEC << " "
            << (double) (c2 - c1) / CLOCKS_PER_SEC << " "
            << (double) (c3 - c2) / CLOCKS_PER_SEC << " "
            << (double) (c4 - c3) / CLOCKS_PER_SEC << " "
            << (double) (c5 - c4) / CLOCKS_PER_SEC << " "
            << "\n";

  // delete [] dbuf;
  // delete [] copy;
  return 0;
}

EDIT: I implemented Alexandru's cache improvements, and that shaved off about 30% of the time on my machine.
EDIT: This implements a recursive sort, so it should work well on Alexandru's 16 core machine.  It also use's Alexandru's last improvement and removes one of the reverse's.  For me, this gave a 20% improvement.
EDIT: Fixed a sign bug that caused inefficiency when there are more than 2 cores.
EDIT: Removed the lambda, so it will compile with older versions of gcc.  It includes the IPP code variation commented out.  I also fixed the documentation for running on 16 cores.  As far as I can tell, this is the fastest implementation.
EDIT: Fixed a bug when STEP is not 8.  Increased the maximum number of threads to 64.  Added some timing info.

Nice. Radix sort is very cache unfriendly. See if you can get better results by changing step (11 was optimal on my laptop).

You have a bug: int cnt[mask] should be int cnt[mask + 1]. For better results use a fixed value int cnt[1 << 16].

I'll try all these solutions later today when I get home.

1.534s!!! I think we have a leader :-D

@static_rtti: Could you try this again?  It has gotten significantly faster than the last time you tried it.  On my machine, it is substantially faster than any other solution.

Sorry, I could'nt try it yesterday because for some reason, the code was gone. I'll try again tonight when I come home.

I get 1.124s on my machine which is actually the fastest solution! I feel bad for not giving you the bounty, but I swear the code wasn't there when I tested a lot of solutions yesterday! Anyways, you have the moral victory, and I will look closely at your code!

@static_rtti Thanks for trying it out.  Do not worry about the bounty.  It turns out that in my last revision, I deleted the old version and forgot to include the new version.  I do feel that the code can be further optimized, but it would be fairly small gains and specific to a particular machine.

I guess this really depends on what you want to do.  If you want to sort a bunch of Gaussians, then this won't help you.  But if you want a bunch of sorted Gaussians, this will.  Even if this misses the problem a bit, I think it will be interesting to compare vs actual sorting routines.
If you want to something to be fast, do less.
Instead of generating a bunch of random samples from the normal distribution and then sorting, you can generate a bunch of samples from the normal distribution in sorted order.
You can use the solution here to generate n uniform random numbers in sorted order.  Then you can use the inverse cdf (scipy.stats.norm.ppf) of the normal distribution to turn the uniform random numbers into numbers from the normal distribution via inverse transform sampling.
import scipy.stats
import random

# slightly modified from linked stackoverflow post
def n_random_numbers_increasing(n):
  """Like sorted(random() for i in range(n))),
  but faster because we avoid sorting."""
  v = 1.0
  while n:
    v *= random.random() ** (1.0 / n)
    yield 1 - v
    n -= 1

def n_normal_samples_increasing(n):
  return map(scipy.stats.norm.ppf, n_random_numbers_increasing(n))

If you want to get your hands dirtier, I'd guess that you might be able to speed up the many inverse cdf calculations by using some kind of iterative method, and using the previous result as your initial guess.  Since the guesses are going to be very close, probably a single iteration will give you great accuracy.

Nice answer, but that would be cheating :) The idea of my question is that while sorting algorithms have been given enormous attention, there is almost no literature on the use of prior knowledge about the data for sorting, even though the few papers that have addressed the issue have reported nice gains. So let's see what's possible!

Try this changing Guvante's solution with this Main(), it starts sorting as soon as 1/4 IO reading is done, it's faster in my test:
    static void Main(string[] args)
    {
        FileStream filestream = new FileStream(@"..\..\..\gaussian.dat", FileMode.Open, FileAccess.Read);
        doubles = new double[][] { new double[count / 4], new double[count / 4], new double[count / 4], new double[count / 4] };
        Thread[] threads = new Thread[4];

        for (int i = 0; i < 4; i++)
        {
            byte[] bytes = new byte[count * 4];
            filestream.Read(bytes, 0, count * 4);

            for (int j = 0; j < count / 4; j++)
            {
                doubles[i][j] = BitConverter.ToDouble(bytes, i * count/4 + j * 8);
            }

            threads[i] = new Thread(ThreadStart);
            waiting[i] = events[i] = new AutoResetEvent(false);
            threads[i].Start(i);
        }

        WaitHandle.WaitAll(waiting);
        double[] left = Merge(doubles[0], doubles[1]);
        double[] right = Merge(doubles[2], doubles[3]);
        double[] result = Merge(left, right);
        Console.ReadKey();
    }
}

8.933s. Slightly faster :)

Since you know the distribution, my idea would be to make k buckets, each with the same expected number of elements (since you know the distribution, you can compute this). Then in O(n) time, sweep the array and put elements into their buckets.
Then concurrently sort the buckets. Suppose you have k buckets, and n elements. A bucket will take (n/k) lg (n/k) time to sort. Now suppose that you have p processors that you can use. Since buckets can be sorted independently, you have a multiplier of ceil(k/p) to deal with. This gives a final runtime of n + ceil(k/p)*(n/k) lg (n/k), which should be a good deal faster than n lg n if you choose k well.

I think this is the best solution.

You don't exactly know the number of elements who will end up in a bucket, so the math is actually wrong. That being said, this is a good answer I think.

@pouejapon: You're right.

This answer sounds really nice.  The problem is -- it isn't really fast.  I implemented this in C99 (see my answer), and it certainly easily beats std::sort(), but it is way slower than Alexandru's radixsort solution.

One low-level optimization idea is to fit two doubles in a SSE register, so each thread would work with two items at a time. This might be complicated to do for some algorithms.
Another thing to do is sorting the array in cache-friendly chunks, then merging the results. Two levels should be used: for example first 4 KB for L1 then 64 KB for L2.
This should be very cache-friendly, since the bucket sort will not go outside the cache, and the final merge will walk memory sequentially.
These days computation is much cheaper than memory accesses. However we have a large number of items, so it's hard to tell which is the array size when dumb cache-aware sort is slower than a low-complexity non-cache-aware version.
But I won't provide an implementation of the above since I would do it in Windows (VC++).

My personal favorite using Intel's Threaded Building Blocks has already been posted, but here's a crude parallel solution using JDK 7 and its new fork/join API:
import java.io.FileInputStream;
import java.nio.channels.FileChannel;
import java.util.Arrays;
import java.util.concurrent.*;
import static java.nio.channels.FileChannel.MapMode.READ_ONLY;
import static java.nio.ByteOrder.LITTLE_ENDIAN;

/**
 *
 * Original Quicksort: https://github.com/pmbauer/parallel/tree/master/src/main/java/pmbauer/parallel
 *
 */
public class ForkJoinQuicksortTask extends RecursiveAction {

    public static void main(String[] args) throws Exception {

        double[] array = new double[Integer.valueOf(args[0])];

        FileChannel fileChannel = new FileInputStream("gaussian.dat").getChannel();
        fileChannel.map(READ_ONLY, 0, fileChannel.size()).order(LITTLE_ENDIAN).asDoubleBuffer().get(array);

        ForkJoinPool mainPool = new ForkJoinPool();

        System.out.println("Starting parallel computation");

        mainPool.invoke(new ForkJoinQuicksortTask(array));
    }

    private static final long serialVersionUID = -642903763239072866L;
    private static final int SERIAL_THRESHOLD = 0x1000;

    private final double a[];
    private final int left, right;

    public ForkJoinQuicksortTask(double[] a) {this(a, 0, a.length - 1);}

    private ForkJoinQuicksortTask(double[] a, int left, int right) {
        this.a = a;
        this.left = left;
        this.right = right;
    }

    @Override
    protected void compute() {
        if (right - left < SERIAL_THRESHOLD) {
            Arrays.sort(a, left, right + 1);
        } else {
            int pivotIndex = partition(a, left, right);
            ForkJoinTask<Void> t1 = null;

            if (left < pivotIndex)
                t1 = new ForkJoinQuicksortTask(a, left, pivotIndex).fork();
            if (pivotIndex + 1 < right)
                new ForkJoinQuicksortTask(a, pivotIndex + 1, right).invoke();

            if (t1 != null)
                t1.join();
        }
    }

    public static int partition(double[] a, int left, int right) {
        // chose middle value of range for our pivot
        double pivotValue = a[left + (right - left) / 2];

        --left;
        ++right;

        while (true) {
            do
                ++left;
            while (a[left] < pivotValue);

            do
                --right;
            while (a[right] > pivotValue);

            if (left < right) {
                double tmp = a[left];
                a[left] = a[right];
                a[right] = tmp;
            } else {
                return right;
            }
        }
    }
}

Important disclaimer: I took the quick-sort adaption for fork/join from: https://github.com/pmbauer/parallel/tree/master/src/main/java/pmbauer/parallel
To run this you need a beta build of JDK 7 (http://jdk7.java.net/download.html).
On my 2.93Ghz Quad core i7 (OS X):
Python reference
time python sort.py 50000000
sorting...

real    1m13.885s
user    1m11.942s
sys     0m1.935s

Java JDK 7 fork/join
time java ForkJoinQuicksortTask 50000000
Starting parallel computation

real    0m2.404s
user    0m10.195s
sys     0m0.347s

I also tried to do some experimenting with parallel reading and converting the bytes to doubles, but I did not see any difference there.
Update:
If anyone wants to experiment with parallel loading of the data, the parallel loading version is below. In theory this could make it go a little faster still, if your IO device has enough parallel capacity (SSDs usually do). There's also some overhead in creating Doubles from bytes, so that could potentially go faster in parallel as well. On my systems (Ubuntu 10.10/Nehalem Quad/Intel X25M SSD, and OS X 10.6/i7 Quad/Samsung SSD) I did not see any real difference.
import static java.nio.ByteOrder.LITTLE_ENDIAN;
import static java.nio.channels.FileChannel.MapMode.READ_ONLY;

import java.io.FileInputStream;
import java.nio.DoubleBuffer;
import java.nio.channels.FileChannel;
import java.util.Arrays;
import java.util.concurrent.ForkJoinPool;
import java.util.concurrent.ForkJoinTask;
import java.util.concurrent.RecursiveAction;

/**
 *
 * Original Quicksort: https://github.com/pmbauer/parallel/tree/master/src/main/java/pmbauer/parallel
 *
 */
public class ForkJoinQuicksortTask extends RecursiveAction {

   public static void main(String[] args) throws Exception {

       ForkJoinPool mainPool = new ForkJoinPool();

       double[] array = new double[Integer.valueOf(args[0])];
       FileChannel fileChannel = new FileInputStream("gaussian.dat").getChannel();
       DoubleBuffer buffer = fileChannel.map(READ_ONLY, 0, fileChannel.size()).order(LITTLE_ENDIAN).asDoubleBuffer();

       mainPool.invoke(new ReadAction(buffer, array, 0, array.length));
       mainPool.invoke(new ForkJoinQuicksortTask(array));
   }

   private static final long serialVersionUID = -642903763239072866L;
   private static final int SERIAL_THRESHOLD = 0x1000;

   private final double a[];
   private final int left, right;

   public ForkJoinQuicksortTask(double[] a) {this(a, 0, a.length - 1);}

   private ForkJoinQuicksortTask(double[] a, int left, int right) {
       this.a = a;
       this.left = left;
       this.right = right;
   }

   @Override
   protected void compute() {
       if (right - left < SERIAL_THRESHOLD) {
           Arrays.sort(a, left, right + 1);
       } else {
           int pivotIndex = partition(a, left, right);
           ForkJoinTask<Void> t1 = null;

           if (left < pivotIndex)
               t1 = new ForkJoinQuicksortTask(a, left, pivotIndex).fork();
           if (pivotIndex + 1 < right)
               new ForkJoinQuicksortTask(a, pivotIndex + 1, right).invoke();

           if (t1 != null)
               t1.join();
       }
   }

   public static int partition(double[] a, int left, int right) {
       // chose middle value of range for our pivot
       double pivotValue = a[left + (right - left) / 2];

       --left;
       ++right;

       while (true) {
           do
               ++left;
           while (a[left] < pivotValue);

           do
               --right;
           while (a[right] > pivotValue);

           if (left < right) {
               double tmp = a[left];
               a[left] = a[right];
               a[right] = tmp;
           } else {
               return right;
           }
       }
   }

}

class ReadAction extends RecursiveAction {

   private static final long serialVersionUID = -3498527500076085483L;

   private final DoubleBuffer buffer;
   private final double[] array;
   private final int low, high;

   public ReadAction(DoubleBuffer buffer, double[] array, int low, int high) {
       this.buffer = buffer;
       this.array = array;
       this.low = low;
       this.high = high;
   }

   @Override
   protected void compute() {
       if (high - low < 100000) {
           buffer.position(low);
           buffer.get(array, low, high-low);
       } else {
           int middle = (low + high) >>> 1;

           invokeAll(new ReadAction(buffer.slice(), array, low, middle),  new ReadAction(buffer.slice(), array, middle, high));
       }
   }
}

Update2:
I executed the code on one of our 12 core dev machines with a slight modification to set a fixed amount of cores. This gave the following results:
Cores  Time
1      7.568s
2      3.903s
3      3.325s
4      2.388s
5      2.227s
6      1.956s
7      1.856s
8      1.827s
9      1.682s
10     1.698s
11     1.620s
12     1.503s

On this system I also tried the Python version which took 1m2.994s and Zjarek's C++ version which took 1.925s (for some reason Zjarek's C++ version seems to run relatively faster on static_rtti's computer).
I also tried what happened if I doubled the file size to 100,000,000 doubles:
Cores  Time
1      15.056s
2      8.116s
3      5.925s
4      4.802s
5      4.430s
6      3.733s
7      3.540s
8      3.228s
9      3.103s
10     2.827s
11     2.784s
12     2.689s

In this case, Zjarek's C++ version took 3.968s. Python just took too long here.
150,000,000 doubles:
Cores  Time
1      23.295s
2      12.391s
3      8.944s
4      6.990s
5      6.216s
6      6.211s
7      5.446s
8      5.155s
9      4.840s
10     4.435s
11     4.248s
12     4.174s

In this case, Zjarek's C++ version was 6.044s. I didn't even attempt Python.
The C++ version is very consistent with its results, where Java swings a little. First it gets a little more efficient when the problem gets larger, but then less efficient again.

This code does not parse the double values correctly for me. Is Java 7 required to correctly parse the values from the file?

Ah, silly me. I forgot to set the endianness again after I locally refactored the IO code from multiple lines to one. Java 7 would normally be needed, unless you added fork/join separately to Java 6 of course.

3.411s on my machine. Not bad, but slower than koumes21's java solution :)

I'll try koumes21's solution here too locally to see what the relative differences are on my system. Anyway, no shame in 'loosing' from koumes21 since it's a much more clever solution. This is just an almost standard quick-sort thrown into a fork/join pool ;)

Here's a linear scan bucket sort implementation. I think it's faster than all current single-threaded implementations except for radix sort. It should have linear expected running time if I'm estimating the cdf accurately enough (I'm using linear interpolation of values I found on the web) and haven't made any mistakes that would cause excessive scanning:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <ctime>

using std::fill;

const double q[] = {
  0.0,
  9.865E-10,
  2.8665150000000003E-7,
  3.167E-5,
  0.001349898,
  0.022750132,
  0.158655254,
  0.5,
  0.8413447460000001,
  0.9772498679999999,
  0.998650102,
  0.99996833,
  0.9999997133485,
  0.9999999990134999,
  1.0,
};
int main(int argc, char** argv) {
  if (argc <= 1)
    return puts("No argument!");
  unsigned count = atoi(argv[1]);
  unsigned count2 = 3 * count;

  bool *ba = new bool[count2 + 1000];
  fill(ba, ba + count2 + 1000, false);
  double *a = new double[count];
  double *c = new double[count2 + 1000];

  FILE *f = fopen("gaussian.dat", "rb");
  if (fread(a, 8, count, f) != count)
    return puts("fread failed!");
  fclose(f);

  int i;
  int j;
  bool s;
  int t;
  double z;
  double p;
  double d1;
  double d2;
  for (i = 0; i < count; i++) {
    s = a[i] < 0;
    t = a[i];
    if (s) t--;
    z = a[i] - t;
    t += 7;
    if (t < 0) {
      t = 0;
      z = 0;
    } else if (t >= 14) {
      t = 13;
      z = 1;
    }
    p = q[t] * (1 - z) + q[t + 1] * z;
    j = count2 * p;
    while (ba[j] && c[j] < a[i]) {
      j++;
    }
    if (!ba[j]) {
      ba[j] = true;
      c[j] = a[i];
    } else {
      d1 = c[j];
      c[j] = a[i];
      j++;
      while (ba[j]) {
        d2 = c[j];
        c[j] = d1;
        d1 = d2;
        j++;
      }
      c[j] = d1;
      ba[j] = true;
    }
  }
  i = 0;
  int max = count2 + 1000;
  for (j = 0; j < max; j++) {
    if (ba[j]) {
      a[i++] = c[j];
    }
  }
  // for (i = 0; i < count; i += 1) {
  //   printf("here %f\n", a[i]);
  // }
  return 0;
}

I'll try this later today when I get home. In the meantime, can I say your code is very ugly? :-D

3.071s! Not bad for a single-threaded solution!

I don't know, why I can't edit my previous post, so here's new version, 0,2 seconds faster (but about 1,5 s faster in CPU time (user)). This solution have 2 programs, first precalculates quantiles for normal distribution for bucket sort, and stores it in table, t[double * scale] = bucket index, where scale is some arbitrary number which makes casting to double possible. Then main program can use this data to put doubles in correct bucket. It has one drawback, if data is not gaussian it will not work correctly (and also there is almost zero chance to work incorrectly for normal distribution), but modification for special case is easy and fast (only number of buckets checks and falling to std::sort()).
Compiling:
g++  => http://pastebin.com/WG7pZEzH helper program
g++ -std=c++0x -O3 -march=native -pthread => http://pastebin.com/T3yzViZP main sorting program

1.621s! I think you're the leader, but I'm rapidly losing track with all these answers :)

Here is another sequential solution. This one uses the fact that the elements are normal distributed, and the I think the idea is generally applicable to get sorting close to linear time.
The algorithm is like this:

Approximate CDF (see phi() function in the implementation)
For all elements compute the approximate position in the sorted array: size * phi(x)
Put elements in a new array close to their final position

In my implementation destination array has some gaps in it so I don't have to shift too many elements when inserting.

Use insertsort to sort the final elements (insertsort is linear if distance to final position is smaller than a constant).

Unfortunately, the hidden constant is quite large and this solution is twice as slow as the radix sort algorithm.

2.470s! Very nice ideas. It doesn't matter that the solution isn't the fastest if the ideas are interesting :)

This is the same as mine, but grouping the phi computations together and the shifts together for better cache performance, right?

@jonderry: I upvoted your solution, now that I understand what it does. Didn't mean to steal your idea. I included your implementation in my (unofficial) set of tests

A version using traditional pthreads. Code for merging copied from Guvante's answer. Compile with g++ -O3 -pthread.
#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <algorithm>

static unsigned int nthreads = 4;
static unsigned int size = 50000000;

typedef struct {
  double *array;
  int size;
} array_t;

void
merge(double *left, int leftsize,
      double *right, int rightsize,
      double *result)
{
  int l = 0, r = 0, insertat = 0;
  while (l < leftsize && r < rightsize) {
    if (left[l] < right[r])
      result[insertat++] = left[l++];
    else
      result[insertat++] = right[r++];
  }

  while (l < leftsize) result[insertat++] = left[l++];
  while (r < rightsize) result[insertat++] = right[r++];
}

void *
run_thread(void *input)
{
  array_t numbers = *(array_t *)input;
  std::sort(numbers.array, numbers.array+numbers.size);
  pthread_exit(NULL);
}

int
main(int argc, char **argv)
{
  double *numbers = (double *) malloc(size * sizeof(double));

  FILE *f = fopen("gaussian.dat", "rb");
  if (fread(numbers, sizeof(double), size, f) != size)
    return printf("Reading gaussian.dat failed");
  fclose(f);

  array_t worksets[nthreads];
  int worksetsize = size / nthreads;
  for (int i = 0; i < nthreads; i++) {
    worksets[i].array=numbers+(i*worksetsize);
    worksets[i].size=worksetsize;
  }

  pthread_attr_t attributes;
  pthread_attr_init(&attributes);
  pthread_attr_setdetachstate(&attributes, PTHREAD_CREATE_JOINABLE);

  pthread_t threads[nthreads];
  for (int i = 0; i < nthreads; i++) {
    pthread_create(&threads[i], &attributes, &run_thread, &worksets[i]);
  }

  for (int i = 0; i < nthreads; i++) {
    pthread_join(threads[i], NULL);
  }

  double *tmp = (double *) malloc(size * sizeof(double));
  merge(numbers, worksetsize, numbers+worksetsize, worksetsize, tmp);
  merge(numbers+(worksetsize*2), worksetsize, numbers+(worksetsize*3), worksetsize, tmp+(size/2));
  merge(tmp, worksetsize*2, tmp+(size/2), worksetsize*2, numbers);

  /*
  printf("Verifying result..\n");
  for (int i = 0; i < size - 1; i++) {
    if (numbers[i] > numbers[i+1])
      printf("Result is not correct\n");
  }
  */

  pthread_attr_destroy(&attributes);
  return 0;
}

On my laptop I get the following results:
real    0m6.660s
user    0m9.449s
sys     0m1.160s

Here is a sequential C99 implementation that tries to really make use of the known distribution.  It basically does a single round of bucket sort using the distribution information, then a few rounds of quicksort on each bucket assuming a uniform distribution within the limits of the bucket and finally a modified selection sort to copy the data back to the original buffer.  The quicksort memorizes the split points, so selection sort only needs to operate on small cunks.  And in spite (because?) of all that complexity, it isn't even really fast.
To make evaluating Φ fast, the values are sampled in a few points and later on only linear interpolation is used.  It actually does not matter if Φ is evaluated exactly, as long as the approximation is strictly monotonic.
The bin sizes are chosen such that the chance of an bin overflow is negligible.  More precisely, with the current parameters, the chance that a dataset of 50000000 elements will cause a bin overflow is 3.65e-09.  (This can be computed using the survival function of the Poisson distribution.)
To compile, please use
gcc -std=c99 -msse3 -O3 -ffinite-math-only

Since there is considerably more computation than in the other solutions, these compiler flags are needed to make it at least reasonably fast.  Without -msse3 the conversions from double to int become really slow.  If your architecture does not support SSE3, these conversions can also be done using the lrint() function.
The code is rather ugly -- not sure if this meets the requirement of being "reasonably readable"...
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>

#define N 50000000
#define BINSIZE 720
#define MAXBINSIZE 880
#define BINCOUNT (N / BINSIZE)
#define SPLITS 64
#define PHI_VALS 513

double phi_vals[PHI_VALS];

int bin_index(double x)
{
    double y = (x + 8.0) * ((PHI_VALS - 1) / 16.0);
    int interval = y;
    y -= interval;
    return (1.0 - y) * phi_vals[interval] + y * phi_vals[interval + 1];
}

double bin_value(int bin)
{
    int left = 0;
    int right = PHI_VALS - 1;
    do
    {
        int centre = (left + right) / 2;
        if (bin < phi_vals[centre])
            right = centre;
        else
            left = centre;
    } while (right - left > 1);
    double frac = (bin - phi_vals[left]) / (phi_vals[right] - phi_vals[left]);
    return (left + frac) * (16.0 / (PHI_VALS - 1)) - 8.0;
}

void gaussian_sort(double *restrict a)
{
    double *b = malloc(BINCOUNT * MAXBINSIZE * sizeof(double));
    double **pos = malloc(BINCOUNT * sizeof(double*));
    for (size_t i = 0; i < BINCOUNT; ++i)
        pos[i] = b + MAXBINSIZE * i;
    for (size_t i = 0; i < N; ++i)
        *pos[bin_index(a[i])]++ = a[i];
    double left_val, right_val = bin_value(0);
    for (size_t bin = 0, i = 0; bin < BINCOUNT; ++bin)
    {
        left_val = right_val;
        right_val = bin_value(bin + 1);
        double *splits[SPLITS + 1];
        splits[0] = b + bin * MAXBINSIZE;
        splits[SPLITS] = pos[bin];
        for (int step = SPLITS; step > 1; step >>= 1)
            for (int left_split = 0; left_split < SPLITS; left_split += step)
            {
                double *left = splits[left_split];
                double *right = splits[left_split + step] - 1;
                double frac = (double)(left_split + (step >> 1)) / SPLITS;
                double pivot = (1.0 - frac) * left_val + frac * right_val;
                while (1)
                {
                    while (*left < pivot && left <= right)
                        ++left;
                    while (*right >= pivot && left < right)
                        --right;
                    if (left >= right)
                        break;
                    double tmp = *left;
                    *left = *right;
                    *right = tmp;
                    ++left;
                    --right;
                }
                splits[left_split + (step >> 1)] = left;
            }
        for (int left_split = 0; left_split < SPLITS; ++left_split)
        {
            double *left = splits[left_split];
            double *right = splits[left_split + 1] - 1;
            while (left <= right)
            {
                double *min = left;
                for (double *tmp = left + 1; tmp <= right; ++tmp)
                    if (*tmp < *min)
                        min = tmp;
                a[i++] = *min;
                *min = *right--;
            }
        }
    }
    free(b);
    free(pos);
}

int main()
{
    double *a = malloc(N * sizeof(double));
    FILE *f = fopen("gaussian.dat", "rb");
    assert(fread(a, sizeof(double), N, f) == N);
    fclose(f);
    for (int i = 0; i < PHI_VALS; ++i)
    {
        double x = (i * (16.0 / PHI_VALS) - 8.0) / sqrt(2.0);
        phi_vals[i] =  (erf(x) + 1.0) * 0.5 * BINCOUNT;
    }
    gaussian_sort(a);
    free(a);
}

4.098s! I had to add -lm to compile it (for erf).

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <memory.h>
#include <algorithm>

// maps [-inf,+inf] to (0,1)
double normcdf(double x) {
        return 0.5 * (1 + erf(x * M_SQRT1_2));
}

int calcbin(double x, int bins) {
        return (int)floor(normcdf(x) * bins);
}

int *docensus(int bins, int n, double *arr) {
        int *hist = calloc(bins, sizeof(int));
        int i;
        for(i = 0; i < n; i++) {
                hist[calcbin(arr[i], bins)]++;
        }
        return hist;
}

void partition(int bins, int *orig_counts, double *arr) {
        int *counts = malloc(bins * sizeof(int));
        memcpy(counts, orig_counts, bins*sizeof(int));
        int *starts = malloc(bins * sizeof(int));
        int b, i;
        starts[0] = 0;
        for(i = 1; i < bins; i++) {
                starts[i] = starts[i-1] + counts[i-1];
        }
        for(b = 0; b < bins; b++) {
                while (counts[b] > 0) {
                        double v = arr[starts[b]];
                        int correctbin;
                        do {
                                correctbin = calcbin(v, bins);
                                int swappos = starts[correctbin];
                                double tmp = arr[swappos];
                                arr[swappos] = v;
                                v = tmp;
                                starts[correctbin]++;
                                counts[correctbin]--;
                        } while (correctbin != b);
                }
        }
        free(counts);
        free(starts);
}

void sortbins(int bins, int *counts, double *arr) {
        int start = 0;
        int b;
        for(b = 0; b < bins; b++) {
                std::sort(arr + start, arr + start + counts[b]);
                start += counts[b];
        }
}

void checksorted(double *arr, int n) {
        int i;
        for(i = 1; i < n; i++) {
                if (arr[i-1] > arr[i]) {
                        printf("out of order at %d: %lf %lf\n", i, arr[i-1], arr[i]);
                        exit(1);
                }
        }
}

int main(int argc, char *argv[]) {
        if (argc == 1 || argv[1] == NULL) {
                printf("Expected data size as argument\n");
                exit(1);
        }
        int n = atoi(argv[1]);
        const int cachesize = 128 * 1024; // a guess
        int bins = (int) (1.1 * n * sizeof(double) / cachesize);
        if (argc > 2) {
                bins = atoi(argv[2]);
        }
        printf("Using %d bins\n", bins);
        FILE *f = fopen("gaussian.dat", "rb");
        if (f == NULL) {
                printf("Couldn't open gaussian.dat\n");
                exit(1);
        }
        double *arr = malloc(n * sizeof(double));
        fread(arr, sizeof(double), n, f);
        fclose(f);

        int *counts = docensus(bins, n, arr);
        partition(bins, counts, arr);
        sortbins(bins, counts, arr);
        checksorted(arr, n);

        return 0;
}

This uses erf() to place each element appropriately into a bin, then sorts each bin.  It keeps the array entirely in-place.
First pass: docensus() counts the
    number of elements in each bin.
Second pass: partition() permutes the array, placing each element into its proper bin
Third pass: sortbins() performs a qsort on each bin.
It's kind of naive, and calls the expensive erf() function twice for every value.  The first and third passes are potentially parallelizable.  The second is highly serial and is probably slowed down by its highly random memory access patterns.  It also might be worthwhile to cache each double's bin number, depending on the CPU power to memory speed ratio.
This program lets you choose the number of bins to use.  Just add a second number to the command line.  I compiled it with gcc -O3, but my machine's so feeble I can't tell you any good performance numbers.
Edit: Poof!  My C program has magically transformed into a C++ program using std::sort!

You can use phi for a faster stdnormal_cdf.

How many bins should I put, approximately?

@Alexandru: I added a piecewise linear approximation to normcdf and only gained about 5% speed.

@static_rtti:  You don't have to put any.  By default, the code chooses the bins count so the average bin size is 10/11 of 128kb.  Too few bins and you don't get the benefit of partitioning.  Too many and the partition phase bogs down due to overflowing the cache.

10.6s! I tried playing a bit with the number of bins, and I got the best results with 5000 (slightly over the default value of 3356). I must say I was expected to see a much better performance for your solution... Maybe it's the fact that you're using qsort instead of the potentially faster std::sort of the C++ solutions?

Have a look at the radix sort implementation by Michael Herf (Radix Tricks). On my machine sorting was 5 times faster compared to the std::sort algorithm in my first answer. The name of the sorting function is RadixSort11.
int main(void)
{
    std::ifstream ifs("C:\\Temp\\gaussian.dat", std::ios::binary | std::ios::in);
    std::vector<float> v;
    v.reserve(50000000);
    double d;
    while (ifs.read(reinterpret_cast<char*>(&d), sizeof(double)))
        v.push_back(static_cast<float>(d));
    std::vector<float> vres(v.size(), 0.0);
    clock_t c0 = clock();
    RadixSort11(&v[0], &vres[0], v.size());
    std::cout << "Finished after: "
              << static_cast<double>(clock() - c0) / CLOCKS_PER_SEC << std::endl;
    return 0;
}

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Box Blur Challenge
The challenge is to implement a box blur in any language. The input will be as a 32-bit uncompressed bitmap, BGRA, preceded by one byte defining the radius of the blur and one byte defining the width of the image. The output should be in the same format as the input (minus the blur radius, of course).
Fastest code wins, and a 50-point bounty to the shortest entry, and another to the "most obfuscated" in my opinion.

Is this 1 dimensional or 2 dimensional?  if 2, how do we determine the image dimensions?

8-bits per channel? Blur is crappy linear averaging rather than taking gamma into account?

Challenge is way too vague as it stands - at the very least you need to specify the blur algorithm much more precisely

@Paul R: No particular algorithm in mind, as long as it blurs. Sorry, I should have specified that. @Keith Randall: By "bitmap" I meant "bitmap image, 2-dimensional". And yes, blur is "crappy averaging", @Peter Taylor.

@Keith Randall: The input will end in EOF in the language of your choice.

@minitech: if the image is 2-D, then we'll need to know the X and Y dimensions, not just the total size (X*Y) which we can get from the input size.

@Keith Randall: Fixed.

@Joey: for greyscale yes, but for RGBA there is more than one way of applying a box blur

Is the input a bitmap (1 bit per pixel), or a greyscale pixmap (with 8? 32? bits per pixel)?

@MtnViewMark: When I say "BGRA" and "32-bit", I mean "full-color with alpha channel" and "32 bits per pixel". Read the question next time.

@minitech: instead of being rude, you should put some effort into improving your question, since it is still way too vague and you haven't answered all the questions raised in the comments.

@Paul R: Please explain how the question is vague. (My point is not to be rude - people should really read questions thouroughly before posting comments.)

How the question is vague: How will you time the code? (Or build it, since you've specified any language.) Any OS/hardware restrictions? Do you actually mean BGRA (which would be somewhat unusual), or was that RGBA? Is the input little- or big-endian? What does the fastest code win, if not the 50-pt bounty?

@minitech: Caleb already beat me to it on most areas where the question is lacking in clarity and definition - I'll just add that you still haven't defined how the blur should be applied to the R, G, B and A channels - perhaps you should read the comments more carefully ?

@Paul R: To quote comment #4: "@Paul R: No particular algorithm in mind, as long as it blurs. Sorry, I should have specified that. @Keith Randall: By "bitmap" I meant "bitmap image, 2-dimensional". And yes, blur is "crappy averaging", @Peter Taylor. – minitech Jun 11 at 19:38 "

@Caleb: Yes, I mean BGRA, that is how the input bytes will be ordered. The fastest code gets marked as the answer, and it will be "fastest" as in "run on my computer". I will download & install the latest versions (or other, if specified) of any interpreters or compilers necessary to run/compile the program. If it's not available online, you should provide the compiler/interpreter in your answer.

Are edges handled by wrapping?

@PeterTaylor: Edges would ideally be handled by transparency, but black or white or wrapping all work too.

eTeX, too many chars
\immediate\openout5\jobname.out\let\x\newlinechar\x0\catcode0=12\let\y
\linepenalty\def\ndef#1{\expandafter\edef\csname\the#1\endcsname}\loop
\lccode0\x\lowercase{\ndef\x{^^@}}\ifnum\x<255\advance\x1\repeat\def\R
#1 #2 #3 {\Sr#1 \def\Y{#2}\x=0 \Rl#3????}\def\Rle#1\fi#2\Rc{\fi\Rl#1}
\def\Ry#1#2{\ndef{\y;\the\x.#1}{\number`#2 }}\def\Rc#1#2#3#4{\ifnum\Y=
\y\Rle#1#2#3#4\fi\advance\y1 \Ry1#1\Ry2#2\Ry3#3\Ry4#4\Rc}\def\Re{\fi\x
-1\immediate\write5{\Bl1;1.!}\end}\def\Rl#1{\edef\X{\the\x}\advance\x1
\y0 \ifx?#1\Re\fi\Rc#1}\def\Bl#1;{\ifnum\X<#1 \Ble\fi\Bc#1;}\def\Ble
#1!{\fi}\def\Bc#1;#2.{\ifnum\Y<#2 \Bce\fi\By#1;#2.1\By#1;#2.2\By
#1;#2.3\By#1;#2.4\expandafter\Bc\number#1\expandafter;\the\numexpr
#2+1.}\def\Bce#1#2#3;{#1\expandafter\Bl\the\numexpr#3+1;1.}\def\Sr#1
{\x-#1\def~{}\Sl#1 \edef~{\def\noexpand\By####1;####2.####3{\noexpand
\csname\noexpand\the\numexpr(\unexpanded\expandafter{~}+0)/(\unexpanded
{\Bt\X}#1 ####1 *\unexpanded{\Bt\Y}#1 ####2 )\endcsname}}~}\def\Sl#1
{\ifnum#1<\x\else\y-#1 \Sc#1 \advance\x1 \Sl#1 \fi}\def\Sc#1 {\ifnum
#1<\y\else\edef~{\unexpanded\expandafter{~}\unexpanded{\expandafter\t
\csname\the}\numexpr####2+\the\x;\noexpand\the\numexpr####1+\the\y
.####3\endcsname}\advance\y1 \Sc#1 \fi}\def\t#1{\ifx\relax#1\else+#1\fi
}\def\Bt#1#2 #3 {\expandafter\use\the\numexpr(\ifnum#2<#3 1+#2\else
#3\fi+\ifnum\numexpr#2+#3>#1 #1-#3\else#2\fi);}\def\use#1;{#1}
\expandafter\R\noexpand

Run as for instance
etex filename 10 3 55555555qqqq55555555qqqq5555qqqqqqqq5555qqqqqqqq

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Make a simple IRC bot.
Task

Make a simple IRC bot which does simple task like greeting user.

Detail
(You may see this RFC document.)
Write a complete program that perform these:

The program inputs nickname and realname from the user.

All nickname and realname which form is [A-Za-z0-9]+ should be allowed, but it isn't  important whether other forms are allowed. For example, if the user wants 1234qwerty as nickname and 5319009 as realname, the program should use those names, but it isn't important whether the program can use ^^☆Rick Astley☆^^ (which includes non-alphabet&decimals) as (nickname or realname) or not, even thought this name may be rejected by the server.

The program inputs serverAddress from the user.

The form of serverAddress is serverIP:serverPort, where serverIP is the IP address of the server and serverPort is the port num of the server.

The program connects to port serverPort at IRC server serverIP, and it should set it's nickname as nickname and realname as realname by sending USER realname 0 * :realname and NICK nickname.

The program inputs a single channelName which doesn't include #, and join to the channel #channelName.

Send a private message Hello, world! to the #channelName channel like this: PRIVMSG #channelName :Hello, world!

Then, the program does these:

A. If someone joins, greeting him by saying Hello, @! to #channelName, where @ is the nickname of him/her.
B. If the program is kicked, try re-joining.
C. If a person is kicked or banned, say POW HAHA.
D. If a person leaves (PART or QUIT), say Goodbye, @!, where @ is the nickname of him/her.
E. If someone says some text including Turing test or turing test, say I'm a human!.
F. If someone says some text including 6*9, quit by sending QUIT :42 to the server.
G. If someone says some text including all your base are belong to us(case-insensitive), quit by sending QUIT :What you say! to the server.
H. If someone says some text preceding with a space and above acts aren't performed by the text, say the sum of charcode of all (except the first space) chars in the text.
If an error is occurred
If an error (such as connection error or invalid nickname) is occurred, you can do what you want to do. For example, if the nickname is already in use, the program may stop, gets a new nickname from the user, or automatically change the nickname and retry.
Restriction
The program

Should not use any internal functions or libraries which is for IRC client/bot making. i.e. something like IRC client libraries
Should prompt what it'll input, before it gets an input. For example, before inputs serverAddress, it should prompt Write the address of server:, serverIP:serverPort > , serverAddress >, or some text that the user can recognize.
Should work on several IRC servers using UTF-8. (i.e. not for single server)
's messages send to server shouldn't be NOTICE. (RFC 2812, 3.3.2 Notice : The difference between NOTICE and PRIVMSG is that automatic replies MUST NEVER be sent in response to a NOTICE message.)
Should send a message where the request came from (#channelName or nickname).
Should PONG when PINGed, of course.

PS
Well, I just made this, because there was no puzzle about IRC bot.
Edit
After reading @dmckee's comment, I changed the rule for choosing winner.

Default bonus = +0

If somehow the program can connect to a channel (even if it's not #channelName), bonus = 1

If the program can perform step 1~5, bonus *= 2.

For each tasks A~H in step 6, if the program implements, bonus *= 1.2.

For each six restrictions, if the program follows that, bonus *= 1.05.

and score = int(codelength/bonus).
Program with lowest score is winner. If two answers have same score, then the answer with highest votes wins.
EDIT Edit : I think 1.2 and 1.05 in step 4 and 5 should be adjusted slightly bigger.. What do you think about this?
The winner is decided, however, I think there was too few entries.
Anyone still may submit the code.. :P

Frankly we have not been very accepting of "There is no metric, just vote.". There are things for you to work with here in terms of setting up an objective metric for this task. See Fake error message for a possible example.

@dmckee : Sorry, I didn't know that votes was quite subjective.. :(

@JiminP: Nice save.

irc bot sounds like a longer golf - I'm in though :)

clarification: "If the program is kicked, try re-joining." - does this mean from the server or the channel?

@lunixbochs From the channel.

no other comers? sadness. playing golf by myself is lonely :)

I think I can't decide winner now, because there isn't many codes... :(

@lunixbochs: as a person who uses IRC waaay too much, kicks are channel only.  If you're not on the server anymore, you got killed.

PHP - 121 points

1396 characters
Follows all steps (including extra) and restrictions

Score
<?php
$bonus = 1;                             //connects to channel
$bonus *= 2;                            //completes steps 1-5
for ($x=0;$x < 8;$x++) $bonus *= 1.2;   //all 8 extra steps
for ($y=0;$y < 6;$y++) $bonus *= 1.05;  //all 6 restrictions
echo (int)(1396 / $bonus);

> 121

Code
<?php function w($t){global$s;echo$t."\n";socket_write($s,$t."\r\n");}function m($c,$t){w("PRIVMSG $c :$t");}@$y=fgets;$z=STDIN;echo"nickname>";$n=trim($y($z));echo"realname>";$r=trim($y($z));echo"serverIP:Port>";$a=explode(":",trim($y($z)));$s=socket_create(2,1,6);socket_connect($s,$a[0],$a[1]);w("NICK $n");w("USER $n 0 * :$r");while(1){if((!($l=trim(socket_read($s,512,1))))||(!preg_match("/^(?:\:(\S+)[ ]+)?([^: ][^ ]*)(?:[ ]+([^:].*?))?(?:[ ]+:(.+))?$/",$l,$m)))continue;echo$l."\n";if(preg_match("/(.+)!(.+)@(.+)/",$m[1],$o))$m[1]=$o;if(!empty($m[3]))$p=explode(" ",$m[3]);else$p=array();if(isset($m[4])) $p[]=$m[4];$p=array($m[1],$m[2],$p);$b=$p[2][0];@$e=$p[0][1];@$f=$p[2][1];if($p[1]=="PING")w("PONG :".$b);elseif($p[1]=="376"){echo"#channel>";$c=trim($y($z));$c=(!empty($c))?$c:"#rintaun";w("JOIN $c");}elseif($p[1]=="JOIN")if($e!=$n)m($b,"Hello, $e!");else m($b,"Hello, world!");elseif($p[1]=="KICK")if($f!=$n)m($b,"POW HAHA");else w("JOIN ".$b);elseif(($p[1]=="PART")&&($e!=$n))m($b,"Goodbye, $e!");elseif(($p[1]=="QUIT")&&($e!=$n))m($c,"Goodbye, $e!");elseif($p[1]=="PRIVMSG")if(preg_match("/[Tt]uring test/",$f))m(($b==$n)?$e:$b,"I'm a human!");elseif(strstr($f,"6*9")){w("QUIT :42");break;}elseif(stristr($f,"all your base are belong to us")){w("QUIT :What you say!");break;}elseif($f[0]==" "){$q=str_split(substr($f,1));$u=0;foreach($q AS $d)$u+=ord($d);m(($b==$n)?$e:$b,$u);}}

Side Note
My answer actually differs from the instructions very slightly; the instructions say to use the realname as the first parameter to the USER command when registering, but I used the nickname instead. This is because the realname is allowed to have spaces, while this parameter is not. It's something of a moot point, though, since switching in the variable for realname would be the exact same character count.

Python - 125 points

304 chars
follows steps 1-5 (works for me on irc.freenode.net:6667. if you try there, remember IDENT might take a while, so give it 20+ seconds to connect and speak)
follows 4 restrictions (assuming the NOTICE restriction which is a bit unclear. the UTF restriction counts because python won't choke if it doesn't need to parse any UTF).

score math (python):
bonus = 1 # connects to channel
bonus *= 2 # steps 1-5
for i in xrange(4):
    bonus *= 1.05 # 4 restrictions

int(304.0/bonus)
> 125

code:
import socket
i=raw_input
u=i('user# ')
n=i('nick# ')
h,p=i('host:port# ').split(':',1)
p=int(p)
c='#'+i('chan# ')
z=0
while 1:
 try:
  def s(m): z.send(m+'\r\n')
  z.recv(9)
 except:
  z=socket.socket();z.connect((h,p));s('USER '+n+' 0 * :'+n);s('NICK '+n);s('JOIN '+c);s('PRIVMSG '+c+' :Hello, world!')

Perl, 66 points

666 characters
all substeps
half of the restrictions

Score
use 5.010;
$bonus = 1;              # connects to channel
$bonus *= 2;             # steps 1 to 5
$bonus *= 1.2 for 1..8;  # substeps A to H
$bonus *= 1.05 for 1..3; # restrictions 3, 4, 6
say int(666 / $bonus);
> 66

Code
(newlines for presentation only, not counted let alone acceptable)
use POE"Component::IRC";$_='sub _start{Zregister,all);Z"connect")}subX001{
Zjoin=>$c);ZY"Hello, world!")}*Xquit=*Xpart=*Xjoin=sub{$g=$_[STATE]=~/t/?G
oodbye:Hello;$_=$_[ARG0];/\w+/;ZY"$g, $&!")};subXkick{$_=$_[ARG2];/\w+/;Z$
&eq$n?"join":Y"POW HAHA")}subXpublic{$_=$_[ARG2];if(/turing test/i){ZY"I\'
m a human!")}elsif(/6\*9/){Zquit,42)}elsif(/all your base are belong to us
/i){Zquit,"What you say!")}elsif(/^ /){$t=-32;$t+=ord for/./g;ZY$t)}}chop(
($n,$r,$s,$c)=<>);$c="#$c";$i=POE::Component::IRC->spawn(nick,$n,ircname,$
r,server,$s)';s/Z/\$i->yield(/g;s/Y/privmsg,\$c,/g;S/X/ irc_/g;eval;POE::S
ession->create(package_states,[main,[grep*$_{CODE},%::]]);POE::Kernel->run

Side Note
The "ALL YOUR BASE ARE BELONG TO US" substep is most probably not worth its character count, but if I dropped it I wouldn't have the nice character count.  Bugger.

+1 for 66 points and 666 characters. Is that intentional or did you get that by accident?

@nyuszika7h I honestly don't remember.

Ruby, 65 points

574 Chars
Steps 1-5, 6(A-G) and restrictions on libs, UTF-8, no NOTICE & PONG
Takes commandline args in the form of server:6667 botnick botuser channel

Score
bonus = 0       # We've done nothing yet
bonus += 1      # Connects to channel
bonus *= 2      # Step 1 - 5

for n in 1..7
    bonus *= 1.2    # A - G
end

for n in 1..4
    bonus *= 1.05   # 4 constraints
end

puts (574/bonus).to_i
> 65

Code
require'socket'
a=ARGV
c="##{a[3]}"
h,p=a[0].split':'
s=TCPSocket.open(h,p)
m="PRIVMSG #{c} :"
s.puts"USER #{a[2]} 0 * :#{a[2]}\nNICK #{a[1]}\nJOIN #{c}\n#{m}Hello, world!"
while l=s.gets
case l
when/\:(.+?)!(.+)JOIN/
s.puts"#{m}Hello, #{$1}!"
when/KICK (.+?) (.+?) \:(.+)/
if $2==a[1]
s.puts"JOIN #{c}"
else
s.puts"#{m}POW HAHA"
end
when/\:(.+?)\!(.+)(PART|QUIT)/
s.puts"#{m}Goodbye, #{$1}!"
when/turing test/i
s.puts"#{m}I'm a human!"
when/6\*9/
s.puts"QUIT :42"
when/all your base are belong to us/i
s.puts"QUIT :What you say!"
when/PING \:(.+)/
puts"PONG #{$1}"
end
end

Ruby, 28 points

249 Chars
Same bot as my other solution, but slightly cheating maybe by applying compression.
http://coolfire.insomnia247.nl/golfbot-inflate.rb (As pasting this code probably won't work too well.)

Score
bonus = 0       # We've done nothing yet
bonus += 1      # Connects to channel
bonus *= 2      # Step 1 - 5

for n in 1..7
    bonus *= 1.2    # A - G
end

for n in 1..4
    bonus *= 1.05   # 4 constraints
end

puts (249/bonus).to_i
>  28

Code
require 'zlib'
b=<<'E'
x�]�_o�0���)��q�-KFBcdf��߃�P�S�B�����(N�Л��s~�� ��(�_$U6��5G�)�r�BB������J�{��� nNhlO�f)QQdJ�g��'�yP�!!����K�ɫ��[Ё�Ə{0�F]ѽ�m�2���GŐP��p`   ��I����E�+�* z�                                                                                )jrmKR�ˮ�%�
#��nQaJ�H��<�ZT���虦T3.�$D('�hw��a�/'��&�_ei�}o��1���M$����H��J�$��������V"���"��'��|A�`<��3L)Y��Z|� e�� ���m�é��ǚ�ڎu��J�����Vq~(ح�
E
eval Zlib::Inflate.new.inflate b

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/