WGSL Utility Functions

fn backIn(t: f32) -> f32 {
  let s = 1.70158;
  let tt = clamp(t, 0.0, 1.0);
  return tt * tt * ((s + 1.0) * tt - s);
}

fn backOut(t: f32) -> f32 {
  let s = 1.70158;
  let tt = clamp(t, 0.0, 1.0);
  let tMinus = tt - 1.0;
  return tMinus * tMinus * ((s + 1.0) * tMinus + s) + 1.0;
}

fn bezierCubic(t: f32, p0: f32, p1: f32, p2: f32, p3: f32) -> vec2<f32> {
  // Clamp t to [0,1]
  let tt = clamp(t, 0.0, 1.0);
  
  // Calculate curve value using cubic Bezier formula
  // P(t) = (1-t)³P₀ + 3(1-t)²tP₁ + 3(1-t)t²P₂ + t³P₃
  let t1 = 1.0 - tt;
  let t1Squared = t1 * t1;
  let t1Cubed = t1Squared * t1;
  let tSquared = tt * tt;
  let tCubed = tSquared * tt;
  
  let value = t1Cubed * p0 + 
             3.0 * t1Squared * tt * p1 + 
             3.0 * t1 * tSquared * p2 + 
             tCubed * p3;
  
  // Calculate derivative for tangent information
  // P'(t) = 3(1-t)²(P₁-P₀) + 6(1-t)t(P₂-P₁) + 3t²(P₃-P₂)
  let derivative = 3.0 * t1Squared * (p1 - p0) +
                  6.0 * t1 * tt * (p2 - p1) +
                  3.0 * tSquared * (p3 - p2);
  
  return vec2<f32>(value, derivative);
}

fn easeIn(t: f32, power: f32) -> f32 {
  return pow(clamp(t, 0.0, 1.0), power);
}

fn easeInOut(t: f32, power: f32) -> f32 {
  let tt = clamp(t, 0.0, 1.0);
  if (tt < 0.5) {
    return 0.5 * pow(2.0 * tt, power);
  } else {
    return 0.5 + 0.5 * (1.0 - pow(2.0 * (1.0 - tt), power));
  }
}

fn easeOut(t: f32, power: f32) -> f32 {
  return 1.0 - pow(1.0 - clamp(t, 0.0, 1.0), power);
}

fn elasticIn(t: f32) -> f32 {
  let tt = clamp(t, 0.0, 1.0);
  return sin(13.0 * 3.14159 * tt) * pow(2.0, 10.0 * (tt - 1.0));
}

fn elasticOut(t: f32) -> f32 {
  let tt = clamp(t, 0.0, 1.0);
  return sin(-13.0 * 3.14159 * (tt + 1.0)) * pow(2.0, -10.0 * tt) + 1.0;
}

fn springPhysics(t: f32, targetPosition: f32, initialPos: f32, initialVel: f32, stiffness: f32, damping: f32, mass: f32) -> vec2<f32> {
  // Ensure positive values for stiffness, damping, and mass
  let k = max(0.0001, stiffness);
  let d = max(0.0, damping);
  let m = max(0.0001, mass);
  
  // Calculate the angular frequency and damping ratio
  let omega = sqrt(k / m);
  let zeta = d / (2.0 * sqrt(k * m));
  
  // Initial displacement from targetPosition position
  let x0 = initialPos - targetPosition;
  let v0 = initialVel;
  
  var position: f32 = 0.0;
  var velocity: f32 = 0.0;
  
  if (zeta < 1.0) {
    // Underdamped case
    let omega_d = omega * sqrt(1.0 - zeta * zeta);
    let A = x0;
    let B = (v0 + zeta * omega * x0) / omega_d;
    
    // Calculate exponential decay term
    let expTerm = exp(-zeta * omega * t);
    
    // Calculate position and velocity
    position = targetPosition + expTerm * (A * cos(omega_d * t) + B * sin(omega_d * t));
    velocity = expTerm * (
                -zeta * omega * A * cos(omega_d * t) - omega_d * A * sin(omega_d * t) +
                -zeta * omega * B * sin(omega_d * t) + omega_d * B * cos(omega_d * t)
               );
  } else if (zeta == 1.0) {
    // Critically damped case
    let A = x0;
    let B = v0 + omega * x0;
    
    // Calculate exponential decay term
    let expTerm = exp(-omega * t);
    
    // Calculate position and velocity
    position = targetPosition + expTerm * (A + B * t);
    velocity = expTerm * (B - omega * (A + B * t));
  } else {
    // Overdamped case
    let omega1 = -omega * (zeta + sqrt(zeta * zeta - 1.0));
    let omega2 = -omega * (zeta - sqrt(zeta * zeta - 1.0));
    
    let A = (v0 - omega2 * x0) / (omega1 - omega2);
    let B = x0 - A;
    
    // Calculate position and velocity
    position = targetPosition + A * exp(omega1 * t) + B * exp(omega2 * t);
    velocity = A * omega1 * exp(omega1 * t) + B * omega2 * exp(omega2 * t);
  }
  
  return vec2<f32>(position, velocity);
}

//! requires hueToRgb
fn hslToRgb(hsl: vec3<f32>) -> vec3<f32> {
  let h = hsl.x;
  let s = hsl.y;
  let l = hsl.z;
  
  if (s == 0.0) {
    return vec3(l); // achromatic
  }
  
  let q = select(l * (1.0 + s), l + s - l * s, l < 0.5);
  let p = 2.0 * l - q;
  
  return vec3(
    hueToRgb(p, q, h + 1.0 / 3.0),
    hueToRgb(p, q, h),
    hueToRgb(p, q, h - 1.0 / 3.0)
  );
}

fn hueToRgb(p: f32, q: f32, t: f32) -> f32 {
  var t_adj = t;
  if (t_adj < 0.0) { t_adj += 1.0; }
  if (t_adj > 1.0) { t_adj -= 1.0; }
  if (t_adj < 1.0 / 6.0) { return p + (q - p) * 6.0 * t_adj; }
  if (t_adj < 1.0 / 2.0) { return q; }
  if (t_adj < 2.0 / 3.0) { return p + (q - p) * (2.0 / 3.0 - t_adj) * 6.0; }
  return p;
}

fn linearToSrgb(color: vec3<f32>) -> vec3<f32> {
  return pow(color, vec3(1.0 / 2.2));
}

fn palette(t: f32, a: vec3<f32>, b: vec3<f32>, c: vec3<f32>, d: vec3<f32>) -> vec3<f32> {
  return a + b * cos(6.28318 * (c * t + d));
}

fn srgbToLinear(color: vec3<f32>) -> vec3<f32> {
  return pow(color, vec3(2.2));
}

fn elasticWave(x: f32, amplitude: f32, frequency: f32, decay: f32, phase: f32) -> f32 {
  let d = max(0.001, decay);
  let decayTerm = exp(-d * x);
  let oscTerm = sin(frequency * x * 6.28318 + phase);
  return amplitude * decayTerm * oscTerm;
}

fn exponentialRamp(x: f32, base: f32, scale: f32, offset: f32) -> vec2<f32> {
  // Ensure base is positive and not 1 (which would make it linear)
  let b = select(base, 2.71828, abs(base - 1.0) < 0.001);
  
  // Calculate the exponential function
  let result = scale * pow(b, x) + offset;
  
  // Calculate the derivative
  let derivative = scale * pow(b, x) * log(b);
  
  return vec2<f32>(result, derivative);
}

