How to create a 3D random gradient out of 3 passed values in a fragment shader?

Question

I need to create an animated smoke-like texture. I can achieve this with the 3D perlin noise using gradients passed from the CPU side, which i did:

But on the current project I cannot pass an array from the normal cpp-code. I'm limited to only writing HLSL shaders (although, all the following stuff is written in GLSL as it's easier to set up). So I thought I need to generate some sort of random values for my gradients inside the fragment shader. While investigating how I can tackle this problem, I figured out that I can actually use hash functions as my pseudo random values. I'm following these articles (the first and the second), so I chose to use PCG hash for my purposes. I managed to generate decently looking value noise with the following code.

#version 420 core

#define MAX_TABLE_SIZE 256
#define MASK (MAX_TABLE_SIZE - 1)

in vec2 TexCoord;

out vec4 FragColor;

uniform float Time;

uint pcg_hash(uint input)
{
    uint state = input * 747796405u + 2891336453u;
    uint word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

// This is taken from here
// https://stackoverflow.com/a/17479300/9778826
float ConvertToFloat(uint n)
{
    uint ieeeMantissa = 0x007FFFFFu; // binary32 mantissa bitmask
    uint ieeeOne = 0x3F800000u;      // 1.0 in IEEE binary32

    n &= ieeeMantissa;
    n |= ieeeOne;

    float f = uintBitsToFloat(n);
    return f - 1.0;
}

float Random1(uint x, uint y)
{
    uint hash = pcg_hash(y ^ pcg_hash(x));
    return ConvertToFloat(hash);
}

float ValueNoise(vec2 p)
{
    int xi = int(p.x);
    uint rx0 = uint(xi & MASK);
    uint rx1 = uint((xi + 1) & MASK);

    int yi = int(p.y);
    uint ry0 = uint(yi & MASK);
    uint ry1 = uint((yi + 1) & MASK);

    float tx = p.x - float(xi);
    float ty = p.y - float(yi);

    float r00 = Random1(rx0, ry0);
    float r10 = Random1(rx1, ry0);
    float r01 = Random1(rx0, ry1);
    float r11 = Random1(rx1, ry1);

    float sx = smoothstep(0, 1, tx);
    float sy = smoothstep(0, 1, ty);

    float lerp0 = mix(r00, r10, sx);
    float lerp1 = mix(r01, r11, sx);

    return mix(lerp0, lerp1, sy);
}

float FractalNoise(vec2 point)
{
    float sum = 0.0;
    float frequency = 0.01;
    float amplitude = 1;
    int nLayers = 5;

    for (int i = 0; i < nLayers; i++)
    {
        float noise = ValueNoise(point * frequency) * amplitude * 0.5;
        sum += noise;
        amplitude *= 0.5;
        frequency *= 2.0;
    }

    return sum;
}

void main()
{
    // Coordinates go from 0.0 to 1.0 both horizontally and vertically
    vec2 Point = TexCoord * 2000;
    float noise = FractalNoise(Point);
    FragColor = vec4(noise, noise, noise, 1.0);
}

What I want, however, is to generate a 3D random gradient (which is actually just a 3D random vector) out of three arguments that I pass, to then feed it into the Perlin noise function. But I don't know how to do it properly. To clarify a bit about these three arguments: see, I need an animated Perlin noise, which means I will need a three component gradient at every joint of the 3D lattice. And the arguments are exactly x, y as well as the time variable but in the strict order. Say, a point (1, 4, 5) produces a gradient (0.1, 0.03, 0.78), but a point (4, 1, 5) should produce a completely different gradient, say, (0.22, 0.95, 0.43). So again, the order matters.

What I came up with (and what I could understand from the articles in question) is that I can hash the arguments sequentially and then use the resulting value as a seed to the same hash function which will now be working as a random number generator. So I wrote this function:

vec3 RandomGradient3(int x, int y, int z)
{
    uint seed = pcg_hash(z ^ pcg_hash(y ^ pcg_hash(x)));
    uint s1 = seed ^ pcg_hash(seed);
    uint s2 = s1 ^ pcg_hash(s1);
    uint s3 = s2 ^ pcg_hash(s2);

    float g1 = ConvertToFloat(s1);
    float g2 = ConvertToFloat(s2);
    float g3 = ConvertToFloat(s3);

    return vec3(g1, g2, g3);
}

And the gradient I then feed to the 3D perlin noise function:

float CalculatePerlin3D(vec2 p)
{
    float z = Time; // a uniform variable passed from the CPU side
    int xi0 = int(floor(p.x)) & MASK;
    int yi0 = int(floor(p.y)) & MASK;
    int zi0 = int(floor(z))   & MASK;

    int xi1 = (xi0 + 1) & MASK;
    int yi1 = (yi0 + 1) & MASK;
    int zi1 = (zi0 + 1) & MASK;

    float tx = p.x - int(floor(p.x));
    float ty = p.y - int(floor(p.y));
    float tz = z   - int(floor(z));

    float u = smoothstep(0, 1, tx);
    float v = smoothstep(0, 1, ty);
    float w = smoothstep(0, 1, tz);

    vec3 c000 = RandomGradient3(xi0, yi0, zi0);
    vec3 c100 = RandomGradient3(xi1, yi0, zi0);
    vec3 c010 = RandomGradient3(xi0, yi1, zi0);
    vec3 c110 = RandomGradient3(xi1, yi1, zi0);
    vec3 c001 = RandomGradient3(xi0, yi0, zi1);
    vec3 c101 = RandomGradient3(xi1, yi0, zi1);
    vec3 c011 = RandomGradient3(xi0, yi1, zi1);
    vec3 c111 = RandomGradient3(xi1, yi1, zi1);

    float x0 = tx, x1 = tx - 1;
    float y0 = ty, y1 = ty - 1;
    float z0 = tz, z1 = tz - 1;

    vec3 p000 = vec3(x0, y0, z0);
    vec3 p100 = vec3(x1, y0, z0);
    vec3 p010 = vec3(x0, y1, z0);
    vec3 p110 = vec3(x1, y1, z0);

    vec3 p001 = vec3(x0, y0, z1);
    vec3 p101 = vec3(x1, y0, z1);
    vec3 p011 = vec3(x0, y1, z1);
    vec3 p111 = vec3(x1, y1, z1);

    float a = mix(dot(c000, p000), dot(c100, p100), u);
    float b = mix(dot(c010, p010), dot(c110, p110), u);
    float c = mix(dot(c001, p001), dot(c101, p101), u);
    float d = mix(dot(c011, p011), dot(c111, p111), u);

    float e = mix(a, b, v);
    float f = mix(c, d, v);

    float noise = mix(e, f, w);
    float unsignedNoise = (noise + 1.0) / 2.0;

    return unsignedNoise;
}

With this RandomGradient3 function, the following noise texture is produced:

So the gradients seem to be correlated, hence the noise is not really random. The question is, how can I properly randomize these s1, s2 and s3 from RandomGradient3? I'm a real beginner in all this random numbers generating stuff and is certainly not a math guy.

The 3D perlin noise function, which I have, seems to be fine because if I feed it with predefined gradients from CPU it produces the expected result.