gistfile1.c

#include <stdint.h>


#define LCG64(N, a, c) \
static uint64_t lcg_ ## N (uint64_t x) { \
  x *= (a); x += (c); \
  return x; \
}

LCG64(C1, 3935559000370003845, 2691343689449507681)
LCG64(C2, 3202034522624059733, 4354685564936845319)
LCG64(C3, 2862933555777941757, 7046029254386353087)

LCG64(knuth, 6364136223846793005, 1442695040888963407)


#define XORSHIFT64(N, a, b, c) \
static uint64_t xorshift_ ## N ## l (uint64_t x) { \
  x ^= x << (a); x ^= x >> (b); x ^= x << (c); \
  return x; \
} \
static uint64_t xorshift_ ## N ## r (uint64_t x) { \
  x ^= x >> (a); x ^= x << (b); x ^= x >> (c); \
  return x; \
}

XORSHIFT64(A1, 21, 35,  4)
XORSHIFT64(A2, 20, 41,  5)
XORSHIFT64(A3, 17, 31,  8)
XORSHIFT64(A4, 11, 29, 14)
XORSHIFT64(A5, 14, 29, 11)
XORSHIFT64(A6, 30, 35, 13)
XORSHIFT64(A7, 21, 37,  4)
XORSHIFT64(A8, 21, 43,  4)
XORSHIFT64(A9, 23, 41, 18)

/* A1-9 & C1-3 are parameters suggested by Numerical Recipes, 3rd ed. */

static uint64_t reverse(uint64_t x) {
  x = ((x >> 1) & 0x5555555555555555) | ((x & 0x5555555555555555) << 1);
  x = ((x >> 2) & 0x3333333333333333) | ((x & 0x3333333333333333) << 2);
  x = ((x >> 4) & 0x0F0F0F0F0F0F0F0F) | ((x & 0x0F0F0F0F0F0F0F0F) << 4);
  x = ((x >> 8) & 0x00FF00FF00FF00FF) | ((x & 0x00FF00FF00FF00FF) << 8);
  x = ((x >>16) & 0x0000FFFF0000FFFF) | ((x & 0x0000FFFF0000FFFF) <<16);
  return (x >> 32) | (x << 32);
}


#define COMMON_SEED 16807

#define RAW164(f, seed) \
static uint64_t roll_ ## f (void) { \
  static uint64_t state = (seed); \
  state = (f)(state); \
  return state; \
}

#define COMP264(f, g, seed) \
static uint64_t roll_ ## f ## _ ## g (void) { \
  static uint64_t state = (seed); \
  state = (f)(state); \
  return (g)(state); \
}

/* Generators used in GNU APL */

/* r415 */ RAW164(lcg_knuth, COMMON_SEED)
/* r439 */ COMP264(lcg_knuth, reverse, COMMON_SEED)

/* Baseline for comparison */

RAW164(lcg_C1, COMMON_SEED)
RAW164(lcg_C2, COMMON_SEED)
RAW164(lcg_C3, COMMON_SEED)

/* Proposed improvements -- klg */

#define KX(n) \
COMP264(lcg_knuth, xorshift_A ## n ## l, COMMON_SEED) \
COMP264(lcg_knuth, xorshift_A ## n ## r, COMMON_SEED)

KX(1)
KX(2)
KX(3)
KX(4)
KX(5)
KX(6)
KX(7)
KX(8)
KX(9)

#define KXN(n) \
  { "k+a" #n "l", roll_lcg_knuth_xorshift_A ## n ## l }, \
  { "k+a" #n "r", roll_lcg_knuth_xorshift_A ## n ## r }

static struct {
  char *name;
  uint64_t (*rng)(void);
} types[] = {
  KXN(1), KXN(2), KXN(3), KXN(4), KXN(5), KXN(6), KXN(7), KXN(8), KXN(9),
  { "lcg1", roll_lcg_C1 },
  { "lcg2", roll_lcg_C2 },
  { "lcg3", roll_lcg_C3 },
  { "knuth", roll_lcg_knuth },
  { "r415", roll_lcg_knuth },
  { "r439", roll_lcg_knuth_reverse },
};


/* Generate test data to pipe to dieharder */

#include <inttypes.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

uint32_t get_low(uint64_t x) {
  return 0xffffffff & x;
}

uint32_t get_high(uint64_t x) {
  return x >> 32;
}

int main(int argc, char **argv) {

  uint64_t (*roll)(void) = NULL;
  uint32_t (*u32)(uint64_t) = &get_low;
  int i = sizeof types / sizeof *types;

  if (argc < 2) {
    fprintf(stderr, "Usage: %s type [--low|--high]\n", *argv);
    fprintf(stderr, "Types:");
    while (i--) {
      fprintf(stderr, " %s%s", types[i].name, i ? "," : "\n");
    }
    return -1;
  }
  while (i--) {
    if (strcmp(argv[1], types[i].name) == 0) {
      roll = types[i].rng;
      break;
    }
  }
  if (roll == NULL) {
    fprintf(stderr, "Unknown type %s.\n", argv[1]);
    return -1;
  }

  if (argc > 2) {
    while (argv[2][0] == '-')
      argv[2]++;
    if (argv[2][0] == 'H' || argv[2][0] == 'h')
      u32 = &get_high;
  }
  fprintf(stderr, "Generating %s 32 bit words with %s\n",
      u32 == &get_high ? "high" : "low", argv[1]);

  while (1) {
    uint32_t x = u32(roll());
    fwrite(&x, sizeof x, 1, stdout);
  }

  return 0;
}