title: Formatting numbers

text:
# Question
Title: Formatting numbers
Today computers use millions of bytes of memory, and decimal or hexadecimal notations get less and less informative.

Most human-readable formatting just display the number with the highest unit prefix. What about the rest?

I devised a notation that I use for myself, and since then, I'm often faced with re-implementing it in various languages. I thought I could share it, and that it make for a good case for code-golfing too.

**The rules**

* The program is given an integer as its input. It can be taken from a decimal or hexadecimal representation in a text file and then parsed, or from the "system innards".
* It should outputs a representation of that number somehow
* The limits are the system's limit. For practical purposes, assume 64 bits unsigned, so that the maximum number is 0xFFFFFFFFFFFFFFFF, and the unit suffixes will be k (kibi), M (mebi), G (gibi), T (tebi), P (pebi), E (exbi).
* zeros should be skipped. Example: 5G + 10k should be written 5G;10k, and not 5G;0M;10k;0
* It frequently happens, e.g. for memory ranges, that a number is just one less than a number with a nicer representation. For instance, 0xffff would be displayed as 63k;1023, though it is 64k minus 1. An exclamation mark denotes this, so that 0xffff should be displayed as !64k.

***Examples***

* 0 : 0
* 1023: !1k
* 1024: 1k
* 1025: 1k;1
* 4096: 4k
* 0x100000 : 1M
* 0xFFFFF : !1M
* 0x20000f : 2M;15
* 0x80003c00 : 2G;15k
* 0xFFFFFFFFFFFFFFFF : !16E
    + For 32 bits users, use 0xFFFFFFFF, which is, as everybody knows, !4G.
    + And 128 bits users can go up to... wait that's big! Anyway, feel free...

# Answer
> 2 votes
## Windows PowerShell, 207 ~~214~~ ~~226~~ ~~230~~ ~~247~~ ~~271~~ ~~319~~ ~~339~~ ~~360~~

```
$v=1PB*1KB
filter t{[Math]::Floor($n/$v)}$(,0*!($n=[Convert]::ToUInt64(($i="$input"),10+6*($i[1]-gt60)))
'EPTGMk'[0..5]|%{if($n%$v+1-eq$v){$v--
"!$(t)$_"
$n=0}if(t){"$(t)$_"
$n%=$v}$v/=1kb}
,$n*!!$n)-join';'

```

Test run:

```
./friendly2.ps1: PASS: '2047' -> '!2k'
./friendly2.ps1: PASS: '4096' -> '4k'
./friendly2.ps1: PASS: '0x80003c00' -> '2G;15k'
./friendly2.ps1: PASS: '2049' -> '2k;1'
./friendly2.ps1: PASS: '0x1' -> '1'
./friendly2.ps1: PASS: '4097' -> '4k;1'
./friendly2.ps1: PASS: '0xFFFFF' -> '!1M'
./friendly2.ps1: PASS: '0xa000' -> '40k'
./friendly2.ps1: PASS: '1022' -> '1022'
./friendly2.ps1: PASS: '0x100000' -> '1M'
./friendly2.ps1: PASS: '2048' -> '2k'
./friendly2.ps1: PASS: '0x20000f' -> '2M;15'
./friendly2.ps1: PASS: '0' -> '0'
./friendly2.ps1: PASS: '0x0' -> '0'
./friendly2.ps1: PASS: '1' -> '1'
./friendly2.ps1: PASS: '0x80aef3' -> '8M;43k;755'
./friendly2.ps1: PASS: '14636699861675008' -> '13P;1023M;27k'
./friendly2.ps1: PASS: '0xFFFFFFFFFFFFFFFe' -> '15E;1023P;1023T;1023G;1023M;1023k;1022'
./friendly2.ps1: PASS: '1023' -> '!1k'
./friendly2.ps1: PASS: '4095' -> '!4k'
./friendly2.ps1: PASS: '0xFFFFFFFFFFFFFFFF' -> '!16E'
./friendly2.ps1: PASS: '1025' -> '1k;1'
./friendly2.ps1: PASS: '1024' -> '1k'

```

# Answer
> 2 votes
## Haskell, ~~335~~ ~~333~~ 332

```
f 0=["0"]
f n=let l=mod n 1024;s=f$div(n+1)1024;in if l==1023 then "!":s else show l:s
g []=[]
g ((_,"0"):s)=g s
g ((_,"!"):s)=let(c,z):t=g s in(c,'!':z):t
g ((c,n):s)=(c,n):g s
h ('_',s)=';':s
h (c,s)=';':(s++[c])
r n=let x=tail$concatMap h$reverse$g$zip"_kMGTPE"$f$read n in if x=="" then "0" else x
main=do x<-getLine;return$r x

```

# Answer
> 1 votes
# JavaScript, ~~461~~ ~~452~~ ~~423~~ ~~269~~ ~~267~~ 247

```
function(n){for(k=1024,M=k*k,G=M*k,T=G*k,P=T*k,E=P*k,y='EPTGMk',r='';n;){for(i=-1,f=0;++i<6;)n==(v=eval(x=y[i]))-1?(r+='!1'+x,n=0,f=1):n>=v?(r+=Math.floor(n/v)+x+';',n%=v,f=1):0;f?0:(r+=n,n=0);}return r?r[l=r.length-1]==';'?r.substring(0,l):r:'0'}

```

# Answer
> 1 votes
## Ruby 1.9, 130 ~~167 181 191 203~~ characters

```
n=eval gets
u=1023
x=n&u<u ?n<1??0:"":(n+=1;?!)
s=[]
6.downto(0){|i|a=1024**i;n/a>0&&s<<[n/a,"EPTGMk"[6-i]]*""&&n=n%a}
puts x+s*?;

```

