1/*
2Redistribution and use in source and binary forms, with or without
3modification, are permitted provided that the following conditions are met:
4
5    * Redistributions of source code must retain the above copyright
6    notice, this list of conditions and the following disclaimer.
7
8    * Redistributions in binary form must reproduce the above copyright
9    notice, this list of conditions and the following disclaimer in the
10    documentation and/or other materials provided with the distribution.
11
12    * Neither the name of "The Computer Language Benchmarks Game" nor the
13    name of "The Computer Language Shootout Benchmarks" nor the names of
14    its contributors may be used to endorse or promote products derived
15    from this software without specific prior written permission.
16
17THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27POSSIBILITY OF SUCH DAMAGE.
28*/
29
30/* The Computer Language Benchmarks Game
31 * http://shootout.alioth.debian.org/
32 *
33 * contributed by The Go Authors.
34 * Based on C program by by Petr Prokhorenkov.
35 */
36
37package main
38
39import (
40	"flag"
41	"os"
42)
43
44var out = make(buffer, 0, 32768)
45
46var n = flag.Int("n", 1000, "length of result")
47
48const Line = 60
49
50func Repeat(alu []byte, n int) {
51	buf := append(alu, alu...)
52	off := 0
53	for n > 0 {
54		m := n
55		if m > Line {
56			m = Line
57		}
58		buf1 := out.NextWrite(m + 1)
59		copy(buf1, buf[off:])
60		buf1[m] = '\n'
61		if off += m; off >= len(alu) {
62			off -= len(alu)
63		}
64		n -= m
65	}
66}
67
68const (
69	IM = 139968
70	IA = 3877
71	IC = 29573
72
73	LookupSize          = 4096
74	LookupScale float64 = LookupSize - 1
75)
76
77var rand uint32 = 42
78
79type Acid struct {
80	sym   byte
81	prob  float64
82	cprob float64
83	next  *Acid
84}
85
86func computeLookup(acid []Acid) *[LookupSize]*Acid {
87	var lookup [LookupSize]*Acid
88	var p float64
89	for i := range acid {
90		p += acid[i].prob
91		acid[i].cprob = p * LookupScale
92		if i > 0 {
93			acid[i-1].next = &acid[i]
94		}
95	}
96	acid[len(acid)-1].cprob = 1.0 * LookupScale
97
98	j := 0
99	for i := range lookup {
100		for acid[j].cprob < float64(i) {
101			j++
102		}
103		lookup[i] = &acid[j]
104	}
105
106	return &lookup
107}
108
109func Random(acid []Acid, n int) {
110	lookup := computeLookup(acid)
111	for n > 0 {
112		m := n
113		if m > Line {
114			m = Line
115		}
116		buf := out.NextWrite(m + 1)
117		f := LookupScale / IM
118		myrand := rand
119		for i := 0; i < m; i++ {
120			myrand = (myrand*IA + IC) % IM
121			r := float64(int(myrand)) * f
122			a := lookup[int(r)]
123			for a.cprob < r {
124				a = a.next
125			}
126			buf[i] = a.sym
127		}
128		rand = myrand
129		buf[m] = '\n'
130		n -= m
131	}
132}
133
134func main() {
135	defer out.Flush()
136
137	flag.Parse()
138
139	iub := []Acid{
140		{prob: 0.27, sym: 'a'},
141		{prob: 0.12, sym: 'c'},
142		{prob: 0.12, sym: 'g'},
143		{prob: 0.27, sym: 't'},
144		{prob: 0.02, sym: 'B'},
145		{prob: 0.02, sym: 'D'},
146		{prob: 0.02, sym: 'H'},
147		{prob: 0.02, sym: 'K'},
148		{prob: 0.02, sym: 'M'},
149		{prob: 0.02, sym: 'N'},
150		{prob: 0.02, sym: 'R'},
151		{prob: 0.02, sym: 'S'},
152		{prob: 0.02, sym: 'V'},
153		{prob: 0.02, sym: 'W'},
154		{prob: 0.02, sym: 'Y'},
155	}
156
157	homosapiens := []Acid{
158		{prob: 0.3029549426680, sym: 'a'},
159		{prob: 0.1979883004921, sym: 'c'},
160		{prob: 0.1975473066391, sym: 'g'},
161		{prob: 0.3015094502008, sym: 't'},
162	}
163
164	alu := []byte(
165		"GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG" +
166			"GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA" +
167			"CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT" +
168			"ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA" +
169			"GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG" +
170			"AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC" +
171			"AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA")
172
173	out.WriteString(">ONE Homo sapiens alu\n")
174	Repeat(alu, 2**n)
175	out.WriteString(">TWO IUB ambiguity codes\n")
176	Random(iub, 3**n)
177	out.WriteString(">THREE Homo sapiens frequency\n")
178	Random(homosapiens, 5**n)
179}
180
181type buffer []byte
182
183func (b *buffer) Flush() {
184	p := *b
185	if len(p) > 0 {
186		os.Stdout.Write(p)
187	}
188	*b = p[0:0]
189}
190
191func (b *buffer) WriteString(s string) {
192	p := b.NextWrite(len(s))
193	copy(p, s)
194}
195
196func (b *buffer) NextWrite(n int) []byte {
197	p := *b
198	if len(p)+n > cap(p) {
199		b.Flush()
200		p = *b
201	}
202	out := p[len(p) : len(p)+n]
203	*b = p[:len(p)+n]
204	return out
205}
206