title: Filter a large file quickly

text:
# Question
Title: Filter a large file quickly
The challenge is to filter a large file quickly.

* Input: Each line has three space separated positive integers.
* Output: All input lines `A` `B`, `T`that satisfy either of the following criterion.

1. There exists another input line `C`, `D`, `U` where `D = A` and `0 <= T - U < 100`.
    2. There exists another input line `C`, `D`, `U` where `B = C` and `0 <= U - T < 100`.

To make a test file use the following python script which will also be used for testing. It will make a 1.3G file. You can of course reduce nolines for testing.

```
import random
nolines = 50000000 # 50 million
for i in xrange(nolines):
    print random.randint(0,nolines-1), random.randint(0,nolines-1), random.randint(0,nolines-1)

```

**Rules.** The fastest code when tested on an input file I make using the above script on my computer wins. The deadline is one week from the time of the first correct entry.

**My Machine** The timings will be run on my machine. This is a standard 8GB RAM ubuntu install on an AMD FX-8350 Eight-Core Processor. This also means I need to be able to run your code.

**Some relevant timing information**

Timings updated to run the following before each test.

```
sync && sudo bash -c 'echo  3 > /proc/sys/vm/drop_caches'

time wc test.file

real    0m26.835s
user    0m18.363s
sys     0m0.495s

time sort -n largefile.file  > /dev/null

real    1m32.344s
user    2m9.530s
sys     0m6.543s

```

**Status of entries**

I run the following line before each test.

```
sync && sudo bash -c 'echo  3 > /proc/sys/vm/drop_caches'

```
* **Perl** (Waiting for bug fix.)
* **Scala** 1 minutes 37 seconds by @James\_pic. (Using scala -J-Xmx6g Filterer largefile.file output.txt)
* **Java**. 1 minute 23 seconds by @Geobits. (Using java -Xmx6g Filter\_26643)
* **C**. 2 minutes 21 seconds by @ScottLeadley.
* **C**. 28 seconds by @James\_pic.
* **Python+pandas**. Maybe there is a simple "groupby" solution?
* **C**. 28 seconds by @KeithRandall.

**The winners are Keith Randall and James\_pic.**

I couldn't tell their running times apart and both are almost as fast as wc!

# Answer
## C, ~ 7 4.1 seconds

Radix sort on T, then walk through the array looking for matches.

It is fast because it is cache friendly. The radix sort reasonably so, and the final walk very. I have to check each row against about 100 others, but they are all adjacent in the cache.

Added: I no longer have to check each row against a scan of 100 other rows. A small table of counts of low order bits of b's in the window suffices to eliminate most instances of this scan.

Now about 1/2 time parsing, 1/3 time sorting, 1/6 time doing the actual matching.

```
#include <stdio.h>
#include <stdlib.h>
#include <sys/mman.h>
#include <fcntl.h>

// B = # of bits per radix pass
// R = # of radix passes
#define B 9
#define R 3
#define M ((1<<B)-1)
#define MAXN 50000000

int count[R][1<<B];

typedef struct {
  int a,b,t,print;
} entry;

entry A[MAXN];
entry C[MAXN];

// Sized to fit well in L1 cache
unsigned char bcount[16384];

int main(int argc, char *argv[]) {
  FILE *f = fopen(argv[1], "r");
  fseek(f, 0, SEEK_END);
  int size = ftell(f);
  fclose(f);

  int fd = open(argv[1], O_RDONLY);
  const char *p = (const char*)mmap(0, size, PROT_READ, MAP_SHARED, fd, 0);
  const char *endp = p + size;

  // parse, insert into array
  int n = 0;
  while(p < endp) {

    // parse line
    int a = 0;
    while(*p != ' ') {
      a *= 10;
      a += *p - '0';
      p++;
    }
    p++;
    int b = 0;
    while(*p != ' ') {
      b *= 10;
      b += *p - '0';
      p++;
    }
    p++;
    int t = 0;
    while(*p != '\n') {
      t *= 10;
      t += *p - '0';
      p++;
    }
    p++;

    // insert it
    if(n == MAXN) {
      printf("too many elements\n");
      exit(1);
    }
    A[n].a = a;
    A[n].b = b;
    A[n].t = t;
    n++;

    // compute counts for radix sort
    count[0][t&M]++;
    count[1][(t>>B)&M]++;
    count[2][t>>2*B]++;
  }

  // accumulate count entries
  for(int r = 0; r < R; r++) {
    for(int i = 0; i < M; i++) {
      count[r][i+1]+=count[r][i];
    }
  }

  // radix sort, 3 rounds
  for(int i = n-1; i >= 0; i--) {
    C[--count[0][A[i].t&M]] = A[i];
  }
  for(int i = n-1; i >= 0; i--) {
    A[--count[1][(C[i].t>>B)&M]] = C[i];
  }
  for(int i = n-1; i >= 0; i--) {
    C[--count[2][A[i].t>>2*B]] = A[i];
  }

  // Walk through array (now sorted by T) and find matches.
  // We maintain a window of T values that might match.
  // To facilitate finding matches within that window, bcount
  // keeps track of a count of how many b's in that window
  // have the given low 14 bits.
  int j = 0;
  for(int i = 0; i < n; i++) {
    int a = C[i].a;
    int t = C[i].t;
    while(C[j].t <= t - 100) {
      int x = C[j].b & 16383;
      if(bcount[x] != 255) bcount[x]--;
      j++;
    }
    if(bcount[a & 16383] > 0) {
      // somewhere in the window is a b that matches the
      // low 14 bits of a.  Find out if there is a full match.
      for(int k = j; k < i; k++) {
        if(a == C[k].b)
          C[k].print = C[i].print = 1;
      }
    }
    int x = C[i].b & 16383;
    if(bcount[x] != 255) bcount[x]++;
  }
  for(int i = 0; i < n; i++) {
    if(C[i].print)
      printf("%d %d %d\n", C[i].a, C[i].b, C[i].t);
  }
}

```
> 10 votes

# Answer
# Scala 2.10 - 0:41

The problem is basically:

```
select * from data x, data x where x.a = y.b and 0 <= x.t - y.t and x.t - y.t < 100

```

Most RDBMSs would notice that the join from `x.a` to `y.b` has the highest specificity, and plan this as a hash join.

So that's what we'll do. We create a hashtable of the data on `a`, hash join it with the same table on `b`, and filter on the difference in `t`.

```
import scala.io.Source
import scala.reflect.ClassTag
import java.io._

object Filterer {
  def roundUpToNextPowerOfTwo(x: Int) = {
    // blatantly stolen from http://bits.stephan-brumme.com/roundUpToNextPowerOfTwo.html
    var y = x - 1
    y |= y >> 1
    y |= y >> 2
    y |= y >> 4
    y |= y >> 8
    y |= y >> 16
    y + 1
  }

  // We hash join the array with itself, a to b, and emit both rows if t is within 100. 50m records should fit into 8GB OK.
  def main(args: Array[String]): Unit = {
    val input = Source.fromFile(args(0), "ASCII").getLines()
    val output = new BufferedWriter(new OutputStreamWriter(new FileOutputStream(args(1)), "US-ASCII"))
    try {
      val data1: Array[Row] = input.map{line =>
        Row(line)
      }.toArray

      /*
       * In theory, data1 and data2 could be created in parallel, but OpenHashMultiMap needs
       * to know its size at creation time, to sidestep the need for rehashing. I could just
       * hard-code 50 million (the size of the data in the challenge), but that seems dishonest.
       */
      val data2 = new OpenHashMultiMap[Int, Row](roundUpToNextPowerOfTwo(data1.size) * 2, -1)
      for (r <- data1) data2.insert(r.a, r) // data2 is hashed by a

      for (row1 <- data1.par) {
        val Row(a, b, t) = row1
        for (Row(c, d, u) <- data2.get(b) if (0 <= u - t) && (u - t < 100)) {
          // The conditions are symmetric, so if row1 matches, so does row2
          output.write(s"$a $b $t\n$c $d $u\n")
        }
      }
    } finally {
      output.close()
    }
  }
}

object Row {
  def apply(data: String): Row = {
    val l = data.length
    var i = 0
    var a = 0
    var b = 0
    var c = 0
    while (data.charAt(i) != ' ') {
      a = a * 10 + (data.charAt(i) - '0')
      i += 1
    }
    i += 1
    while (data.charAt(i) != ' ') {
      b = b * 10 + (data.charAt(i) - '0')
      i += 1
    }
    i += 1
    while (i < l) {
      c = c * 10 + (data.charAt(i) - '0')
      i += 1
    }
    Row(a, b, c)
  }
}

final case class Row(a: Int, b: Int, t: Int)

/*
 * None of the standard Java or Scala collections are particularly efficient as large MultiMaps,
 * so we write our own. We use open hashing with quadratic probing.
 */
class OpenHashMultiMap[@specialized(Int) K: ClassTag, V: ClassTag](capacity: Int, default: K) {
  require((capacity & (capacity - 1)) == 0) // Power of 2 capacity
  private val keys = Array.fill(capacity)(default)
  private val values = new Array[V](capacity)
  private val mask = capacity - 1

  private def hash(k: K) = {
    // Hash mingling - Int has a particularly poor hash
    k.hashCode * 428916315
  }

  def insert(k: K, v: V) = {
    var found = false
    var loc = hash(k) & mask
    var inc = 0
    while (inc <= capacity && !found) {
      loc = (loc + inc) & mask
      inc += 1
      found = keys(loc) == default
    }
    keys(loc) = k
    values(loc) = v
  }

  def get(key: K) = new Traversable[V] {
    override def foreach[U](f: V => U) = {
      var break = false
      var loc = hash(key) & mask
      var inc = 0
      while (inc <= capacity && !break) {
        loc = (loc + inc) & mask
        inc += 1
        val k = keys(loc)
        if (key == k) f(values(loc))
        else if (k == default) break = true
      }
    }
  }
}

```

Compile with:

```
scalac Filterer.scala

```

And run with:

```
scala -J-server -J-XX:+AggressiveOpts -J-Xms6g -J-Xmx6g Filterer input_file.dat output_file.txt

```

On my machine, this runs in 2 minutes 27.

Still, it might be interesting to try the approach from @Lembik's answer, but in a faster language. It corresponds to something like a merge join on `t`. On paper, it should be slower, but it has better cache locality, which may push it ahead.

## Update

I've managed to shave a large chunk of time off with a surprisingly small change - a better hash mingler. The hash map is very sensitive to hash clumping, so this change brings it down to 1:45 on my machine.

Most of that time is spent reading the data into an array.

I'm curious as to why my data reading code is so much slower than @Geobits. It takes my code 70 seconds to read the data in - longer than @Geobits entire program, once the `Thread.start` bug is fixed. I'm tempted to steal @Geobits approach to reading data, but I'm not sure how the Stack Exchange gods would feel about that.

## Update 2

I've made further improvements, this time to the data reader. Using pattern matching and monad operations inside the loop was hurting performance, so I've simplified it. I think `scala.io.Source` is the next bottleneck to tackle.

It's down to 1:26 on my machine now.

## Update 3

Got rid of `probe` from OpenHashMultiMap. The code is now more java-ish, and runs in 1:15.

## Update 4

I'm now using an FSM to parse input. Run time is down to 0:41
> 7 votes

# Answer
# C - 12 seconds

I decided to port my Scala answer to C, to see how much more performance I could get.

It's more or less the same approach (build an open hash table on `a`), except that I skip the step where I build the initial array, and iterate from the hash table directly (for some reason I could never get this approach to perform in Scala - I suspect JVM inlining was to blame).

I didn't bother with threads, as it's a pain to do portably.

The code is:

```
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>

// Should be 37% occupied with 50m entries
#define TABLE_SIZE 0x8000000
#define MASK (TABLE_SIZE - 1)
#define BUFFER_SIZE 16384
#define END_OF_FILE (-1)
#define DEFAULT_VALUE (-1)

typedef struct Row {
  int32_t a;
  int32_t b;
  int32_t t;
} Row;

int32_t hash(int32_t a) {
  return a * 428916315;
}

void insert(Row * table, Row row) {
  long loc = hash(row.a) & MASK; // Entries are hashed on a
  long inc = 0;
  while (inc <= TABLE_SIZE) {
    loc = (loc + inc) & MASK;
    inc++;
    if (table[loc].a == DEFAULT_VALUE) {
      table[loc] = row;
      break;
    }
  }
}

int readChar(FILE * input, char * buffer, int * pos, int * limit) {
  if (*limit < *pos) {
    return buffer[(*limit)++];
  } else {
    *limit = 0;
    *pos = fread(buffer, sizeof(char), BUFFER_SIZE, input);
    if (*limit < *pos) {
      return buffer[(*limit)++];
    } else return END_OF_FILE;
  }
}

void readAll(char * fileName, Row * table) {
  char* buffer = (char*) malloc(sizeof(char) * BUFFER_SIZE);
  int limit = 0;
  int pos = 0;

  FILE * input = fopen(fileName, "rb");

  int lastRead;
  Row currentRow;
  uint32_t * currentElement = &(currentRow.a);

  // As with the Scala version, we read rows with an FSM. We can
  // roll up some of the code using the `currentElement` pointer
  while (1) {
    switch(lastRead = readChar(input, buffer, &pos, &limit)) {
      case END_OF_FILE:
        fclose(input);
        return;
      case ' ':
        if (currentElement == &(currentRow.a)) currentElement = &(currentRow.b);
        else currentElement = &(currentRow.t);
        break;
      case '\n':
        insert(table, currentRow);
        currentRow.a = 0;
        currentRow.b = 0;
        currentRow.t = 0;
        currentElement = &(currentRow.a);
        break;
      default:
        *currentElement = *currentElement * 10 + (lastRead - '0');
        break;
    }
  }
  //printf("Read %d", lastRead);
}

int main() {
  Row* table = (Row*) malloc(sizeof(Row) * TABLE_SIZE);
  memset(table, 255, sizeof(Row) * TABLE_SIZE);

  readAll("test.file", table);

  // We'll iterate through our hash table inline - passing a callback
  // is trickier in C than in Scala, so we just don't bother
  for (size_t i = 0; i < TABLE_SIZE; i++) {
    Row * this = table + i;
    if (this->a != DEFAULT_VALUE) {
      // Lookup entries `that`, where `that.a == this.b`
      long loc = hash(this->b) & MASK;
      long inc = 0;
      while (inc <= TABLE_SIZE) {
        loc = (loc + inc) & MASK;
        inc++;
        Row * that = table + loc;
        if ((this->b == that->a) && (0 <= that->t - this->t) && (that->t - this->t < 100)) {
          // Conditions are symmetric, so we output both rows
          printf("%d %d %d\n", this->a, this->b, this->t);
          printf("%d %d %d\n", that->a, that->b, that->t);
        }
        else if (that->b == DEFAULT_VALUE) break;
      }
    }
  }

  free(table);
  return 0;
}

```

Compile with:

```
gcc -std=c99 -O3 -m64 filter.c

```

And run with:

```
./a.out

```

The location of the test file is hard-coded as "test.file".

Once again, reading the data takes up most of the time (just under 9 seconds). Matching takes the rest of the time.

Again, it'll be interesting to see how this fares against Scott Leadley's answer, which uses the same language but a different strategy. Scott is joining on T, which in principle means he'll have more to join against, but again, joining on T gives better cache locality.
> 7 votes

# Answer
# Java : 1m54s

**(On my i7)**

Since every match will be within 100 `t` of its mate, I decided to bucket the inputs by `t`. There is a bucket for each 100, so to check a number, it only has to check against +/-1 buckets.

On average, each bucket only contains 100 entries, so it doesn't take long to scan a few buckets for each. Well over half the time is spent reading and bucketing, matching takes only 40 seconds or so.

*Note:* Depending on your JVM setup, you may need to increase the heap size. This also assumes file name of `test.file`. Simply change it on line 24 if that's not the case.

```
import java.io.BufferedReader;
import java.io.FileReader;
import java.util.ArrayList;
import java.util.StringTokenizer;

public class Filter_26643 {

    final static int numThreads = 8;
    final static int numInputs = 50000000;
    final static int bucketSize = 100;
    final static int numBuckets = numInputs/bucketSize;
    ArrayList<ArrayList<int[]>> buckets;

    public static void main(String[] args) {
        new Filter_26643().run();
    }

    void run(){
        try{
            buckets = new ArrayList<ArrayList<int[]>>(numBuckets);
            for(int i=0;i<numBuckets;i++)
                buckets.add(new ArrayList<int[]>(bucketSize*2));

            BufferedReader reader = new BufferedReader(new FileReader("test.file"));
            int c=0,e[];
            while(c++<numInputs){
                StringTokenizer tokenizer = new StringTokenizer(reader.readLine());
                e = new int[] {
                                Integer.parseInt(tokenizer.nextToken()),
                                Integer.parseInt(tokenizer.nextToken()),
                                Integer.parseInt(tokenizer.nextToken())
                                };
                buckets.get(e[2]/100).add(e);
            }
            reader.close();

            MatchThread[] threads = new MatchThread[numThreads];
            for(int i=0;i<numThreads;i++){
                threads[i] = new MatchThread(i);
                threads[i].start();
            }
            for(int i=0;i<numThreads;i++)
                threads[i].join();

        } catch(Exception e){
            e.printStackTrace();
        }
    }

    class MatchThread extends Thread{
        int index;

        public MatchThread(int index){
            this.index = index;
        }

        @Override
        public void run() {
            for(int i=index;i<numBuckets;i+=numThreads){
                int max = i+2 >= numBuckets ? numBuckets : i+2;
                int min = i-1 < 0 ? i : i-1;
                for(int[] entry : buckets.get(i)){
                    outer:
                    for(int j=min;j<max;j++){
                        ArrayList<int[]> bucket = buckets.get(j);
                        for(int[] other : bucket){
                            if(((entry[0]==other[1] && entry[2]-other[2]<100 && entry[2]>=other[2]) ||
                                (entry[1]==other[0] && other[2]-entry[2]<100 && other[2]>=entry[2]))
                                && entry != other){
                                 System.out.println(entry[0] + " " + entry[1] + " " + entry[2]);
                                 break outer;
                            }
                        }

                    }
                }
            }
        }
    }
}

```
> 6 votes

# Answer
# perl, 17m46s on a core i7 w/ 8GB mem

First, we use sort -n -k3 to get the most important field in order, taking advantage of the built-in parallelism on modern versions of sort(1). Then, since perl is greatly hampered by the fact that a simple scalar takes on the order of 80 bytes each (50 million * 3 * 80 is too much - at least 12GB), we slurp the output in to a 50 million * 12 byte array (12 bytes per line, each line contains 3 integers that can be represented as a 32 bit integer). Then we fire off 8 threads covering (roughly) 1/8th of the data each (+ some overlap).

```
#!perl

use strict;
use warnings;

# find lines s.t. $lines[$M]->{a} == $lines[$N]->{b} and
#                 0 <= $lines[$M]->{t} - $lines[$N]->{t} < 100
# OR              $lines[$M]->{b} == $lines[$N]->{a} and
#                 0 <= $lines[$N]->{t} - $lines[$M]->{t} < 100

my $infile = shift;
open(my $fh, "sort -n -k3 $infile |") || die "open sort pipe: $@";

my @lines;
my $bytes_per_int = 4;
my $bytes_per_line = $bytes_per_int * 3;
my $nlines = 50_000_000;
my $buf = "\0" x ($nlines * $bytes_per_line);
my $ln = 0;
my $nprocs = 8;
my $last_group_start = 0;
my $this_group_start;
my $group = $nlines / $nprocs;
my @pids;
while(<$fh>) {
  my ($A, $B, $T) = split/\s+/;
  substr($buf, $ln * $bytes_per_line, $bytes_per_line, pack "L3", ($A, $B, $T));
  if( defined $this_group_start ) {
    if( $T - $last_group_start >= $group + 100 ) {
      if(my $pid = fork()) {
        push @pids, $pid;
        $last_group_start = $this_group_start;
        undef $this_group_start;
      } else {
#warn "checking $last_group_start - $ln...\n";
        for(my $l=$last_group_start; $l<=$ln; ++$l) {
          my $lpos = $l * $bytes_per_line;
          my ($A, $B, $T) = unpack "L3", substr($buf, $lpos, $bytes_per_line);
          my ($lA, $lB);
          my $lT = $T;
          for(my $lb=$l; $lb>=$last_group_start && $T - $lT <= 100; $lb--, $lpos -= $bytes_per_line) {
            ($lA, $lB, $lT) = unpack "L3", substr($buf, $lpos, $bytes_per_line);
            if($A == $lB || $B == $lA) {
              #print "($last_group_start) $A $B $T matches $lA $lB $lT\n";
              print "$lA $lB $lT\n$A $B $T\n";
            }
          }
        }
        exit;
      }
    }
  } elsif( !defined $this_group_start && $T - $last_group_start >= $group ) {
    $this_group_start = $ln;
  }
  $ln++;
}

waitpid $_, 0 for @pids;

```

Unsorted output:

```
8455767 30937130 50130
20468509 8455767 50175
47249523 17051933 111141
17051933 34508661 111215
39504040 36752393 196668
42758015 39504040 196685
25072294 28422439 329284
35458609 25072294 329375
45340163 42711710 6480186
39315845 45340163 6480248
1435779 49643646 12704996
38229692 1435779 12705039
18487099 24556657 6665821
24556657 28498505 6665884
6330540 35363455 18877328
22500774 6330540 18877347
10236123 22026399 598647
39941282 10236123 598717
45756517 24831687 6726642
34578158 45756517 6726670
29385533 7181838 621179
7181838 29036551 621189
40647929 11895227 25075557
11895227 1900895 25075652
17921258 42642822 18935923
40140275 17921258 18935949
44573044 38139831 12899467
38139831 1321655 12899468
11223983 1788656 12920946
1788656 21905607 12921040
1357565 8148234 801402
8148234 46556089 801498
30929735 303373 19105532
31258424 30929735 19105543
34899776 9929507 6990057
9929507 49221343 6990078
49779853 43951357 25306335
41120244 49779853 25306424
6177313 41551055 25343755
24462722 6177313 25343804
16392217 32915797 31472388
32915797 19696674 31472479
6834305 36264354 25440771
44983650 6834305 25440800
26559923 47360227 19356637
47360227 49749757 19356700
33018256 36233269 37654651
36233269 5459333 37654671
6932997 23123567 25502355
23123567 7882426 25502356
5878434 43421728 25510707
43421728 40827189 25510765
38695636 33504665 1099515
13504170 38695636 1099605
32832720 40188845 37689854
8335398 32832720 37689927
35858995 41917651 1130028
41917651 28797444 1130096
47102665 6796460 43806189
6796460 6113288 43806229
21248273 5422675 43819677
48011830 21248273 43819728
32187324 39177373 25624030
39177373 42539402 25624102
41722647 14351373 25626925
14351373 45070518 25627013
22298566 25860163 37862683
2273777 22298566 37862692
10617763 32776583 7561272
35581425 10617763 7561279
18526541 18709244 31960780
18709244 32777622 31960867
36976439 24222624 31973215
24222624 9534777 31973262
25751007 11612449 38066826
43652333 25751007 38066923
8303520 2615666 7633297
2615666 29961938 7633357
22317573 31811902 31982722
14298221 22317573 31982819
43089781 7653813 44154683
8732322 43089781 44154769
24227311 43800700 13711475
40906680 24227311 13711539
48061947 30109196 7660402
43993467 48061947 7660488
29580639 5292950 38140285
5292950 21293538 38140356
17646232 47737600 32058831
47737600 42934248 32058836
13262640 23462343 1617194
23462343 1901587 1617259
5150775 7046596 44270140
7046596 22819218 44270181
17749796 34924638 32171251
8386063 17749796 32171346
30095973 12202864 38257881
12202864 42679593 38257912
10353022 40646034 26158412
40646034 36237182 26158412
8416485 16245525 32223010
16245525 32420032 32223081
20420340 1371966 7893319
1371966 2031617 7893335
2864137 20279212 26199008
29145409 2864137 26199080
29141766 19729396 44433106
44115780 29141766 44433141
6513924 34515379 32283579
12686666 6513924 32283636
20116056 49736865 44464394
49736865 47918939 44464416
38212450 3465543 32302772
3465543 39217131 32302873
12019664 37367876 44485630
3639658 12019664 44485639
18053021 1279896 7973955
2220749 18053021 7974031
19701732 12984505 1857435
24625926 19701732 1857528
9876789 34881917 26285125
27687743 9876789 26285134
5696632 6064263 44534580
34888313 5696632 44534629
14865531 46418593 38457138
5929897 14865531 38457191
44378135 4051962 38485208
4051962 10804515 38485308
11865822 21793388 14142622
7760360 11865822 14142649
32333570 24478420 44702533
24478420 23749609 44702588
29098286 25015092 44723985
32171647 29098286 44723985
20522503 20522503 2127735
20522503 20522503 2127735
22597975 20938239 8260902
20938239 48618802 8260905
8310032 34659671 2153994
34659671 25406149 2154075
49085033 5708432 26644257
5708432 32265692 26644305
18751513 18226037 32726402
18226037 33885794 32726424
45877488 23211339 20566948
23211339 26209405 20567002
48554034 25770643 38853402
9683274 48554034 38853467
9770420 14556349 2309265
27255587 9770420 2309324
32926392 16744099 44954824
24840989 32926392 44954840
29066838 49434549 26755357
49434549 12635292 26755407
21927714 32352409 20626921
32352409 15895076 20626932
7422009 23559357 14550898
32743911 7422009 14550982
38816601 5850890 26851026
5850890 32996623 26851107
42148171 47021378 26872907
47021378 32628418 26872908
9850929 10501741 32998960
10501741 24899993 32999043
27491904 4393602 33033499
4393602 17712085 33033570
37978226 42226216 39114479
42226216 2511412 39114525
42859989 49908919 45241083
48131208 42859989 45241088
39753103 30674979 14807321
30674979 45637890 14807371
30154199 11988643 2641926
11988643 11241926 2641976
7191871 13518594 45370275
13518594 45354921 45370344
54745 19711137 8871851
24814115 54745 8871937
38770495 34574748 2756244
41962321 38770495 2756337
26229406 39306415 21057327
10735951 26229406 21057347
46704290 11506122 39359422
18181795 46704290 39359481
38796645 28410469 45452212
28410469 13478996 45452222
412456 27727741 39466147
27727741 19639136 39466226
24470627 13030982 21266756
13030982 21713410 21266825
6058593 23139172 27435254
19236012 6058593 27435303
14457750 39190113 39701131
30253141 14457750 39701227
26898421 39016446 45812750
40952330 26898421 45812829
18647206 27663400 45817956
27663400 21728474 45817989
5559358 41319001 33664547
41319001 37210583 33664636
29066692 30653068 39759813
30653068 38963132 39759856
12086617 49971187 3232640
49971187 32302154 3232649
12008399 13656671 3239395
43088998 12008399 3239439
10061612 38594475 39804389
38594475 6327106 39804405
16703379 21150436 39851057
21150436 34093320 39851136
1035486 4199407 3314170
26974438 1035486 3314196
21869320 14532221 33851404
15208937 21869320 33851473
38840190 4742355 3402401
4742355 46055869 3402462
34432016 8734566 39966972
27614117 34432016 39967002
9988172 49209666 46063078
49209666 29374155 46063087
3208946 47030309 21722002
47030309 39809983 21722030
10928661 46423741 3496854
46423741 29486710 3496862
42464855 22978174 46154827
22978174 3814497 46154901
47090840 16768393 46169667
39523858 47090840 46169714
28186104 11618234 34024001
11618234 33711158 34024019
45471813 37332848 3585557
37332848 4607526 3585600
14885742 38990612 15863749
38990612 3710491 15863779
42391514 33643913 22005928
33643913 32254640 22006022
4299590 19482026 34202327
19482026 35838894 34202406
24298776 16276160 3858885
16276160 3198758 3858958
29322567 12536696 40433239
12536696 26083938 40433317
16080151 9648322 22221443
9648322 43846385 22221458
999302 19218350 10078183
10296062 999302 10078189
40544377 34492433 34463953
19908418 40544377 34463993
10765321 45143043 34542584
39154522 10765321 34542646
48642526 31097951 4104790
2940654 48642526 4104887
26972730 47422139 46846889
39228577 26972730 46846901
13788696 11503551 34728076
11503551 9151627 34728130
8676030 30463644 10406398
15204754 8676030 10406405
42984277 41087708 34805119
48741576 42984277 34805143
29634598 2151247 22699609
12264074 29634598 22699614
47525963 48470003 16667878
48470003 4566846 16667953
9725907 43325112 4498307
26465445 9725907 4498368
306967 11708860 10633595
11708860 31017081 10633669
39420965 46595640 41089015
46595640 41260374 41089048
29232745 39705052 16754836
4739295 29232745 16754840
35246405 42811088 41273637
48986699 35246405 41273719
2398239 36985098 35181790
36985098 7460784 35181841
18955749 23678549 35221035
47264406 18955749 35221129
18105816 26003002 17044057
26003002 17467477 17044087
14430126 46039962 47492180
46039962 29118827 47492275
30329324 40926812 41425850
43304610 30329324 41425912
34966996 36567528 17095113
3967517 34966996 17095144
42829171 42530474 23209891
25923738 42829171 23209967
28187681 26297990 35474412
48986691 28187681 35474475
5707126 41598794 17298139
40466899 5707126 17298188
28838696 30725820 5142797
30725820 35360418 5142798
44642019 42570370 17339657
42570370 19022469 17339727
42193681 8389736 17386517
48906013 42193681 17386586
42303185 30337820 41795129
30337820 42473956 41795170
30935782 8441903 17515229
41549758 30935782 17515275
41239019 10011768 23619001
10011768 25386353 23619062
494288 13341166 29815779
49113152 494288 29815876
7106674 26227442 29833029
47459682 7106674 29833047
17246497 35389391 17628365
35389391 34005133 17628371
23347674 48243185 17792799
48243185 22907892 17792836
21852744 1662414 36088704
8040124 21852744 36088775
32384657 27122374 36100767
24980361 32384657 36100782
31016207 26300043 42222489
26300043 36869529 42222544
17178756 44315094 42223989
44315094 11222466 42224042
34139317 39164101 36197907
39164101 27563542 36197947
31638631 22215137 17999735
22215137 10771707 17999769
30257199 32883043 24127009
32883043 179099 24127047
47774058 17451960 30283073
44583527 47774058 30283162
13816647 12695130 24145102
12695130 42284941 24145188
42749234 20004242 5893793
20004242 38129713 5893819
22210359 22178109 18109989
22178109 112961 18110049
42509645 28599506 42508465
28599506 3722411 42508513
34412629 22547405 48610262
22547405 16664124 48610296
2330283 32267749 24256113
35915758 2330283 24256157
44560231 49353986 12101694
6471293 44560231 12101780
23289721 8186827 30407293
10624448 23289721 30407389
12329357 35765163 30560085
4511908 12329357 30560158
31332240 39704929 12269193
39704929 47770487 12269249
22286152 22082044 36734758
22082044 25076919 36734833
47381309 9459604 36735886
9459604 31071680 36735890
43832763 45342283 30707519
45342283 26992816 30707602
2883029 18642608 42989696
14697025 2883029 42989793
15149987 40746227 24700535
40746227 34776566 24700549
2387554 49015265 43057085
49015265 21103141 43057139
23057202 13308993 30982514
34596334 23057202 30982553
44598498 31714790 43285828
18170064 44598498 43285841
38273701 11976319 31179763
15344094 38273701 31179764
3651338 27427037 37188945
12876654 3651338 37189007
10081580 3418061 37221143
3418061 38353019 37221143
172544 18699860 37295343
824744 172544 37295372
13914 8890169 37303853
8890169 14008003 37303898
18716557 29456130 49605004
29456130 16390535 49605083
15398102 22446674 43711290
22446674 38760679 43711383

```

I'm sure this would be an order of magnitude faster in C, but I probably won't take the time to do so.
> 2 votes

# Answer
# C# - 30 Seconds

A different approach from most if I read right - I don't use any hash based structures.

~~I tend to get no results, not sure if this a statistical anomaly, or error in my reasoning.~~ Fixed, comparison for binary sort was flawed.