Passes all testcases from Joey's answer.

# Answer
> 1 votes
## Scala, 295 characters

```
object F extends App{type S=String;def f(l:BigInt,s:S,p:S,g:Int,q:Int):S={val j=l&1023;return if(l+g>0)f(l>>10,if((j+g<1024&&j+g-q>0)||l-q<0)";"+"!"*(g-q)+(j+g-q)+p(0)+s else s,p.tail,(j.toInt+g)/1024,0)else s};val x=readLine.split('x');print(f(BigInt(x.last,4+x.size*6),""," kMGTPE",1,1).tail)}

```

With some line breaks to aid readability (though not much):

```
object F extends App
{
    type S=String
    def f(l:BigInt,s:S,p:S,g:Int,q:Int):S=
    {
        val j=l&1023
        return if(l+g>0)f(l>>10,if((j+g<1024&&j+g-q>0)||l-q<0)";"+"!"*(g-q)+(j+g-q)+p(0)+s else s,p.tail,(j.toInt+g)/1024,0)else s
    }
    val x=readLine.split('x')
    print(f(BigInt(x.last,4+x.size*6),""," kMGTPE",1,1).tail)
}

```

# Answer
> 0 votes
## C#, 302 ~~323~~

Direct conversion of my PowerShell solution.

```
using System;class F{static void
Main(){var
s=Console.ReadLine();ulong
n=Convert.ToUInt64(s,s.IndexOf('x')==1?16:10),v=1L<<60,b;s=n<1?"0":"";foreach(var
x in "EPTGMk"){b=0;if(n%v==v-1){s+="!"+(n/--v)+x;n=0;}if(n/v>0){s=s+(n/v)+x;n%=v;b=1;}if(b>0&&n>0)s+=';';v/=1024;}if(n>0)s+=n;Console.Write(s);}}

```

# Answer
> 0 votes
## Python, 138 chars

```
N=input()
M=1023
s='0'[N:]
if N&M==M:N+=1;s='!'
print s+';'.join('%d'%b+k for b,k in zip([N>>60-10*i&M for i in range(7)],'EPTGMk ')if b)

```

---
Tags: code-golf, number
---


title: Find Possible Word Rectangles

text:
# Question
Title: Find Possible Word Rectangles
Johnny is trying to create crossword puzzles, but he is having difficulty making words match each other.

He has come up with several simple word rectangles: that is, groups of words that form a rectangle where all horizontal and vertical paths form a word.

```
//2x2
PA
AM

//2x3
GOB
ORE

//3x3
BAG
AGO
RED

//3x4
MACE
AGES
WEES

```

However, to make a good puzzle, he needs some word rectangles that are somewhat larger than 3x4. Instead of agonizing over the arrangement letters for hours, Johnny would prefer to have a program that does this for him--and in as few characters as possible, because long blocks of code are extremely intimidating for casual programmers like Johnny.

---

Given

* a text file dictionary where words are separated by newlines in alphabetical order,
* input specifying the number of rows and columns in the words rectangle (which can be provided however is most convenient in your programming language of choice)

generate at least one word rectangle. If it is not possible to generate a word rectangle with the given lexicon and dimensions, the program need not have defined behaviour. It is not necessary for the program to be able to generate rectangles which contain more than 64 letters or have dimensions that exceed 8 in either direction. The program should be able to complete in a reasonable amount of time, say, in thirty minutes or less.

---

EDIT: If you are doing an NxN rectangle, you are permitted to use a smaller dictionary file that only contains words that have a length of N letters.

# Answer
# Haskell, 586 characters

```
import Data.List
import qualified Data.Vector as V
import System
data P=P[String](V.Vector P)
e=P[]$V.replicate 128 e
(a:z)∈(P w v)|a>' '=z∈(V.!)v(fromEnum a);_∈(P w _)=w
p∫w=q p w where q(P u v)=P(w:u).V.accum q v.x;x(a:z)=[(fromEnum a,z)];x _=[]
l%n=foldl'(∫)e.filter((==n).length)$l
d§(p,q:r)=map(\w->p++w:r)$q∈d;_§(p,_)=[p]
(d¶e)i=filter(not.any(null.(∈e))).map transpose.(d§).splitAt i
s i n d e m|i==n=["":m]|1<3=(d¶e)i m>>=(e¶d)i>>=s(i+1)n d e
p[b,a,n]w=take n$s 0(a`max`b)(w%a)(w%b)$replicate b$replicate a ' '
main=do{g<-map read`fmap`getArgs;interact$unlines.concat.p g.lines}

```

Invoked by supplying 3 arguments: number of rows, number of columns, number of solutions; and the word list is accepted on `stdin`:

```
$> ghc -O3 2554-WordRect.hs
[1 of 1] Compiling Main             ( 2554-WordRect.hs, 2554-WordRect.o )
Linking 2554-WordRect ...

$> time ./2554-WordRect 7 7 1 < 2554-words.txt

zosters
overlet
seriema
trimmer
element
remends
startsy

real    0m22.381s
user    0m22.094s
sys     0m0.223s

```

As you can see 7×7 runs relatively fast. Still timing 8×8 and 7×6....

It would be 9 characters shorter to remove the number of solutions argument and just produce all solutions, but then it becomes impossible to time.

---

* Edit: (585 → 455) replaced custom data structure with a simple map of prefix string to possible replacements; oddly, this is a bit slower, perhaps because `Map String a` is slower than a hand-built tree of `Map Char a`...
* Edit: (455 → 586) *Bigger?!?!!* This version does more optimization of the search space, using both the techniques of my original solution, and of the python and awk solutions. Further, the custom data structure, based on `Vector` is much faster than using a simple `Map`. Why do this? Because I'm think a solution that is closer to the goal of 8x8 in under ½ hour is preferable to a shorter solution.
> 5 votes

