#ifndef NOISELIB_H
#define NOISELIB_H

#include <stdint.h>
#include <string.h>
#include <math.h>
#include "xoroLib.h"
#include "md5.h"

//typedefs
//a single raw noisemap
typedef struct
{
    uint8_t lookupHash[256 + 1]; //shuffled bytes
    float xOffset, yOffset, zOffset; //offsets of the current noisemap. x y and z
    float amplitude; //amplitude of the noisemap
    float lacunarity; //frequency of the noisemap
} PerlinNoise;

//a bunch of stacked noisemaps
typedef struct
{
    int octcnt; //number of octaves
    PerlinNoise octaves[4];
} OctaveNoise;

typedef struct
{
    //name of the noisemap, including the minecraft:* part
    //this is used for the xor constants
    char* name;

    //these can be derived directly from name via md5.
    //we will be hardcoding them for now, inside main
    uint64_t lowMask;
    uint64_t highMask;

    //first octave is octaveNinimum
    int firstOctave;

    //this is the first amplitude. pretty sure everything is just scaled by this. seems
    //to be mostly useless, as it is 1 except for one place
    float firstAmplitude;

    //amplitudes
    int ampCount;
    float* amplitudes;
} NoiseParameters;

//two octave noisemaps that get stacked
typedef struct
{
    float amplitude;
    OctaveNoise octA;
    OctaveNoise octB;
} DoublePerlinNoise;

//==================HELPER_FUNCTIONS==================

float lerp(float part, float from, float to)
{
    return from + part * (to - from);
}

//modify NoiseParameters to have a valid mask
void md5NoiseCalc(NoiseParameters* np)
{
    np->highMask = 0;
    np->lowMask = 0;

    uint8_t hash[16];
    MD5_CTX ctx;
    md5_init(&ctx);
    md5_update(&ctx, np->name, strlen(np->name));
    md5_final(&ctx, hash);

    for(int i = 0; i < 8; i++)
    {
        np->lowMask |= (uint64_t)hash[i] << (8 * (7 - i));
        np->highMask |= (uint64_t)hash[i + 8] << (8 * (7 - i));
    }
    printf("%.16lx %.16lx\n", np->lowMask, np->highMask);
}

#define FASTLERP 0

#if !FASTLERP
    const float LERP_LOOKUP[16][3] = {
        { 1.0f,  1.0f,  0.0f}, //  0:  a + b
        {-1.0f,  1.0f,  0.0f}, //  1: -a + b
        { 1.0f, -1.0f,  0.0f}, //  2:  a - b
        {-1.0f, -1.0f,  0.0f}, //  3: -a - b
        { 1.0f,  0.0f,  1.0f}, //  4:  a + c
        {-1.0f,  0.0f,  1.0f}, //  5: -a + c
        { 1.0f,  0.0f, -1.0f}, //  6:  a - c
        {-1.0f,  0.0f, -1.0f}, //  7: -a - c
        { 0.0f,  1.0f,  1.0f}, //  8:  b + c
        { 0.0f, -1.0f,  1.0f}, //  9: -b + c
        { 0.0f,  1.0f, -1.0f}, // 10:  b - c
        { 0.0f, -1.0f, -1.0f}, // 11: -b - c
        { 1.0f,  1.0f,  0.0f}, // 12:  a + b 
        { 0.0f, -1.0f,  1.0f}, // 13: -b + c 
        {-1.0f,  1.0f,  0.0f}, // 14: -a + b 
        { 0.0f, -1.0f, -1.0f}  // 15: -b - c 
    };

    float indexedLerp(uint8_t idx, float a, float b, float c)
    {
        idx &= 0xF;
        const float* m = LERP_LOOKUP[idx];

        return (a * m[0]) + (b * m[1]) + (c * m[2]);
    }
#endif

#if FASTLERP
    float indexedLerp(uint8_t idx, float a, float b, float c)
    {
        idx &= 0xF;

        float sum = 0.0;

        const uint32_t bitFieldA = 0xbf5550ff;
        const uint32_t bitFieldB = 0x55f3ff0f;
        const uint32_t bitFieldC = 0x733faff0;

        sum += a * ((bitFieldA >> idx) & 1) * ((int32_t)(2 * ((bitFieldA >> (idx + 16)) & 1)) - 1);
        sum += b * ((bitFieldB >> idx) & 1) * ((int32_t)(2 * ((bitFieldB >> (idx + 16)) & 1)) - 1);
        sum += c * ((bitFieldC >> idx) & 1) * ((int32_t)(2 * ((bitFieldC >> (idx + 16)) & 1)) - 1);

        return sum;
    }
#endif

//===================START_INIT=================