```
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.IO;
using System.Linq;
using System.Threading;
using System.Threading.Tasks;

namespace FilterFile
{
    class Program
    {
        const int COUNT = 50000000;

        static string inputFile = "data" + COUNT + ".txt";
        static string outputFile = "results.txt";

        static void Main(string[] args)
        {
            Console.WriteLine("Prepping Test");
            if (args.Length > 0) inputFile = args[0];
            if (args.Length > 1) outputFile = args[1];

            if (!File.Exists(inputFile))
            {
                Console.WriteLine(inputFile);

                File.WriteAllLines(inputFile,
                                     GenerateData(COUNT)
                                     .Select(r => string.Format("{0} {1} {2}", r.A, r.B, r.C)));
            }

            File.Delete("results.txt");

            Console.WriteLine("Starting Test \n\n");

            using (Timer.Create("Total Time"))
            {
                Row[] sortedA, sortedB;
                //http://codegolf.stackexchange.com/questions/26643/filter-a-large-file-quickly
                using (Timer.Create("Reading Data"))
                    FillData(out sortedA, out sortedB);

                using (Timer.Create("Parallel Sort A"))
                    ParallelSort.QuicksortParallel(sortedA);
                using (Timer.Create("Parallel Sort B"))
                    ParallelSort.QuicksortParallel(sortedB, (x, y) => x.B - y.B);

                object rLock = new object();
                List<Row> results = new List<Row>();

                var comparison = Comparer<Row>.Create((B, A) => B.B - A.A);
                using (Timer.Create("Compute Results"))
                    Parallel.ForEach(sortedA, row =>
                    //foreach (var row in sortedA)
                    {
                        var i = Array.BinarySearch(sortedB, row, comparison);
                        if (i < 0) return;

                        Row other;
                        bool solved = false;
                        for (var tempI = i; tempI < sortedB.Length && row.A == (other = sortedB[tempI]).B; tempI++)
                        {
                            var diff = row.C - other.C;
                            if (diff >= 0 && diff < 100)
                            {
                                lock (rLock) results.Add(row);
                                return;
                            }
                        }

                        for (var tempI = i - 1; tempI >= 0 && row.A == (other = sortedB[tempI]).B; tempI--)
                        {
                            var diff = row.C - other.C;
                            if (diff >= 0 && diff < 100)
                            {
                                lock (rLock) results.Add(row);
                                return;
                            }
                        }
                    });

                using (Timer.Create("Save Results"))
                {
                    File.WriteAllLines(outputFile, results.Select(r => r.ToString()));
                }
            }
        }

        private static void FillData(out Row[] sortedA, out Row[] sortedB)
        {
            var tempA = new Row[COUNT];
            var tempB = tempA;//new Row[COUNT];

            const int PARTITION_SIZE = 1 << 22;

            ReadAndSort(tempA, tempB, PARTITION_SIZE);

            sortedA = tempA;
            sortedB = new Row[COUNT];
            Array.Copy(sortedA, sortedB, COUNT);
            /*using (Timer.Create("MergeA"))
            {
                int destIndex = 0;
                int[][] partitions = Enumerable.Range(0, COUNT / PARTITION_SIZE + 1)
                    .Select(i => new[] { i * PARTITION_SIZE, Math.Min(i * PARTITION_SIZE + PARTITION_SIZE, COUNT) - 1 })
                    .ToArray();

                for (int i = 0; i < COUNT; i++)
                {
                    foreach (var partition in partitions)
                    {
                        while (partition[0] <= partition[1] && tempA[partition[0]].A == i)
                        {
                            sortedA[destIndex++] = tempA[partition[0]++];
                        }
                    }
                }
            }*/

            /*//Verify Paritioning Works
            var results = new List<Tuple<Row, int>> { Tuple.Create(tempA[0], 0) };
            for (int i = 1; i < tempA.Length; i++)
            {
                var r = tempA[i];
                if (r.A < tempA[i-1].A)
                    results.Add(Tuple.Create(r, i % PARTITION_SIZE));
            }
            results.ForEach(t => Console.WriteLine(t.Item1 + " " + t.Item2));*/
        }

        private static void ReadAndSort(Row[] tempA, Row[] tempB, int PARTITION_SIZE)
        {
            List<Task> tasks = new List<Task>();

            using (var stream = File.OpenRead(inputFile))
            {
                int b;
                int tempMember = 0;
                int memberIndex = 0;
                int elementIndex = 0;

                using (Timer.Create("Read From Disk"))
                    while ((b = stream.ReadByte()) >= 0)
                    {
                        switch (b)
                        {
                            case (byte)'\r':
                            case (byte)' ':
                                switch (memberIndex)
                                {
                                    case 0: tempA[elementIndex].A = tempMember; memberIndex = 1; break;
                                    case 1: tempA[elementIndex].B = tempMember; memberIndex = 2; break;
                                    case 2: tempA[elementIndex].C = tempMember; memberIndex = 0; break;
                                }
                                tempMember = 0;
                                break;
                            case (byte)'\n':
                                /*if (elementIndex % PARTITION_SIZE == 0 && elementIndex > 0)
                                {
                                    var copiedIndex = elementIndex;
                                    tasks.Add(Task.Run(() =>
                                    {
                                        var startIndex = copiedIndex - PARTITION_SIZE;
                                        Array.Copy(tempA, startIndex, tempB, startIndex, PARTITION_SIZE);
                                        ParallelSort.QuicksortSequentialInPlace(tempA, startIndex, copiedIndex - 1);
                                        ParallelSort.QuicksortSequentialInPlace(tempB, startIndex, copiedIndex - 1, (x, y) => x.B - y.B);
                                    }));
                                }*/
                                elementIndex++;
                                break;
                            default:
                                tempMember = tempMember * 10 + b - '0';
                                break;
                        }
                    }

                /* tasks.Add(Task.Run(() =>
                 {
                     elementIndex--;  //forget about the last \n
                     var startIndex = (elementIndex / PARTITION_SIZE) * PARTITION_SIZE;
                     Array.Copy(tempA, startIndex, tempB, startIndex, elementIndex - startIndex + 1);
                     ParallelSort.QuicksortParallelInPlace(tempA, startIndex, elementIndex);
                     ParallelSort.QuicksortSequentialInPlace(tempB, startIndex, elementIndex, (x, y) => x.B - y.B);
                 }));

                 using (Timer.Create("WaitForSortingToFinish"))
                     Task.WaitAll(tasks.ToArray());*/
            }
        }

        static Random rand = new Random();

        public struct Row : IComparable<Row>
        {
            public int A;
            public int B;
            public int C;
            public static Row RandomRow(int count)
            {
                return new Row { A = rand.Next(count), B = rand.Next(count), C = rand.Next(count) };
            }

            public int CompareTo(Row other)
            {
                return A - other.A;
            }

            public override string ToString()
            {
                return string.Format("{0} {1} {2}", A, B, C);
            }
        }

        public static Row[] GenerateData(int count)
        {
            var data = new Row[count];
            for (int i = 0; i < count; i++)
                data[i] = Row.RandomRow(count);
            return data;
        }

        public static Row[] GenerateSplitData(int count)
        {
            var data = new Row[count];
            for (int i = 0; i < count; i++)
                data[i] = Row.RandomRow(count);
            return data;
        }

        public class Timer : IDisposable
        {
            string message;
            Stopwatch sw;
            public static Timer Create(string message)
            {
                Console.WriteLine("Started: " + message);
                var t = new Timer();
                t.message = message;
                t.sw = Stopwatch.StartNew();
                return t;
            }
            public void Dispose()
            {
                Console.WriteLine("Finished: " + message + " in " + sw.ElapsedMilliseconds + "ms");
            }
        }

        // <summary>
        /// Parallel quicksort algorithm.
        /// </summary>
        public class ParallelSort
        {
            const int SEQUENTIAL_THRESHOLD = 4096;
            #region Public Static Methods

            /// <summary>
            /// Sequential quicksort.
            /// </summary>
            /// <typeparam name="T"></typeparam>
            /// <param name="arr"></param>
            public static void QuicksortSequential<T>(T[] arr) where T : IComparable<T>
            {
                QuicksortSequentialInPlace(arr, 0, arr.Length - 1);
            }

            /// <summary>
            /// Parallel quicksort
            /// </summary>
            /// <typeparam name="T"></typeparam>
            /// <param name="arr"></param>
            public static void QuicksortParallel<T>(T[] arr) where T : IComparable<T>
            {
                QuicksortParallelInPlace(arr, 0, arr.Length - 1);
            }

            #endregion

            #region Private Static Methods

            public static void QuicksortSequentialInPlace<T>(T[] arr, int left, int right)
                where T : IComparable<T>
            {
                if (right > left)
                {
                    int pivot = Partition(arr, left, right);
                    QuicksortSequentialInPlace(arr, left, pivot - 1);
                    QuicksortSequentialInPlace(arr, pivot + 1, right);
                }
            }

            public static void QuicksortParallelInPlace<T>(T[] arr, int left, int right)
                where T : IComparable<T>
            {
                if (right > left)
                {
                    if (right - left < SEQUENTIAL_THRESHOLD)
                        QuicksortSequentialInPlace(arr, left, right);
                    else
                    {
                        int pivot = Partition(arr, left, right);
                        Parallel.Invoke(() => QuicksortParallelInPlace(arr, left, pivot - 1),
                                        () => QuicksortParallelInPlace(arr, pivot + 1, right));
                    }
                }
            }

            private static void Swap<T>(T[] arr, int i, int j)
            {
                T tmp = arr[i];
                arr[i] = arr[j];
                arr[j] = tmp;
            }

            private static int Partition<T>(T[] arr, int low, int high)
                where T : IComparable<T>
            {
                // Simple partitioning implementation
                int pivotPos = (high + low) / 2;
                T pivot = arr[pivotPos];
                Swap(arr, low, pivotPos);

                int left = low;
                for (int i = low + 1; i <= high; i++)
                {
                    if (arr[i].CompareTo(pivot) < 0)
                    {
                        left++;
                        Swap(arr, i, left);
                    }
                }

                Swap(arr, low, left);
                return left;
            }

            #endregion

            #region Public Static Methods

            /// <summary>
            /// Sequential quicksort.
            /// </summary>
            /// <typeparam name="T"></typeparam>
            /// <param name="arr"></param>
            public static void QuicksortSequential<T>(T[] arr, Func<T, T, int> comparer)
            {
                QuicksortSequentialInPlace(arr, 0, arr.Length - 1, comparer);
            }

            /// <summary>
            /// Parallel quicksort
            /// </summary>
            /// <typeparam name="T"></typeparam>
            /// <param name="arr"></param>
            public static void QuicksortParallel<T>(T[] arr, Func<T, T, int> comparer)
            {
                QuicksortParallelInPlace(arr, 0, arr.Length - 1, comparer);
            }

            #endregion

            #region Private Static Methods

            public static void QuicksortSequentialInPlace<T>(T[] arr, int left, int right, Func<T, T, int> comparer)
            {
                if (right > left)
                {
                    int pivot = Partition(arr, left, right, comparer);
                    QuicksortSequentialInPlace(arr, left, pivot - 1, comparer);
                    QuicksortSequentialInPlace(arr, pivot + 1, right, comparer);
                }
            }

            public static void QuicksortParallelInPlace<T>(T[] arr, int left, int right, Func<T, T, int> comparer)
            {
                if (right > left)
                {
                    if (right - left < SEQUENTIAL_THRESHOLD)
                    {
                        QuicksortSequentialInPlace(arr, left, right, comparer);
                    }
                    else
                    {
                        int pivot = Partition(arr, left, right, comparer);
                        Parallel.Invoke(() => QuicksortParallelInPlace(arr, left, pivot - 1, comparer),
                                        () => QuicksortParallelInPlace(arr, pivot + 1, right, comparer));
                    }
                }
            }

            private static int Partition<T>(T[] arr, int low, int high, Func<T, T, int> comparer)
            {
                // Simple partitioning implementation
                int pivotPos = (high + low) / 2;
                T pivot = arr[pivotPos];
                Swap(arr, low, pivotPos);

                int left = low;
                for (int i = low + 1; i <= high; i++)
                {
                    if (comparer(arr[i], pivot) < 0)
                    {
                        left++;
                        Swap(arr, i, left);
                    }
                }

                Swap(arr, low, left);
                return left;
            }
            #endregion
        }
    }
}

```
> 2 votes

# Answer
### C

Brutal, brute force, in-your-face ugly C. On a do-over I would choose just about any other compiled language.