# Answer
# Python, 232 chars

```
x,y=input()
H=[]
P={}
for w in open('words.txt'):
 l=len(w)-2
 if l==x:H+=[w]
 if l==y:
  for i in range(y+1):P[w[:i]]=1
def B(s):
 if(x+2)*y-len(s):[B(s+w)for w in H if all((s+w)[i::x+2]in P for i in range(x))]
 else:print s
B('')

```

It can only handle up to 6x6 in the 1/2 hour limit, though.
> 5 votes

# Answer
# Java (1065 bytes)

```
import java.util.*;public class W{public static void main(String[]a){new
W(Integer.parseInt(a[0]),Integer.parseInt(a[1]));}W(int w,int h){M
H=new M(),V=new M();String L;int i,j,l,m,n=w*h,p[]=new int[n];long
I,J,K,M,C=31;long[]G=new long[h],T=new long[w],W[]=new long[n][],X;try{Scanner
S=new Scanner(new java.io.File("words.txt"));while(0<1){L=S.nextLine();l=L.length();for(i=0;i>>l<1;i++){K=0;for(j=0;j<l;j++)K+=(i>>j&1)*(L.charAt(j)-96L)<<5*j;if(l==w)H.put(K,H.g(K)+1);if(l==h)V.put(K,V.g(K)+1);}}}catch(Exception
E){}while(n-->0){j=1;if(W[n]==null){M=1L<<62;for(i=w*h;i-->0;){m=i/w;l=i%w*5;if((G[m]>>l&C)<1){X=new
long[27];I=K=0;for(;K++<26;){J=H.g(G[m]+(K<<l))*V.g(T[i%w]+(K<<5*m));X[(int)K]=K-32*J;I+=J;}if(I<1)j=0;if(I<M){M=I;p[n]=i;W[n]=X;}}}}X=W[n];Arrays.sort(X);M=X[0]*j;X[0]=0;K=M&C;i=p[n]%w;j=p[n]/w;l=5*i;m=5*j;G[j]&=~(C<<l);G[j]+=K<<l;T[i]&=~(C<<m);T[i]+=K<<m;if(M>=0){W[n]=null;n+=2;}}for(long
A:G){L="";for(i=0;i<w;)L+=(char)(96+(C&A>>5*i++));System.out.println(L);}}class
M extends HashMap<Long,Long>{long g(Long s){return get(s)!=null?get(s):0;}}}

```

A long way from being the shortest, but I think it's the closest to meeting the timing constraints. I saved 14 bytes by assuming that the input file has been filtered to words of the right lengths; on my netbook, if you feed the whole `words.txt` then it spends the first minute preprocessing it, discarding most of what it produces, and then takes a mere 20 or so seconds to solve 7x7. On my desktop it does the whole thing in under 15 seconds, giving:

```
rascals
areolae
serrate
coroner
alanine
latents
seeress

```

I've let it run for over 50 hours without finding a solution to 8x7 or 8x8. 8-letter words seem to be a critical boundary for this problem - it just hovers around half-full without making much progress.

The approach used is full pivoting and a heuristic based on the number of possible horizontal completions times the number of possible vertical completions. E.g. if we have intermediate grid

```
*ean*
algae
*ar**
*ier*
*nee*

```

then we give the top-left corner a heuristic value of `count(aean*)count(aa***) + count(bean*)count(ba***) + ... + count(zean*)count(za***)`. Of all the cells we pick the one with the smallest heuristic value (i.e. hardest to satisfy), and then work though the letters in descending order of the amount they contributed to the heuristic value of that cell (i.e. starting with the most likely to succeed).
> 3 votes

# Answer
## F#

Backtracking solution, but I'm going to optimize the search space later.

```
open System

(*-NOTES
    W=7 H=3
    abcdefg<-- searching from wordsW
    abcdefg
    abcdefg
    ^
    partial filtering from wordsH
  *)

let prefix (s : char[]) (a : char[]) =
  a.[0..s.Length - 1] = s

let filterPrefix (s : char[]) =
  Array.filter (prefix s)

let transpose (s : char[][]) =
  [|
    for y = 0 to s.[0].Length - 1 do
      yield [|
        for x = 0 to s.Length - 1 do
          yield s.[x].[y]
      |]
  |]

[<EntryPoint>]
let main (args : String[]) =
  let filename, width, height = "C:\Users\AAA\Desktop\words.txt", 3, 3
  let lines = System.IO.File.ReadAllLines filename |> Array.map (fun x -> x.ToCharArray())
  let wordsW = Array.filter (Array.length >> ((=) width)) lines
  let wordsH = Array.filter (Array.length >> ((=) height)) lines

  let isValid (partial : char[][]) =
    if partial.Length = 0 then
      true
    else
      seq {
        for part in transpose partial do
          yield Seq.exists (prefix part) wordsH
      }
      |> Seq.reduce (&&)

  let slns = ref []
  let rec back (sub : char[][]) =
    if isValid sub then
      if sub.Length = height then
        slns := sub :: !slns
      else
        for word in wordsW do
          back <| Array.append sub [| word |]

  back [| |]
  printfn "%A" !slns
  0

```
> 2 votes

# Answer
## awk, 283

(may need to add 14 for parameter input flags)

Call with e.g. `awk -v x=2 -v y=2`...
Find first match and print it (283 chars):

