text: /*  The Computer Language Benchmarks Game
   *  http://shootout.alioth.debian.org/
   * 
   *  contributed by Petr Prokhorenkov
   */
  
  #include <stdio.h>
  #include <stdlib.h>
  #include <string.h>
  
  //SEPXXX
  #define fputs_unlocked (void)
  #define fwrite_unlocked (void)
  #define fputc_unlocked (void)
  
  #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);
  }
  
  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;
      }
  }
  
  /*
   * 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;
  }
  
  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);
  }
  
  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;
  }
  
  