xorshift

アルゴリズムの勉強とかしてると、試しに書いてみたコードを走らせるためにたくさんのサンプルデータが必要になる。
で、標準関数のrand()を使って要素数1000個とか10000個とかの配列をでっち上げたりしていたわけですが、ここで問題にぶち当たりました。
C言語のrand()は0以上RAND_MAX以下の擬似乱数を発生させるわけですが、mingw-gccだとこのRAND_MAXの値が0x7fff(32767)しかありません。
これではもっと大きい数を扱いたい場合に色々と不便なので、「もっとよい乱数生成関数はないかいな」と探してみたところ、xorshiftというのがなかなか良さそうです。

参考：良い乱数・悪い乱数

なんでもやたら速いうえに周期が128bitもあるそうなので、uint64_tな範囲でも問題なく扱えそう。

というわけで書いてみました。

/* xorshift.h*/
#ifndef XOR_SHIFT_H
#define XOR_SHIFT_H

#include <stdint.h>

uint64_t *srand_xs64(void);
uint64_t rand_xs64(uint64_t *x);
void close_xs64(uint64_t *x);

#endif

/* xor128.c*/
#include <stdlib.h>
#include <stdint.h>
#include <time.h>
#include "xorshift.h"

uint64_t *srand_xs64(void)
{
    uint64_t *x = (uint64_t *)malloc(sizeof(uint64_t) * 4);
    if (!x)
        return NULL;

    uint64_t s;
    do {
        s = (uint64_t)(time(NULL));
    } while (!s);

    *x = (s << 32) | s;
    *(x + 1) = (*x << 8) | ((*x & 0xff00000000000000) >> 56);
    *(x + 2) = (*(x + 1) << 8) | ((*(x + 1) & 0xff00000000000000) >> 56);
    *(x + 3) = (*(x + 2) << 8) | ((*(x + 2) & 0xff00000000000000) >> 56);

    return x;
}

uint64_t rand_xs64(uint64_t *x)
{ 
    uint64_t t = (*x ^ (*x << 11));
    *x = *(x + 1);
    *(x + 1) = *(x + 2);
    *(x + 2) = *(x + 3);
    *(x + 3) = (*(x + 3) ^ (*(x + 3) >> 19)) ^ (t ^ (t >> 8));
    return *(x + 3);
}

void close_xs64(uint64_t *x)
{
    if (x) {
        free(x);
        x = NULL;
    }
}

#include <stdio.h>
#include <inttypes.h>
#include "xorshift.h"

int main(void)
{
    uint64_t *seed = srand_xs64();
    if (!seed)
        return -1;

    for (int i = 0; i < 100000)
        printf("%"PRIu64"\n", rand_xs64(seed));

    close_xs64(seed);
    return 0;
}

とりあえず10万個程度の乱数を発生させてみてもダブりがひとつも出ないし(当たり前なんだろうけど)、スピードも速いみたい。
素晴らしい!