#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 a dot product. 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 wrapDotProduct(float XoY, float k) {
  float lambertian = XoY;
  float half_lambertian = pow(max(1e-4, (XoY + 0.5f) / (1.0f + 0.5f)), 2);
  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);
}

#endif  // __MATH_INC