test/bench: import new fasta C reference, update Go, optimizations

OLD fasta -n 25000000
	gcc -O2 fasta.c	7.59u 0.06s 7.74r
	gc fasta	9.54u 0.15s 9.84r
	gc_B fasta	9.48u 0.10s 9.62r

NEW fasta -n 25000000
	gcc -O2 fasta.c	2.59u 0.02s 2.66r
	gc fasta	3.00u 0.03s 3.09r
	gc_B fasta	2.72u 0.03s 2.81r

R=r
CC=golang-dev
https://golang.org/cl/1054041

2 files changed, 297 insertions, 222 deletions

diff --git a/test/bench/fasta.c b/test/bench/fasta.c
index 65f4d3d35d..78a8490d71 100644
--- a/test/bench/fasta.c
+++ b/test/bench/fasta.c
@@ -28,146 +28,190 @@ POSSIBILITY OF SUCH DAMAGE.
 */
 
 /*
- * http://shootout.alioth.debian.org/u32q/benchmark.php?test=fasta&lang=gcc&id=4
-*/
-/* The Computer Language Benchmarks Game
- * http://shootout.alioth.debian.org/
- * Contributed by Joern Inge Vestgaarden
- * Modified by Jorge Peixoto de Morais Neto
+ * http://shootout.alioth.debian.org/u32/program.php?test=fasta&lang=gcc&id=3
+ */
+
+/*  The Computer Language Benchmarks Game
+ *  http://shootout.alioth.debian.org/
+ *
+ *  contributed by Petr Prokhorenkov
  */
 
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
-#include <err.h>
-
-#define WIDTH 60
-#define MIN(a,b) ((a) <= (b) ? (a) : (b))
-#define NELEMENTS(x) (sizeof (x) / sizeof ((x)[0]))
-
-typedef struct {
-    float p;
-    char c;
-} aminoacid_t;
-
-static inline float myrandom (float max) {
-    unsigned long const IM = 139968;
-    unsigned long const IA = 3877;
-    unsigned long const IC = 29573;
-    static unsigned long last = 42;
-    last = (last * IA + IC) % IM;
-    /*Integer to float conversions are faster if the integer is signed*/
-    return max * (long) last / IM;
+
+// not available on OS X 
+#define fwrite_unlocked fwrite
+#define fputc_unlocked fputc
+#define fputs_unlocked fputs
+
+#define ARRAY_SIZE(a) (sizeof(a)/sizeof(a[0]))
+#define unlikely(x) __builtin_expect((x), 0)
+
+#define IM 139968
+#define IA 3877
+#define IC 29573
+
+#define LINE_LEN 60
+#define LOOKUP_SIZE 4096
+#define LOOKUP_SCALE ((float)(LOOKUP_SIZE - 1))
+
+typedef unsigned random_t;
+
+void
+random_init(random_t *random) {
+    *random = 42;
+}
+
+// Special version with result rescaled to LOOKUP_SCALE.
+static inline
+float
+random_next_lookup(random_t *random) {
+    *random = (*random*IA + IC)%IM;
+
+    return (*random)*(LOOKUP_SCALE/IM);
 }
 
-static inline void accumulate_probabilities (aminoacid_t *genelist, size_t len) {
-    float cp = 0.0;
-    size_t i;
-    for (i = 0; i < len; i++) {
-        cp += genelist[i].p;
-        genelist[i].p = cp;
+struct amino_acid {
+   char sym;
+   float prob;
+   float cprob_lookup;
+};
+
+void
+repeat(const char *alu, const char *title, int n) {
+    int len = strlen(alu);
+    char buffer[len + LINE_LEN];
+    int pos = 0;
+
+    memcpy(buffer, alu, len);
+    memcpy(buffer + len, alu, LINE_LEN);
+
+    fputs_unlocked(title, stdout);
+    while (n > 0) {
+        int bytes = n > LINE_LEN ? LINE_LEN : n;
+
+        fwrite_unlocked(buffer + pos, bytes, 1, stdout);
+        pos += bytes;
+        if (pos > len) {
+            pos -= len;
+        }
+        fputc_unlocked('\n', stdout);
+        n -= bytes;
     }
 }
 
-/* This function prints the characters of the string s. When it */
-/* reaches the end of the string, it goes back to the beginning */
-/* It stops when the total number of characters printed is count. */
-/* Between each WIDTH consecutive characters it prints a newline */
-/* This function assumes that WIDTH <= strlen (s) + 1 */
-static void repeat_fasta (char const *s, size_t count) {
-    size_t pos = 0;
-    size_t len = strlen (s);
-    char *s2 = malloc (len + WIDTH);
-    memcpy (s2, s, len);
-    memcpy (s2 + len, s, WIDTH);
-    do {
-     	size_t line = MIN(WIDTH, count);
-     	fwrite (s2 + pos,1,line,stdout);
-     	putchar_unlocked ('\n');
-     	pos += line;
-     	if (pos >= len) pos -= len;
-     	count -= line;
-    } while (count);
-    free (s2);
+/*
+ * Lookup table contains mapping from real values to cumulative
+ * probabilities. Careful selection of table size allows lookup
+ * virtually in constant time.
+ *
+ * All cumulative probabilities are rescaled to LOOKUP_SCALE,
+ * this allows to save one multiplication operation on each iteration
+ * in randomize().
+ */
+
+void *
+fill_lookup(struct amino_acid **lookup, struct amino_acid *amino_acid, int amino_acid_size) {
+    float p = 0;
+    int i, j;
+
+    for (i = 0; i < amino_acid_size; i++) {
+        p += amino_acid[i].prob;
+        amino_acid[i].cprob_lookup = p*LOOKUP_SCALE;
+    }
+
+    // Prevent rounding error.
+    amino_acid[amino_acid_size - 1].cprob_lookup = LOOKUP_SIZE - 1;
+
+    for (i = 0, j = 0; i < LOOKUP_SIZE; i++) {
+        while (amino_acid[j].cprob_lookup < i) {
+            j++;
+        }
+        lookup[i] = &amino_acid[j];
+    }
+
+    return 0;
 }
 
-/* This function takes a pointer to the first element of an array */
-/* Each element of the array is a struct with a character and */
-/* a float number p between 0 and 1. */
-/* The function generates a random float number r and */
-/* finds the first array element such that p >= r. */
-/* This is a weighted random selection. */
-/* The function then prints the character of the array element. */
-/* This is done count times. */
-/* Between each WIDTH consecutive characters, the function prints a newline */
-static void random_fasta (aminoacid_t const *genelist, size_t count) {
-    do {
-	size_t line = MIN(WIDTH, count);
-	size_t pos = 0;
-	char buf[WIDTH + 1];
-	do {
-	    float r = myrandom (1.0);
-	    size_t i = 0;
-	    while (genelist[i].p < r)
-		++i; /* Linear search */
-	    buf[pos++] = genelist[i].c;
-	} while (pos < line);
-	buf[line] = '\n';
-	fwrite (buf, 1, line + 1, stdout);
-	count -= line;
-    } while (count);
+void
+randomize(struct amino_acid *amino_acid, int amino_acid_size,
+        const char *title, int n, random_t *rand) {
+    struct amino_acid *lookup[LOOKUP_SIZE];
+    char line_buffer[LINE_LEN + 1];
+    int i, j;
+
+    line_buffer[LINE_LEN] = '\n';
+
+    fill_lookup(lookup, amino_acid, amino_acid_size);
+
+    fputs_unlocked(title, stdout);
+
+    for (i = 0, j = 0; i < n; i++, j++) {
+        if (j == LINE_LEN) {
+            fwrite_unlocked(line_buffer, LINE_LEN + 1, 1, stdout);
+            j = 0;
+        }
+
+        float r = random_next_lookup(rand);
+        struct amino_acid *u = lookup[(short)r];
+        while (unlikely(u->cprob_lookup < r)) {
+            ++u;
+        }
+        line_buffer[j] = u->sym;
+    }
+    line_buffer[j] = '\n';
+    fwrite_unlocked(line_buffer, j + 1, 1, stdout);
 }
 
-int main (int argc, char **argv) {
-    size_t n;
-    if (argc > 1) {
-	char const *arg = argv[1];
- 	char *tail;
- 	n = strtoul (arg, &tail, 0);
- 	if (tail == arg)
-	    errx (1, "Could not convert \"%s\" to an unsigned long integer", arg);
-    } else n = 1000;
-
-    static aminoacid_t iub[] = {
-	{ 0.27, 'a' },
-	{ 0.12, 'c' },
-	{ 0.12, 'g' },
-	{ 0.27, 't' },
-	{ 0.02, 'B' },
-	{ 0.02, 'D' },
-	{ 0.02, 'H' },
-	{ 0.02, 'K' },
-	{ 0.02, 'M' },
-	{ 0.02, 'N' },
-	{ 0.02, 'R' },
-	{ 0.02, 'S' },
-	{ 0.02, 'V' },
-	{ 0.02, 'W' },
-	{ 0.02, 'Y' }};
-
-    static aminoacid_t homosapiens[] = {
-	{ 0.3029549426680, 'a' },
-	{ 0.1979883004921, 'c' },
-	{ 0.1975473066391, 'g' },
-	{ 0.3015094502008, 't' }};
-
-    accumulate_probabilities (iub, NELEMENTS(iub));
-    accumulate_probabilities (homosapiens, NELEMENTS(homosapiens));
-
-    static char const *const alu ="\
-GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG\
-GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA\
-CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT\
-ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA\
-GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG\
-AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC\
-AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA";
-
-    fputs (">ONE Homo sapiens alu\n", stdout);
-    repeat_fasta (alu, 2 * n);
-    fputs (">TWO IUB ambiguity codes\n", stdout);
-    random_fasta (iub, 3 * n);
-    fputs (">THREE Homo sapiens frequency\n", stdout);
-    random_fasta (homosapiens, 5 * n);
+struct amino_acid amino_acid[] = {
+   { 'a', 0.27 },
+   { 'c', 0.12 },
+   { 'g', 0.12 },
+   { 't', 0.27 },
+
+   { 'B', 0.02 },
+   { 'D', 0.02 },
+   { 'H', 0.02 },
+   { 'K', 0.02 },
+   { 'M', 0.02 },
+   { 'N', 0.02 },
+   { 'R', 0.02 },
+   { 'S', 0.02 },
+   { 'V', 0.02 },
+   { 'W', 0.02 },
+   { 'Y', 0.02 },
+};
+
+struct amino_acid homo_sapiens[] = {
+   { 'a', 0.3029549426680 },
+   { 'c', 0.1979883004921 },
+   { 'g', 0.1975473066391 },
+   { 't', 0.3015094502008 },
+};
+
+static const char alu[] =
+   "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTG"
+   "GGAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGA"
+   "GACCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAA"
+   "AATACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAAT"
+   "CCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAAC"
+   "CCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTG"
+   "CACTCCAGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA";
+
+int
+main(int argc, const char **argv) {
+    int n = argc > 1 ? atoi( argv[1] ) : 512;
+    random_t rand;
+
+    random_init(&rand);
+
+    repeat(alu, ">ONE Homo sapiens alu\n", n*2);
+    randomize(amino_acid, ARRAY_SIZE(amino_acid),
+            ">TWO IUB ambiguity codes\n", n*3, &rand);
+    randomize(homo_sapiens, ARRAY_SIZE(homo_sapiens),
+            ">THREE Homo sapiens frequency\n", n*5, &rand);
+
     return 0;
-}
+}
\ No newline at end of file
diff --git a/test/bench/fasta.go b/test/bench/fasta.go
index f79ff680fb..470bdb3285 100644
--- a/test/bench/fasta.go
+++ b/test/bench/fasta.go
@@ -31,135 +31,137 @@ POSSIBILITY OF SUCH DAMAGE.
  * http://shootout.alioth.debian.org/
  *
  * contributed by The Go Authors.
- * Based on C program by Joern Inge Vestgaarden
- * and Jorge Peixoto de Morais Neto.
+ * Based on C program by by Petr Prokhorenkov.
  */
 
 package main
 
 import (
-	"bufio"
+	"bytes"
 	"flag"
 	"os"
 )
 
-var out *bufio.Writer
+var out = make(buffer, 0, 32768)
 
 var n = flag.Int("n", 1000, "length of result")
 
-const WIDTH = 60 // Fold lines after WIDTH bytes
+const Line = 60
 
-func min(a, b int) int {
-	if a < b {
-		return a
+func Repeat(alu []byte, n int) {
+	buf := bytes.Add(alu, alu)
+	off := 0
+	for n > 0 {
+		m := n
+		if m > Line {
+			m = Line
+		}
+		buf1 := out.NextWrite(m + 1)
+		copy(buf1, buf[off:])
+		buf1[m] = '\n'
+		if off += m; off >= len(alu) {
+			off -= len(alu)
+		}
+		n -= m
 	}
-	return b
 }
 
-type AminoAcid struct {
-	p float
-	c byte
-}
+const (
+	IM = 139968
+	IA = 3877
+	IC = 29573
 
-func AccumulateProbabilities(genelist []AminoAcid) {
-	for i := 1; i < len(genelist); i++ {
-		genelist[i].p += genelist[i-1].p
-	}
+	LookupSize  = 4096
+	LookupScale float64 = LookupSize - 1
+)
+
+var rand uint32 = 42
+
+type Acid struct {
+	sym   byte
+	prob  float64
+	cprob float64
+	next  *Acid
 }
 
-// RepeatFasta prints the characters of the byte slice s. When it
-// reaches the end of the slice, it goes back to the beginning.
-// It stops after generating count characters.
-// After each WIDTH characters it prints a newline.
-// It assumes that WIDTH <= len(s) + 1.
-func RepeatFasta(s []byte, count int) {
-	pos := 0
-	s2 := make([]byte, len(s)+WIDTH)
-	copy(s2, s)
-	copy(s2[len(s):], s)
-	for count > 0 {
-		line := min(WIDTH, count)
-		out.Write(s2[pos : pos+line])
-		out.WriteByte('\n')
-		pos += line
-		if pos >= len(s) {
-			pos -= len(s)
+func computeLookup(acid []Acid) *[LookupSize]*Acid {
+	var lookup [LookupSize]*Acid
+	var p float64
+	for i := range acid {
+		p += acid[i].prob
+		acid[i].cprob = p * LookupScale
+		if i > 0 {
+			acid[i-1].next = &acid[i]
 		}
-		count -= line
 	}
-}
+	acid[len(acid)-1].cprob = 1.0 * LookupScale
 
-var lastrandom uint32 = 42
+	j := 0
+	for i := range lookup {
+		for acid[j].cprob < float64(i) {
+			j++
+		}
+		lookup[i] = &acid[j]
+	}
 
-const (
-	IM = 139968
-	IA = 3877
-	IC = 29573
-)
+	return &lookup
+}
 
-// Each element of genelist is a struct with a character and
-// a floating point number p between 0 and 1.
-// RandomFasta generates a random float r and
-// finds the first element such that p >= r.
-// This is a weighted random selection.
-// RandomFasta then prints the character of the array element.
-// This sequence is repeated count times.
-// Between each WIDTH consecutive characters, the function prints a newline.
-func RandomFasta(genelist []AminoAcid, count int) {
-	buf := make([]byte, WIDTH+1)
-	for count > 0 {
-		line := min(WIDTH, count)
-		for pos := 0; pos < line; pos++ {
-			lastrandom = (lastrandom*IA + IC) % IM
-			// Integer to float conversions are faster if the integer is signed.
-			r := float(int32(lastrandom)) / IM
-			for _, v := range genelist {
-				if v.p >= r {
-					buf[pos] = v.c
-					break
-				}
+func Random(acid []Acid, n int) {
+	lookup := computeLookup(acid)
+	for n > 0 {
+		m := n
+		if m > Line {
+			m = Line
+		}
+		buf := out.NextWrite(m + 1)
+		f := LookupScale / IM
+		myrand := rand
+		for i := 0; i < m; i++ {
+			myrand = (myrand*IA + IC) % IM
+			r := float64(int(myrand)) * f
+			a := lookup[int(r)]
+			for a.cprob < r {
+				a = a.next
 			}
+			buf[i] = a.sym
 		}
-		buf[line] = '\n'
-		out.Write(buf[0 : line+1])
-		count -= line
+		rand = myrand
+		buf[m] = '\n'
+		n -= m
 	}
 }
 
 func main() {
-	out = bufio.NewWriter(os.Stdout)
 	defer out.Flush()
 
 	flag.Parse()
 
-	iub := []AminoAcid{
-		AminoAcid{0.27, 'a'},
-		AminoAcid{0.12, 'c'},
-		AminoAcid{0.12, 'g'},
-		AminoAcid{0.27, 't'},
-		AminoAcid{0.02, 'B'},
-		AminoAcid{0.02, 'D'},
-		AminoAcid{0.02, 'H'},
-		AminoAcid{0.02, 'K'},
-		AminoAcid{0.02, 'M'},
-		AminoAcid{0.02, 'N'},
-		AminoAcid{0.02, 'R'},
-		AminoAcid{0.02, 'S'},
-		AminoAcid{0.02, 'V'},
-		AminoAcid{0.02, 'W'},
-		AminoAcid{0.02, 'Y'},
+	iub := []Acid{
+		Acid{prob: 0.27, sym: 'a'},
+		Acid{prob: 0.12, sym: 'c'},
+		Acid{prob: 0.12, sym: 'g'},
+		Acid{prob: 0.27, sym: 't'},
+		Acid{prob: 0.02, sym: 'B'},
+		Acid{prob: 0.02, sym: 'D'},
+		Acid{prob: 0.02, sym: 'H'},
+		Acid{prob: 0.02, sym: 'K'},
+		Acid{prob: 0.02, sym: 'M'},
+		Acid{prob: 0.02, sym: 'N'},
+		Acid{prob: 0.02, sym: 'R'},
+		Acid{prob: 0.02, sym: 'S'},
+		Acid{prob: 0.02, sym: 'V'},
+		Acid{prob: 0.02, sym: 'W'},
+		Acid{prob: 0.02, sym: 'Y'},
 	}
 
-	homosapiens := []AminoAcid{
-		AminoAcid{0.3029549426680, 'a'},
-		AminoAcid{0.1979883004921, 'c'},
-		AminoAcid{0.1975473066391, 'g'},
-		AminoAcid{0.3015094502008, 't'},
+	homosapiens := []Acid{
+		Acid{prob: 0.3029549426680, sym: 'a'},
+		Acid{prob: 0.1979883004921, sym: 'c'},
+		Acid{prob: 0.1975473066391, sym: 'g'},
+		Acid{prob: 0.3015094502008, sym: 't'},
 	}
 
-	AccumulateProbabilities(iub)
-	AccumulateProbabilities(homosapiens)
-
 	alu := []byte(
 		"GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG" +
 			"GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA" +
@@ -170,9 +172,38 @@ func main() {
 			"AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA")
 
 	out.WriteString(">ONE Homo sapiens alu\n")
-	RepeatFasta(alu, 2**n)
+	Repeat(alu, 2**n)
 	out.WriteString(">TWO IUB ambiguity codes\n")
-	RandomFasta(iub, 3**n)
+	Random(iub, 3**n)
 	out.WriteString(">THREE Homo sapiens frequency\n")
-	RandomFasta(homosapiens, 5**n)
+	Random(homosapiens, 5**n)
+}
+
+
+type buffer []byte
+
+func (b *buffer) Flush() {
+	p := *b
+	if len(p) > 0 {
+		os.Stdout.Write(p)
+	}
+	*b = p[0:0]
+}
+
+func (b *buffer) WriteString(s string) {
+	p := b.NextWrite(len(s))
+	for i := 0; i < len(s); i++ {
+		p[i] = s[i]
+	}
+}
+
+func (b *buffer) NextWrite(n int) []byte {
+	p := *b
+	if len(p)+n > cap(p) {
+		b.Flush()
+		p = *b
+	}
+	out := p[len(p) : len(p)+n]
+	*b = p[0 : len(p)+n]
+	return out
 }