#ifndef __MATH_INC
#define __MATH_INC

#define PI              3.14159265358979323846264f
#define TAU             (2.0f * PI)
#define HALF_PI         (PI * 0.5f)
#define RCP_PI          (1.0f / PI)
#define RCP_TAU         (1.0f / TAU)
#define PHI             1.618033989f
#define RCP_PHI         0.618033989f
#define SQRT_2          1.414213562f
#define SQRT_2_RCP      0.707106781f
#define RCP_SQRT_2      0.707106781f
#define RCP_SQRT_3      0.577350269f
#define TWO_OVER_THREE  0.6666666666666666f
#define SQRT_3          1.73205081f
#define SQRT_3_OVER_2   0.8660254037844386f
#define EULERS_CONSTANT 2.718281828f

float sin_noise_3d(float3 uvw) {
  return sin(uvw[0]) * sin(uvw[1]) * sin(uvw[2]);
}

float sin_noise_3d_fbm(float3 uvw, uint octaves, float k, float strength) {
  float result = 0;
  float factor = 1.0f;
  for (uint i = 0; i < octaves; i++) {
    result += sin_noise_3d(uvw) * factor * strength;
    uvw *= k;
    factor /= k;
  }
  return result;
}

// Wrap NoL. Assume it's already clamped.
// At k=0, you get standard lambertian shading.
// At k=0.5, you get half-lambertian shading.
// At k=1.0, you get flat shading.
// k must be on [0, 1].
// Energy preserving, within some small bound.
float wrapNoL(float NoL, float k) {
  float lambertian = NoL;
  float tmp = (NoL + 0.5f) / (1.0f + 0.5f);
  float half_lambertian = tmp * tmp;
  float flat = RCP_PI;

  if (k < 0.5) {
    return lerp(lambertian, half_lambertian, k * 2.0f);
  } else {
    return lerp(half_lambertian, flat, k * 2.0f - 1.0f);
  }
}

// Vector rotation with a quaternion
// https://blog.molecular-matters.com/2013/05/24/a-faster-quaternion-vector-multiplication/
float3 rotate_vector(float3 v, float4 q)
{
  float3 t = 2.0 * cross(q.xyz, v);
  return v + q.w * t + cross(q.xyz, t);
}

// Cartesian to cube hexagonal coordinates.
// Based on this: https://backdrifting.net/post/064_hex_grids
float3 cart_to_hex(float2 cart) {
  float p = cart.x;
  float q = dot(cart, float2(0.5f, SQRT_3_OVER_2));
  float r = dot(cart, float2(0.5f, -SQRT_3_OVER_2));

  return float3(p, q, r) * 0.5f;
}

float2 hex_to_cart(float3 cart) {
  return float2(
      cart[0] + (cart[1] + cart[2]) * 0.5f,
      (cart[1] - cart[2]) * SQRT_3_OVER_2);
}


float3 round_hex(float3 hex_coord) {
    float3 rounded = round(hex_coord);
    float3 diff = abs(rounded - hex_coord);
    float error = rounded.x - rounded.y - rounded.z;

    float3 max_mask = float3(
        diff.x > diff.y && diff.x > diff.z,
        diff.y > diff.x && diff.y > diff.z,
        diff.z > diff.x && diff.z > diff.y);

    rounded += error * float3(-1, 1, 1) * max_mask;
    return rounded;
}

// Reoriented normal mapping
// https://blog.selfshadow.com/publications/blending-in-detail/
// Inputs are in tangent space.
float3 blendNormalsHill12(float3 n0, float3 n1) {
  n0.z += 1.0;
  n1.xy = -n1.xy;
  
  return normalize(n0 * dot(n0, n1) - n1 * n0.z);
}

// https://en.wikipedia.org/wiki/Relative_luminance
float luminance(float3 rgb) {
  return dot(float3(0.2126, 0.7152, 0.0722), rgb);
}

float4 alpha_blend(float4 front, float4 back) {
  return float4(front.rgb * front.a + back.rgb * (1 - front.a), front.a + back.a * (1 - front.a));
}

// Cheap procedural 3D hash -> [0,1]^3. Based on the "hashwithoutsine" family.
float3 hash33_fast(float3 p) {
  p = frac(p * float3(0.1031, 0.1030, 0.0973));
  p += dot(p, p.yxz + 33.33);
  return frac((p.xxy + p.yxx) * p.zyx);
}

// Procedural value noise in [0,1]^3 — trilinear interpolation of hashed corners.
float3 value_noise3(float3 p) {
  float3 i = floor(p);
  float3 f = frac(p);
  float3 u = f * f * (3.0 - 2.0 * f);

  return lerp(
    lerp(lerp(hash33_fast(i + float3(0, 0, 0)), hash33_fast(i + float3(1, 0, 0)), u.x),
         lerp(hash33_fast(i + float3(0, 1, 0)), hash33_fast(i + float3(1, 1, 0)), u.x), u.y),
    lerp(lerp(hash33_fast(i + float3(0, 0, 1)), hash33_fast(i + float3(1, 0, 1)), u.x),
         lerp(hash33_fast(i + float3(0, 1, 1)), hash33_fast(i + float3(1, 1, 1)), u.x), u.y),
    u.z);
}

// Domain warping using procedural value noise.
float3 domain_warp_procedural(float3 uvw, float strength,
    uint octaves, float lacunarity, float gain) {
  float3 noise = 0;
  float g = 1;

  for (uint ii = 0; ii < octaves; ++ii) {
    noise += value_noise3(uvw + noise * strength) * g;
    uvw *= lacunarity;
    g *= gain;
  }

  noise *= (1 - gain) / (1 - pow(gain, octaves));

  return noise;
}

// Domain warping using a 3D noise texture. Texture should have an EV of
// 0.5.
float3 domain_warp_3d_tex(texture3D noise_tex, float3 uvw, float strength,
    uint octaves, float lacunarity, float gain) {
  float3 noise = 0;
  float g = 1;

  float3 uvw_dx = ddx(uvw);
  float3 uvw_dy = ddy(uvw);

  for (uint ii = 0; ii < octaves; ++ii) {
    noise += noise_tex.SampleGrad(aniso4_trilinear_repeat_s, uvw + noise * strength, uvw_dx, uvw_dy).rgb * g;
    uvw *= lacunarity;
    g *= gain;
    uvw_dx *= lacunarity;
    uvw_dy *= lacunarity;
  }

  // Normalize: geometric series 1 + r + ... + r^{n-1} = (1 - r^n) / (1 - r)
  noise *= (1 - gain) / (1 - pow(gain, octaves));

  return noise;
}

// Return distance to the nearest cell edge in a Voronoi pattern.
float voronoi_edge_distance(float3 x) {
  float3 p = floor(x);
  float3 f = frac(x);
  float d1 = 1e6;
  float d2 = 1e6;
  float3 p1 = 0;
  float3 p2 = 0;

  for (int k = -1; k <= 1; k++) {
    for (int j = -1; j <= 1; j++) {
      for (int i = -1; i <= 1; i++) {
        float3 b = float3(i, j, k);
        float3 r = b + hash33_fast(p + b) - f;
        float d = dot(r, r);
        if (d < d1) {
          d2 = d1;
          p2 = p1;
          d1 = d;
          p1 = r;
        } else if (d < d2) {
          d2 = d;
          p2 = r;
        }
      }
    }
  }
  return (d2 - d1) / (2.0 * max(1e-4, length(p2 - p1)));
}

#endif  // __MATH_INC