fn logisticCurve(x: f32, midpoint: f32, steepness: f32, minValue: f32, maxValue: f32) -> vec2<f32> {
  // Scale factor for steepness
  let k = max(0.001, steepness);
  
  // Shift x relative to midpoint
  let z = -k * (x - midpoint);
  
  // Calculate the exponent
  let expTerm = exp(z);
  
  // Calculate the logistic function value
  let logistic = 1.0 / (1.0 + expTerm);
  
  // Scale to min-max range
  let range = maxValue - minValue;
  let value = minValue + range * logistic;
  
  // Calculate the derivative
  let derivative = range * k * expTerm / ((1.0 + expTerm) * (1.0 + expTerm));
  
  return vec2<f32>(value, derivative);
}

fn rotate2D(v: vec2<f32>, angle: f32) -> vec2<f32> {
  let c = cos(angle);
  let s = sin(angle);
  return vec2(v.x * c - v.y * s, v.x * s + v.y * c);
}

fn smoothStep(edge0: f32, edge1: f32, x: f32) -> f32 {
  let t = clamp((x - edge0) / (edge1 - edge0), 0.0, 1.0);
  return t * t * (3.0 - 2.0 * t);
}

fn smoothStepVec2(edge0: vec2f, edge1: vec2f, x: vec2f) -> vec2f {
  let t = clamp((x - edge0) / (edge1 - edge0), vec2f(0.0), vec2f(1.0));
  return t * t * (3.0 - 2.0 * t);
}

fn stepSequence(x: f32, steps: f32, smoothing: f32, minValue: f32, maxValue: f32) -> vec2<f32> {
  // Ensure at least 1 step and positive smoothing
  let numSteps = max(1.0, floor(steps));
  let smoothFactor = max(0.0, smoothing);
  
  // Normalize x to 0-1 range
  let normalizedX = fract(x);
  
  // Calculate the size of each step
  let stepSize = 1.0 / numSteps;
  
  // Calculate the current step (0 to numSteps-1)
  let currentStep = floor(normalizedX * numSteps);
  let nextStep = fract(currentStep + 1.0);
  
  // Calculate progress within the current step
  let stepProgress = fract(normalizedX * numSteps);
  
  // Calculate the progress values for current and next steps
  let currentStepValue = currentStep / (numSteps - 1.0);
  
  // Prepare next step value, handle the last step case
  var nextStepValue: f32 = 0.0;
  if (currentStep >= numSteps - 1.0) {
    nextStepValue = 1.0;
  } else {
    nextStepValue = nextStep / (numSteps - 1.0);
  }
  
  // Apply smoothing between steps if needed
  var result: f32 = 0.0;
  
  if (smoothFactor > 0.0 && stepProgress > (1.0 - smoothFactor) && numSteps > 1.0) {
    // Calculate smoothing factor
    let t = (stepProgress - (1.0 - smoothFactor)) / smoothFactor;
    
    // Smoothstep for better transition
    let smoothT = t * t * (3.0 - 2.0 * t);
    
    // Interpolate between current and next step
    result = mix(currentStepValue, nextStepValue, smoothT);
  } else {
    result = currentStepValue;
  }
  
  // Scale to min-max range
  let range = maxValue - minValue;
  let finalResult = minValue + result * range;
  
  return vec2<f32>(finalResult, currentStep);
}

fn taylorInvSqrt4(r: vec4f) -> vec4f { 
    return 1.79284291400159 - 0.85373472095314 * r; 
}

//! requires noise2D
fn fbm2D(p: vec2<f32>, octaves: i32) -> f32 {
  var value = 0.0;
  var amplitude = 0.5;
  var frequency = 1.0;
  for (var i = 0; i < octaves; i++) {
    value += amplitude * noise2D(p * frequency);
    amplitude *= 0.5;
    frequency *= 2.0;
  }
  return value;
}

//! requires noise3D
  fn fbm3D(x: vec3f, seed: f32) -> f32 {
    var p = x + seed;
    var f = 0.0;
    var w = 0.5;
    for (var i = 0; i < 5; i++) {
      f += w * noise3D(p);
      p *= 2.0;
      w *= 0.5;
    }
    return f;
  }

fn hash1D(p: f32) -> f32 {
  // Convert to integer and apply bit manipulation
  let x = bitcast<u32>(p + 123.456789);
  var h = x;
  
  // Wang hash function
  h = (h ^ 61u) ^ (h >> 16u);
  h = h + (h << 3u);
  h = h ^ (h >> 4u);
  h = h * 0x27d4eb2du;
  h = h ^ (h >> 15u);
  
  // Convert back to float and normalize
  return f32(h) / 4294967296.0;
}

fn hash22(p: vec2<f32>) -> vec2<f32> {
  var p3 = fract(vec3<f32>(p.xyx) * vec3<f32>(0.1031, 0.1030, 0.0973));
  p3 += dot(p3, p3.yzx + 33.33);
  return fract((p3.xx + p3.yz) * p3.zy);
}

fn hash31(p: vec3<f32>) -> f32 {
  var p3 = fract(p * vec3<f32>(0.1031, 0.1030, 0.0973));
  p3 += dot(p3, p3.yxz + 33.33);
  return fract((p3.x + p3.y) * p3.z);
}

fn hash3D(p: vec3<f32>) -> vec3<f32> {
  var p3 = fract(p * vec3<f32>(0.1031, 0.1030, 0.0973));
  p3 += dot(p3, p3.yxz + 33.33);
  return fract((p3.xxy + p3.yxx) * p3.zyx) * 2.0 - 1.0;
}

fn mod289(x: vec4f) -> vec4f { 
    return x - floor(x * (1. / 289.)) * 289.; 
}

//! requires hash22
fn noise2D(p: vec2<f32>) -> f32 {
  let i = floor(p);
  let f = fract(p);
  let u = f * f * (3.0 - 2.0 * f);
  return mix(
    mix(dot(hash22(i + vec2<f32>(0.0, 0.0)), f - vec2<f32>(0.0, 0.0)),
        dot(hash22(i + vec2<f32>(1.0, 0.0)), f - vec2<f32>(1.0, 0.0)), u.x),
    mix(dot(hash22(i + vec2<f32>(0.0, 1.0)), f - vec2<f32>(0.0, 1.0)),
        dot(hash22(i + vec2<f32>(1.0, 1.0)), f - vec2<f32>(1.0, 1.0)), u.x), u.y);
}

//! requires hash31
fn noise3D(x: vec3<f32>) -> f32 {
  let p = floor(x);
  let f = fract(x);
  
  return mix(
    mix(
      mix(hash31(p), 
          hash31(p + vec3<f32>(1.0, 0.0, 0.0)), 
          f.x),
      mix(hash31(p + vec3<f32>(0.0, 1.0, 0.0)), 
          hash31(p + vec3<f32>(1.0, 1.0, 0.0)), 
          f.x),
      f.y),
    mix(
      mix(hash31(p + vec3<f32>(0.0, 0.0, 1.0)), 
          hash31(p + vec3<f32>(1.0, 0.0, 1.0)), 
          f.x),
      mix(hash31(p + vec3<f32>(0.0, 1.0, 1.0)), 
          hash31(p + vec3<f32>(1.0, 1.0, 1.0)), 
          f.x),
      f.y),
    f.z);
}

fn pcg(n: u32) -> u32 {
    var h = n * 747796405u + 2891336453u;
    h = ((h >> ((h >> 28u) + 4u)) ^ h) * 277803737u;
    return (h >> 22u) ^ h;
}

fn pcg2d(p: vec2u) -> vec2u {
    var v = p * 1664525u + 1013904223u;
    v.x += v.y * 1664525u; v.y += v.x * 1664525u;
    v ^= v >> vec2u(16u);
    v.x += v.y * 1664525u; v.y += v.x * 1664525u;
    v ^= v >> vec2u(16u);
    return v;
}