```
{if(x==L=length)w[++n]=$0;if(y==L)for(;L>0;L--)W[substr($0,1,L)]++}END{for(i[Y=1]++;i[j=1]<=n;){b[Y]=w[i[Y]];while(j<=x){s="";for(k=1;k<=Y;k++)s=s substr(b[k],j,1);W[s]?0:j=x;j++}if(W[s])if(Y-y)i[++Y]=0;else{for(k=1;k<=Y;k++)print b[k];exit}i[Y]++;while(Y>1&&i[Y]>n)i[--Y]++}print N}

```

Find number of matches (245 chars, much slower):

```
{if(x==L=length)w[++n]=$0;if(y==L)for(;L>0;L--)W[substr($0,1,L)]++}END{for(i[Y=1]++;i[j=1]<=n;){b[Y]=w[i[Y]];while(j<=x){s="";for(k=1;k<=Y;k++)s=s substr(b[k],j,1);W[s]?0:j=x;j++}W[s]?Y-y?i[++Y]=0:N++:0;i[Y]++;while(Y>1&&i[Y]>n)i[--Y]++}print N}

```

For both programs (of course more so for the solutions count), the running time far exceeds 30 minutes for some values of x and y.

Just as a matter of interest, here is the word count for each word length:

```
 2     85
 3    908
 4   3686
 5   8258
 6  14374
 7  21727
 8  26447
 9  16658
10   9199
11   5296
12   3166
13   1960
14   1023
15    557
16    261
17    132
18     48
19     16
20      5
21      3

```
> 1 votes

---
Tags: code-golf, crossword
---


title: Pascal's Pyramid and Higher Dimensions

text:
# Question
Title: Pascal's Pyramid and Higher Dimensions
Pascal's Triangle is a familiar mathematical construct with many interesting properties. It is constructed by starting with a 1 on top, and generating the numbers in the next row from the sum of the number to the left and right in the row above (with implied 0's around the boundary):

```
    1
   1 1
  1 2 1
 1 3 3 1
1 4 6 4 1

```

The Triangle can be expanded into N dimensions by using the same rules. For example, the 3-D Pascal's Triangle would be a pyramid with each row being a cross-section of the pyramid, i.e., a triangle. Each triangle's entry is generated by adding the three numbers in the "row" above forming a triangle directly above that entry. An example of the first four "rows" is shown below:

```
1

 1
1 1

  1
 2 2
1 2 1

   1
  3 3
 3 6 3
1 3 3 1

```

The 4-D Pascal's Triangle would have a 3-D Pascal's Triangle for the cross section, and so forth.

Your task is to write a program that takes as input the dimension of the Pascal's Triange, and a row index to display. The row index starts at 0, so row 0 is `1` for every dimension. Your program should abide by the following specifications:

* Take input N = dimension number, R = row number (0-based). The input can be either command-line arguments or from stdin. N \>= 2, R \>= 0
* Output the R'th row in the N-dimensional Pascal's Triangle, by displaying each row or triangle on a separate line as in the above example.
* The format can vary as long as it is 'triangle-like' and obvious where each row and entry in the row belongs with respect to other rows.

Here is some sample output (I believe these are correct):

```
>./pascal 2 5
1 5 10 10 5 1
>./pascal 3 4
        1
      4   4
    6  12   6
  4  12  12   4
1   4   6   4   1
>./pascal 4 3
1

 3
3 3

  3
 6 6
3 6 3

   1
  3 3
 3 6 3
1 3 3 1
>./pascal 5 2
1

2

 2
2 2

1

 2
2 2

  1
 2 2
1 2 1

```

Bonus: What is the significance of each "row" in an N-dimensional Pascal's triangle? I can think of at least three interesting facts.

# Answer
## Ruby, 181

```
def P(a)a.min<0?0:a.max<1?1:(b=a*1;t=k=0;b.map{b[k]-=1;k+=1;t+=P(b)};t)end
def d(n,a)print n>1?((0..a[-1]).map{|i|d(n-1,a+[i])};"\n"):"#{P(a)} "end
d(ARGV[0].to_i,[ARGV[1].to_i])

```

The first approach using the following recursive formula

```
P(...,-1,...) = 0
P(0,0,...,0) = 1
P(a,b,c,....) = P(a-1,b,c,...) + P(a-1,b-1,c,...) + P(a-1,b-1,c-1,...) + ...

```

and then printing row `P(R,...)`.

The output is triangle-like - for `pascal 4 2` it looks like

```
1

3
3 3

3
6 6
3 6 3

1
3 3
3 6 3
1 3 3 1

```
> 5 votes

---
Tags: code-golf, ascii-art
---


title: Code 39 Barcode Scanner

text:
# Question
Title: Code 39 Barcode Scanner
Code 39, developed in 1974, is one of the most commonly used symbologies, or types, of barcodes, although it is the UPC/EAN system that is most often seen in retail sales. Code 39 barcodes can encode uppercase letters, numbers, and some symbols and are trivial to print from computer software using a special font. This led to their widespread commercial and industrial use (e.g. company ID badges, asset tracking, factory automation).

Create the shortest program or function to read a Code 39 barcode in any orientation from a 512x512 pixel grayscale image; the barcode might not be aligned horizontally or vertically.

* Your program must accept a standard image file format and produce the data encoded in the barcode as its standard output or return value (not including any start/stop character).
* No image contains more than one valid Code 39 barcode, and no barcode encodes a space character (ASCII 32).
* If no valid Code 39 barcode is shown in the image, the program must output a single question mark (`?`).

I have prepared a JavaScript reference implementation and test suite of images in PNG format, both with valid barcodes and without. The reference implementation, which fails only 3 of 46 test cases in most recent Web browsers, is intended to show one possible decoding algorithm, not to strictly conform to the above specification.

A valid submission passes at least 80% of these tests (37/46) and takes no more than one minute to do so for each image on a reasonably fast CPU (e.g. 2.6 GHz quad-core). My reference implementation passes 93% of tests and processes each image within 10 seconds (on my desktop PC running Google Chrome).

