3 * http://shootout.alioth.debian.org/
15 double gen_random (double max) {
16 static long last = 42;
17 return max * (last = (last * IA + IC) % IM) / IM;
25 /* Weighted selection from alphabet */
27 void makeCumulative (struct aminoacids * genelist, int count) {
31 for (i=0; i < count; i++) {
37 char selectRandom (const struct aminoacids * genelist, int count) {
38 double r = gen_random (1);
41 if (r < genelist[0].p) return genelist[0].c;
48 if (r < genelist[i].p) hi = i; else lo = i;
50 return genelist[hi].c;
53 /* Generate and write FASTA format */
55 #define LINE_LENGTH (60)
57 void makeRandomFasta (const char * id, const char * desc, const struct
58 aminoacids * genelist, int count, int n) {
62 printf (">%s %s\n", id, desc);
64 for (; todo > 0; todo -= LINE_LENGTH) {
65 char pick[LINE_LENGTH+1];
66 if (todo < LINE_LENGTH) m = todo; else m = LINE_LENGTH;
67 for (i=0; i < m; i++) pick[i] = selectRandom (genelist, count);
73 void makeRepeatFasta (const char * id, const char * desc, const char *
76 int todo = n, k = 0, kn = strlen (s);
79 ss = (char *) malloc (kn + 1);
82 printf (">%s %s\n", id, desc);
84 for (; todo > 0; todo -= LINE_LENGTH) {
85 if (todo < LINE_LENGTH) m = todo; else m = LINE_LENGTH;
102 /* Main -- define alphabets, make 3 fragments */
104 struct aminoacids iub[] = {
123 #define IUB_LEN (sizeof (iub) / sizeof (struct aminoacids))
125 struct aminoacids homosapiens[] = {
126 { 'a', 0.3029549426680 },
127 { 'c', 0.1979883004921 },
128 { 'g', 0.1975473066391 },
129 { 't', 0.3015094502008 },
132 #define HOMOSAPIENS_LEN (sizeof (homosapiens) / sizeof (struct aminoacids))
135 "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG"
136 "GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA"
137 "CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT"
138 "ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA"
139 "GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG"
140 "AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC"
141 "AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA";
143 int main (int argc, char * argv[]) {
146 if (argc > 1) sscanf (argv[1], "%d", &n);
148 makeCumulative (iub, IUB_LEN);
149 makeCumulative (homosapiens, HOMOSAPIENS_LEN);
151 makeRepeatFasta ("ONE", "Homo sapiens alu", alu, n*2);
152 makeRandomFasta ("TWO", "IUB ambiguity codes", iub, IUB_LEN, n*3);
153 makeRandomFasta ("THREE", "Homo sapiens frequency", homosapiens,
154 HOMOSAPIENS_LEN, n*5);