#define V3TEST

/*******************************************************************************

Velox 3b

This is a high-quality non-cryptographic pseudorandom number generator.
It is pretty fast, simple and generates a fine stream.

Usage something like
    v3state ctx;
    v3srand(&ctx, time(0));
    randomness = v3rand(&ctx);

--------------------------------------------------------------------------------

Guaranteed period for 128-bit blocks is 2^128.[1]
Seeds can't begin to overlap before generating at least 2^128 blocks.[2]
The expected average period for 128-bit blocks is 2^255.[3]

[1] Proof that the minimum period for 128-bit output blocks is 2^128:

The structure of this algorithm is
    output(x = f(x + (ctr = ctr + 1)));

Let's say the output blocks go into an array out[].
If there is a cycle, then at points `a` and `b` in time,
    out[a] == out[b].
Now, the next blocks are as follows:
    out[a+1] == f(out[a]+a+1)
and
    out[b+1] == f(out[b]+b+1).
As f is a 128-bit bijective function, we get that a == b (mod 2^128). This can
only happen when the difference between a and b is a multiple of 2^128.

[2] Proof of non-overlapping of streams by different seeds:

Any particular fixed number of state update iterations gives a bijection of the
initial block. The seed is placed in the initial block. With two different
seeds, two streams can thus never generate the same 128-bit block at the same
stream position. If they generate the same 128-bit block at different positions,
the next block will necesarily be different unless the counters also equal,
which can only happen after generating at least 2^128 blocks, as explained in
the previous proof.

[3] Explanation of expected period 2^255:

The possibility for cycling happens every 2^128 blocks. 2^128 rounds of the
block mix is expected to be a pseudorandom permutation. The average period for
an n-bit (n not small) pseudorandom permutation is about 2^(n-1).

This gives expected total average period 2^128 * 2^(128-1) == 2^255.


2011-03-14 significant structural update
2010-11-23 corrected the "expected period" number
2009-11-20 guaranteed non-overlapping seeds
           added rough proofs of properties
2009-11-10 mainly cosmetic update
2009-11-09 initial release

This file is released into the public domain through CC0. No rights reserved.
Legal stuff: http://creativecommons.org/publicdomain/zero/1.0/

*******************************************************************************/

#include <stdint.h>

// -----------------------------------------------------------------------------

typedef struct { uint32_t v[4], ctr[4], pos; } v3state;

inline uint32_t v3rand(v3state *ctx);
void v3srand(v3state *ctx, uint32_t seed);

#define rol32(a, b) (((a)<<(b))|((a)>>(32-(b))))
inline uint32_t v3rand(v3state *ctx)
{
    if(ctx->pos == 0)
    {
        int i;

        ctx->v[0] = rol32(ctx->v[0] + ctx->v[3], 21);
        ctx->v[1] = rol32(ctx->v[1], 12) + ctx->v[2];
        ctx->v[2] = ctx->v[2] ^ ctx->v[0];
        ctx->v[3] = ctx->v[3] ^ ctx->v[1];
        ctx->v[0] = rol32(ctx->v[0] + ctx->v[3], 19);
        ctx->v[1] = rol32(ctx->v[1], 24) + ctx->v[2];
        ctx->v[2] = ctx->v[2] ^ ctx->v[0];
        ctx->v[3] = ctx->v[3] ^ ctx->v[1];
        ctx->v[0] = rol32(ctx->v[0] + ctx->v[3],  7);
        ctx->v[1] = rol32(ctx->v[1], 12) + ctx->v[2];
        ctx->v[2] = ctx->v[2] ^ ctx->v[0];
        ctx->v[3] = ctx->v[3] ^ ctx->v[1];
        ctx->v[0] = rol32(ctx->v[0] + ctx->v[3], 27);
        ctx->v[1] = rol32(ctx->v[1], 17) + ctx->v[2];
        ctx->v[2] = ctx->v[2] ^ ctx->v[0];
        ctx->v[3] = ctx->v[3] ^ ctx->v[1];

        for(i=0; i<4; i++) ctx->v[i] += ctx->ctr[i];
        for(i=0; i<4; i++) if(++ctx->ctr[i]) break;

        ctx->pos = 4;
    }
    return ctx->v[--ctx->pos];
}

void v3srand(v3state *ctx, uint32_t seed)
{
    int i;
    for(i=0; i<4; i++) ctx->v[i] = ctx->ctr[i] = i*0x9e3779b9;
    ctx->v[0] = seed;
    ctx->pos = 0;
    for(i=0; i<16; i++) v3rand(ctx);
}

#ifdef V3TEST

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

int main()
{
    unsigned int i;

    v3state ctx;
    v3srand(&ctx, 0);

    uint32_t x = 0;
    for(i=0; i<1<<20; i++)
        x ^= v3rand(&ctx);
    printf("Test value: %08X -- %s\n", x, x==0x3EE2E740 ? "ok!" : "FAILED!");

    /*
    FILE *k = fopen("v3-256M.bin", "wb");
    for(i=0; i<1024*1024*256/4; i++)
    {
        uint32_t num = v3rand(&ctx);
        fwrite(&num, 4, 1, k);
    }
    fclose(k);
    */

    return 0;
}

#endif