src/test/bench/shootout-fasta-redux.rs

   1 // The Computer Language Benchmarks Game
   2 // http://benchmarksgame.alioth.debian.org/
   3 //
   4 // contributed by the Rust Project Developers
   5
   6 // Copyright (c) 2013-2014 The Rust Project Developers
   7 //
   8 // All rights reserved.
   9 //
  10 // Redistribution and use in source and binary forms, with or without
  11 // modification, are permitted provided that the following conditions
  12 // are met:
  13 //
  14 // - Redistributions of source code must retain the above copyright
  15 //   notice, this list of conditions and the following disclaimer.
  16 //
  17 // - Redistributions in binary form must reproduce the above copyright
  18 //   notice, this list of conditions and the following disclaimer in
  19 //   the documentation and/or other materials provided with the
  20 //   distribution.
  21 //
  22 // - Neither the name of "The Computer Language Benchmarks Game" nor
  23 //   the name of "The Computer Language Shootout Benchmarks" nor the
  24 //   names of its contributors may be used to endorse or promote
  25 //   products derived from this software without specific prior
  26 //   written permission.
  27 //
  28 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  29 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  30 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  31 // FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
  32 // COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
  33 // INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
  34 // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  35 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  36 // HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
  37 // STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  38 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
  39 // OF THE POSSIBILITY OF SUCH DAMAGE.
  40
  41 use std::cmp::min;
  42 use std::io::{self, Write};
  43 use std::sync::{Arc, Mutex};
  44 use std::thread;
  45
  46
  47 const LINE_LEN: usize = 60;
  48
  49 const BLOCK_LINES: usize = 512;
  50 const BLOCK_THOROUGHPUT: usize = LINE_LEN * BLOCK_LINES;
  51 const BLOCK_LEN: usize = BLOCK_THOROUGHPUT + BLOCK_LINES;
  52
  53 const STDIN_BUF: usize = (LINE_LEN + 1) * 1024;
  54 const LOOKUP_SIZE: usize = 4 * 1024;
  55 const LOOKUP_SCALE: f32 = (LOOKUP_SIZE - 1) as f32;
  56
  57 const ALU: &'static [u8] =
  58     b"GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG\
  59       GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA\
  60       CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT\
  61       ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA\
  62       GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG\
  63       AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC\
  64       AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA";
  65
  66 const IUB: &'static [(u8, f32)] =
  67     &[(b'a', 0.27), (b'c', 0.12), (b'g', 0.12),
  68       (b't', 0.27), (b'B', 0.02), (b'D', 0.02),
  69       (b'H', 0.02), (b'K', 0.02), (b'M', 0.02),
  70       (b'N', 0.02), (b'R', 0.02), (b'S', 0.02),
  71       (b'V', 0.02), (b'W', 0.02), (b'Y', 0.02)];
  72
  73 const HOMOSAPIENS: &'static [(u8, f32)] =
  74     &[(b'a', 0.3029549426680),
  75       (b'c', 0.1979883004921),
  76       (b'g', 0.1975473066391),
  77       (b't', 0.3015094502008)];
  78
  79 // We need a specific Rng,
  80 // so implement this manually
  81
  82 const MODULUS: u32 = 139968;
  83 const MULTIPLIER: u32 = 3877;
  84 const ADDITIVE: u32 = 29573;
  85
  86 // Why doesn't rust already have this?
  87 // Algorithm directly taken from Wikipedia
  88 fn powmod(mut base: u64, mut exponent: u32, modulus: u64) -> u64 {
  89     let mut ret = 1;
  90     base %= modulus;
  91
  92     while exponent > 0 {
  93         if exponent & 1 == 1 {
  94            ret *= base;
  95            ret %= modulus;
  96         }
  97         exponent >>= 1;
  98         base *= base;
  99         base %= modulus;
 100     }
 101
 102     ret
 103 }
 104
 105 // Just a typical LCRNG
 106 pub struct Rng {
 107     last: u32
 108 }
 109
 110 impl Rng {
 111     pub fn new() -> Rng {
 112         Rng { last: 42 }
 113     }
 114
 115     pub fn max_value() -> u32 {
 116         MODULUS - 1
 117     }
 118
 119     pub fn normalize(p: f32) -> u32 {
 120         (p * MODULUS as f32).floor() as u32
 121     }
 122
 123     pub fn gen(&mut self) -> u32 {
 124         self.last = (self.last * MULTIPLIER + ADDITIVE) % MODULUS;
 125         self.last
 126     }
 127
 128     // This allows us to fast-forward the RNG,
 129     // allowing us to run it in parallel.
 130     pub fn future(&self, n: u32) -> Rng {
 131         let a = MULTIPLIER as u64;
 132         let b = ADDITIVE as u64;
 133         let m = MODULUS as u64;
 134
 135         //                          (a^n - 1) mod (a-1) m
 136         // x_k = ((a^n x_0 mod m) + --------------------- b) mod m
 137         //                                   a - 1
 138         //
 139         // Since (a - 1) divides (a^n - 1) mod (a-1) m,
 140         // the subtraction does not overflow and thus can be non-modular.
 141         //
 142         let new_seed =
 143             (powmod(a, n, m) * self.last as u64) % m +
 144             (powmod(a, n, (a-1) * m) - 1) / (a-1) * b;
 145
 146         Rng { last: (new_seed % m) as u32 }
 147     }
 148 }
 149
 150
 151 // This will end up keeping track of threads, like
 152 // in the other multithreaded Rust version, in
 153 // order to keep writes in order.
 154 //
 155 // This is stolen from another multithreaded Rust
 156 // implementation, although that implementation
 157 // was not able to parallelize the RNG itself.
 158 struct BlockSubmitter<W: io::Write> {
 159     writer: W,
 160     pub waiting_on: usize,
 161 }
 162
 163 impl<W: io::Write> BlockSubmitter<W> {
 164     fn submit(&mut self, data: &[u8], block_num: usize) -> Option<io::Result<()>> {
 165         if block_num == self.waiting_on {
 166             self.waiting_on += 1;
 167             Some(self.submit_async(data))
 168         }
 169         else {
 170             None
 171         }
 172     }
 173
 174     fn submit_async(&mut self, data: &[u8]) -> io::Result<()> {
 175         self.writer.write_all(data)
 176     }
 177 }
 178
 179
 180 // For repeating strings as output
 181 fn fasta_static<W: io::Write>(
 182     writer: &mut W,
 183     header: &[u8],
 184     data: &[u8],
 185     mut n: usize
 186 ) -> io::Result<()>
 187 {
 188     // The aim here is to print a short(ish) string cyclically
 189     // with line breaks as appropriate.
 190     //
 191     // The secret technique is to repeat the string such that
 192     // any wanted line is a single offset in the string.
 193     //
 194     // This technique is stolen from the Haskell version.
 195
 196     try!(writer.write_all(header));
 197
 198     // Maximum offset is data.len(),
 199     // Maximum read len is LINE_LEN
 200     let stream = data.iter().cloned().cycle();
 201     let mut extended: Vec<u8> = stream.take(data.len() + LINE_LEN + 1).collect();
 202
 203     let mut offset = 0;
 204     while n > 0 {
 205         let write_len = min(LINE_LEN, n);
 206         let end = offset + write_len;
 207         n -= write_len;
 208
 209         let tmp = extended[end];
 210         extended[end] = b'\n';
 211         try!(writer.write_all(&extended[offset..end + 1]));
 212         extended[end] = tmp;
 213
 214         offset = end;
 215         offset %= data.len();
 216     }
 217
 218     Ok(())
 219 }
 220
 221
 222 // For RNG streams as output
 223 fn fasta<W: io::Write + Send + 'static>(
 224     submitter: &Arc<Mutex<BlockSubmitter<W>>>,
 225     header: &[u8],
 226     table: &'static [(u8, f32)],
 227     rng: &mut Rng,
 228     n: usize
 229 ) -> io::Result<()>
 230 {
 231     // Here the lookup table is part of the algorithm and needs the
 232     // original probabilities (scaled with the LOOKUP_SCALE), because
 233     // Isaac says so :-)
 234     fn sum_and_scale(a: &'static [(u8, f32)]) -> Vec<(u8, f32)> {
 235         let mut p = 0f32;
 236         let mut result: Vec<(u8, f32)> = a.iter().map(|e| {
 237             p += e.1;
 238             (e.0, p * LOOKUP_SCALE)
 239         }).collect();
 240         let result_len = result.len();
 241         result[result_len - 1].1 = LOOKUP_SCALE;
 242         result
 243     }
 244
 245     fn make_lookup(a: &[(u8, f32)]) -> [(u8, f32); LOOKUP_SIZE] {
 246         let mut lookup = [(0, 0f32); LOOKUP_SIZE];
 247         let mut j = 0;
 248         for (i, slot) in lookup.iter_mut().enumerate() {
 249             while a[j].1 < (i as f32) {
 250                 j += 1;
 251             }
 252             *slot = a[j];
 253         }
 254         lookup
 255     }
 256
 257     {
 258         try!(submitter.lock().unwrap().submit_async(header));
 259     }
 260
 261     let lookup_table = Arc::new(make_lookup(&sum_and_scale(table)));
 262
 263     let thread_count = 4;
 264     let mut threads = Vec::new();
 265     for block_num in (0..thread_count) {
 266         let offset = BLOCK_THOROUGHPUT * block_num;
 267
 268         let local_submitter = submitter.clone();
 269         let local_lookup_table = lookup_table.clone();
 270         let local_rng = rng.future(offset as u32);
 271
 272         threads.push(thread::spawn(move || {
 273             gen_block(
 274                 local_submitter,
 275                 local_lookup_table,
 276                 local_rng,
 277                 n.saturating_sub(offset),
 278                 block_num,
 279                 thread_count
 280             )
 281         }));
 282     }
 283
 284     for thread in threads {
 285         try!(thread.join().unwrap());
 286     }
 287
 288     *rng = rng.future(n as u32);
 289
 290     Ok(())
 291 }
 292
 293 // A very optimized writer.
 294 // I have a feeling a simpler version wouldn't slow
 295 // things down too much, though, since the RNG
 296 // is the really heavy hitter.
 297 fn gen_block<W: io::Write>(
 298     submitter: Arc<Mutex<BlockSubmitter<W>>>,
 299     lookup_table: Arc<[(u8, f32)]>,
 300     mut rng: Rng,
 301     mut length: usize,
 302     mut block_num: usize,
 303     block_stride: usize,
 304 ) -> io::Result<()>
 305 {
 306     // Include newlines in block
 307     length += length / LINE_LEN;
 308     let block: &mut [u8] = &mut [b'\n'; BLOCK_LEN];
 309
 310     while length > 0 {
 311         {
 312             let gen_into = &mut block[..min(length, BLOCK_LEN)];
 313
 314             // Write random numbers, skipping newlines
 315             for (i, byte) in gen_into.iter_mut().enumerate() {
 316                 if (i + 1) % (LINE_LEN + 1) != 0 {
 317                     let p = rng.gen() as f32 * (LOOKUP_SCALE / MODULUS as f32);
 318                     *byte = lookup_table[p as usize..LOOKUP_SIZE].iter().find(
 319                         |le| le.1 >= p).unwrap().0;
 320                 }
 321             }
 322         }
 323
 324         let write_out = {
 325             if length >= BLOCK_LEN               { &mut *block }
 326             else if length % (LINE_LEN + 1) == 0 { &mut block[..length] }
 327             else                                 { &mut block[..length + 1] }
 328         };
 329
 330         *write_out.last_mut().unwrap() = b'\n';
 331         loop {
 332             // Make sure to release lock before calling `yield_now`
 333             let res = { submitter.lock().unwrap().submit(write_out, block_num) };
 334
 335             match res {
 336                 Some(result) => { try!(result); break; }
 337                 None => std::thread::yield_now()
 338             }
 339         }
 340         block_num += block_stride;
 341         rng = rng.future((BLOCK_THOROUGHPUT * (block_stride - 1)) as u32);
 342         length = length.saturating_sub(BLOCK_LEN * (block_stride - 1));
 343
 344         length = length.saturating_sub(BLOCK_LEN);
 345     }
 346
 347     Ok(())
 348 }
 349
 350 fn run<W: io::Write + Send + 'static>(writer: W) -> io::Result<()> {
 351     let n = std::env::args_os().nth(1)
 352         .and_then(|s| s.into_string().ok())
 353         .and_then(|n| n.parse().ok())
 354         .unwrap_or(1000);
 355
 356     let rng = &mut Rng::new();
 357
 358     // Use automatic buffering for the static version...
 359     let mut writer = io::BufWriter::with_capacity(STDIN_BUF, writer);
 360     try!(fasta_static(&mut writer, b">ONE Homo sapiens alu\n", ALU, n * 2));
 361
 362     // ...but the dynamic version does its own buffering already
 363     let writer = try!(writer.into_inner());
 364     let submitter = Arc::new(Mutex::new(BlockSubmitter { writer: writer, waiting_on: 0 }));
 365
 366     { submitter.lock().unwrap().waiting_on = 0; }
 367     try!(fasta(&submitter, b">TWO IUB ambiguity codes\n", &IUB, rng, n * 3));
 368     { submitter.lock().unwrap().waiting_on = 0; }
 369     try!(fasta(&submitter, b">THREE Homo sapiens frequency\n", &HOMOSAPIENS, rng, n * 5));
 370
 371     Ok(())
 372 }
 373
 374 fn main() {
 375     run(io::stdout()).unwrap()
 376 }