//generate a randomly shuffled vector, and x y and z offsets
void xPerlinInit(PerlinNoise *noise, Xoroshiro *xr)
{
    int i = 0;
    noise->xOffset = xNextDouble(xr) * 256.0;

    //the y offset is a randomly generated number between 0 and 256
    //this controls the height of the noisemap, and also makes it 3d
    noise->yOffset = xNextDouble(xr) * 256.0;
    noise->zOffset = xNextDouble(xr) * 256.0;

    noise->amplitude = 1.0;
    noise->lacunarity = 1.0;

    //get the pointer to noise values so we can easily work with them
    uint8_t *lookupHash = noise->lookupHash;

    //init 256 identity permutation
    for (i = 0; i < 256; i++)
    {
        lookupHash[i] = i;
    }

    //shuffle
    for (i = 0; i < 256; i++)
    {
        int j = xNextInt(xr, 256 - i) + i;
        uint8_t n = lookupHash[i];
        lookupHash[i] = lookupHash[j];
        lookupHash[j] = n;
    }

    lookupHash[256 + 1] = lookupHash[0];
}

int xOctaveInit(OctaveNoise *noise, Xoroshiro *xr, PerlinNoise *octaves, NoiseParameters* np, int maxOctave)
{
    //this function initializes stacked octaves

    //bunch of constants that shall be xor'd with xrng state
    const uint64_t md5_octave_n[][2] = {
        {0xb198de63a8012672, 0x7b84cad43ef7b5a8}, // md5 "octave_-12"
        {0x0fd787bfbc403ec3, 0x74a4a31ca21b48b8}, // md5 "octave_-11"
        {0x36d326eed40efeb2, 0x5be9ce18223c636a}, // md5 "octave_-10"
        {0x082fe255f8be6631, 0x4e96119e22dedc81}, // md5 "octave_-9"
        {0x0ef68ec68504005e, 0x48b6bf93a2789640}, // md5 "octave_-8"
        {0xf11268128982754f, 0x257a1d670430b0aa}, // md5 "octave_-7"
        {0xe51c98ce7d1de664, 0x5f9478a733040c45}, // md5 "octave_-6"
        {0x6d7b49e7e429850a, 0x2e3063c622a24777}, // md5 "octave_-5"
        {0xbd90d5377ba1b762, 0xc07317d419a7548d}, // md5 "octave_-4"
        {0x53d39c6752dac858, 0xbcd1c5a80ab65b3e}, // md5 "octave_-3"
        {0xb4a24d7a84e7677b, 0x023ff9668e89b5c4}, // md5 "octave_-2"
        {0xdffa22b534c5f608, 0xb9b67517d3665ca9}, // md5 "octave_-1"
        {0xd50708086cef4d7c, 0x6e1651ecc7f43309}, // md5 "octave_0"
    };

    //scalings
    const float lacuna_ini[] = { // -minimumOctave = 3..12
        1, .5, .25, 1./8, 1./16, 1./32, 1./64, 1./128, 1./256, 1./512, 1./1024,
        1./2048, 1./4096,
    };

    static const double persist_ini[] = { // len = 4..9
        0, 1, 2./3, 4./7, 8./15, 16./31, 32./63, 64./127, 128./255, 256./511,
    };

    float lacuna = lacuna_ini[-np->firstOctave];
    float persist = persist_ini[np->ampCount];

    uint64_t xlo = xNextLong(xr);
    uint64_t xhi = xNextLong(xr);
    int i;

    //static const float amplitudes[] = {1.0, 1.0, 1.0, 1.0};

    for (i = 0; i < np->ampCount && i != maxOctave; i++)
    {
        Xoroshiro pxr;
        pxr.lo = xlo ^ md5_octave_n[12 + np->firstOctave + i][0];
        pxr.hi = xhi ^ md5_octave_n[12 + np->firstOctave + i][1];
        xPerlinInit(&octaves[i], &pxr);
        octaves[i].amplitude = np->amplitudes[i] * persist;
        octaves[i].lacunarity = lacuna;

        lacuna *= 2.0;
        persist *= 0.5;
    }

    //copy local octaves

    noise->octcnt = i;
    for(i = 0; i < noise->octcnt; i++) noise->octaves[i] = octaves[i];
    return i;
}

void xDoublePerlinInit(DoublePerlinNoise *noise, Xoroshiro *xr, NoiseParameters* np, int maxOctave)
{
    //this splits the octaves so na + nb = maxOctave BUT na gets more if its odd.
    //this is so we evenly spread the octaves
    int na = (maxOctave + 1) >> 1;
    int nb = maxOctave - na;

    //init the first octave noise in float perlin noise. na is octave count
    //append octaves initalized to octaves
    xOctaveInit(&noise->octA, xr, noise->octA.octaves, np, na);

    //second octave noise
    xOctaveInit(&noise->octB, xr, noise->octB.octaves, np, nb);

    const float amp_ini[] = { // (5 ./ 3) * len / (len + 1), len = 2..9
        0, 5./6, 10./9, 15./12, 20./15, 25./18, 30./21, 35./24, 40./27, 45./30,
    };

    //used because 4 octaves
    noise->amplitude = amp_ini[np->ampCount];
}