fn pcg3d(p: vec3u) -> vec3u {
    var v = p * 1664525u + 1013904223u;
    v.x += v.y*v.z; v.y += v.z*v.x; v.z += v.x*v.y;
    v ^= v >> vec3u(16u);
    v.x += v.y*v.z; v.y += v.z*v.x; v.z += v.x*v.y;
    return v;
}

fn pcg4d(p: vec4u) -> vec4u {
    var v = p * 1664525u + 1013904223u;
    v.x += v.y*v.w; v.y += v.z*v.x; v.z += v.x*v.y; v.w += v.y*v.z;
    v ^= v >> vec4u(16u);
    v.x += v.y*v.w; v.y += v.z*v.x; v.z += v.x*v.y; v.w += v.y*v.z;
    return v;
}

fn perlinNoise2_permute4(x: vec4f) -> vec4f { 
    return ((x * 34. + 1.) * x) % vec4f(289.); 
}

fn perlinNoise2_fade2(t: vec2f) -> vec2f { 
    return t * t * t * (t * (t * 6. - 15.) + 10.); 
}

fn perlinNoise2(P: vec2f) -> f32 {
    var Pi: vec4f = floor(P.xyxy) + vec4f(0., 0., 1., 1.);
    let Pf = fract(P.xyxy) - vec4f(0., 0., 1., 1.);
    Pi = Pi % vec4f(289.); // To avoid truncation effects in permutation
    let ix = Pi.xzxz;
    let iy = Pi.yyww;
    let fx = Pf.xzxz;
    let fy = Pf.yyww;
    let i = perlinNoise2_permute4(perlinNoise2_permute4(ix) + iy);
    var gx: vec4f = 2. * fract(i * 0.0243902439) - 1.; // 1/41 = 0.024...
    let gy = abs(gx) - 0.5;
    let tx = floor(gx + 0.5);
    gx = gx - tx;
    var g00: vec2f = vec2f(gx.x, gy.x);
    var g10: vec2f = vec2f(gx.y, gy.y);
    var g01: vec2f = vec2f(gx.z, gy.z);
    var g11: vec2f = vec2f(gx.w, gy.w);
    let norm = 1.79284291400159 - 0.85373472095314 *
        vec4f(dot(g00, g00), dot(g01, g01), dot(g10, g10), dot(g11, g11));
    g00 = g00 * norm.x;
    g01 = g01 * norm.y;
    g10 = g10 * norm.z;
    g11 = g11 * norm.w;
    let n00 = dot(g00, vec2f(fx.x, fy.x));
    let n10 = dot(g10, vec2f(fx.y, fy.y));
    let n01 = dot(g01, vec2f(fx.z, fy.z));
    let n11 = dot(g11, vec2f(fx.w, fy.w));
    let fade_xy = perlinNoise2_fade2(Pf.xy);
    let n_x = mix(vec2f(n00, n01), vec2f(n10, n11), vec2f(fade_xy.x));
    let n_xy = mix(n_x.x, n_x.y, fade_xy.y);
    return 2.3 * n_xy;
}

//! requires taylorInvSqrt4
fn perlinNoise3_permute4(x: vec4f) -> vec4f { 
    return ((x * 34. + 1.) * x) % vec4f(289.); 
}

fn perlinNoise3_fade3(t: vec3f) -> vec3f { 
    return t * t * t * (t * (t * 6. - 15.) + 10.); 
}

fn perlinNoise3(P: vec3f) -> f32 {
    var Pi0 : vec3f = floor(P); // Integer part for indexing
    var Pi1 : vec3f = Pi0 + vec3f(1.); // Integer part + 1
    Pi0 = Pi0 % vec3f(289.);
    Pi1 = Pi1 % vec3f(289.);
    let Pf0 = fract(P); // Fractional part for interpolation
    let Pf1 = Pf0 - vec3f(1.); // Fractional part - 1.
    let ix = vec4f(Pi0.x, Pi1.x, Pi0.x, Pi1.x);
    let iy = vec4f(Pi0.yy, Pi1.yy);
    let iz0 = Pi0.zzzz;
    let iz1 = Pi1.zzzz;

    let ixy = perlinNoise3_permute4(perlinNoise3_permute4(ix) + iy);
    let ixy0 = perlinNoise3_permute4(ixy + iz0);
    let ixy1 = perlinNoise3_permute4(ixy + iz1);

    var gx0: vec4f = ixy0 / 7.;
    var gy0: vec4f = fract(floor(gx0) / 7.) - 0.5;
    gx0 = fract(gx0);
    var gz0: vec4f = vec4f(0.5) - abs(gx0) - abs(gy0);
    var sz0: vec4f = step(gz0, vec4f(0.));
    gx0 = gx0 + sz0 * (step(vec4f(0.), gx0) - 0.5);
    gy0 = gy0 + sz0 * (step(vec4f(0.), gy0) - 0.5);

    var gx1: vec4f = ixy1 / 7.;
    var gy1: vec4f = fract(floor(gx1) / 7.) - 0.5;
    gx1 = fract(gx1);
    var gz1: vec4f = vec4f(0.5) - abs(gx1) - abs(gy1);
    var sz1: vec4f = step(gz1, vec4f(0.));
    gx1 = gx1 - sz1 * (step(vec4f(0.), gx1) - 0.5);
    gy1 = gy1 - sz1 * (step(vec4f(0.), gy1) - 0.5);

    var g000: vec3f = vec3f(gx0.x, gy0.x, gz0.x);
    var g100: vec3f = vec3f(gx0.y, gy0.y, gz0.y);
    var g010: vec3f = vec3f(gx0.z, gy0.z, gz0.z);
    var g110: vec3f = vec3f(gx0.w, gy0.w, gz0.w);
    var g001: vec3f = vec3f(gx1.x, gy1.x, gz1.x);
    var g101: vec3f = vec3f(gx1.y, gy1.y, gz1.y);
    var g011: vec3f = vec3f(gx1.z, gy1.z, gz1.z);
    var g111: vec3f = vec3f(gx1.w, gy1.w, gz1.w);

    let norm0 = taylorInvSqrt4(
        vec4f(dot(g000, g000), dot(g010, g010), dot(g100, g100), dot(g110, g110)));
    g000 = g000 * norm0.x;
    g010 = g010 * norm0.y;
    g100 = g100 * norm0.z;
    g110 = g110 * norm0.w;
    let norm1 = taylorInvSqrt4(
        vec4f(dot(g001, g001), dot(g011, g011), dot(g101, g101), dot(g111, g111)));
    g001 = g001 * norm1.x;
    g011 = g011 * norm1.y;
    g101 = g101 * norm1.z;
    g111 = g111 * norm1.w;

    let n000 = dot(g000, Pf0);
    let n100 = dot(g100, vec3f(Pf1.x, Pf0.yz));
    let n010 = dot(g010, vec3f(Pf0.x, Pf1.y, Pf0.z));
    let n110 = dot(g110, vec3f(Pf1.xy, Pf0.z));
    let n001 = dot(g001, vec3f(Pf0.xy, Pf1.z));
    let n101 = dot(g101, vec3f(Pf1.x, Pf0.y, Pf1.z));
    let n011 = dot(g011, vec3f(Pf0.x, Pf1.yz));
    let n111 = dot(g111, Pf1);

    var fade_xyz: vec3f = perlinNoise3_fade3(Pf0);
    let temp = vec4f(f32(fade_xyz.z)); // simplify after chrome bug fix
    let n_z = mix(vec4f(n000, n100, n010, n110), vec4f(n001, n101, n011, n111), temp);
    let n_yz = mix(n_z.xy, n_z.zw, vec2f(f32(fade_xyz.y))); // simplify after chrome bug fix
    let n_xyz = mix(n_yz.x, n_yz.y, fade_xyz.x);
    return 2.2 * n_xyz;
}

//! requires mod289
fn perm4(x: vec4f) -> vec4f { 
    return mod289(((x * 34.) + 1.) * x); 
}