*(This question was proposed on Meta on May 28, 2011.)*

# Answer
## Python, 899 chars

```
import sys,random
raw_input()
X,Y=map(int,raw_input().split())
input()
I=[' x'[v<'~']for v in sys.stdin.read()]
M={196:' ',168:'$',148:'*',388:'.',52:'0',97:'2',49:'4',112:'6',292:'8',73:'B',25:'D',88:'F',268:'H',28:'J',67:'L',19:'N',82:'P',262:'R',22:'T',193:'V',145:'X',208:\
'Z',42:'%',138:'+',133:'-',162:'/',289:'1',352:'3',304:'5',37:'7',100:'9',265:'A',328:'C',280:'E',13:'G',76:'I',259:'K',322:'M',274:'O',7:'Q',70:'S',385:'U',448:'W'\
,400:'Y'}
N=500
for w in' '*30000:
 a,b,c,d=eval('random.random(),'*4);A=''.join(I[int((a+(c-a)*i/N)*X)+X*int((b+(d-b)*i/N)*Y)]for i in range(N)).lstrip();T=A.count(' x')+1;K=A.count('x')/T;L=A.count\
(' ')/T;s='';z=c=0
 while A:
  z*=2;y=A.find(' ')
  if y<0:y=len(A)
  z+=y>K;A=A[y:]
  z*=2;y=A.find('x')
  if y<0:y=len(A)
  z+=y>L;A=A[y:];c+=2
  if c>9:
   if z/2in M:s+=M[z/2];z=c=0
   else:break
 if s and'*'==s[0]and'*'==s[-1]and'*'!=s:print s[1:-1];break

```

This code takes a pnm format image as input, so I normally run it like:

```
pngtopnm s01.png | ./barcode.py

```

The code itself just picks lots of random scanlines and tries to match the black and white runs on that scanline to the code39 patterns. It is randomized so it can fail to find barcodes on occasion. (I get about a 20% false negative failure rate on the test images.) When it fails it takes about a minute to run, when it succeeds it often does so much quicker than that. I've never seen a false positive.
> 6 votes

---
Tags: code-golf, restricted-time
---


title: Expanding Expressions

text:
# Question
Title: Expanding Expressions
Given an expression of one-letter variables (`[a-z]`), operators (`*, +, &`) and parenthesis, expand it using the following axioms:

```
a * b != b * a
a * b * c = (a * b) * c = a * (b * c)
a + b = b + a
a + b + c = (a + b) + c = a + (b + c)
a * (b + c) = a * b + a * c
a & b = b & a
a & (b + c) = a & b + a & c

   | Comm | Assoc | Dist
 * | NO   | YES   |
 + | YES  | YES   | *
 & | YES  | YES   | + *

```

The user will input an expression, the syntax of the input expression is called "input form". It has the following grammar:

```
inputform ::= expr
var ::= [a-z] // one lowercase alphabet
expr ::= add
       | add & expr
add ::= mult
      | mult + add
mult ::= prim
       | prim * mult
       | prim mult
prim ::= var
       | ( expr )

```

Before That means, the order of operations is * + &, `a + b * c & d + e = (a + (b * c)) & (d + e)` Furthermore, the operator `*` can be omitted: a b (c + d) = ab(c + d) = a * b * (c + d)

Whitespace is stripped out before parsing.

Examples:

```
(a + b) * (c + d)
= (a + b)(c + d)
= (a + b)c + (a + b)d
= ac + bc + ad + bd
(a & b)(c & d)
= ac & ad & bc & bd
(a & b) + (c & d)
= a + c & a + d & b + c & b + d
((a & b) + c)(d + e)
= ((a & b) + c)d + ((a & b) + c)e (I'm choosing the reduction order that is shortest, but you don't need to)
= ((a & b)d + cd) + ((a & b)e + ce)
= ((ad & bd) + cd) + ((ae & be) + ce)
= (ad + cd & bd + cd) + (ae + ce & be + ce)
= ad + cd + ae + ce & ad + cd + be + ce & bd + cd + ae + ce & bd + cd + be + ce

```

Due to commutativity, order of some terms do not matter.

Your program will read an expression in input form, and expand the expression fully, and output the expanded expression in input form, with one space separating operators, and no spaces for multiplication. (`a + bc` instead of `a+b * c` or `a + b * c` or `a + b c`)

The fully expanded form can be written without any parens, for example, `a + b & a + c` is fully expanded, because it has no parens, and `a(b + c)` is not.

Here is an example interactive session (notice the whitespaces in input)

```
$ expand
> a(b + (c&d))
ab + ac & ab + ad
> x y * (wz)
xywz
> x-y+1
Syntax error
> c(a +b
Syntax error

```

# Answer
## Ruby, 698 648

```
O=[A='&',B='+',C='*']
M='Syntax error'
def u(x)c=x.shift;(c.nil?)?S[0]=~/[a-z]/?(S.shift):S[0]=='('?(S.shift;o=u(O*1);S[0]==')'?(S.shift;o):abort(M)):abort(M):(o=[c];loop{o<<u(x*1);break if S[0]!=c;S.shift};o)end
def d(e,a,b)a!=e[0]?e:(i=e.index{|x|x[0]==b};(i.nil?)?e:y([b]+e[i][1..-1].map{|k|z=e*1;z[i]=k;z}))end
def y(e)e.map!{|x|(x.is_a?Array)?(y(x)):(x)};o=[e[0]];e[1..-1].map{|x|x[0]==o[0]?(o+=x[1..-1]):(o<<x)};d(d(d(o.size==2?o[1]:o,C,B),C,A),B,A)end
def j(e)(e.is_a?Array)?(e.map!{|x|j(x)};o=e.shift;e*(o=='*'?'':" #{o} ")):e;end
S=gets.gsub(/ /,'').gsub(/([a-z)])(?=[a-z(])/,'\1*\2').split(//)
k=u(O*1)
puts S!=["\n"]?(M):(j y k)

```

