xorshiftで高速なランダムデータを出力

そこそこの精度の手頃なランダムデータが欲しい場合
xorshiftでランダムデータを吐き出すプログラムを用意しておくと便利

以下のプログラムは毎回異なるランダムデータを標準出力に吐き出し続けてくれる

c

/* xorshift.c */

#include <stdint.h>
#include <stdio.h>

#define N (1024*1024)

static uint32_t _data[N];

static uint32_t x = 123456789;
static uint32_t y = 362436069;
static uint32_t z = 521288629;
static uint32_t w = 88675123;

uint32_t xor128(void) {
  uint32_t t;

  t = x ^ (x << 11);
  x = y; y = z; z = w;
  return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
}

void init_xor128()
{
  FILE* fp;
  fp = fopen("/dev/urandom","r");
  fread(_data, sizeof(uint32_t),4,fp);
  fclose(fp);
  x = _data[0] ^ x;
  y = _data[1] ^ y;
  z = _data[2] ^ z;
  w = _data[3] ^ w;
}

int main()
{
  int i;

  init_xor128();

  while(1)
  {
    for(i=0;i<N;++i)
    {
      _data[i] = xor128();
    }
    fwrite(_data, sizeof(uint32_t), N, stdout);
  }
  return 0;
}

早速ビルドして実行してみると
/dev/urandom を利用した場合と比較して相当早くなっている事がわかる。

bash

#コンパイル
$ gcc -O3 -o qrand xorshift.c

#xorshiftを利用した場合
$ ./qrand | dd of=/dev/null bs=1M count=1000
996+4 records in
996+4 records out
1046478848 bytes (1.0 GB) copied, 1.87409 s, 558 MB/s

# /dev/urandomを利用した場合
$ dd if=/dev/urandom of=/dev/null bs=1M count=100
100+0 records in
100+0 records out
104857600 bytes (105 MB) copied, 9.46438 s, 11.1 MB/s

方式	速度
xorshift	558MB/s
/dev/urandom	11MB/s

HDDの初期化とかに使いやすいんじゃないでしょうか？

リロール

速度が気になる人は関数呼び出しコストとか気になるよね！
アンロールもやりたいよね！

/* xorshift.c */

#include <stdint.h>
#include <stdio.h>

#define N (1024*1024)

static uint32_t _data[N];

int main()
{
	int i,t;
	uint32_t x,y,z,w;

	FILE* fp;
	fp = fopen("/dev/urandom","r");
	fread(_data, sizeof(uint32_t),4,fp);
	fclose(fp);
	x = _data[0] ^ 123456789;
	y = _data[1] ^ 362436069;
	z = _data[2] ^ 521288629;
	w = _data[3] ^ 88675123;

	while(1)
	{
		for(i=0;i<N;i+=4)
		{
			t = x ^ (x << 11);
			x = y; y = z; z = w;
			w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
			_data[i+0] = w;

			t = x ^ (x << 11);
			x = y; y = z; z = w;
			w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
			_data[i+1] = w;

			t = x ^ (x << 11);
			x = y; y = z; z = w;
			w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
			_data[i+2] = w;

			t = x ^ (x << 11);
			x = y; y = z; z = w;
			w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
			_data[i+3] = w;
		}
		fwrite(_data, sizeof(uint32_t), N, stdout);
	}
	return 0;
}

$ gcc -O3 -o qrand xorshift.c

$ ./qrand | dd of=/dev/null bs=1M count=3000
1827+1173 records in
1827+1173 records out
2530217984 bytes (2.5 GB) copied, 3.06039 s, 827 MB/s

方式	速度
修正前	558MB/s
修正後	827MB/s

ちょっと早くなったよ！

64bit化

こちらの記事を参考に64bit化版を試しました
https://qiita.com/Nabetani/items/f8357e736f989633a2c0

#include <stdint.h>
#include <stdio.h>

static inline uint64_t rotl(const uint64_t x, int k) {
	return (x << k) | (x >> (64 - k));
}

static inline uint64_t next64(uint64_t* s){
	const uint64_t result = rotl(s[1] * 5, 7) * 9;

	const uint64_t t = s[1] << 17;

	s[2] ^= s[0];
	s[3] ^= s[1];
	s[1] ^= s[2];
	s[0] ^= s[3];

	s[2] ^= t;

	s[3] = rotl(s[3], 45);

	return result;
}

#define N (1024*8)
static uint64_t s[4];
static uint64_t _buf[N];

static void init(){
	FILE* fp;
	fp = fopen("/dev/urandom","r");
	size_t size = fread(s, sizeof(uint64_t), 4, fp);
	fclose(fp);
}

static inline void update_buffer(){
  for(int i=0;i<N;i+=8){
    _buf[i+0] = next64(s);
    _buf[i+1] = next64(s);
    _buf[i+2] = next64(s);
    _buf[i+3] = next64(s);
    _buf[i+4] = next64(s);
    _buf[i+5] = next64(s);
    _buf[i+6] = next64(s);
    _buf[i+7] = next64(s);
  }
}

void main()
{
  init();
  while(1){
    update_buffer();
    fwrite(_buf, sizeof(uint64_t), N, stdout);
  }
}

単純に倍速化した感じに

方式	速度
修正前	820MB/s
修正後	1.5GB/s