fn rand11(n: f32) -> f32 { 
    return fract(sin(n) * 43758.5453123); 
}

fn rand22(n: vec2f) -> f32 { 
    return fract(sin(dot(n, vec2f(12.9898, 4.1414))) * 43758.5453); 
}

fn snoise2D_permute3(x: vec3<f32>) -> vec3<f32> {
    return ((x * 34.0 + 1.0) * x) % vec3<f32>(289.0);
}

fn snoise2D(v: vec2<f32>) -> f32 {
    let C = vec4<f32>(0.211324865405187, 0.366025403784439, -0.577350269189626, 0.024390243902439);
    var i = floor(v + dot(v, C.yy));
    let x0 = v - i + dot(i, C.xx);
    
    var i1: vec2<f32>;
    if (x0.x > x0.y) {
        i1 = vec2<f32>(1.0, 0.0);
    } else {
        i1 = vec2<f32>(0.0, 1.0);
    }
    
    var x12 = x0.xyxy + C.xxzz;
    x12 = vec4<f32>(x12.xy - i1, x12.zw);
    
    i = i % vec2<f32>(289.0);
    let p = snoise2D_permute3(snoise2D_permute3(i.y + vec3<f32>(0.0, i1.y, 1.0)) + i.x + vec3<f32>(0.0, i1.x, 1.0));
    
    var m = max(0.5 - vec3<f32>(dot(x0, x0), dot(x12.xy, x12.xy), dot(x12.zw, x12.zw)), vec3<f32>(0.0));
    m = m * m;
    m = m * m;
    
    let x = 2.0 * fract(p * C.www) - 1.0;
    let h = abs(x) - 0.5;
    let ox = floor(x + 0.5);
    let a0 = x - ox;
    
    m *= 1.79284291400159 - 0.85373472095314 * (a0 * a0 + h * h);
    
    var g: vec3<f32>;
    g.x = a0.x * x0.x + h.x * x0.y;
    g.y = a0.y * x12.x + h.y * x12.y;
    g.z = a0.z * x12.z + h.z * x12.w;
    
    return 130.0 * dot(m, g);
}

fn snoise3D_permute4(x: vec4<f32>) -> vec4<f32> {
    return ((x * 34.0 + 1.0) * x) % vec4<f32>(289.0);
}

fn snoise3D_taylorInvSqrt4(r: vec4<f32>) -> vec4<f32> {
    return 1.79284291400159 - 0.85373472095314 * r;
}

fn snoise3D(v: vec3<f32>) -> f32 {
    let C = vec2<f32>(1.0 / 6.0, 1.0 / 3.0);
    let D = vec4<f32>(0.0, 0.5, 1.0, 2.0);
    
    // First corner
    var i = floor(v + dot(v, C.yyy));
    let x0 = v - i + dot(i, C.xxx);
    
    // Other corners
    let g = step(x0.yzx, x0.xyz);
    let l = 1.0 - g;
    let i1 = min(g.xyz, l.zxy);
    let i2 = max(g.xyz, l.zxy);
    
    // x0 = x0 - 0. + 0.0 * C
    let x1 = x0 - i1 + 1.0 * C.xxx;
    let x2 = x0 - i2 + 2.0 * C.xxx;
    let x3 = x0 - 1.0 + 3.0 * C.xxx;
    
    // Permutations
    i = i % vec3<f32>(289.0);
    let p = snoise3D_permute4(snoise3D_permute4(snoise3D_permute4(
        i.z + vec4<f32>(0.0, i1.z, i2.z, 1.0)) +
        i.y + vec4<f32>(0.0, i1.y, i2.y, 1.0)) +
        i.x + vec4<f32>(0.0, i1.x, i2.x, 1.0));
    
    // Gradients
    // (N*N points uniformly over a square, mapped onto an octahedron.)
    let n_ = 1.0 / 7.0; // N=7
    let ns = n_ * D.wyz - D.xzx;
    
    let j = p - 49.0 * floor(p * ns.z * ns.z); // mod(p,N*N)
    
    let x_ = floor(j * ns.z);
    let y_ = floor(j - 7.0 * x_); // mod(j,N)
    
    let x = x_ * ns.x + ns.yyyy;
    let y = y_ * ns.x + ns.yyyy;
    let h = 1.0 - abs(x) - abs(y);
    
    let b0 = vec4<f32>(x.xy, y.xy);
    let b1 = vec4<f32>(x.zw, y.zw);
    
    let s0 = floor(b0) * 2.0 + 1.0;
    let s1 = floor(b1) * 2.0 + 1.0;
    let sh = -step(h, vec4<f32>(0.0));
    
    let a0 = b0.xzyw + s0.xzyw * sh.xxyy;
    let a1 = b1.xzyw + s1.xzyw * sh.zzww;
    
    var p0 = vec3<f32>(a0.xy, h.x);
    var p1 = vec3<f32>(a0.zw, h.y);
    var p2 = vec3<f32>(a1.xy, h.z);
    var p3 = vec3<f32>(a1.zw, h.w);
    
    // Normalise gradients
    let norm = snoise3D_taylorInvSqrt4(vec4<f32>(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3)));
    p0 *= norm.x;
    p1 *= norm.y;
    p2 *= norm.z;
    p3 *= norm.w;
    
    // Mix final noise value
    var m = max(0.6 - vec4<f32>(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3)), vec4<f32>(0.0));
    m = m * m;
    return 42.0 * dot(m * m, vec4<f32>(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3)));
}

//! requires mod289 taylorInvSqrt4
fn snoise_permute4(x: vec4<f32>) -> vec4<f32> {
    return mod289(((x * 34.0) + 10.0) * x);
}

fn snoise_mod289f(x: f32) -> f32 {
    return x - floor(x * (1.0 / 289.0)) * 289.0;
}

fn snoise_permute(x: f32) -> f32 {
    return snoise_mod289f(((x * 34.0) + 10.0) * x);
}

fn snoise_taylorInvSqrt(r: f32) -> f32 {
    return 1.79284291400159 - 0.85373472095314 * r;
}

fn snoise_grad4(j: f32, ip: vec4<f32>) -> vec4<f32> {
    let ones = vec4(1.0, 1.0, 1.0, -1.0);

    var p: vec4<f32>;
    var s: vec4<f32>;

    p = vec4(floor(fract(vec3(j) * ip.xyz) * 7.0) * ip.z - 1.0, p.w);
    p.w = 1.5 - dot(abs(p.xyz), ones.xyz);
    s = select(vec4(0.0), vec4(1.0), p < vec4(0.0));
    p = vec4(p.xyz + (s.xyz * 2.0 - 1.0) * s.www, p.w);

    return p;
}

