Commit f8f83e80 authored by Russ Cox's avatar Russ Cox

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
parent 7d7ebd2f
......@@ -28,146 +28,190 @@ POSSIBILITY OF SUCH DAMAGE.
*/
/*
* http://shootout.alioth.debian.org/u32q/benchmark.php?test=fasta&lang=gcc&id=4
*/
* 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 Joern Inge Vestgaarden
* Modified by Jorge Peixoto de Morais Neto
*
* 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;
}
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;
// 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);
}
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
......@@ -31,134 +31,136 @@ 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
}
return b
}
type AminoAcid struct {
p float
c byte
}
func AccumulateProbabilities(genelist []AminoAcid) {
for i := 1; i < len(genelist); i++ {
genelist[i].p += genelist[i-1].p
buf1 := out.NextWrite(m + 1)
copy(buf1, buf[off:])
buf1[m] = '\n'
if off += m; off >= len(alu) {
off -= len(alu)
}
}
// 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)
}
count -= line
n -= m
}
}
var lastrandom uint32 = 42
const (
IM = 139968
IA = 3877
IC = 29573
LookupSize = 4096
LookupScale float64 = LookupSize - 1
)
// 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
}
}
}
buf[line] = '\n'
out.Write(buf[0 : line+1])
count -= line
var rand uint32 = 42
type Acid struct {
sym byte
prob float64
cprob float64
next *Acid
}
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]
}
}
acid[len(acid)-1].cprob = 1.0 * LookupScale
j := 0
for i := range lookup {
for acid[j].cprob < float64(i) {
j++
}
lookup[i] = &acid[j]
}
return &lookup
}
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
}
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'},
}
homosapiens := []AminoAcid{
AminoAcid{0.3029549426680, 'a'},
AminoAcid{0.1979883004921, 'c'},
AminoAcid{0.1975473066391, 'g'},
AminoAcid{0.3015094502008, 't'},
}
AccumulateProbabilities(iub)
AccumulateProbabilities(homosapiens)
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 := []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'},
}
alu := []byte(
"GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG" +
......@@ -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
}
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment