MRG32k3a

Javascript

  function MRG32k3a() {
    return (function(args) {
      // Copyright (c) 1998, 2002 Pierre L'Ecuyer, DIRO, Université de Montréal.
      // http://www.iro.umontreal.ca/~lecuyer/
      var m1 = 4294967087;
      var m2 = 4294944443;
      var s10 = 12345,
          s11 = 12345,
          s12 = 123,
          s20 = 12345,
          s21 = 12345,
          s22 = 123;
 
      if (args.length === 0) {
        args = [+new Date()];
      }
      var mash = Mash();
      for (var i = 0; i < args.length; i++) {
        s10 += mash(args[i]) * 0x100000000; // 2 ^ 32
        s11 += mash(args[i]) * 0x100000000;
        s12 += mash(args[i]) * 0x100000000;
        s20 += mash(args[i]) * 0x100000000;
        s21 += mash(args[i]) * 0x100000000;
        s22 += mash(args[i]) * 0x100000000;
      }
      s10 %= m1;
      s11 %= m1;
      s12 %= m1;
      s20 %= m2;
      s21 %= m2;
      s22 %= m2;
      mash = null;
 
      var uint32 = function() {
        var m1 = 4294967087;
        var m2 = 4294944443;
        var a12 = 1403580;
        var a13n = 810728;
        var a21 = 527612;
        var a23n = 1370589;
 
        var k, p1, p2;
 
        /* Component 1 */
        p1 = a12 * s11 - a13n * s10;
        k = p1 / m1 | 0;
        p1 -= k * m1;
        if (p1 < 0) p1 += m1;
        s10 = s11;
        s11 = s12;
        s12 = p1;
 
        /* Component 2 */
        p2 = a21 * s22 - a23n * s20;
        k = p2 / m2 | 0;
        p2 -= k * m2;
        if (p2 < 0) p2 += m2;
        s20 = s21;
        s21 = s22;
        s22 = p2;
 
        /* Combination */
        if (p1 <= p2) return p1 - p2 + m1;
        else return p1 - p2;
      };
 
      var random = function() {
        return uint32() * 2.3283064365386963e-10; // 2^-32
      };
      random.uint32 = uint32;
      random.fract53 = function() {
        return random() +
          (uint32() & 0x1fffff) * 1.1102230246251565e-16; // 2^-53
      };
      random.version = 'MRG32k3a 0.9';
      random.args = args;
 
      return random;
    } (Array.prototype.slice.call(arguments)));
  };

C

#include <stdio.h>
#include <stdint.h>
#include "mash.h"
 
// Copyright (c) 2002 Pierre L'Ecuyer, DIRO, Université de Montréal.
//  http://www.iro.umontreal.ca/~lecuyer/
#define norm 2.328306549295728e-10
#define m1    4294967087.0
#define m2    4294944443.0
#define a12      1403580.0
#define a13n      810728.0
#define a21       527612.0
#define a23n     1370589.0
 
static double s10 = 12345, s11 = 12345, s12 = 123,
  s20 = 12345, s21 = 12345, s22 = 123;
 
double MRG32k3a(void) {
  long k;
  double p1, p2;
  /* Component 1 */
  p1 = a12 * s11 - a13n * s10;
  k = p1 / m1;
  p1 -= k * m1;
  if (p1 < 0.0)
    p1 += m1;
  s10 = s11;
  s11 = s12;
  s12 = p1;
  /* Component 2 */
  p2 = a21 * s22 - a23n * s20;
  k = p2 / m2;
  p2 -= k * m2;
  if (p2 < 0.0)
    p2 += m2;
  s20 = s21;
  s21 = s22;
  s22 = p2;
  /* Combination */
  if (p1 <= p2)
    return ((p1 - p2 + m1) * norm);
  else
    return ((p1 - p2) * norm);
}
 
int main(int argc, unsigned char *argv[]) {
  double buffer[256];
  int i;
 
  for (i = 1; i < argc; i++) {
    s10 += mash(argv[i]) * 0x100000000ULL;   // 2^32
    s11 += mash(argv[i]) * 0x100000000ULL;   // 2^32
    s12 += mash(argv[i]) * 0x100000000ULL;   // 2^32
    s20 += mash(argv[i]) * 0x100000000ULL;   // 2^32
    s21 += mash(argv[i]) * 0x100000000ULL;   // 2^32
    s22 += mash(argv[i]) * 0x100000000ULL;   // 2^32
  }
 
  s10 = (uint64_t) s10 % (uint32_t) m1;
  s11 = (uint64_t) s11 % (uint32_t) m1;
  s12 = (uint64_t) s12 % (uint32_t) m1;
  s20 = (uint64_t) s20 % (uint32_t) m2;
  s21 = (uint64_t) s21 % (uint32_t) m2;
  s22 = (uint64_t) s22 % (uint32_t) m2;
 
  while (1) {
    for (i = 0; i < 256; i++) {
      buffer[i] = MRG32k3a();
    }
    fwrite(buffer, sizeof buffer, 1, stdout);
  }
}