void initClimateSeed(
    DoublePerlinNoise* dpn,
    NoiseParameters* np, uint64_t seed, int maxOctave
    )
{
    Xoroshiro pxr;
    xSetSeed(&pxr, seed);
    uint64_t xlo = xNextLong(&pxr);
    uint64_t xhi = xNextLong(&pxr);

    
    // md5 "minecraft:calcite"
    pxr.lo = xlo ^ np->lowMask;
    pxr.hi = xhi ^ np->highMask;

    //md5 minecraft:powder_snow
    //pxr.lo = xlo ^ 0x7520c0b89fe4562f;
    //pxr.hi = xhi ^ 0x62ed8c50c818845d;

    xDoublePerlinInit(dpn, &pxr, np, maxOctave);
}

//===========================END_INIT=====================

//===========================SAMPLING==============================

//sample single perlin noisemap. this takes a noisemap and x z coordinates
//within the range 0.0 .. 256.0 on the y level 0
float samplePerlin(const PerlinNoise *noise, float x, float y, float z)
{
    uint8_t xInt, yInt, zInt;
    float t1, t2, t3;

    //shift the noisemap by xoffset and zoffset
    x += noise->xOffset;
    y += noise->yOffset;
    z += noise->zOffset;

    //get the integer cell we're in
    float xIntTemp = floor(x);
    float yIntTemp = floor(y);
    float zIntTemp = floor(z);

    //restrict x and z to be only the fractional part. 0.0 - 1.0
    x -= xIntTemp;
    y -= yIntTemp;
    z -= zIntTemp;
    xInt = (int) xIntTemp;
    yInt = (int) yIntTemp;
    zInt = (int) zIntTemp;

    //smoothstep function. potentially will remove
    t1 = x*x*x * (x * (x*6.0-15.0) + 10.0);
    t2 = y*y*y * (y * (y*6.0-15.0) + 10.0);
    t3 = z*z*z * (z * (z*6.0-15.0) + 10.0);

    const uint8_t *lookupHash = noise->lookupHash;

    //get psuedorandom values from a lookup table, to "randomize"
    //the 8 corners of the unit cube we're in

    //it looks like this because its optimized
    uint8_t a1 = lookupHash[xInt]   + yInt;
    uint8_t b1 = lookupHash[(uint8_t)xInt+1] + yInt;

    uint8_t a2 = lookupHash[a1]   + zInt;
    uint8_t b2 = lookupHash[b1]   + zInt;
    uint8_t a3 = lookupHash[(uint8_t)a1+1] + zInt;
    uint8_t b3 = lookupHash[(uint8_t)b1+1] + zInt;

    //computes the gradients, and calculate the dot products across all 8 vectors
    float l1 = indexedLerp(lookupHash[a2],   x,   y,   z);
    float l2 = indexedLerp(lookupHash[b2],   x-1, y,   z);
    float l3 = indexedLerp(lookupHash[a3],   x,   y-1, z);
    float l4 = indexedLerp(lookupHash[b3],   x-1, y-1, z);
    float l5 = indexedLerp(lookupHash[(uint8_t)a2+1], x,   y,   z-1);
    float l6 = indexedLerp(lookupHash[(uint8_t)b2+1], x-1, y,   z-1);
    float l7 = indexedLerp(lookupHash[(uint8_t)a3+1], x,   y-1, z-1);
    float l8 = indexedLerp(lookupHash[(uint8_t)b3+1], x-1, y-1, z-1);

    //linearlly interpolate down to 2d(on the x axis)
    l1 = lerp(t1, l1, l2);
    l3 = lerp(t1, l3, l4);
    l5 = lerp(t1, l5, l6);
    l7 = lerp(t1, l7, l8);

    //linearlly interpolate down to 1d(on the y axis)
    l1 = lerp(t2, l1, l3);
    l5 = lerp(t2, l5, l7);

    //linearlly interpolate down to 0d(on the z axis)
    return lerp(t3, l1, l5);
}

//sample stacked perlin noisemaps(one OctaveNoise)
float sampleOctave(const OctaveNoise *noise, int octaveMax, float x, float y, float z)
{
    float v = 0;
    for (int i = 0; i < noise->octcnt && i < octaveMax; i++)
    {
        const PerlinNoise *p = &noise->octaves[i];
        float lf = p->lacunarity;
        float ax = x * lf;
        float az = z * lf;
        float pv = samplePerlin(p, ax, y, az);
        v += p->amplitude * pv;
    }
    return v;
}

//sample together two stacked octave noisemaps to get the final continentalness
//noisemap at a layer
float sampleDoublePerlin(const DoublePerlinNoise *noise,
        int octaveMax, float x, float y, float z)
{
    //biomes are actually sampled at 1:4. look here if scales are being weird
    //x /= 4;
    //z /= 4;
    float v = 0.0;

    //DRY this DRY that whatever
    int na = (octaveMax + 1) >> 1;
    int nb = octaveMax - na;

    const float f = 337.0 / 331.0;
    v += sampleOctave(&noise->octA, na, x, y, z);
    v += sampleOctave(&noise->octB, nb, x*f, y, z*f);

    return v * noise->amplitude;
}

//END_SAMPLING

#endif // NOISELIB_H