let
plus = bool32_sum;
bits = nat_to_bool32;
in
{
// The 'rand' functions hash one or more arbitrary numbers onto one or more
// uniformly distributed random numbers in the range 0 <= x < 1.
// These random numbers are meant to be high quality (no visual artifacts for
// typical graphics operations), scale independent, and the same results across
// different GPUs.

rand1x1 n = do
    // after Alexis Naveros
    local x = float_to_bool32 n;
    x := x `plus` (x `rshift` 11);
    x := x `xor` (x `lshift` 7);
    x := x `plus` (x `rshift` 15);
    x := x `xor` (x `lshift` 5);
    x := x `plus` (x `rshift` 12);
    x := x `xor` (x `lshift` 9);
    x := x `and` bits 0x007FFFFF `or` bits 0x3F800000;
    in bool32_to_float(x) - 1;

rand2x1 xy = do
    // after Alexis Naveros
    local [x,y] = float_to_bool32 xy;
    x := x `plus` (x `rshift` 11);
    x := x `xor` (x `lshift` 7);
    x := x `plus` y;
    x := x `xor` (x `lshift` 6);
    x := x `plus` (x `rshift` 15);
    x := x `xor` (x `lshift` 5);
    x := x `plus` (x `rshift` 12);
    x := x `xor` (x `lshift` 9);
    x := x `and` bits 0x007FFFFF `or` bits 0x3F800000;
    in bool32_to_float(x) - 1;

rand3x1 xyz = do
    // after Alexis Naveros
    local [x,y,z] = float_to_bool32 xyz;
    x := x `plus` (x `rshift` 11);
    x := x `xor` (x `lshift` 7);
    x := x `plus` y;
    x := x `xor` (x `lshift` 3);
    x := x `plus` (z `xor` (x `rshift` 14));
    x := x `xor` (x `lshift` 6);
    x := x `plus` (x `rshift` 15);
    x := x `xor` (x `lshift` 5);
    x := x `plus` (x `rshift` 12);
    x := x `xor` (x `lshift` 9);
    x := x `and` bits 0x007FFFFF `or` bits 0x3F800000;
    in bool32_to_float(x) - 1;

rand4x1 [a,b,c,d] =
    // TODO: Optimize, in the style of rand3x1. Or using SIMD.
    rand2x1[rand2x1[a,b],rand2x1[c,d]];

rand2x2 xy = do
    // Marc Olano, 2010, GPU Random Numbers via the Tiny Encryption Algorithm
    // https://www.csee.umbc.edu/~olano/papers/GPUTEA.pdf
    // https://umbcgaim.wordpress.com/2010/07/01/gpu-random-numbers/
    local v = float_to_bool32 xy;
    // two rounds of TEA hashing
    v.[X] := v.[X] `plus` xor[
                v.[Y] `lshift` 4 `plus` bits 0xA341316C,
                v.[Y] `plus` bits 0x9E3779B9,
                v.[Y] `rshift` 5 `plus` bits 0xC8013EA4
            ];
    v.[Y] := v.[Y] `plus` xor[
                v.[X] `lshift` 4 `plus` bits 0xAD90777D,
                v.[X] `plus` bits 0x9E3779B9,
                v.[X] `rshift` 5 `plus` bits 0x7E95761E
            ];
    v.[X] := v.[X] `plus` xor[
                v.[Y] `lshift` 4 `plus` bits 0xA341316C,
                v.[Y] `plus` bits 0x3C6EF372,
                v.[Y] `rshift` 5 `plus` bits 0xC8013EA4
            ];
    v.[Y] := v.[Y] `plus` xor[
                v.[X] `lshift` 4 `plus` bits 0xAD90777D,
                v.[X] `plus` bits 0x3C6EF372,
                v.[X] `rshift` 5 `plus` bits 0x7E95761E
            ];
    // Normalize to the range 1.0 <= v < 2.0:
    v := v `and` [bits 0x007FFFFF,  // strip sign bit and exponent
                  bits 0x007FFFFF]
           `or` [bits 0x3F800000,   // set sign and exponent to 0
                 bits 0x3F800000];
    in bool32_to_float(v) - 1;

}