This solution consists of a basic recursive descent parser (method `u`) which builds an syntax tree which is then flattened (`y`) and transformed (`d`) according to the distribution axioms. The output (`j`) simply joins the operands again.

Unfortunately I broke the functionality when I later tried to reduce the code by restructuring. Nevertheless I think there is still some chance for golfing.
> 3 votes

# Answer
### Scala - 1049 characters

This is not fully golfed, but it gets increasingly incomprehensible as I continue, and since the examples do not all appear to be correct, I want to leave the code at least slightly sane so it can be modified.

To be run as a script file. Empty input ends the program.

```
type S=String
case class R(rs: List[R], s: S)
def n(s: S)=List(" + "," & ").indexOf(s)
def pr(r:R,v:Int=2):S = if (r.rs.isEmpty) r.s else {
  val x = r.rs.map(pr(_,n(r.s))).mkString(r.s)
  if (n(r.s)>v) "("+x+")" else x
}
object P extends scala.util.parsing.combinator.RegexParsers {
  val v="[a-z]"r
  def e:Parser[R]=repsep(a,"&") ^^ {R(_," & ")}
  def a=repsep(m,"+") ^^ {R(_," + ")}
  def m=repsep(p,"\\*?"r) ^^ {R(_,"")}
  def p=v ^^ {R(Nil,_)} | "("~>e<~")" ^^ {r=>R(r::Nil,"")}
}
def up(r: R):R = r match {
  case R(x::Nil,_) => up(x)
  case _ => R(r.rs.map(up),r.s)
}
def asc(r: R) = R(r.rs.flatMap(i=>if (i.s==r.s) i.rs else List(i)),r.s)
def dis(r: R):R = {
  val q = asc(r.rs.find(i=>n(i.s)>n(r.s)).map(i => {
    val (a,b)=r.rs.span(_ ne i)
    R(i.rs.map(j=>R(a ::: j :: b.tail,r.s)),i.s)
  }).getOrElse(R(r.rs.map(dis),r.s)))
  if (q==r) r else dis(q)
}
def parse(s: S) = P.parseAll(P.e,s).map(x=>pr(dis(up(x)))).getOrElse("Syntax error")
Iterator.continually{print("> ");readLine}.
takeWhile(_.length>0).foreach(s=>println(parse(s)))

```
> 1 votes

# Answer
# Haskell, ~~472~~ 437 characters

```
import Text.Parsec
data E=V Char|O Int E E
r o=fmap$foldr1(O o)
c=char
s=sepBy1
e=r 0$s(r 1$s (r 2$s(fmap V lower<|>between(c '(')(c ')')e)(optional$c '*'))(c '+'))(c '&')
g(O a x(O b y z))|b<a=g$O b(g$O a x y)(g$O a x z)
g(O a (O b x y)z)|b<a=g$O b(g$O a x z)(g$O a y z)
g e=e
d(V c)=[c]
d(O o x y)=d x++[" & "," + ",""]!!o++d y
x=either(const "Syntax error")d.fmap g.parse e"".filter(`notElem`" \t")
main=interact$unlines.map x.lines

```
> 0 votes

---
Tags: code-golf
---


# Question
Title: 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!

# Answer
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.
> 14 votes

# Answer
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.
> 13 votes

# Answer
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());

  }
}

```
> 10 votes

# Answer
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();
        }
    }
}

```
> 9 votes

# Answer
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.
> 8 votes

# Answer
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.
> 8 votes

# Answer
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

```
> 6 votes

# Answer
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.
> 5 votes

# Answer
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.
> 5 votes

# Answer
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.
> 4 votes

# Answer
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.
> 4 votes

# Answer
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.
> 2 votes

# Answer
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();
    }
}

```
> 2 votes

# Answer
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.
> 2 votes

# Answer
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++).
> 2 votes

# Answer
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.
> 2 votes

# Answer
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;
}

```
> 2 votes

# Answer
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
> 2 votes

# Answer
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 votes

# Answer
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

```
> 1 votes

# Answer
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);
}

```
> 1 votes

# Answer
```
#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!
> 1 votes

# Answer
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;
}

```
> 1 votes

---
Tags: sorting, fastest-code, restricted-memory
---


title: Daily horoscope for zodiac signs

text:
# Question
Title: Daily horoscope for zodiac signs
Code a daily horoscope for the 12 zodiac signs, it needs to be self-sustained (no external ressources), quite randomized but fixed by zodiac sign for each day.

It needs to extract three "variables" from 1 to 3 :

1. love
2. work
3. money

Each one ranked from 1 to 3.

Ex. of variables that can be used (vars fixed for a day or more, can be extended) :

* Year
* Month
* Day
* Day of week
* Math.PI
* Zodiac signs order (from 1 to 12)

# Answer
> 2 votes
## Python, 106 chars

```
import time,zlib
for i in range(12):n=zlib.crc32(time.strftime('%%Y%%j%d'%i));print 1+n%3,1+n/3%3,1+n/9%3