```
/*
Filter a file based on these rules:

Input:
    - each item is an ordered list of three integers ( A B T )
    - each line represents an item
    - each line is formated as <number> <w> <number> <w> <number>
    - <w> is whitespace (a single blank in the challenge)
    - <number> is an integer in the range 0..49_999_999
    - the first number on a line is A, second B, third T

Output a given item ( A B T ) if:
    1 - there exists an item ( C D U ) such that 0 <= T-U < 100 and D == A
    OR
    2 - there exists an item ( C D U ) such that 0 <= U-T < 100 and B == C

CLARIFICATION:
An item should be output only once, even if there is more than one match.

We're sorting on T, we know the number of Ts to be sorted and the Ts are random.
Trade space for speed and create a lookup table that can handle collisions
(AKA hash table).
*/

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <stdbool.h>
#include <pthread.h>
#include <assert.h>

#define NTHREADS    (16)
#define BINSPERTHREAD   (1*1000*1000)
bool    oneThread = false;

typedef struct {
    pthread_t   tid;
    long        begin;
    long        end;
} threadState;

void *initTID() {
    return NULL;
}

#define MAXITEMS    (50*1000*1000)
//  items on the boundary are not included in the search
#define SEARCHBOUNDARY  (100)

void usage(char *name) {
    fprintf(stderr, "usage: %s [-n 1..%d]\n", name, MAXITEMS);
}

typedef struct item {
    long    A;
    long    B;
    long    T;
    bool    unprinted;
    struct item *b;         // b(ackward to previous item)
    struct item *f;         // f(orward to next item)
    struct item *BINb;          // backward to previous bin
    struct item *BINf;          // forward to next bin
#ifdef DEVTEST
    long    lineNumber;
#endif
} item;
#ifdef DEVTEST
bool    printVerbose = false;
#endif

//  Why global variables? Because large MAXITEMS overflow the stack.
long    maxItems;           // entries allocated in list & lookup
item    *list;
long    listN;              // number of entries (max index + 1)
item    **lookup;
long    lookupN;            // number of entries (max index + 1)

/*
input -
    n       - index of current bin
    list        - global
    lookup      - global
    lookupN     - global
side-effects -
    list[] (.unprinted)
    stdout
*/
static inline void *walkThisBin(long n) {
    item    *p;
    item    *searchHead;
    item    *searchTail;
    item    *currentItem;
    long    i;

    //  for all items in bin
    for ( currentItem = lookup[n]; currentItem != lookup[n]->BINf;
        currentItem = currentItem->f)
    {
    /*
        merged forward&backward search
    */
    searchHead = currentItem;
    //  step to index min((T+100-1),lookupN-1), find largest U<T+100
    i = ((n+SEARCHBOUNDARY-1) < lookupN) ?
        n+SEARCHBOUNDARY-1 :
        lookupN-1;
    //  find largest i such that U-T<100 (U is lookup[i]->T)
    //  degenerate case is i == n
    for(p=lookup[i];
        !((p->T-searchHead->T)<SEARCHBOUNDARY);
        p=lookup[i]) {
        i--;
    }
    searchTail = p->BINf;       // boundary, not included in search
    p = currentItem;
    do {
        if (searchHead->B == p->A) {
        //  matches are symmetric
        if (searchHead->unprinted) {
            printf("%ld %ld %ld\n", searchHead->A, searchHead->B,
            searchHead->T);
            searchHead->unprinted = false;
        }
        if (p->unprinted) {
            printf("%ld %ld %ld\n", p->A, p->B, p->T);
            p->unprinted = false;
        }
        }
        p = p->f;
    } while (p!=searchTail);
    }
    return NULL;
}

/*
Must handle out-of-range indexes for lookup.

input -
    n       - index of current bin
    list        - global
    lookup      - global
    lookupN     - global
side-effects -
    list (.unprinted)
    stdout
*/

static inline void *walkTheseBins(void *tState) {
    long    startIndex = ((threadState *)tState)->begin;
    long    finishIndex = ((threadState *)tState)->end;
    long    n;

    startIndex = (startIndex<0) ? 0 : startIndex;
    finishIndex = (finishIndex>lookupN-1) ? lookupN-1 : finishIndex;
    for (n=startIndex; n<=finishIndex; n++) {
    walkThisBin(n);
    }
    return NULL;
}

int main(int argc, char *argv[]) {
#ifdef DEVTEST
item    *head;
item    *tail;
long    count = 0;
#endif
    //  subroutines? subroutines? we don't need no stinkin' subroutines
    //  this is all the scoping you're going to need
    //                      ... truuuuust me
    /*
    Allocate list[] and lookup[]. Set maxItems.

    input -
        argc
        argv
    side-effects -
        list
        lookup
        maxItems
        malloc()
        DEVTEST stuff
    */
    {
    int c;          // option character

    maxItems = MAXITEMS;
    while ((c = getopt(argc, argv, ":n:sv")) != -1) {
        switch(c) {
#ifdef DEVTEST
        case 'v':
        //  print some reassuring messages
        printVerbose = true;
        break;
#else
        case 'v':
        fprintf(stderr, "unknown option -%c\n", optopt);
        usage(argv[0]);
        exit(1);
        break;
#endif
        case 'n':
        if (sscanf(optarg, "%ld", &maxItems) != 1) {
            fprintf(stderr, "-n argument \"%s\" unscannable\n", optarg);
            usage(argv[0]);
            exit(1);
        }
        break;
        case 's':
        //  use only one thread?
        oneThread = true;
        break;
        case ':':           // -s needs an argument
        usage(argv[0]);
        exit(1);
        break;
        case '?':           // not a valid option
        fprintf(stderr, "unknown option -%c\n", optopt);
        usage(argv[0]);
        exit(1);
        break;
        }
    }
    if ((maxItems<1) || (maxItems>MAXITEMS)) {
        fprintf(stderr, "-s argument \"%ld\" out of range\n", maxItems);
        usage(argv[0]);
        exit(1);
    }
    list = (item *) malloc(sizeof(item) * maxItems);
    if (list == NULL) {
        fprintf(stderr, "ERROR: list = malloc() failure\n");
        exit(1);
    }
    lookup = (item **) malloc(sizeof(item *) * maxItems);
    if (lookup == NULL) {
        fprintf(stderr, "ERROR: lookup = malloc() failure\n");
        exit(1);
    }
    }

    /*
    Convert STDIN into an array of items.

    input -
        list
        lookup
        maxItems
    side-effects -
        list
        lookup
        listN
        stdin
    */
    {
    long    largestT = 0;
    item    x;

    for (listN=0; scanf("%ld%ld%ld", &x.A, &x.B, &x.T)==3; listN++) {
        if (listN == maxItems) {
        fprintf(stderr, "ERROR: > %ld input items read\n", maxItems);
        exit(1);
        }
        x.b = x.f = NULL;
        x.unprinted = true;
        x.BINb = x.BINf = NULL;
        largestT = (x.T>largestT) ? x.T : largestT;
#ifdef DEVTEST
        x.lineNumber = listN + 1;
#endif
        list[listN] = x;
    }
    if (!feof(stdin)) {
        fprintf(stderr, "ERROR: ferror() = %d\n", ferror(stdin));
        exit(1);
    }
    //  Be paranoid. Because cores are obnoxious.
    if (largestT>=maxItems) {
        fprintf(stderr, "ERROR: T:%ld > %ld \n", largestT, maxItems-1);
        exit(1);
    }
    }
#ifdef DEVTEST
(printVerbose) && printf("in: %ld\n", listN);
#endif
    //  Short-circuit on 0 items. Simplifies things like finding the head.
    if  (listN == 0) {
    exit(0);
    }

    /*
    Populate the lookup table. Build a doubly linked list through it.

    input -
        list
        lookup
        listN
    side-effects -
        list[]
        lookup[]
        lookupN
        DEVTEST stuff
    */
    {
    long    n;

    /*
        Populate the lookup table. The lookup table is an array-of-lists.
    The lists are LIFO. This is the most primitive of hashes, where the
    key, item.T, is used as the index into the lookup table.
    */
    for (n=0; n<maxItems; n++) {
        lookup[n] = NULL;
    }
    for (n=0; n<listN; n++) {
        long    t = list[n].T;

        if (lookup[t] == NULL) {
        lookup[t] = &(list[n]);
        } else {
        // collision
        list[n].f = lookup[t];  // forward pointer assigned
        lookup[t] = &(list[n]);
        }
    }
    /*
        Collapse lookup to squeeze out NULL references. This breaks
    the linear mapping between T value & lookup index, but worth it for
    simpler search logic. Build a doubly linked list of bins.
    */
    item    *previousBin = NULL;    // last non-NULL lookup entry
    lookupN = 0;
    for (n=0; n<maxItems; n++) {
        if (lookup[n] != NULL) {
        lookup[lookupN] = lookup[n];
        lookup[lookupN]->BINb = previousBin;
        if (previousBin) {
            previousBin->BINf = lookup[lookupN];
        }
        previousBin = lookup[lookupN];
        lookupN++;
        }
    }
    previousBin->BINf = NULL;

    /*
        Build a doubly linked list. The forward pointers already exist
    within each lookup table bin.
    */
    item    *p;
    item    *binHead;
    item    *previous;

    //  create a loop in each bin
    for (n=0; n<lookupN; n++) {
#ifdef DEVTEST
count++;
#endif
        binHead = lookup[n];
        for (p=binHead; p->f; p=p->f) {
        p->f->b = p;
#ifdef DEVTEST
count++;
#endif
        }
        p->f = binHead;
        binHead->b = p;
    }
    //  break the loops and connect them tail-to-head
#ifdef DEVTEST
head = lookup[0];
#endif
    previous = NULL;
    for (n=0; n<lookupN; n++) {
        binHead = lookup[n];
        p = binHead->b;     // p => tail of this bin list
        binHead->b = previous;  // connect bin head to list
        if (previous) {     // connect list to bin head
        previous->f = binHead;
        }
        previous = p;
    }
    previous->f = NULL;
#ifdef DEVTEST
tail = previous;
#endif
    }

#ifdef DEVTEST
if (printVerbose) {
    printf("out: %ld\n", count);

    //  run through the list forwards
    item    *p;
    count = 0;
    for (p=head; p; p=p->f) {
    count++;
    }
    printf("forwards: %ld\n", count);
    //  run through the list backwards
    count = 0;
    for (p=tail; p; p=p->b) {
    count++;
    }
    printf("backwards: %ld\n", count);
    /*
        //  print the list
        for (p=head; p; p=p->f) {
        printf("%ld %ld %ld\n", p->A, p->B, p->T);
        }
    */
}
#endif

    /*
    Find matches & print.

    (authoritative statement)
    Print item ( A B T ) if:
    1 - there exists an item ( C D U ) such that 0 <= T-U < 100 and D == A
        OR
    2 - there exists an item ( C D U ) such that 0 <= U-T < 100 and B == C

    TBD
    - threading

    input -
        lookupN
    side-effects -
        lots hidden in walkTheseBins(), all thread-local or thread-safe
    */
    {
    volatile threadState    tState[NTHREADS]; // use as cicular buffer
    long                h;  // cicular buffer head
    long                n;

    if (oneThread) {
        tState[0].begin = 0;
        tState[0].end = lookupN-1;
        walkTheseBins((void *)tState);
    } else {
        //  every slot has a thread to wait for
        for (h=0; h<NTHREADS; h++) {
        assert( pthread_create(&(tState[h].tid), NULL, initTID, NULL) == 0);
        }
        h = 0;
        for (n=0; n<lookupN+BINSPERTHREAD; n+=BINSPERTHREAD) {
        pthread_join(tState[h].tid, NULL);
        tState[h].begin = n;
        tState[h].end = n + BINSPERTHREAD - 1;
        assert( pthread_create(&(tState[h].tid), NULL, walkTheseBins, (void *)(tState+h)) == 0);
        h = (h + 1) % NTHREADS;
        }
        //  wait for any remaining threads
        for (h=0; h<NTHREADS; h++) {
        pthread_join(tState[h].tid, NULL); // may have already join'ed some
        }
    }
    }

    return 0;
}

```

Compile with "gcc -m64 -pthreads -O". Expects input on stdin. Runs multi-threaded by default. Use "-s" option to use only one thread.
> 1 votes

# Answer
I finally got a chance to build a physical Ubuntu 14.04 system similar to Lembik's and do a post-mortem on my solution to this puzzle. In my choice of importance:

1. The true master is James\_pic, because he didn't optimize prematurely.
    * he had a plan
    * he executed the plan at a high level of abstraction (Scala) and refined it there
    * he refined it further in C
    * he didn't over-refine it (see next point)
2. File system I/O time is probably the lower-bound on elapsed time for the target system.
    * Lembik alludes to this, i.e. "The winners ... both are almost as fast as wc!"
3. Some of the reasons my original solution ~~sucked~~ was non-performant:
    * Locality of reference is the dominant factor on the target system.
    * Sorting on A or B is a good idea when doing a hash sort. Sorting on T adds complexity (and cache-hostile indirection) to the hash sort, at least the way I did it.
    * Scanf() is a pig.
    * Massive bandwidth (disk-\>memory-\>cache) changes where the bottlenecks are. The target system doesn't have massive bandwidth. (see next point)
4. Rapid development works best if done in the target environment.
    * Duh! But, I was stuck with Solaris/SPARC originally and couldn't play otherwise.
    * It's difficult to eliminate caching effects in a virtualized and SAN environment.
    * Linux VMs typically have the same issues.
5. A little bit of math helps.
    * Fetching one tuple directly from the hash table cuts the probability of indirect references to ~37% (~1/e).
    * Fetching two tuples directly from the hash table would have cut references to the overflow table to ~10%. It wasn't necessary.