fn snoise(v: vec4<f32>) -> f32 {
    let C: vec4<f32> = vec4(
        0.138196601125011, // (5 - sqrt(5))/20  G4
        0.276393202250021, // 2 * G4
        0.414589803375032, // 3 * G4
        -0.447213595499958 // -1 + 4 * G4
    );

    // (sqrt(5) - 1)/4 = F4, used once below
    let F4: f32 = 0.309016994374947451;

    // First corner
    var i: vec4<f32> = floor(v + dot(v, vec4(F4)));
    var x0: vec4<f32> = v - i + dot(i, C.xxxx);

    // Other corners

    // Rank sorting originally contributed by Bill Licea-Kane, AMD (formerly ATI)
    var i0: vec4<f32>;
    var isX: vec3<f32> = step(x0.yzw, x0.xxx);
    var isYZ: vec3<f32> = step(x0.zww, x0.yyz);

    //  i0.x = dot( isX, vec3( 1.0 ) );
    i0.x = isX.x + isX.y + isX.z;
    var minusX = 1.0 - isX;
    i0.y = minusX.x;
    i0.z = minusX.y;
    i0.w = minusX.z;

    //  i0.y += dot( isYZ.xy, vec2( 1.0 ) );
    i0.y += isYZ.x + isYZ.y;
    var minusY = 1.0 - isYZ;
    i0.z += minusY.x;
    i0.w += minusY.y;
    i0.z += isYZ.z;
    i0.w += minusY.z;

    // i0 now contains the unique values 0,1,2,3 in each channel
    var i3: vec4<f32> = vec4<f32>(clamp(i0, vec4(0.0), vec4(1.0)));
    var i2: vec4<f32> = clamp(i0 - vec4(1.0), vec4(0.0), vec4(1.0));
    var i1: vec4<f32> = clamp(i0 - vec4(2.0), vec4(0.0), vec4(1.0));

    //  x0 = x0 - 0.0 + 0.0 * C.xxxx
    //  x1 = x0 - i1  + 1.0 * C.xxxx
    //  x2 = x0 - i2  + 2.0 * C.xxxx
    //  x3 = x0 - i3  + 3.0 * C.xxxx
    //  x4 = x0 - 1.0 + 4.0 * C.xxxx
    var x1: vec4<f32> = x0 - i1 + C.xxxx;
    var x2: vec4<f32> = x0 - i2 + C.yyyy;
    var x3: vec4<f32> = x0 - i3 + C.zzzz;
    var x4: vec4<f32> = x0 + C.wwww;

    // Permutations
    i = mod289(i);
    var j0: f32 = snoise_permute(snoise_permute(snoise_permute(snoise_permute(i.w) + i.z) + i.y) + i.x);
    var j1: vec4<f32> = snoise_permute4(snoise_permute4(snoise_permute4(snoise_permute4(i.w + vec4(i1.w, i2.w, i3.w, 1.0)) + i.z + vec4(i1.z, i2.z, i3.z, 1.0)) + i.y + vec4(i1.y, i2.y, i3.y, 1.0)) + i.x + vec4(i1.x, i2.x, i3.x, 1.0));

    // Gradients: 7x7x6 points over a cube, mapped onto a 4-cross polytope
    // 7*7*6 = 294, which is close to the ring size 17*17 = 289.
    var ip: vec4<f32> = vec4(1.0 / 294.0, 1.0 / 49.0, 1.0 / 7.0, 0.0);

    var p0: vec4<f32> = snoise_grad4(j0, ip);
    var p1: vec4<f32> = snoise_grad4(j1.x, ip);
    var p2: vec4<f32> = snoise_grad4(j1.y, ip);
    var p3: vec4<f32> = snoise_grad4(j1.z, ip);
    var p4: vec4<f32> = snoise_grad4(j1.w, ip);

    // Normalise gradients
    var norm: vec4<f32> = taylorInvSqrt4(vec4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3)));
    p0 *= norm.x;
    p1 *= norm.y;
    p2 *= norm.z;
    p3 *= norm.w;
    p4 *= snoise_taylorInvSqrt(dot(p4, p4));

    // Mix contributions from the five corners
    var m0: vec3<f32> = max(0.6 - vec3(dot(x0, x0), dot(x1, x1), dot(x2, x2)), vec3(0.0));
    var m1: vec2<f32> = max(0.6 - vec2(dot(x3, x3), dot(x4, x4)), vec2(0.0));
    m0 = m0 * m0;
    m1 = m1 * m1;
    return 49.0 * (dot(m0 * m0, vec3(dot(p0, x0), dot(p1, x1), dot(p2, x2))) + dot(m1 * m1, vec2(dot(p3, x3), dot(p4, x4))));
}

//! requires rand11Sin
fn noise(p: f32) -> f32 {
    let fl = floor(p);
    return mix(rand11(fl), rand11(fl + 1.), fract(p));
}

//! requires rand22Sin smoothStepVec2
fn noise2(n: vec2f) -> f32 {
    let d = vec2f(0., 1.);
    let b = floor(n);
    let f = smoothStepVec2(vec2f(0.), vec2f(1.), fract(n));
    return mix(mix(rand22(b), rand22(b + d.yx), f.x), mix(rand22(b + d.xy), rand22(b + d.yy), f.x), f.y);
}

//! requires perm4
fn noise3(p: vec3f) -> f32 {
    let a = floor(p);
    var d: vec3f = p - a;
    d = d * d * (3. - 2. * d);

    let b = a.xxyy + vec4f(0., 1., 0., 1.);
    let k1 = perm4(b.xyxy);
    let k2 = perm4(k1.xyxy + b.zzww);

    let c = k2 + a.zzzz;
    let k3 = perm4(c);
    let k4 = perm4(c + 1.);

    let o1 = fract(k3 * (1. / 41.));
    let o2 = fract(k4 * (1. / 41.));

    let o3 = o2 * d.z + o1 * (1. - d.z);
    let o4 = o3.yw * d.x + o3.xz * (1. - d.x);

    return o4.y * d.y + o4.x * (1. - d.y);
}

//! requires noise3D
fn warpNoise3D(x: vec3<f32>, seedVal: f32) -> vec3<f32> {
  var p = x + seedVal;
  var nx = 0.0;
  var ny = 0.0;
  var nz = 0.0;
  var w = 0.5;
  
  for (var i = 0; i < 3; i++) {
    nx += w * noise3D(p);
    ny += w * noise3D(p + vec3<f32>(13.5, 41.3, 17.8));
    nz += w * noise3D(p + vec3<f32>(31.2, 23.7, 11.9));
    p *= 2.0;
    w *= 0.5;
  }
  
  return vec3<f32>(nx, ny, nz) * 2.0 - 1.0;
}

fn xxhash32(n: u32) -> u32 {
    var h32 = n + 374761393u;
    h32 = 668265263u * ((h32 << 17) | (h32 >> (32 - 17)));
    h32 = 2246822519u * (h32 ^ (h32 >> 15));
    h32 = 3266489917u * (h32 ^ (h32 >> 13));
    return h32^(h32 >> 16);
}

fn xxhash322d(p: vec2u) -> u32 {
    let p2 = 2246822519u; let p3 = 3266489917u;
    let p4 = 668265263u; let p5 = 374761393u;
    var h32 = p.y + p5 + p.x * p3;
    h32 = p4 * ((h32 << 17) | (h32 >> (32 - 17)));
    h32 = p2 * (h32^(h32 >> 15));
    h32 = p3 * (h32^(h32 >> 13));
    return h32^(h32 >> 16);
}

fn xxhash323d(p: vec3u) -> u32 {
    let p2 = 2246822519u; let p3 = 3266489917u;
    let p4 = 668265263u; let p5 = 374761393u;
    var h32 =  p.z + p5 + p.x*p3;
    h32 = p4 * ((h32 << 17) | (h32 >> (32 - 17)));
    h32 += p.y * p3;
    h32 = p4 * ((h32 << 17) | (h32 >> (32 - 17)));
    h32 = p2 * (h32^(h32 >> 15));
    h32 = p3 * (h32^(h32 >> 13));
    return h32^(h32 >> 16);
}

fn xxhash324d(p: vec4u) -> u32 {
    let p2 = 2246822519u; let p3 = 3266489917u;
    let p4 = 668265263u; let p5 = 374761393u;
    var h32 = p.w + p5 + p.x * p3;
    h32 = p4 * ((h32 << 17) | (h32 >> (32 - 17)));
    h32 += p.y * p3;
    h32 = p4 * ((h32 << 17) | (h32 >> (32 - 17)));
    h32 += p.z  * p3;
    h32 = p4 * ((h32 << 17) | (h32 >> (32 - 17)));
    h32 = p2 * (h32^(h32 >> 15));
    h32 = p3 * (h32^(h32 >> 13));
    return h32 ^ (h32 >> 16);
}