```

Calculates the CRC of a string made up of the current year, day of year, and zodiac sign. Grabs some bits of the result for the horoscope.

# Answer
> 0 votes
From SO : "Add year+month\*11+day\*13+sign, feed this in as a seed to your local PRNG. After that you can either generate one "next random number" mod 27, extracting your 3 3s from there, or run it three times mod 3."

---
Tags: code-golf
---


title: Jumping and Running

text:
# Question
Title: Jumping and Running
Matthew likes solving puzzles. Whenever he manages to solve one he skips around happily. Recently he really *needs* to do this as a meteor shower has opened craters and holes in the ground in which he wouldn't like to fall.

You are given a part of landscape that Matthew wants to cross, hopefully arriving healthy at the end. The ground is given in meters, with each meter being either normal ground or a hole. When he is running he manages to cross one meter per step; the alternative is jumping which crosses four meters per step. Matthew starts at the far left on the first ground meter and wants to get to the last one (not beyond it, though – just imagine an endless hole beyond the last meter given in the landscape).

**Input**

Input is given as a single line on standard input, terminated by a line break. The line consists of either dashes (`-`) or underscores (`_`), representing a ground or hole meter, respectively. A sample input could be:

```
----__--___---

```

The given landscape is at least one and at most 30 meters long and always starts with ground.

**Output**

Output is given on standard output and represents a series of movement commands to Matthew, either *run* (`R`) or *jump* (`J`). As noted above, a *run* command causes Matthew to run one meter while jumping carries him forward exactly four meters. For the example given above the following movement is possible:

```
RRJRJRR

```

which looks approximately as follows:

If there is no safe path through the landscape, then a single exclamation mark (`!`) should be printed.

**Sample inputs**

```
--------
----__--___---
-_______
-_-_-_-_-_-
-

```

**Sample outputs**

```
JRRR
RRJRJRR
!
!

```

*(the last output is blank as no movement is necessary, but I guess, Markdown cannot parse this)*

**Note**

Only a single possible path is necessary, so the program output does not have to conform exactly to the sample outputs. As long as a solution is given if it exists and every movement command moves to ground and the last meter is reached eventually, the output is valid.

Additional output on standard error is ignored.

**Winning condition**

Shortest code wins, as is customary in golf. In case of a tie, the earlier solution wins.

**Test cases**

There are two tests scripts, containing identical test cases:

Invocation is in both cases: `<test script> <my program> [arguments]`, e.g. `./test ruby jumprun.rb` or `./test.ps1 ./jumprun.exe`.

**Another note**

This task was part of a golf contest held at my university during 2011-W24. The scores and languages of our contestants were as follows:

* 104 – Haskell
* 131 – Haskell
* 154 – C
* 170 – C
* 275 – VB.NET
* 286 – Common Lisp

Our own solutions were

* 92 – Ruby
* 124 – PowerShell

# Answer
## Perl, 53 characters

```
s/-...(?=(-|-...)*-$)/J/g;y/-/R/;/_/?$_="!":s/.$//

```

Run this with `perl -p jumpnrun.pl`. I've counted 3 characters for the `-p` option, which is the length difference between `perl jumpnrun.pl` and `perl -p jumpnrun.pl`.

I'm not that fluent in Perl, so I'm pretty sure this can be shortened further. This uses a regexp similar to Howard's solution.
> 7 votes

# Answer
### Ruby, 93 90 79 70 characters

I thought a regex solution would be quite fine and compact (let the matcher do the work). Unfortunately all the edge-cases and special treatments made this one such long - at least I didn't touch the 100 ;-).

```
puts gets.gsub(/-...(?=(-|-...)*-$)/,?J).tr(?-,?R)=~/^([JR]*)R$/?$1:?!

```

It passes all testcases of the provided script.

Saved several characters in comparison to the previous script (now a single call to `gsub` is sufficient).

*Edit 1:* Changed `puts z!=?-??!:''` to `z!=?-&&$><<?!` after the test script allowed no output for test case 1.

*Edit 2:* The previous version was

```
z=gets.chop
z.chars{z.sub!(/^(-|-...)((-|-...)*-)$/){$><<($1==?-??R:?J);$2}}
z!=?-&&$><<?!

```

My original idea was to replace the characters by using a look-behind and look-ahead strategy like this: The pattern was `^(?<=[RJ]*)(-|-...)(?=(-|-...)*-$)` and I then would replace '-' with 'R' and otherwise with 'J'. Unfortunately Ruby does not allow variable-length look-behind and another capturing group for the first part made the code even longer.

So then I did the iterative approach: if I can start with a step or jump `^(-|-...)` followed by series of other steps or jumps until the last platform `(-|-...)*-$` then I print the corresponding letter, remove the first one/four characters and start over again. On can even tune the RJ vs. JR priority by switching the choices inside the expression (currently it prefers RJ).

*Edit 3:* Splitting the single subtitution

```
puts (z=gets.chop.gsub(/(-...|-)(?=(-|-...)*-$)/){$1==?-??R:?J})=~/_/??!:z.chop

```

into two

```
puts (z=gets.chop.gsub(/-...(?=(-|-...)*-$)/,?J).tr(?-,?R))=~/_/??!:z.chop

```

gave another few chars. Finally I managed to get rid of this end-of-line issue but at a cost: the fail-detection costs some more characters.
> 3 votes

# Answer
# Perl - 71 60

```
$_=<>;y/-/R/;s/R...(?=(R(...)?)*R$)/J/g;print/_/?"!":s/.$//r

