#pragma once

#include <array>
#include <cstdint>
#include <cmath>
#include <numeric>
#include <algorithm>

// Classic seeded Perlin noise.
class Noise {
public:
    explicit Noise(uint64_t seed = 0) {
        std::array<uint8_t, 256> arr{};
        std::iota(arr.begin(), arr.end(), 0);
        // Fisher-Yates with xorshift64
        uint64_t s = seed ^ 0x9e3779b97f4a7c15ULL;
        for (int i = 255; i > 0; --i) {
            s ^= s << 13; s ^= s >> 7; s ^= s << 17;
            int j = static_cast<int>(s % static_cast<uint64_t>(i + 1));
            std::swap(arr[i], arr[j]);
        }
        for (int i = 0; i < 256; ++i) p_[i] = p_[i + 256] = arr[i];
    }

    // Returns value in [-1, 1]
    float noise2(float x, float y) const {
        const int X = static_cast<int>(std::floor(x)) & 255;
        const int Y = static_cast<int>(std::floor(y)) & 255;
        x -= std::floor(x);
        y -= std::floor(y);
        const float u = fade(x), v = fade(y);
        const int A = p_[X] + Y, B = p_[X + 1] + Y;
        return lerp(
            lerp(grad(p_[A],     x,   y  ), grad(p_[B],   x-1, y  ), u),
            lerp(grad(p_[A+1],   x,   y-1), grad(p_[B+1], x-1, y-1), u),
            v);
    }

    // Fractal Brownian Motion, result in approximately [-1, 1]
    float octaves2(float x, float y, int n,
                   float persistence = 0.5f, float lacunarity = 2.0f) const {
        float val = 0, amp = 1, freq = 1, maxv = 0;
        for (int i = 0; i < n; ++i) {
            val  += noise2(x * freq, y * freq) * amp;
            maxv += amp;
            amp  *= persistence;
            freq *= lacunarity;
        }
        return val / maxv;
    }

    // 3-D Perlin noise, returns value in [-1, 1]
    float noise3(float x, float y, float z) const {
        const int X = static_cast<int>(std::floor(x)) & 255;
        const int Y = static_cast<int>(std::floor(y)) & 255;
        const int Z = static_cast<int>(std::floor(z)) & 255;
        x -= std::floor(x); y -= std::floor(y); z -= std::floor(z);
        const float u = fade(x), v = fade(y), w = fade(z);
        const int A  = p_[X]     + Y, AA = p_[A]     + Z, AB = p_[A + 1] + Z;
        const int B  = p_[X + 1] + Y, BA = p_[B]     + Z, BB = p_[B + 1] + Z;
        return lerp(
            lerp(lerp(grad3(p_[AA  ], x,   y,   z  ), grad3(p_[BA  ], x-1, y,   z  ), u),
                 lerp(grad3(p_[AB  ], x,   y-1, z  ), grad3(p_[BB  ], x-1, y-1, z  ), u), v),
            lerp(lerp(grad3(p_[AA+1], x,   y,   z-1), grad3(p_[BA+1], x-1, y,   z-1), u),
                 lerp(grad3(p_[AB+1], x,   y-1, z-1), grad3(p_[BB+1], x-1, y-1, z-1), u), v), w);
    }

    float octaves3(float x, float y, float z, int n,
                   float persistence = 0.5f, float lacunarity = 2.0f) const {
        float val = 0, amp = 1, freq = 1, maxv = 0;
        for (int i = 0; i < n; ++i) {
            val  += noise3(x * freq, y * freq, z * freq) * amp;
            maxv += amp;
            amp  *= persistence;
            freq *= lacunarity;
        }
        return val / maxv;
    }

private:
    std::array<int, 512> p_{};

    static float fade(float t) { return t*t*t*(t*(t*6.f-15.f)+10.f); }
    static float lerp(float a, float b, float t) { return a + t*(b-a); }
    static float grad(int h, float x, float y) {
        h &= 7;
        const float u = (h < 4) ? x : y;
        const float v = (h < 4) ? y : x;
        return ((h & 1) ? -u : u) + ((h & 2) ? -v : v);
    }
    static float grad3(int h, float x, float y, float z) {
        h &= 15;
        const float u = (h < 8) ? x : y;
        const float v = (h < 4) ? y : ((h == 12 || h == 14) ? x : z);
        return ((h & 1) ? -u : u) + ((h & 2) ? -v : v);
    }
};