fn sdfBend(position: vec3<f32>, angle: f32, axis: vec3<f32>, center: vec3<f32>) -> vec3<f32> {
  // Normalize the bend axis
  let axisNorm = normalize(axis);
  
  // Translate position relative to bend center
  let localPos = position - center;
  
  // Find perpendicular vectors to the bend axis to define the bend plane
  var perpVec1: vec3<f32>;
  if (abs(axisNorm.y) < 0.999) {
    perpVec1 = normalize(cross(vec3<f32>(0.0, 1.0, 0.0), axisNorm));
  } else {
    perpVec1 = normalize(cross(vec3<f32>(1.0, 0.0, 0.0), axisNorm));
  }
  let perpVec2 = normalize(cross(axisNorm, perpVec1));
  
  // Project the position onto the perpendicular vectors
  let proj1 = dot(localPos, perpVec1);
  let proj2 = dot(localPos, perpVec2);
  let axisProj = dot(localPos, axisNorm);
  
  // Calculate radius for the bend
  let radius = proj1;
  
  // Calculate the angle based on the distance along the bend direction
  let bendAngle = proj2 * angle;
  
  // Calculate the bent position using polar coordinates
  let c = cos(bendAngle);
  let s = sin(bendAngle);
  
  // Apply the transformation
  let bentPos = center + 
                axisNorm * axisProj +
                perpVec1 * (c * radius) +
                perpVec2 * (s * radius);
  
  return bentPos;
}

//! requires hash3D
fn sdfDisplace(position: vec3<f32>, amount: f32, frequency: f32, seed: f32) -> vec3<f32> {
  let noisePos = position * frequency + seed;
  let displacement = hash3D(noisePos) * amount;
  
  return position + displacement;
}

//! requires warpNoise3D
fn sdfDomainRepeat(position: vec3<f32>, cellSize: vec3<f32>, warpAmount: f32, warpScale: f32, seed: f32) -> vec3<f32> {
  // Calculate warping for position
  let warp = warpNoise3D(position * warpScale, seed) * warpAmount;
  
  // Apply warping to position
  let warpedPos = position + warp;
  
  // Calculate repetition
  return warpedPos - cellSize * round(warpedPos / cellSize);
}

fn sdfFiniteRepeat(position: vec3<f32>, spacing: vec3<f32>, count: vec3<f32>) -> vec3<f32> {
  let id = clamp(round(position / spacing), -count * 0.5, count * 0.5);
  return position - spacing * id;
}

fn sdfInfiniteRepeat(position: vec3<f32>, spacing: vec3<f32>) -> vec3<f32> {
  return position - spacing * round(position / spacing);
}

fn sdfTaper(position: vec3<f32>, amount: f32, axis: vec3<f32>, height: f32, offset: f32) -> vec3<f32> {
  let axisNorm = normalize(axis);
  
  // Project position onto the taper axis
  let axisPos = dot(position, axisNorm) - offset;
  
  // Calculate taper factor based on position along axis
  let t = clamp(axisPos / height, 0.0, 1.0);
  let taperFactor = 1.0 - amount * t;
  
  // Apply taper to the perpendicular components
  let axisComponent = axisPos * axisNorm;
  let perpComponent = position - dot(position, axisNorm) * axisNorm;
  
  return axisComponent + perpComponent * taperFactor;
}

fn sdfTwist(position: vec3<f32>, angle: f32, axis: vec3<f32>) -> vec3<f32> {
  // Normalize the axis
  let axisNorm = normalize(axis);
  
  // Project position onto the twist axis
  let proj = dot(position, axisNorm);
  
  // Calculate twist angle based on projection along axis
  let twistAngle = proj * angle;
  
  // Get sin and cos of the twist angle
  let s = sin(twistAngle);
  let c = cos(twistAngle);
  
  // Calculate vector from axis (the part that will be rotated)
  let axisProj = proj * axisNorm;
  let fromAxis = position - axisProj;
  
  // Find a perpendicular vector for the rotation
  let basis1 = normalize(fromAxis);
  let basis2 = cross(axisNorm, basis1);
  
  // Rotate using the basis vectors
  let rotated = axisProj + 
                basis1 * length(fromAxis) * c + 
                basis2 * length(fromAxis) * s;
  
  return rotated;
}

fn sdfChamferIntersect(distanceA: f32, distanceB: f32, radius: f32) -> f32 {
  return max(max(distanceA, distanceB), (distanceA + radius + distanceB) * 0.5);
}

fn sdfChamferSubtract(distanceA: f32, distanceB: f32, radius: f32) -> f32 {
  return max(max(distanceA, -distanceB), (distanceA + radius - distanceB) * 0.5);
}

fn sdfChamferUnion(distanceA: f32, distanceB: f32, radius: f32) -> f32 {
  return min(min(distanceA, distanceB), (distanceA - radius + distanceB) * 0.5);
}

fn sdfOpIntersect(distanceA: f32, distanceB: f32) -> f32 {
  return max(distanceA, distanceB);
}

fn sdfOpSubtract(distanceA: f32, distanceB: f32) -> f32 {
  return max(distanceA, -distanceB);
}

fn sdfOpUnion(distanceA: f32, distanceB: f32) -> f32 {
  return min(distanceA, distanceB);
}

fn sdfSmoothIntersect(distanceA: f32, distanceB: f32, smoothing: f32) -> f32 {
  let h = clamp(0.5 - 0.5 * (distanceB - distanceA) / smoothing, 0.0, 1.0);
  return mix(distanceB, distanceA, h) + smoothing * h * (1.0 - h);
}

fn sdfSmoothSubtract(distanceA: f32, distanceB: f32, smoothing: f32) -> f32 {
  let h = clamp(0.5 - 0.5 * (distanceB + distanceA) / smoothing, 0.0, 1.0);
  return mix(distanceA, -distanceB, h) + smoothing * h * (1.0 - h);
}

fn sdfSmoothUnion(distanceA: f32, distanceB: f32, smoothing: f32) -> f32 {
  let h = clamp(0.5 + 0.5 * (distanceB - distanceA) / smoothing, 0.0, 1.0);
  return mix(distanceB, distanceA, h) - smoothing * h * (1.0 - h);
}

fn sdfBox(p: vec2<f32>, b: vec2<f32>) -> f32 {
  let d = abs(p) - b;
  return length(max(d, vec2(0.0))) + min(max(d.x, d.y), 0.0);
}

fn sdfBoxFrame(position: vec3<f32>, size: vec3<f32>, thickness: f32) -> f32 {
  let q = abs(position) - size;
  let w = abs(q + thickness) - thickness;
  return min(min(
    length(max(vec3<f32>(q.x, w.y, w.z), vec3<f32>(0.0))) + min(max(q.x, max(w.y, w.z)), 0.0),
    length(max(vec3<f32>(w.x, q.y, w.z), vec3<f32>(0.0))) + min(max(w.x, max(q.y, w.z)), 0.0)),
    length(max(vec3<f32>(w.x, w.y, q.z), vec3<f32>(0.0))) + min(max(w.x, max(w.y, q.z)), 0.0));
}

fn sdfCappedTorus(position: vec3<f32>, majorRadius: f32, minorRadius: f32, angle: f32) -> f32 {
  let sc = vec2<f32>(sin(angle), cos(angle));
  let q = vec3<f32>(abs(position.x), position.y, position.z);
  let k = select(
    length(q.xy), 
    dot(q.xy, sc),
    sc.y * q.x > sc.x * q.y
  );
  return sqrt(dot(q, q) + 
              majorRadius * majorRadius - 
              2.0 * majorRadius * k) - 
              minorRadius;
}