```

Now passes all of the testcases. :) Turns out I was removing the last character too soon... and half of my original regex was entirely redundant.

$\_=$ARGV\[0\];y/-/R/;s/(R...(?=R))(R\*(?=R))/J$2/g;chop;print/*/?"!":$*,$/

Yet another regex solution, passes the 5 testcases in the post.

Could be shortened by running as a one-liner with `-E` and `say` instead of `print`, but then perl tries to interpret the input as a switch... (`Unrecognized switch: -_-_-_-_-_-`)
> 1 votes

# Answer
## Haskell, 90 characters:

My first solution -- long, but works in linear time, using dynamic programming. :) **150 characters**:

```
x!y="!"
q '-'=max
q c=(!)
s i=r where r=zipWith3 q i(0&'R')$3&'J';n&c="":replicate n"!"++map(c#)r
c#"!"="!"
c#s=c:s
main=interact$reverse.last.s.init

```

The second solution -- much slower (exponential time), but much shorter: **90 characters**

```
s"-\n"=""
s('-':t)=max('R'#s t)$'J'#s(drop 3 t)
s _="!"
c#"!"="!"
c#s=c:s
main=interact s

```
> 1 votes

---
Tags: code-golf, path-finding
---


title: The strange life of a beehive

text:
# Question
Title: The strange life of a beehive
Researchers recently discovered an interesting bee colony that lives in an infinite field of honeycomb:

Each cell can house a bee or not. In fact, the lives of those creatures appear to be a bit ... chaotic. It could be calculated that a colony always starts with the following pattern:

<sub>(Bee drawn by Emmanuel Boutet on Wikimedia Commons. This honeycomb-and-bees image thus is released under CC-By-SA. *grumbles*)</sub>

After that the bee's life cycles are divided into so-called generations. Each generation old bees die and new ones hatch and it primarily depends on their cell's neighbors:

* If a bee has less than two neighbors it dies due to loneliness.
* If a bee has more than three neighbors it dies due to overcrowding.
* If a cell has two, three or four living bees in neighboring cells a new bee hatches there in the next generation.

Dying bees don't die until the end of a generation so they still affect surrounding cells that might hatch bees in the next generation.

Now that we know how such a colony works, we can simulate it through any number of generations.

**Input**

Input is a single number **N**, given on standard input, terminated by a line break. 0 ≤ **N** ≤ 150. This is the number of generations to simulate.

**Output**

Output is a single number, on standard output and optionally followed by a single line break, which represents the number of living bees after **N** generations.

Additional output on standard error is ignored.

**Sample inputs**

```
0
5
42
100

```

**Sample outputs**

```
6
44
1029
5296

```

**Winning condition**

Shortest code wins, as is customary in golf. In case of a tie, the earlier solution wins.

**Test cases**

There are two tests scripts, containing identical test cases:

Invocation is in both cases: `<test script> <my program> [arguments]`, e.g. `./test ruby beehive.rb` or `./test.ps1 ./beehive.exe`.

I know there are only 22 tests instead of 151 (mainly because solutions are often quite slow). Please refrain from embedding the exact test cases instead of solving the task. These scripts are a convenience for you to test whether a change still causes the program to behave correctly; not that you can adapt your code to the specific test cases.

**Another note**

This task was part of a golf contest held at my university during 2011-W24. The scores and languages of our contestants were as follows:

* 336 – C
* 363 – C
* 387 – C
* 389 – Haskell
* 455 – C

Our own solution was

* 230 – Ruby

# Answer
> 9 votes
### Ruby, 181 163 153 146 characters

```
h=[0]*4e4
[0,-200,201,202,2,3].map{|i|h[i]=1}
gets.to_i.times{h=h.map{[g=1,200,201].map{|x|g+=h[x]+h[-x]};g>5-h.rotate![-1]||g<3?0:1}}
p h.count 1

```

This implementation follows a standard approach using an array `h` (dimensions `200` x `200` flattened), where each element is either `0` (no bee) or `1` (bee included). The array `[0,-200,201,202,2,3]` describes the initial positions of bees (relative to any initial cell).

Input and output as specified above, passes all test cases defined.

*Edit 1:* changed back to a wrapping solution instead of the "additional space"-version (which was shorter in an intermediate version but now is several characters longer).

*Edit 2:* removed variable `b` completely.

*Edit 3:* warning: this edit made the program horribly slow. I thus reduced the dimensions to 200 each which still is sufficient for up to 150 iterations. Instead of indexing the array by a variable we constantly rotate the array forward. Not really good design, but now we are considerably below the 150.

# Answer
> 7 votes
## Python, 152 characters

```
P=[0,2,3,1j,1+1j,1-1j]
for i in' '*input():Q=[p+d for d in(1,-1,1j,-1j,1j-1,1-1j)for p in P];P=set(p for p in Q if 1<Q.count(p)<5-(p in P))
print len(P)

```

This solution keeps track of bee locations with a set of complex numbers. It is pretty slow because the inner loop is quadratic in the number of bees. I've tested up to 50 and it works.

# Answer
> 5 votes
## Python, 171 169 158 chars

```
w=300
s=w*w
x=[0]*297
h=[1,1,0]+x+[1,0,1,1]+x+[1]+x*s
exec('h=[1<sum(h[(i+x)%s]for x in[-1,-w-1,-w,1,w+1,w])<5-h[i]for i in range(s)];'*input())
print sum(h)

```

I model the world as a 300\*300=900000 1D array (`h` is actually bigger, but the end is not used), where a bee is a 1 and empty is 0. Size of 300 is fine because at *most* the growth would be of 2 in each dimension for each generation, and there is no more than 150 generations.

Here is a slightly ungolfed and commented version:

```
w=300 # width and height of the world
s=w*w
# create initial grid
l=[1,1]+[0]*298+[1,0,1,1]+[0]*297+[1]+[0]*s

for k in range(input()):
  h=l[:]

  for i in range(s):

    # for each point, compute the number of neighbors
    n=sum(map(lambda x:h[(i+x)%s],[-1,-w-1,-w,1,w+1,w]))

    # if that number verifies the conditions, put 1 here, if not put 0
    l[i]=1<n<5-h[i]

print sum(l)

```

---
Tags: code-golf, hexagonal-grid, simulation, cellular-automata
---