*/
/*
- * 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
* 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" +
"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
}