fn sdfCapsule(position: vec3<f32>, radius: f32, height: f32) -> f32 {
  let d = abs(length(position.xz)) - radius;
  let p = vec2<f32>(d, abs(position.y) - height * 0.5);
  return length(max(p, vec2<f32>(0.0))) + min(max(p.x, p.y), 0.0) - radius;
}

fn sdfCircle(p: vec2<f32>, r: f32) -> f32 {
  return length(p) - r;
}

fn sdfCone(position: vec3<f32>, radius: f32, height: f32) -> f32 {
  let q = vec2<f32>(length(position.xz), position.y);
  let h = height;
  let r = radius;
  
  // Calculate distance
  let d1 = -q.y - h;
  let d2 = max(q.x * h - q.y * r, q.y * h + q.x * r);
  
  return length(max(vec2<f32>(d1, d2), vec2<f32>(0.0))) + min(max(d1, d2), 0.0);
}

fn sdfCylinder(position: vec3<f32>, radius: f32, height: f32) -> f32 {
  let d = vec2<f32>(length(position.xz), abs(position.y)) - vec2<f32>(radius, height * 0.5);
  return min(max(d.x, d.y), 0.0) + length(max(d, vec2<f32>(0.0)));
}

fn sdfEllipsoid(position: vec3<f32>, radius: vec3<f32>) -> f32 {
  let k0 = length(position / radius);
  let k1 = length(position / (radius * radius));
  return k0 * (k0 - 1.0) / k1;
}

fn sdfGyroid(position: vec3<f32>, scale: f32, thickness: f32) -> f32 {
  let p = position * scale;
  return (abs(dot(sin(p), cos(p.zxy))) - thickness) / scale;
}

fn sdfHexagonalPrism(position: vec3<f32>, radius: f32, height: f32) -> f32 {
  // Project into 2D
  var p = abs(position);
  let k = vec3<f32>(-0.866025404, 0.5, 0.577350269);
  
  // Hexagon in xy-plane
  p = vec3<f32>(p.x + p.y * k.x, p.y * k.y, p.z);
  p = vec3<f32>(p.x - min(p.x, p.y), p.y, p.z);
  let d = vec2<f32>(length(vec2<f32>(p.x, p.y - radius * k.z)) - radius, abs(p.z) - height * 0.5);
  
  return min(max(d.x, d.y), 0.0) + length(max(d, vec2<f32>(0.0)));
}

fn sdfIcosahedron(position: vec3<f32>, size: f32) -> f32 {
  var p = position;
  let s = size;
  
  // Constants for icosahedron
  let phi = 1.618033988749895;
  let a = s;
  let b = s * phi;
  
  // Compute distance to icosahedron
  p = abs(p / s);
  let d = p.x * p.y * p.z;
  let m = max(max(p.x, p.y), p.z);
  let n = min(min(p.x, p.y), p.z);
  let mid = p.x + p.y + p.z - m - n;
  
  // Calculate the signed distance
  let q = select(mid, d, m < phi * n);
  return (length(p) - phi) * s;
}

fn sdfIntersection(d1: f32, d2: f32) -> f32 {
  return max(d1, d2);
}

fn sdfJulia(position: vec3<f32>, c: vec4<f32>, iterations: f32, bailout: f32) -> f32 {
  var z = vec4<f32>(position, 0.0);
  var dz = vec4<f32>(1.0, 0.0, 0.0, 0.0);
  var m = dot(z, z);
  var i = 0;
  
  // Julia set iteration
  for (i = 0; i < i32(iterations) && m < bailout * bailout; i += 1) {
    // Quaternion multiplication for dz = 2.0 * z * dz
    dz = 2.0 * vec4<f32>(
      z.x * dz.x - z.y * dz.y - z.z * dz.z - z.w * dz.w,
      z.x * dz.y + z.y * dz.x + z.z * dz.w - z.w * dz.z,
      z.x * dz.z - z.y * dz.w + z.z * dz.x + z.w * dz.y,
      z.x * dz.w + z.y * dz.z - z.z * dz.y + z.w * dz.x
    );
    
    // Quaternion multiplication for z = z * z + c
    z = vec4<f32>(
      z.x * z.x - z.y * z.y - z.z * z.z - z.w * z.w,
      z.x * z.y + z.y * z.x + z.z * z.w - z.w * z.z,
      z.x * z.z - z.y * z.w + z.z * z.x + z.w * z.y,
      z.x * z.w + z.y * z.z - z.z * z.y + z.w * z.x
    ) + c;
    
    m = dot(z, z);
  }
  
  // Compute the distance
  let dist = 0.5 * log(m) * sqrt(m) / length(dz);
  return dist;
}

fn sdfOctahedron(position: vec3<f32>, size: f32) -> f32 {
  let p = abs(position);
  let m = p.x + p.y + p.z - size;
  
  // Calculate the distance
  var q: vec3<f32>;
  if (3.0 * p.x < m) {
    q = p;
  } else if (3.0 * p.y < m) {
    q = vec3<f32>(p.x, p.z, p.y);
  } else if (3.0 * p.z < m) {
    q = vec3<f32>(p.x, p.y, p.z);
  } else {
    q = p;
  }
  
  let k = clamp(0.5 * (q.z - q.y + size), 0.0, size);
  return length(vec3<f32>(q.x, q.y - size + k, q.z - k));
}

fn sdfPlane(position: vec3<f32>, normal: vec3<f32>) -> f32 {
  let n = normalize(normal);
  return dot(position, n);
}

fn sdfPyramid(position: vec3<f32>, size: f32, height: f32) -> f32 {
  // Normalize position
  var p = position;
  let h = height;
  let m2 = h * h + size * size;
  
  // Project into 2D
  let q = abs(p);
  p.y -= h;
  p.y = max(p.y, 0.0);
  
  // Distance calculation
  var d: f32;
  if (max(q.x, q.z) < size) {
    d = length(vec2<f32>(length(p.xz), p.y)) - sqrt(m2);
  } else {
    d = length(vec2<f32>(length(max(abs(p.xz) - vec2<f32>(size), vec2<f32>(0.0))), p.y));
  }
  
  // Account for position below base
  d = select(d, length(p) - sqrt(m2), p.y < 0.0);
  
  return d;
}

fn sdfRhombus(position: vec3<f32>, dimensions: vec3<f32>, sharpness: f32) -> f32 {
  var p = abs(position);
  let b = dimensions;
  let e = sharpness;
  
  // Calculate distance to rhombus
  p = p - b;
  let q = abs(p.x + p.y + p.z) + e;
  let h = max(vec3<f32>(q) - vec3<f32>(e), vec3<f32>(0.0));
  
  return min(max(p.x, max(p.y, p.z)), 0.0) + length(h);
}

fn sdfRoundBox(position: vec3<f32>, size: vec3<f32>, radius: f32) -> f32 {
  let q = abs(position) - size;
  return length(max(q, vec3<f32>(0.0))) + 
         min(max(q.x, max(q.y, q.z)), 0.0) - 
         radius;
}

fn sdfRoundedCone(position: vec3<f32>, radius1: f32, radius2: f32, height: f32, roundness: f32) -> f32 {
  // Calculate distances
  let p = position;
  let r1 = radius1 - roundness;
  let r2 = radius2 - roundness;
  let h = height;
  
  // Squared distance from axis
  let q = length(p.xz);
  
  // Project into 2D space
  let k1 = (r2 - r1) / h;
  let k2 = h / (r1 - r2);
  let projected = vec2<f32>(q - r1 + r1 * (p.y / h) * (r1 - r2) / r1, p.y - h);
  let ca = p.y * k1 - q;
  let cb = p.y - r1 * k2 + q * k2;
  
  var s: f32;
  if (ca < 0.0 && projected.y < 0.0) {
    s = length(projected) - roundness;
  } else if (ca > 0.0 && cb < 0.0) {
    s = -ca - roundness;
  } else {
    s = length(vec2<f32>(max(ca, 0.0), max(projected.y, 0.0))) - roundness;
  }
  
  return s;
}