6. The 32-bit memory model (gcc -m32) is a distraction.
    * Sometimes a small win for unthreaded programs, sometimes a small loss.
    * Sometimes a significant loss for threaded programs.
    * If 32-bit is a significant win (and the target isn't an embedded controller), it's probably cheaper to refresh the hardware.
    * Take the extra registers and larger address space and don't look back.
7. Scanf() is a pig, but using stdio is not hopeless.
    * Most of scanf()'s overhead appears to be in the format-driven parsing and string to integer conversion.
    * Replacing sscanf() with:
        + strtok() + atoi() is ~2x faster (see table below)
        + strtol() is ~3x faster
        + a custom, local strtol() is ~6.5x faster
        + replacing strtol() with a local solution puts it on par with "wc"
        + a FSM using getc\_unlocked() is almost as fast as Keith Randall's minimalist mmap() solution
        + the results of my experiments while re-implementing in C \[in CSV, since Stack Exchange apparently doesn't do tables\]:
        ```

        "solution (64-bit unless noted)","disposition of input","user","system","elapsed"
        "dd if=? of=/dev/null bs=1024k","","0.010","1.107","26.47"
        "wc {LANG=C}","","4.921","0.752","26.38"
        "","","","",""
        "fscanf()","discard","13.130","0.490","26.43"
        "fgets(), no integer conversion","discard","1.636","0.468","26.42"
        "fgets() + sscanf()","discard","16.173","0.498","26.48"
        "fgets() + strtok(), no integer conversion","discard","4.659","0.481","26.48"
        "fgets() + strtok() + atoi()","discard","8.929","0.490","26.49"
        "fgets() + strtol()","discard","6.009","0.483","26.50"
        "fgets() + custom-strtol()","discard","3.842","0.474","26.43"
        "fgets() + custom-strtol()","sort (load hash) while reading","7.118","1.207","26.70"
        "fgets() + custom-strtol()","sort, match & print","10.096","1.357","28.40"
        "fgets() + custom-strtol(), 32-bit","sort, match & print","10.065","1.159","28.38"
        "","","","",""
        "james_pic's solution","sort, match & print","9.764","1.113","28.21"

        ```

Rather than boring you with yet another FSM parser, the solution below uses fgets() and a local strtol() replacement \[look for s2i()\].

A reference implementation in Ruby:

```
#!/usr/bin/ruby2.0
# only tested against ruby v1.9 & v2.0
=begin
Filter a file based on these rules:
Input:
  - each line is a set of three integers
  - each line is formatted as <number> <w> <number> <w> <number>
    - <w> is whitespace (a single blank in the challenge)
    - <number> is an integer in the range 1..50_000_000
Output a given tuple ( A B T ) if:
  - there exists a tuple ( C D U ) 0 <= T - U < 100 and D == A
    OR
  - there exists a tuple ( C D U ) 0 <= U - T < 100 and B == C

Typical use:
  filter.rb test.input | sort | uniq > test.output
=end
list = Array.new
lookupB = Hash.new { |hash, key| hash[key] = Array.new }
ARGF.each_with_index do |line, index|
  abt = line.split.map { |s| s.to_i }
  list << abt
  lookupB[abt[1]] << index
end
for abt in list do
  for i in Array( lookupB[abt[0]] ) do
    delta = abt[2] - list[i][2]     # T - U
    if (0<=delta) && (delta<100)
      puts "#{abt.join(' ')}"
      puts "#{list[i].join(' ')}"
    end
  end
end

```

It's a dog, ~50x slower than a C solution, but perl is just as slow and less concise.
The C solution:

```

#include <stdlib.h>
#include <stdio.h>
#include <ctype.h>
//      Throw caution, and error checking, to the winds.
// #include <assert.h>

#define RANGEMIN        (1)
#define RANGEMAX        (50*1000*1000)
#define SEARCHBOUNDARY  (100)
typedef struct {
    int             A;
    int             B;
    int             T;
} tuple_t;
typedef struct bin {
    tuple_t         slot;
    struct bin     *next;       // NULL=>0 items, self=>1 item, other=>overflow
} bin_t;
#define LISTSIZE        (RANGEMAX)
tuple_t         list[LISTSIZE];
#define HASH(x)         (x-1)
#define LOOKUPSIZE      (LISTSIZE)
bin_t           lookup[LOOKUPSIZE];
bin_t           overflow[LISTSIZE];
int             overflowNext = 0;

// based on strtol()
static inline int s2i(char *s, char **r)
{
    char            c;
    int             l = 0;

    do {
        c = *s++;
    } while (!isdigit(c));
    do {
        l = l * 10 + (c - '0');
        c = *s++;
    } while (isdigit(c));
    *r = s - 1;
    return l;
}

static inline void lookupInsert(tuple_t x)
{
    bin_t          *p = lookup + HASH(x.B);

    if (p->next) {
        overflow[overflowNext].slot = x;
        overflow[overflowNext].next = (p->next == p) ? p : p->next;
        p->next = overflow + overflowNext;
        overflowNext++;
    } else {
        p->slot = x;
        p->next = p;
    }
}

static void printOverflow(bin_t * head, bin_t * tail)
{
    if (head->next != tail) {
        printOverflow(head->next, tail);
    }
    printf("%d %d %d\n", head->slot.A, head->slot.B, head->slot.T);
}

static inline void dumpLookupSortedOnB()
{
    bin_t          *p;

    for (p = lookup; p < (lookup + LOOKUPSIZE); p++) {
        if (p->next) {
            printf("%d %d %d\n", p->slot.A, p->slot.B, p->slot.T);
            if (p != p->next) {
                printOverflow(p->next, p);
            }
        }
    }
}

static inline void printIfMatch(tuple_t abt, tuple_t cdu)
{
    int             A, B, T;
    int             C, D, U;

    A = abt.A;
    D = cdu.B;
    if (D == A) {
        T = abt.T;
        U = cdu.T;
        if ((0 <= (T - U)) && ((T - U) < SEARCHBOUNDARY)) {
            B = abt.B;
            C = cdu.A;
            printf("%d %d %d\n", A, B, T);
            printf("%d %d %d\n", C, D, U);
        }
    }
}

static inline void printMatches(int n)
{
    tuple_t        *p;

    for (p = list; p < (list + n); p++) {
        bin_t          *b = lookup + HASH(p->A);

        if (b->next) {
            bin_t          *q;

            printIfMatch(*p, b->slot);
            for (q = b->next; q != b; q = q->next) {
                printIfMatch(*p, q->slot);
            }
        }
    }
}

static inline void overflowTattle(int n)
{
    fprintf(stderr, "%d/%d items in overflow\n", overflowNext, n);
}

int main(int argc, char *argv[])
{
    int             n;

    // initialize lookup[]
    {
        bin_t          *p = lookup;

        for (n = 0; n < LOOKUPSIZE; n++) {
            p->next = NULL;
            p++;
        }
    }
    // read all tuples into list[] and insert into lookup[] & overflow[]
    {
        char            line[64];
        char           *lp;
        tuple_t        *p = list;

        for (n = 0; fgets(line, sizeof(line), stdin); n++) {
            p->A = s2i(line, &lp);
            p->B = s2i(lp, &lp);
            p->T = s2i(lp, &lp);
            lookupInsert(*p);
            p++;
        }
    }
    printMatches(n);
    exit(0);
}

```

Compile with "gcc -O3 -std=c99 -Wall -m64".
> 1 votes

---
Tags: fastest-code, file-system
---


title: Calculate the block entropy

text:
# Question
Title: Calculate the block entropy
I once needed to write a function that calculates the block entropy of a given symbol series for a given block size and was surprised how short the result was. Thus I am challenging you to codegolf such a function. I am not telling you what I did for now (and in which language), but I will in a week or so, if nobody came up with the same or better ideas first.

**Definition of the block entropy:**

Given a symbol sequence A = A\_1, …, A\_n and a block size m:

* A block of size m is a segment of m consecutive elements of the symbol sequence, i.e., A\_i, …, A\_(i+m−1) for any appropriate i.
* If x is a symbol sequence of size m, N(x) denotes the number of blocks of A which are identical to x.
* p(x) is the probability that a block from A is identical to a symbol sequence x of size m, i.e., p(x) = N(x)/(n−m+1)
* Finally, the block entropy for block size m of A is the average of −log(p(x)) over all blocks x of size m in A or (which is equivalent) the sum of −p(x)·log(p(x)) over every x of size m occurring in A. (You can choose any reasonable logarithm you want.)

**Restricions and clarifications:**

* Your function should take the symbol sequence A as well as the block size m as an argument.
* You may assume that the symbols are represented as zero-based integers or in another convenient format.
* Your program should be capable of taking any reasonable argument in theory and in reality should be able to calculate the example case (see below) on a standard computer.
* Built-in functions and libraries are allowed, as long as they do not perform big portions of the procedure in one call, i.e., extracting all blocks of size m from A, counting the number of occurrences of a given block x or calculating the entropies from a sequence of p values – you have to do those things yourself.

**Test:**

```
[2, 3, 4, 1, 2, 3, 0, 0, 3, 2, 3, 0, 2, 2, 4, 4, 4, 1, 1, 1, 0, 4, 1,
2, 2, 4, 0, 1, 2, 3, 0, 2, 3, 2, 3, 2, 0, 1, 3, 4, 4, 0, 2, 1, 4, 3,
0, 2, 4, 1, 0, 4, 0, 0, 2, 2, 0, 2, 3, 0, 0, 4, 4, 2, 3, 1, 3, 1, 1,
3, 1, 3, 1, 0, 0, 2, 2, 4, 0, 3, 2, 2, 3, 0, 3, 3, 0, 0, 4, 4, 1, 0,
2, 3, 0, 0, 1, 4, 4, 3]

```

The first block entropies of this sequence are (for the natural logarithm):

* m = 1: 1.599
* m = 2: 3.065
* m = 3: 4.067
* m = 4: 4.412
* m = 5: 4.535
* m = 6: 4.554

# Answer
> 6 votes
## Mathematica - ~~81~~ ~~78~~ ~~75~~ ~~72~~ ~~67~~ ~~65~~ ~~62~~ 56 bytes

I haven't golfed anything in Mathematica before, so I suppose there's room for improvement. This one doesn't quite conform to the rules due to the `Partition` and `Tally` functions, but it's quite neat so I thought I'd post it anyway.

```
f=N@Tr[-Log[p=#2/Length@b&@@@Tally[b=##~Partition~1]]p]&

```

This works with any set of symbols, and can be used like

```
sequence = {2, 3, 4, 1, 2, 3, 0, 0, 3, 2, 3, 0, 2, 2, 4, 4, 4, 1, 1,
   1, 0, 4, 1, 2, 2, 4, 0, 1, 2, 3, 0, 2, 3, 2, 3, 2, 0, 1, 3, 4, 4,
   0, 2, 1, 4, 3, 0, 2, 4, 1, 0, 4, 0, 0, 2, 2, 0, 2, 3, 0, 0, 4, 4,
   2, 3, 1, 3, 1, 1, 3, 1, 3, 1, 0, 0, 2, 2, 4, 0, 3, 2, 2, 3, 0, 3,
   3, 0, 0, 4, 4, 1, 0, 2, 3, 0, 0, 1, 4, 4, 3};
f[sequence, 3]

> 4.06663

```

Here is a somewhat ungolfed version:

```
f[sequence_, m_] := (
    blocks = Partition[sequence, m, 1];
    probabilities = Apply[#2/Length[blocks] &, Tally[blocks], {1}];
    N[Tr[-Log[probabilities]*probabilities]]
)

```

It will probably run faster if I apply `N` directly to the result of `Tally`.

By the way, Mathematica does actually have an `Entropy` function, that reduces this to **28 bytes**, but that's definitely against the rules.

```
f=N@Entropy@Partition[##,1]&

```

On the other hand, here is a **128 byte** version that reimplements `Partition` and `Tally`:

```
f=N@Tr[-Log[p=#2/n&@@@({#[[i;;i+#2-1]],1}~Table~{i,1,(n=Length@#-#2+1)}//.{p___,{s_,x_},q___,{s_,y_},r___}:>{p,{s,x+y},q,r})]p]&

```

Ungolfed:

```
f[sequence_, m_] := (
    n = Length[sequence]-m+1; (*number of blocks*)
    blocks = Table[{Take[sequence, {i, i+m-1}], 1},
                   {i, 1, n}];
    blocks = b //. {p___, {s_, x_}, q___, {s_, y_}, r___} :> {p,{s,x+y},q,r};
    probabilities = Apply[#2/n &, blocks, {1}];
    N[Tr[-Log[probabilities]*probabilities]]
)

```

# Answer
> 4 votes
# Perl, 140 bytes

The following Perl script defines a function `E` that takes the symbol sequence, followed by segment size as arguments.

```
sub E{$m=pop;$E=0;%h=();$"=',';$_=",@_,";for$i(0..@_-$m){next
if$h{$s=",@_[$i..$i+$m-1],"}++;$E-=($p=s|(?=$s)||g/(@_-$m+1))*log$p;}return$E}

```

## Ungolfed version with test

```
sub E { # E for "entropy"
    # E takes the sequence and segment size as arguments
    # and returns the calculated entropy.
    $m = pop;    # get segment size (last argument)
    $E = 0;      # initialize entropy
    %h = ();     # hash that remembers already calculated segments
    $" = ',';#"  # comma is used as separator
    $_ = ",@_,"; # $_ takes sequence as string, with comma as delimiters
    for $i (0 .. @_-$m) {
        $s = ",@_[$i..$i+$m-1],"; # segment
        next if$h{$s}++;          # check, if this segment is already calculated
        $p = s|(?=\Q$s\E)||g / (@_ - $m + 1); # calculate probability
             # N(x) is calculated using the substitution operator
             # with a zero-width look-ahead pattern
             # (golfed version without "\Q...\E", see below)
        $E -= $p * log($p); # update entropy
    }
    return $E
}

# Test

my @A = (
    2, 3, 4, 1, 2, 3, 0, 0, 3, 2, 3, 0, 2, 2, 4, 4, 4, 1, 1, 1, 0, 4, 1,
    2, 2, 4, 0, 1, 2, 3, 0, 2, 3, 2, 3, 2, 0, 1, 3, 4, 4, 0, 2, 1, 4, 3,
    0, 2, 4, 1, 0, 4, 0, 0, 2, 2, 0, 2, 3, 0, 0, 4, 4, 2, 3, 1, 3, 1, 1,
    3, 1, 3, 1, 0, 0, 2, 2, 4, 0, 3, 2, 2, 3, 0, 3, 3, 0, 0, 4, 4, 1, 0,
    2, 3, 0, 0, 1, 4, 4, 3
);

print "m = $_: ", E(@A, $_), "\n" for 1 .. @A;

```

**Result:**

```
m = 1: 1.59938036027528
m = 2: 3.06545141203611
m = 3: 4.06663334311518
m = 4: 4.41210802885304
m = 5: 4.53546705894451
m = 6: 4.55387689160055
m = 7: 4.54329478227001
m = 8: 4.53259949315326
m = 9: 4.52178857704904
...
m = 97: 1.38629436111989
m = 98: 1.09861228866811
m = 99: 0.693147180559945
m = 100: 0

```

**Symbols:**

The symbols are not restricted to integers, because pattern matching based on strings is used. The string representation of a symbol must not contain the comma, because it is uses as delimiter. Of course, different symbols must have different string representations.

In the golfed version, the string representation of the symbols should not contain specials characters of patterns. The additional four bytes `\Q`...`\E` are not needed for numbers.

# Answer
> 2 votes
# Python ~~127 152B~~ 138B

```
import math
def E(A,m):N=len(A)-m+1;R=range(N);return sum(math.log(float(N)/b) for b in [sum(A[i:i+m]==A[j:j+m] for i in R) for j in R])/N

```

Adjusted to not break the rules any more and have a slightly cuter algorithm. Adjusted to be smaller

### Older version:

```
import math
def E(A,m):
 N=len(A)-m+1
 B=[A[i:i+m] for i in range(N)]
 return sum([math.log(float(N)/B.count(b)) for b in B])/N

```

My first ever Python script! See it in action: http://pythonfiddle.com/entropy

# Answer
> 1 votes
# Python with Numpy, ~~146~~ 143 Bytes

As promised, here is my own solution. It requires an input of non-negative integers:

```
from numpy import*
def e(A,m):
    B=zeros(m*[max(A)+1]);j=0
    while~len(A)<-j-m:B[tuple(A[j:j+m])]+=1;j+=1
    return -sum(x*log(x)for x in B[B>0]/j)

```

The disadvantage is that this bursts your memory for a large `m` or `max(A)`.

Here is the mostly ungolfed and commented version:

```
from numpy import *
def e(A,m):
    B = zeros(m*[max(A)+1])          # Generate (max(A)+1)^m-Array of zeroes for counting.
    for j in range(len(A)-m+1):
        B[tuple(A[j:j+m])] += 1      # Do the counting by directly using the array slice
                                     # for indexing.
    C = B[B>0]/(len(A)-m+1)          # Flatten array, take only non-zero entries,
                                     # divide for probability.
    return -sum(x*log(x) for x in C) # Calculate entropy

```

# Answer
> 1 votes
## MATLAB

```
function E =BlockEntropy01(Series,Window,Base )

%-----------------------------------------------------------
% Calculates BLOCK ENTROPY of Series
% Series: a Vector of numbers
% Base: 2 or 10 (affects logarithm of the Calculation)
% for 2 we use log2, for 10 log10
% Windows: Length of the "Sliding" BLOCK
% E: Entropy
%-----------------------------------------------------------
% For the ENTROPY Calculations
% http://matlabdatamining.blogspot.gr/2006/....
% 11/introduction-to-entropy.html
% BlogSpot: Will Dwinnell
%-----------------------------------------------------------
% For the BLOCK ENTROPY
% http://codegolf.stackexchange.com/...
% questions/24316/calculate-the-block-entropy
%-----------------------------------------------------------
% Test (Base=10)
% Series=[2, 3, 4, 1, 2, 3, 0, 0, 3, 2, 3, 0, ....
%     2, 2, 4, 4, 4, 1, 1, 1, 0, 4, 1,2, 2, 4, 0, ....
%     1, 2, 3, 0, 2, 3, 2, 3, 2, 0, 1, 3, 4, 4, 0, ....
%     2, 1, 4, 3,0, 2, 4, 1, 0, 4, 0, 0, 2, 2, 0, ....
%     2, 3, 0, 0, 4, 4, 2, 3, 1, 3, 1, 1,3, 1, 3, 1, ....
%     0, 0, 2, 2, 4, 0, 3, 2, 2, 3, 0, 3, 3, 0, 0, 4, ...
%     4, 1, 0,2, 3, 0, 0, 1, 4, 4, 3]';
%
% Results
%
% Window=1: 1.599
% Window=2: 3.065
% Window=3: 4.067
% Window=4: 4.412
% Window=5: 4.535
% Window=6: 4.554
%-----------------------------------------------------------
n=length(Series);
D=zeros(n,Window); % Pre Allocate Memory
for k=1:Window;    D(:,k)=circshift(Series,1-k);end
D=D(1:end-Window+1,:); % Truncate Last Part
%
% Repace each Row with a "SYMBOL"
% in this Case a Number ...............
[K l]=size(D);
for k=1:K; MyData(k)=polyval(D(k,:),Base);end
clear D
%-----------------------------------------------------------
% ENTROPY CALCULATIONS on MyData
% following  Will Dwinnell
%-----------------------------------------------------------
UniqueMyData = unique(MyData);
nUniqueMyData = length(UniqueMyData);
FreqMyData = zeros(nUniqueMyData,1); % Initialization
for i = 1:nUniqueMyData
    FreqMyData(i) = ....
        sum(double(MyData == UniqueMyData(i)));
end
% Calculate sample class probabilities
P = FreqMyData / sum(FreqMyData);
% Calculate entropy in bits
% Note: floating point underflow is never an issue since we are
%   dealing only with the observed alphabet
if Base==10
    E= -sum(P .* log(P));
elseif BASE==2
    E= -sum(P .* log2(P));
else
end
end

WITH TEST SCRIPT
%-----------------------------------------------------------
Series=[2, 3, 4, 1, 2, 3, 0, 0, 3, 2, 3, 0, ....
    2, 2, 4, 4, 4, 1, 1, 1, 0, 4, 1,2, 2, 4, 0, ....
    1, 2, 3, 0, 2, 3, 2, 3, 2, 0, 1, 3, 4, 4, 0, ....
    2, 1, 4, 3,0, 2, 4, 1, 0, 4, 0, 0, 2, 2, 0, ....
    2, 3, 0, 0, 4, 4, 2, 3, 1, 3, 1, 1,3, 1, 3, 1, ....
    0, 0, 2, 2, 4, 0, 3, 2, 2, 3, 0, 3, 3, 0, 0, 4, ...
    4, 1, 0,2, 3, 0, 0, 1, 4, 4, 3]';
Base=10;
%-----------------------------------------------------------
for Window=1:6
    E =BlockEntropy01(Series,Window,Base )
end

```

---
Tags: code-golf
---


title: Recover a bzip2 file

text:
# Question
Title: Recover a bzip2 file
The popular `.bz2` compression format is used to store all sorts of things. One of the more interesting features is that it consists of several independently decompressible blocks, allowing recovery of partial segments if the file is damaged. The `bzip2recover` program, included with `bzip2`, allows recovery of blocks from a `bzip2` archive. In this challenge, you will implement a simple clone of `bzip2recover`.

The basic idea is to scan a bitstream for a particular magic number, and output blocks delimited by this magic number to separate files. Read on for gory details.

---

Technically, bzip2 files are a *bitstream*, with every group of 8 bits being stuffed into a byte. The first bit is placed in the most significant bit of a byte, with the eighth bit going in the least significant bit of a byte. Files begin with the header "BZh", and a digit "N" from 1 to 9 indicating the compression level (9 is the most common). Then the bitstream begins.

The bitstream is split up into blocks. Each block encodes N\*100KB of uncompressed text (e.g. 900KB of text for compression level 9). Blocks can begin basically anywhere in the bitstream, and do not necessarily begin on 8-bit boundaries. Blocks begin with a magic number, which is set to the 48-bit sequence 0x314159265359 (in ASCII, "1AY&SY"). For example, the sequence of bytes "73 14 15 92 65 35 90" contains the magic number, offset by 4 bits within the byte.

Your goal will be to take a bzip2 file on stdin, and output each detected BZip2 block to a separate decompressible file. Specifically, your program does the following:

* Read and skip the bzip2 header ("BZh" + N)
* Scan the bitstream for the magic sequence
* Copy the bitstream between two consecutive magic sequences to a new file, with a new BZip2 header and block header. Note that this may require realigning the bitstream to start at a byte boundary, and padding the last byte.

You don't have to care about parsing the blocks, and you can assume that a magic sequence always starts a block (i.e. that it isn't in the middle of a block).

How you name the output files is up to you. The only requirement is that they be clearly sequential somehow (so, e.g. "a b c ... z aa ab" or "1 2 3 4 5 .. 9 10 11 12" are both acceptable).

---

For a test case, you may use this file from the IPSC 2014 problem-solving contest. The file is a bzip2 file with four blocks. The file as a whole actually decompresses to another bzip2 file with 162 complete blocks. Therefore, you may use these block counts to validate that your solution is getting all of the blocks. (Note that the file actually decompresses even further into more bzip2 files, but you probably don't want to decompress it all the way).

---

This is code golf, with a twist: your score is the number of bytes in the bz2-compressed version of your source code. (I'm aware this may make submissions *larger*, but these are the rules!). Since bzip2 compressors may differ in performance, provide a base64-encoded version of the compressed file to prove it's size. (You may therefore be able to obtain a better score by optimizing your program for bzip2, though I don't necessarily expect anyone to do this). Standard loopholes apply.

# Answer
# Python 3 -- 375(534 chars)

It can process the example file `d2.in` in about 10 seconds, using Pypy3. With CPython 3, it will cost 60 seconds.

It's interesting and weird that some old golf tricks would make the score worse! For example, replacing `<space>+<space>` by `+` would increase the score from 375 to 384; replacing all `range` by `R` will also increase the score to 384.

## Golfed

```
import sys
from collections import deque
def T(B):
 for i in range(0,len(B),8):
  x=0
  for j in range(min(8,len(B)-i)): x |= B[i + j] << (7-j)
  yield x
V=sys.argv
C=F=0
B=bytearray()
I=open(V[1],'rb')
H=I.read(4)
P=deque([(x>>(7-j)) & 1 for x in b'1AY&SY' for j in range(8)])
Q=deque([0]*48)
for x in I.read():
 for j in range(8):
  b=(x>>(7-j)) & 1
  B.append(b)
  Q.popleft()
  Q.append(b)
  if Q!=P: continue
  C+=1
  M,B=B[:-48],B[-48:]
  if F:F.write(bytes(T(M)))
  F=open(V[2] + str(C),'wb')
  F.write(H)
F.write(bytes(T(B)))

```

## Base64

```
QlpoOTFBWSZTWSoMWQ0AAMLfgAAQYf901zlibSo/t//kMAFaCIaCQANTTTKMnlGg8oeRpPaoNDDAAAAAaAAAAABKaSGmmpplD0jQaAABo00ZNqLfKsdHBd3QgkQ2wY3RfYpMMjPpsyit5Ee6+k6E6ObzlFT2vmgbSKPkni7ZgDkcwTRiLFBmNCuEX/HtOC8KJ+Rg51SkMCiAGXZpXbvaK+DzMLvaEhAM8reCB8HFAsNYfRghS9+i5/IpLrDzGZVhsbVOpQzHIIVZ7zCCgDkcElOMvTuvfw2bruKz2WQuvEg99gDa8T3Kw2GdyMNjEVFVVQj4sU1MaOTUcwQROwt86LfUZAOJE6HqnB60CFO38NJmpifx+trcQ8x5alqRiEGoKC2wJpPkCXpT7fbtFDbLCLQaplK23qMTC5gqGwCghjN3zEtwyKhCz3kixqBr7gya6yPG03rMOGY9hlmV3tS0W/SdC2taUpCoNSEaEnUXckU4UJAqDFkN

```
> 3 votes

---
Tags: code-golf, compression, bitwise
---


title: 'Number' Crunching Benchmark

text:
# Question
Title: 'Number' Crunching Benchmark
I just want to do some benchmarking of different programming languages. So here is a tip: You can use as many characters as you want. In the end there is only the time your program needed.

So, here is the Task:

# Task

After searching for a good task that would require a lot of computation, I came up with the problem `Is a word in a given Type-1 grammar?`

For people not knowing what this is: Formal Grammar Chomsky-Hierarchy

Maybe there is a better algorithm but I came only up with brute forcing.

So, write a code that brute forces the Grammars as fast as possible.

# Rules

* Parallel execution is allowed
* Delete your cache after every tested word
* Add your complete compiler command for optimized result
* You can use any algorithm and language you want, as long as the job gets done.
* Add to your answer what files are read/how to set the files that should be read.

# I Guarantee

* Input is always well formatted
* Variables are always capital letters, Terminals always small ones
* Words only consist of Terminals
* No empty Word

# Input

For equal footing we need to define what the input looks like. The input is a file. I write multiple grammars and words in one file. The structure can be described with following Pseudocode:

```
foreach(Grammar G) do
    write('{');
    foreach (Variable V in G) do
        write(T + ',');
    od
    deleteLast();
    write('}');
    linebreak();
    write('{');
    foreach (Terminal T in G) do
        write(V + ',');
    od
    deleteLast();
    write('}');
    linebreak();
    write('{');
    foreach(Production P in G) do
        write(P + ',');
    od
    write('}');
    linebreak;
    write(StartSymbol(G));
    linebreak;
    foreach(Word W for G) do
        write(W + ',');
    od
    deleteLast();
    linebreak;
od
deleteLastLinebreak();

```

An Example

# Output

The output has to be put in a file as well. The formatting of the output does not matter. It should be recognizable which result belongs to which word.

An Example

# Example

Here is an example written in C# by myself.

Main Class

Data Classes

# Test System

The system the code will run on is following:

* Arch Linux
* RunTime you need (if available; Please add in answer)
* Intel i7-4700MQ 3.4GHz
* 16GB RAM

# How do I test

After the code is compiled I'll run a bash script that runs your program a certain amount of times and calculates the average and median run time.

Good Luck and may the best code win.

~~My Benchmark: Not completed after 8-9 hours of runtime.~~

My Benchmark (after changing some words):

* mean: 39.5157
* median: 39.448687194

PS: Here is the originally planned input, which took ages with my program. Grammars are slightly changed in the new input and words are way shorter.

# Answer
> 1 votes
As the code is rather long, it can be found here (Java 8): http://pastebin.com/Q5Buuqfx

How I did it:

Instead of brute-force creating the possible word combinations and then try to find a match I instead did the reverse and tried to collapse the target word to find the starting point. As a result the path the algorithm has to traverse is much less branched, and the words from the example can all be verified in about 35 ms (on my machine).

As a side-effect of the reverse-search I do not need the first two lines of each task (variables, terminals).

In addition I utilize the Java 8 streaming technology, which allows to easily activate parallel processing with just a single change in one line of code.

Usage: `java Grammar <inputFile>`

**Edit:** Did some code cleanup and further speed tests and it ironically turned out, that the algorithm is so fast, that setting up the parallel stream takes longer than the actual task. I now use a single-CPU `stream()` instead and the complete example input set is processed in ~35 ms on my machine instead of ~80 ms if I use `parallelStream()`.

**Edit2:** Further testing reveals some strange effects with parallelStream and certain types of Collections. A TreeSet seems to work best in both single and parallel processing. I did some tests on your extended input set and it seems to work fine, except for the task with the looping rules (CD-\>DC,DC-\>CD). Detecting and removing those would cause the search to verify instantly.

---
Tags: fastest-code
---


title: Three-pointer! But what kind?

text:
# Question
Title: Three-pointer! But what kind?
From http://en.wikipedia.org/wiki/Triangle:

---

Write a program that takes three 2d coordinate tuples (Cartesian), and classifies what shape these three points describe.

In almost all cases these points will describe a triangle of varying types. In some degenerate cases, the points will either describe a singular point or a straight line. The program will determine which of the following tags apply to the described shape:

* Point (3 points are co-incident)
* Line (3 points lie on a straight line - no more than 2 points may be co-incident)
* Equilateral (3 sides equal, 3 angles equal)
* Isosceles (2 sides equal, 2 angles equal)
* Scalene (0 sides equal, 0 angles equal)
* Right (1 angle exactly π/2 (or 90°))
* Oblique (0 angles exactly π/2 (or 90°))
* Obtuse (1 angle \> π/2 (or 90°))
* Acute (3 angles \< π/2 (or 90°))

Note that for some described shapes, more than one of the above tags will apply. For example, any right-angled will also either be isosceles or scalene.

### Input

* The program may read the 3 input coordinates from STDIN, command-line, environment variables or whatever method is convenient for your language of choice.
* The input co-ordinates my be formatted however is convenient for your language of choice. It can be assumed that all input numbers are well-formed with respect to the datatypes you end up using.
* Nothing can be assumed about the ordering of the input coordinates.

### Output

* The program will output to STDOUT, dialog box or whatever display method is convenient for your language of choice.
* The output will display all the tags applicable to the shape described by the input coordinates.
* Tags may be output in any order.

### Other Rules

* Your language's trigonometric libraries/APIs are allowed, but any APIs that specifically calculate triangle types are banned.
* When determining equality of angles or lengths of sides, you will likely end up comparing floating-point values. Two such values are to be considered "equal" if one is within 1% of the other.
* Standard “loopholes” which are no longer funny
* This is , so the shortest answer in bytes wins.

### Examples

```
Input                   Output
(1,2) (1,2) (1,2)       Point
(1,2) (3,4) (5,6)       Line
(0,0) (1,1) (2,0)       Isosceles Right
(0,0) (2,1) (10,1)      Scalene Oblique Obtuse

```

# Answer
> 10 votes
# C (451 bytes)

Uses only squared lengths and slopes.

```
p[2],q[2],r[2];z(c){char*y[]={"Line","Point","Isosceles ","Equilateral ","Scalene ","Right","Oblique ","Acute","Obtuse"};printf(y[c]);}d(int*a,int*b){int c=*a++-*b++,e=*a-*b;return c*c+e*e;}main(){scanf("%d%d%d%d%d%d",p,p+1,q,q+1,r,r+1);int a=d(p,q),b=d(q,r),c=d(r,p),e=!a+!b+!c,f=(a==b)+(b==c)+(c==a),g=a>b&&b>c?a:b>c?b:c,h=g^a?g^b?a+b:c+a:b+c;e?z(e/2):(1[q]-1[p])*(*r-*q)^(1[r]-1[q])*(*q-*p)?f?z(2+f/2),f-1&&z(2):z(4),h^g?z(6),z(7+(h<g)):z(5):z(0);}

```

Ungolfed (and ternary operator replaced with if/else):

```
int p[2],q[2],r[2];

void print(c){
    char *y[]={"Line","Point","Isosceles ","Equilateral ","Scalene ","Right","Oblique ","Acute","Obtuse"};
    printf(y[c]);
}
squared_distance(int *a,int *b){
    int c = *a++ - *b++, e = *a - *b;
    return c*c+e*e;
}
main(){
    scanf("%d%d%d%d%d%d",p,p+1,q,q+1,r,r+1); // read in coordinates
    int a = squared_distance(p,q),b = squared_distance(q,r),c = squared_distance(r,p),
    e=!a+!b+!c, // number of sides of length 0
    f=(a==b)+(b==c)+(c==a), // number of equal-length pairs
    g = a > b && b > c ? a : (b > c ? b : c), // longest side
    h = g != a ? g != b ? a + b : c + a : b + c; // sum of squares of length of other two sides
    if(e)
        print(e/2); // 1 side of len 0: line, 3 sides: point
    // comparing slopes PQ and QR
    else if((q[1]-p[1])*(*r-*q) != (r[1]-q[1])*(*q-*p)){ // not line
        if(f){
            print(2+f/2); // 1 pair of equal length sides: isosceles, 3: equilateral
            if(f-1) print(2); // equilateral therefore also isosceles
        }else print(4); // 0: scalene
        if(h!=g){ // a^2+b^2!=c^2: not right
            print(6); // oblique
            print(7+(h<g)); // a^2+b^2<c^2:obtuse, acute otherwise
        }else print(5); // right
    }else
        print(0); // line
}

```

Input (through stdin) Format: x y x y x y

ex. 0 0 1 1 2 0 for Isosceles Right

# Answer
> 3 votes
# C,333

```
z,c,r,b,s,i,t[14],g[14];
main(){
  for(;i<14;i++){
    g[i]=r=t[(i+2)%6]-t[i%6];r*=r;t[i|1]+=r;
    i<6&&scanf("%d",t+i);
    i>7&&(b<t[i]&&(b=t[i]),s+=t[i],z+=!t[i],c+=t[i]==t[i-2]);
  }

  if(g[6]*g[9]==g[8]*g[7])puts(z==6?"point":"line");else
    printf(b*2==s?"right ":"oblique %s",b*2>s?"obtuse ":"acute "),puts(c>3?c>5?"equilateral":"isosceles":"scalene");
}

```

I left the whitespace in for the moment. This works but could probably do with some tidying up and golfing. Similar math to `@es1024`'s answer, but uses a loop and arrays. Input format `x y x y x y`

Variables

`t[]` stores both the input and the squares of the lengths. By the end of the program it looks like the table below (increasing the number of iterations of the loop would lead to indefinite repitition of the squared lengths.) At the start of the loop squares of lengths (garbage as not all the data is available) are needlessly stored in cells 1,3 and 5, but are promptly overwritten by `scanf.` Useful data is written to `z,b,c` ahd `s` when `i`=9,11,13 (`t[i]` and `t[i-2]` are accessed.)

```
01 23 45 67 89 1011 1213
aa bb cc  a  b    c    a
xy xy xy  L  L    L    L

```

`g[]` was added late on to hold the dx and dy values needed for slope calculation. The only cells used are 6 through 9 (0 through 5 are unreliable as not all the data is available when they are written.) If `g[6]/g[7]==g[8]/g[9]` the slopes of 2 lines are equal and the triangle is just a line (or a point.) The equation is rearranged in the program to avoid division.

`r`is an intermediate value used for squaring

`z`counts the number of sides of length zero. It has an offset of +3 because the loop reads 3 blank cells in `t[]`.

`c`counts the number of sides that are of identical length. It also has an offset of +3. Note that side `a` is written to `t[]` twice in order to be able to check a=b,b=c,c=a.

`b`is the biggest length of a side, squared. `s`is the sum of the squares of all sides.

Note that comparing side lengths A^2+B^2+C^2 with 2\*B^2 is the same as comparing A^2+C^2 with B^2 (just subtract B^2 from both sides.) Thus if B^2=A^2+C^2 it is a right triangle. if B^2 is greater it is obtuse, if smaller it is acute.

# Answer
> 2 votes
# Golfscript (175 bytes)

```
~..|,({.)2$([\;\]@(;]{{~}/@- 2?@@- 2?+}%$.{2-1??100*}/-+abs 1<{;"Line"}{.[]|,((["Isosceles ""Scalene "]=\~-+.!["Oblique ""Right "]=\.!\0>-)["Acute ""Obtuse "]=}if}{;"Point "}if

```

You can test it here (test set included).

**Input format:**

```
"[x y][x y][x y]"

```

**Commented version:**

```
~                       # evaluates input string
..|,(                   # pushes the number of unique coords - 1
{
  .)2$([\;\]@(;]        # makes all 3 possible pairings of coords
  {{~}/@- 2?@@- 2?+}%$  # gets all squares of side lengths
  .{2-1??100*}/-+abs 1< # 0 if triangle, 1 if line
  {;"Line"}
  {
     .[]|,((["Isosceles ""Scalene "]=\   # removes duplicate side-squares,
                                         #   and use the count to determine
                                         #   if isosceles or scalene (no
                                         #   equilaterals will exist)
     ~-+.!["Oblique ""Right "]=\         # compute a^2 + b^2 - c^2. Use to
                                         #   determine if oblique or right.
                                         #   c = max side length
     .!\0>-)["Acute ""Obtuse "]=         # use same value to determine if
                                         #   acute, obtuse, or right
  }
  if
}
{;"Point "}
if

```

**NOTE:**

The reason my code contains no "equilateral" output is because:

* OP said "all input numbers are well-formed with respect to the datatypes you end up using"
* Golfscript does not have floating point numbers - not inherently anyway
* It's impossible (in a 2-dimensional grid) to have an equilateral triangle with integer coordinates, as proven here.

# Answer
> 0 votes
# Mathematica (~~313~~ 307 characters)

Golfed:

```
f@p_:=(P=Print;R=RotateLeft;L=Length;U=Union;If[L@U@p==1,P@"Point",If[Det[Join[#,{1}]&/@p]==0,P@"Line",v=p-R@p;a=MapThread[VectorAngle[#,#2]&,{-v,R@v}];u=L@U[Norm/@v];If[u==1,P@"Equilateral",If[u==2,P@"Isosceles",P@"Scalene"]];If[MemberQ[a,Pi/2],P@"Right",P@"Oblique";If[Max@a>Pi/2,P@"Obtuse",P@"Acute"]]]])

```

Ungolfed:

```
f@p_ := (
  P = Print;    (* make aliases for functions used more than once *)
  R = RotateLeft;
  L = Length;
  U = Union;
  If[L@U@p == 1,    (* if all points identical *)
   P@"Point",
   If[Det[Join[#, {1}] & /@ p] == 0,    (* if area is zero *)
    P@"Line",
    v = p - R@p;    (* cyclic vectors *)
    a = MapThread[VectorAngle[#, #2] &, {-v, R@v}];    (* interior angles *)
    u = L@U[Norm /@ v];    (* number of unique side lengths *)
    If[u == 1,
     P@"Equilateral",
     If[u == 2,
      P@"Isosceles",
      P@"Scalene"
      ]
     ];
    If[MemberQ[a, Pi/2],
     P@"Right",
     P@"Oblique";
     If[Max@a > Pi/2,
      P@"Obtuse",
      P@"Acute"
      ]
     ]
    ]
   ]
  )

```

The input format is a list of points, upon which the function is called:

```
points = {{x1,y1},{x2,y2},{x3,y3}};
f@points

```

---
Tags: code-golf, math, geometry
---