fn sdfRoundedCylinder(position: vec3<f32>, radius: f32, height: f32, roundness: f32) -> f32 {
  // Calculate distances
  let radiusOffset = radius - roundness;
  let heightOffset = height * 0.5 - roundness;
  
  // Generate rounded cylinder
  let d = vec2<f32>(length(position.xz) - radiusOffset, abs(position.y) - heightOffset);
  return min(max(d.x, d.y), 0.0) + length(max(d, vec2<f32>(0.0))) - roundness;
}

fn sdfSphere(position: vec3<f32>, radius: f32) -> f32 {
  return length(position) - radius;
}

fn sdfSubtraction(d1: f32, d2: f32) -> f32 {
  return max(-d1, d2);
}

fn sdfTetrahedron(position: vec3<f32>, size: f32) -> f32 {
  var p = position;
  let s = size;
  
  // Set initial values
  let signVal = sign(p.x + p.y + p.z);
  p.x = abs(p.x);
  p.y = abs(p.y);
  p.z = abs(p.z);
  
  // Calculate the distance
  if (p.x < p.y) {
    let t = p.x;
    p.x = p.y;
    p.y = t;
  }
  if (p.x < p.z) {
    let t = p.x;
    p.x = p.z;
    p.z = t;
  }
  if (p.y < p.z) {
    let t = p.y;
    p.y = p.z;
    p.z = t;
  }
  
  let k = clamp((p.x + p.z - p.y) * 0.5, 0.0, p.z);
  return signVal * (length(vec3<f32>(p.x, p.y - s, p.z - k)) - s);
}

fn sdfTorus(position: vec3<f32>, majorRadius: f32, minorRadius: f32) -> f32 {
  let q = vec2<f32>(length(position.xz) - majorRadius, position.y);
  return length(q) - minorRadius;
}

fn sdfTriangularPrism(position: vec3<f32>, radius: f32, height: f32) -> f32 {
  var q = abs(position);
  
  // Triangle distance in xy-plane
  let k = sqrt(3.0);
  q.x = abs(q.x - q.y * k * 0.5);
  q.y = q.y * 0.866025404 + q.x * 0.5;
  
  // Combine with z distance
  let d1 = vec2<f32>(q.x - radius, q.y);
  let d2 = vec2<f32>(q.y - radius, q.x);
  let d = min(d1, d2);
  
  // Account for height
  let h = height * 0.5;
  let dz = q.z - h;
  let dz2 = max(dz, 0.0);
  
  return length(max(vec2<f32>(max(d.x, 0.0), dz2), vec2<f32>(0.0))) + min(max(d.x, dz), 0.0);
}

fn sdfUnion(d1: f32, d2: f32) -> f32 {
  return min(d1, d2);
}

fn sdfCylindricalRepeat(position: vec3<f32>, angleRepeat: f32, heightRepeat: f32, axis: vec3<f32>) -> vec3<f32> {
  let n = normalize(axis);
  
  // Project onto axis for height
  let h = dot(position, n);
  let radial = position - h * n;
  
  // Repeat in height
  let newH = h - heightRepeat * round(h / heightRepeat);
  
  // Repeat in angle
  let radius = length(radial);
  if (radius < 0.001) {
    return newH * n;
  }
  
  let angle = atan2(radial.y, radial.x);
  let sectorAngle = 6.28318530718 / angleRepeat;
  let newAngle = angle - sectorAngle * round(angle / sectorAngle);
  
  let newRadial = radius * vec2<f32>(cos(newAngle), sin(newAngle));
  
  // This assumes axis is along Z - for general axis, need proper basis vectors
  return newH * n + vec3<f32>(newRadial.x, newRadial.y, 0.0);
}

fn sdfMirror(position: vec3<f32>, normal: vec3<f32>, offset: f32) -> vec3<f32> {
  let n = normalize(normal);
  let d = dot(position, n) - offset;
  return position - 2.0 * max(0.0, d) * n;
}

fn sdfPolarRepeat(position: vec3<f32>, count: f32, axis: vec3<f32>) -> vec3<f32> {
  let n = normalize(axis);
  
  // Project position onto axis
  let axisProj = dot(position, n) * n;
  let radial = position - axisProj;
  
  // Get angle in the plane perpendicular to axis
  let radius = length(radial);
  if (radius < 0.001) {
    return position;
  }
  
  let angle = atan2(radial.y, radial.x);
  let sectorAngle = 6.28318530718 / count;
  let snappedAngle = round(angle / sectorAngle) * sectorAngle;
  
  // Reconstruct position with snapped angle
  let newRadial = radius * vec2<f32>(cos(snappedAngle), sin(snappedAngle));
  
  // This assumes axis is along Z - for general axis, need proper basis vectors
  return axisProj + vec3<f32>(newRadial.x, newRadial.y, 0.0);
}

fn sdfRotate(position: vec3<f32>, angles: vec3<f32>, pivot: vec3<f32>) -> vec3<f32> {
  // First translate to origin relative to pivot point
  let centered = position - pivot;
  
  // Create rotation matrices (inverse rotation = negative angles)
  let cx = cos(-angles.x);
  let sx = sin(-angles.x);
  let cy = cos(-angles.y);
  let sy = sin(-angles.y);
  let cz = cos(-angles.z);
  let sz = sin(-angles.z);
  
  // Rotate around X axis
  let rx = vec3<f32>(
    centered.x,
    centered.y * cx - centered.z * sx,
    centered.y * sx + centered.z * cx
  );
  
  // Rotate around Y axis
  let ry = vec3<f32>(
    rx.x * cy + rx.z * sy,
    rx.y,
    -rx.x * sy + rx.z * cy
  );
  
  // Rotate around Z axis
  let rz = vec3<f32>(
    ry.x * cz - ry.y * sz,
    ry.x * sz + ry.y * cz,
    ry.z
  );
  
  // Translate back from pivot point
  return rz + pivot;
}

fn sdfScale(position: vec3<f32>, scale: vec3<f32>) -> vec3<f32> {
  return position / scale;
}

fn sdfSphericalRepeat(position: vec3<f32>, phiRepeat: f32, thetaRepeat: f32) -> vec3<f32> {
  let radius = length(position);
  if (radius < 0.001) {
    return position;
  }
  
  // Convert to spherical coordinates
  let theta = acos(clamp(position.z / radius, -1.0, 1.0));
  let phi = atan2(position.y, position.x);
  
  // Repeat in spherical coordinates
  let phiSector = 6.28318530718 / phiRepeat;
  let thetaSector = 3.14159265359 / thetaRepeat;
  
  let newPhi = phi - phiSector * round(phi / phiSector);
  let newTheta = theta - thetaSector * round(theta / thetaSector);
  
  // Convert back to Cartesian
  let sinTheta = sin(newTheta);
  return radius * vec3<f32>(
    sinTheta * cos(newPhi),
    sinTheta * sin(newPhi),
    cos(newTheta)
  );
}

fn sdfTranslate(position: vec3<f32>, offset: vec3<f32>) -> vec3<f32> {
  return position - offset;
}

fn sdfToSmoothSolid(signedDistance: f32, threshold: f32, smoothing: f32) -> f32 {
  return 1.0 - smoothstep(threshold - smoothing, threshold + smoothing, signedDistance);
}

fn sdfToSmoothStroke(signedDistance: f32, thickness: f32, center: f32, smoothing: f32) -> f32 {
  let distance = abs(signedDistance - center);
  return 1.0 - smoothstep(thickness * 0.5 - smoothing, thickness * 0.5 + smoothing, distance);
}