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

fn sdEllipsoid(p: vec3<f32>, r: vec3<f32>) -> f32 {
  let k0 = length(p / r);
  let k1 = length(p / (r * r));
  return k0 * (k0 - 1.) / k1;
}

fn sdBox(p: vec3<f32>, b: vec3<f32>) -> f32 {
  let q = abs(p) - b;
  return length(max(q, vec3<f32>(0.))) + min(max(q.x, max(q.y, q.z)), 0.);
}

fn sdRoundBox(p: vec3<f32>, b: vec3<f32>, r: f32) -> f32 {
  let q = abs(p) - b;
  return length(max(q, vec3<f32>(0.))) + min(max(q.x,max(q.y, q.z)), 0.) - r;
}

fn sdBoxFrame(p: vec3<f32>, b: vec3<f32>, e: f32) -> f32 {
  let q = abs(p) - b;
  let w = abs(q + e) - e;
  return min(min(
      length(max(vec3<f32>(q.x, w.y, w.z), vec3<f32>(0.))) + min(max(q.x, max(w.y, w.z)), 0.),
      length(max(vec3<f32>(w.x, q.y, w.z), vec3<f32>(0.))) + min(max(w.x, max(q.y, w.z)), 0.)),
      length(max(vec3<f32>(w.x, w.y, q.z), vec3<f32>(0.))) + min(max(w.x, max(w.y, q.z)), 0.));
}

fn sdGyroid(p: vec3<f32>, h: f32) -> f32 {
  return abs(dot(sin(p), cos(p.zxy))) - h;
}

fn sdTorus(p: vec3<f32>, R: f32, r: f32) -> f32 {
  let q = vec2<f32>(length(p.xz) - R, p.y);
  return length(q) - r;
}

fn sdCappedTorus(p: vec3<f32>, R: f32, r: f32, sincos: vec2<f32>) -> f32 {
  let q = vec3<f32>(abs(p.x), p.y, p.z);
  let k = select(dot(q.xy, sincos), length(q.xy), sincos.y * q.x > sincos.x * q.y);
  return sqrt(dot(q, q) + R * R - 2. * R * k) - r;
}

fn sdLink(p: vec3<f32>, R: f32, r: f32, le: f32) -> f32 {
  let q = vec3<f32>(p.x, max(abs(p.y) - le, 0.), p.z);
  return length(vec2<f32>(length(q.xy) - R, q.z)) - r;
}

fn sdCapsule(p: vec3<f32>, a: vec3<f32>, b: vec3<f32>, r: f32) -> f32 {
  let pa = p - a;
  let ba = b - a;
  let h = clamp(dot(pa, ba) / dot(ba, ba), 0., 1.);
  return length(pa - ba * h) - r;
}

fn sdVerticalCapsule(p: vec3<f32>, h: f32, r: f32) -> f32 {
  let q = vec3<f32>(p.x, p.y - clamp(p.y, 0., h), p.z);
  return length(q) - r;
}

fn sdCylinder(p: vec3<f32>, a: vec3<f32>, b: vec3<f32>, r: f32) -> f32 {
  let ba = b - a;
  let pa = p - a;
  let baba = dot(ba, ba);
  let paba = dot(pa, ba);
  let x = length(pa * baba - ba * paba) - r * baba;
  let y = abs(paba - baba * 0.5) - baba * 0.5;
  let x2 = x * x;
  let y2 = y * y * baba;
  let d = x2 * step(0., x) + y2 * step(0., y);
  let d2 = select(-min(x2, y2), d, max(x, y) < 0.);
  return sign(d2) * sqrt(abs(d2)) / baba;
}

fn sdVerticalCylinder(p: vec3<f32>, h: f32, r: f32) -> f32 {
  let d = abs(vec2<f32>(length(p.xz), p.y)) - vec2<f32>(r, h);
  return min(max(d.x, d.y), 0.) + length(max(d, vec2<f32>(0.)));
}

fn sdRoundedCylinder(p: vec3<f32>, h: f32, r: f32, re: f32) -> f32 {
  let d = vec2<f32>(length(p.xz) - 2. * r + re, abs(p.y) - h);
  return min(max(d.x, d.y), 0.) + length(max(d, vec2<f32>(0.))) - re;
}

fn sdInfiniteCylinder(p: vec3<f32>, c: vec3<f32>) -> f32 {
  return length(p.xz - c.xy) - c.z;
}

fn sdCone(p: vec3<f32>, h: f32, sincos: vec2<f32>) -> f32 {
  // Alternatively pass q instead of (sin(alpha), cos(alpha))
  let q = h * vec2<f32>(sincos.x / sincos.y, -1.);

  let w = vec2<f32>(length(p.xz), p.y);
  let a = w - q * clamp(dot(w,q) / dot(q,q), 0., 1.);
  let b = w - q * vec2<f32>(clamp(w.x / q.x, 0., 1.), 1.);
  let k = sign(q.y);
  let d = min(dot(a, a), dot(b, b));
  let s = max(k * (w.x * q.y - w.y * q.x), k * (w.y - q.y));
  return sqrt(d) * sign(s);
}

fn sdConeBound(p: vec3<f32>, h: f32, sincos: vec2<f32>) -> f32 {
  return max(dot(sincos.yx, vec2<f32>(length(p.xz), p.y)), -h - p.y);
}

fn sdInfiniteCone(p: vec3<f32>, sincos: vec2<f32>) -> f32 {
  let q = vec2<f32>(length(p.xz), -p.y);
  let d = length(q - sincos * max(dot(q, sincos), 0.));
  return d * select(1., -1., q.x * sincos.y - q.y * sincos.x > 0.0);
}

fn sdCappedVerticalCone(p: vec3<f32>, h: f32, r1: f32, r2: f32) -> f32 {
  let q = vec2<f32>(length(p.xz), p.y);
  let k1 = vec2<f32>(r2, h);
  let k2 = vec2<f32>(r2 - r1, 2. * h);
  let ca = vec2<f32>(q.x - min(q.x, select(r1, r2, q.y < 0.)), abs(q.y) - h);
  let cb = q - k1 + k2 * clamp(dot(k1 - q, k2) / dot(k2, k2), 0., 1.);
  let s = select(-1., 1., cb.x < 0. && ca.y < 0.);
  return s * sqrt(min(dot(ca, ca), dot(cb, cb)));
}

fn sdCappedCone(p: vec3<f32>, a: vec3<f32>, b: vec3<f32>, ra: f32, rb: f32) -> f32 {
  let rba = rb - ra;
  let baba = dot(b - a, b - a);
  let papa = dot(p - a, p - a);
  let paba = dot(p - a, b - a) / baba;
  let x = sqrt(papa - paba * paba * baba);
  let cax = max(0.0, x - select(ra, rb, paba < 0.5));
  let cay = abs(paba - 0.5) - 0.5;
  let k = rba * rba + baba;
  let f = clamp((rba * (x - ra) + paba * baba) / k, 0.0, 1.0);
  let cbx = x - ra - f * rba;
  let cby = paba - f;
  let s = select(-1., 1., cbx < 0.0 && cay < 0.0);
  return s * sqrt(min(cax * cax + cay * cay * baba, cbx * cbx + cby * cby * baba));
}

fn sdRoundVerticalCone(p: vec3<f32>, h: f32, r1: f32, r2: f32) -> f32 {
  let q = vec2<f32>(length(p.xz), p.y);
  let b = (r1 - r2) / h;
  let a = sqrt(1. - b * b);
  let k = dot(q, vec2<f32>(-b, a));
  if (k < 0.) { return length(q) - r1; }
  if (k > a * h) { return length(q - vec2<f32>(0., h)) - r2; }
  return dot(q, vec2<f32>(a, b)) - r1;
}

fn sdRoundCone(p: vec3<f32>, a: vec3<f32>, b: vec3<f32>, r1: f32, r2: f32) -> f32 {
  let ba = b - a;
  let l2 = dot(ba, ba);
  let rr = r1 - r2;
  let a2 = l2 - rr * rr;
  let il2 = 1. / l2;

  let pa = p - a;
  let y = dot(pa, ba);
  let z = y - l2;
  let w = pa * l2 - ba * y;
  let x2 = dot(w, w);
  let y2 = y * y * l2;
  let z2 = z * z * l2;

  let k = sign(rr) * rr * rr * x2;
  if (sign(z) * a2 * z2 > k) { return sqrt(x2 + z2) * il2 - r2; }
  if (sign(y) * a2 * y2 < k) { return sqrt(x2 + y2) * il2 - r1; }
  return (sqrt(x2 * a2 * il2) + y * rr) * il2 - r1;
}

fn sdSolidAngle(p: vec3<f32>, sincos: vec2<f32>, r: f32) -> f32 {
  let q = vec2<f32>(length(p.xz), p.y);
  let l = length(q) - r;
  let m = length(q - sincos * clamp(dot(q, sincos), 0., r));
  return max(l, m * sign(sincos.y * q.x - sincos.x * q.y));
}

fn sdPlane(p: vec3<f32>, n: vec3<f32>, h: f32) -> f32 {
  // n must be normalized
  return dot(p, n) + h;
}

fn sdOctahedron(p: vec3<f32>, s: f32) -> f32 {
  var q: vec3<f32> = abs(p);
  let m = q.x + q.y + q.z - s;
  if (3. * q.x < m) {q = q.xyz;}
  else {if (3. * q.y < m) {q = q.yzx;}
        else {if (3. * q.z < m) {q = q.zxy;}
              else {return m * 0.57735027;}}}
  let k = clamp(0.5 * (q.z - q.y + s), 0., s);
  return length(vec3<f32>(q.x, q.y - s + k, q.z - k));
}

fn sdOctahedronBound(p: vec3<f32>, s: f32) -> f32 {
  let q = abs(p);
  return (q.x + q.y + q.z - s) * 0.57735027;
}

fn sdPyramid(p: vec3<f32>, h: f32) -> f32 {
  let m2 = h * h + 0.25;
  var xz: vec2<f32> = abs(p.xz);
  xz = select(xz.yx, xz, xz[1] > xz[0]);
  xz = xz - vec2<f32>(0.5);

  let q = vec3<f32>(xz[1], h * p.y - 0.5 * xz[0], h * xz[0] + 0.5 * p.y);
  let s = max(-q.x, 0.);
  let t = clamp((q.y - 0.5 * xz[1]) / (m2 + 0.25), 0., 1.);

  let a = m2 * (q.x + s) * (q.x + s) + q.y * q.y;
  let b = m2 * (q.x + 0.5 * t) * (q.x + 0.5 * t) + (q.y - m2 * t) * (q.y - m2 * t);

  let d2 = min(a, b) * step(min(q.y, -q.x * m2 - q.y * 0.5), 0.);
  return sqrt((d2 + q.z * q.z) / m2) * sign(max(q.z, -p.y));
}

fn sdHexPrism(p: vec3<f32>, h: vec2<f32>) -> f32 {
  let k = vec3<f32>(-0.8660254, 0.5, 0.57735);
  let a = abs(p);
  let v = a.xy - 2. * min(dot(k.xy, a.xy), 0.) * k.xy;
  let d1 = length(v - vec2<f32>(clamp(v.x, -k.z * h.x, k.z * h.x), h.x)) * sign(v.y - h.x);
  let d2 = a.z - h.y;
  return min(max(d1, d2), 0.) + length(max(vec2<f32>(d1, d2), vec2<f32>(0.)));
}

fn sdTriPrism(p: vec3<f32>, h: vec2<f32>) -> f32 {
  let q = abs(p);
  return max(q.z - h.y, max(q.x * 0.866025 + p.y * 0.5, -p.y) - h.x * 0.5);
}

fn fmat4init(a1: f32, a2: f32, a3: f32, a4: f32, b1: f32, b2: f32, b3: f32, b4: f32,
             c1: f32, c2: f32, c3: f32, c4: f32, d1: f32, d2: f32, d3: f32, d4: f32) -> mat4x4<f32> {
  return mat4x4<f32>(vec4<f32>(a1, a2, a3, a4), vec4<f32>(b1, b2, b3, b4),
                     vec4<f32>(c1, c2, c3, c4), vec4<f32>(d1, d2, d3, d4));
}

fn bunny(p: vec3<f32>) -> f32 {
  if (dot(p, p) > 1.) { return length(p) - .8; }
  let q = vec4<f32>(p, 1.);
  let f00=sin(fmat4init(-3.02,1.95,-3.42,-.6,3.08,.85,-2.25,-.24,-.29,1.16,-3.74,2.89,-.71,4.5,-3.24,-3.5)*q);
  let f01=sin(fmat4init(-.4,-3.61,3.23,-.14,-.36,3.64,-3.91,2.66,2.9,-.54,-2.75,2.71,7.02,-5.41,-1.12,-7.41)*q);
  let f02=sin(fmat4init(-1.77,-1.28,-4.29,-3.2,-3.49,-2.81,-.64,2.79,3.15,2.14,-3.85,1.83,-2.07,4.49,5.33,-2.17)*q);
  let f03=sin(fmat4init(-.49,.68,3.05,.42,-2.87,.78,3.78,-3.41,-2.65,.33,.07,-.64,-3.24,-5.9,1.14,-4.71)*q);
  let f10=sin(fmat4init(-.34,.06,-.59,-.76,.1,-.19,-.12,.44,.64,-.02,-.26,.15,-.16,.21,.91,.15)*f00+
      fmat4init(.01,.54,-.77,.11,.06,-.14,.43,.51,-.18,.08,.39,.2,.33,-.49,-.1,.19)*f01+
      fmat4init(.27,.22,.43,.53,.18,-.17,.23,-.64,-.14,.02,-.1,.16,-.13,-.06,-.04,-.36)*f02+
      fmat4init(-.13,.29,-.29,.08,1.13,.02,-.83,.32,-.32,.04,-.31,-.16,.14,-.03,-.2,.39)*f03+
      vec4<f32>(.73,-4.28,-1.56,-1.8))+f00;
  let f11=sin(fmat4init(-1.11,.55,-.12,-1.00,.16,.15,-.3,.31,-.01,.01,.31,-.42,-.29,.38,-.04,.71)*f00+
      fmat4init(.96,-.02,.86,.52,-.14,.6,.44,.43,.02,-.15,-.49,-.05,-.06,-.25,-.03,-.22)*f01+
      fmat4init(.52,.44,-.05,-.11,-.56,-.1,-.61,-.4,-.04,.55,.32,-.07,-.02,.28,.26,-.49)*f02+
      fmat4init(.02,-.32,.06,-.17,-.59,.00,-.24,.6,-.06,.13,-.21,-.27,-.12,-.14,.58,-.55)*f03+
      vec4<f32>(-2.24,-3.48,-.8,1.41))+f01;
  let f12=sin(fmat4init(.44,-.06,-.79,-.46,.05,-.6,.3,.36,.35,.12,.02,.12,.4,-.26,.63,-.21)*f00+
      fmat4init(-.48,.43,-.73,-.4,.11,-.01,.71,.05,-.25,.25,-.28,-.2,.32,-.02,-.84,.16)*f01+
      fmat4init(.39,-.07,.9,.36,-.38,-.27,-1.86,-.39,.48,-.2,-.05,.1,-.00,-.21,.29,.63)*f02+
      fmat4init(.46,-.32,.06,.09,.72,-.47,.81,.78,.9,.02,-.21,.08,-.16,.22,.32,-.13)*f03+
      vec4<f32>(3.38,1.2,.84,1.41))+f02;
  let f13=sin(fmat4init(-.41,-.24,-.71,-.25,-.24,-.75,-.09,.02,-.27,-.42,.02,.03,-.01,.51,-.12,-1.24)*f00+
      fmat4init(.64,.31,-1.36,.61,-.34,.11,.14,.79,.22,-.16,-.29,-.70,.02,-.37,.49,.39)*f01+
      fmat4init(.79,.47,.54,-.47,-1.13,-.35,-1.03,-.22,-.67,-.26,.1,.21,-.07,-.73,-.11,.72)*f02+
      fmat4init(.43,-.23,.13,.09,1.38,-.63,1.57,-.2,.39,-.14,.42,.13,-.57,-.08,-.21,.21)*f03+
      vec4<f32>(-.34,-3.28,.43,-.52))+f03;
  let f20=sin(fmat4init(-.72,.23,-.89,.52,.38,.19,-.16,-.88,.26,-.37,.09,.63,.29,-.72,.3,-.95)*f10+
      fmat4init(-.22,-.51,-.42,-.73,-.32,.00,-1.03,1.17,-.2,-.03,-.13,-.16,-.41,.09,.36,-.84)*f11+
      fmat4init(-.21,.01,.33,.47,.05,.2,-.44,-1.04,.13,.12,-.13,.31,.01,-.34,.41,-.34)*f12+
      fmat4init(-.13,-.06,-.39,-.22,.48,.25,.24,-.97,-.34,.14,.42,-.00,-.44,.05,.09,-.95)*f13+
      vec4<f32>(.48,.87,-.87,-2.06))/1.4+f10;
  let f21=sin(fmat4init(-.27,.29,-.21,.15,.34,-.23,.85,-.09,-1.15,-.24,-.05,-.25,-.12,-.73,-.17,-.37)*f10+
      fmat4init(-1.11,.35,-.93,-.06,-.79,-.03,-.46,-.37,.6,-.37,-.14,.45,-.03,-.21,.02,.59)*f11+
      fmat4init(-.92,-.17,-.58,-.18,.58,.6,.83,-1.04,-.8,-.16,.23,-.11,.08,.16,.76,.61)*f12+
      fmat4init(.29,.45,.3,.39,-.91,.66,-.35,-.35,.21,.16,-.54,-.63,1.1,-.38,.2,.15)*f13+
      vec4<f32>(-1.72,-.14,1.92,2.08))/1.4+f11;
  let f22=sin(fmat4init(1.00,.66,1.3,-.51,.88,.25,-.67,.03,-.68,-.08,-.12,-.14,.46,1.15,.38,-.1)*f10+
      fmat4init(.51,-.57,.41,-.09,.68,-.5,-.04,-1.01,.2,.44,-.6,.46,-.09,-.37,-1.3,.04)*f11+
      fmat4init(.14,.29,-.45,-.06,-.65,.33,-.37,-.95,.71,-.07,1.00,-.6,-1.68,-.2,-.00,-.7)*f12+
      fmat4init(-.31,.69,.56,.13,.95,.36,.56,.59,-.63,.52,-.3,.17,1.23,.72,.95,.75)*f13+
      vec4<f32>(-.9,-3.26,-.44,-3.11))/1.4+f12;
  let f23=sin(fmat4init(.51,-.98,-.28,.16,-.22,-.17,-1.03,.22,.7,-.15,.12,.43,.78,.67,-.85,-.25)*f10+
      fmat4init(.81,.6,-.89,.61,-1.03,-.33,.6,-.11,-.06,.01,-.02,-.44,.73,.69,1.02,.62)*f11+
      fmat4init(-.1,.52,.8,-.65,.4,-.75,.47,1.56,.03,.05,.08,.31,-.03,.22,-1.63,.07)*f12+
      fmat4init(-.18,-.07,-1.22,.48,-.01,.56,.07,.15,.24,.25,-.09,-.54,.23,-.08,.2,.36)*f13+
      vec4<f32>(-1.11,-4.28,1.02,-.23))/1.4+f13;
  return dot(f20,vec4<f32>(.09,.12,-.07,-.03))+dot(f21,vec4<f32>(-.04,.07,-.08,.05))+
      dot(f22,vec4<f32>(-.01,.06,-.02,.07))+dot(f23,vec4<f32>(-.05,.07,.03,.04))- 0.16;
}

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

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

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

fn opUnionChamfer(d1: f32, d2: f32, r: f32) -> f32 {
  return min(min(d1, d2), (d1 - r + d2) * 0.5);
}

fn opSubtractionChamfer(d1: f32, d2: f32, r: f32) -> f32{
  return max(max(d1, -d2), (d1 + r - d2) * 0.5);
}

fn opIntersectionChamfer(d1: f32, d2: f32, r: f32) -> f32 {
  return max(max(d1, d2), (d1 + r + d2) * 0.5);
}

fn opUnionBlend(d1: f32, d2: f32, k: f32) -> f32 {
  let h = clamp(0.5 + 0.5 * (d2 - d1) / k, 0., 1.);
  return mix(d2, d1, h) - k * h * (1. - h);
}

fn opSubtractionBlend(d1: f32, d2: f32, k: f32) -> f32 {
  let h = clamp(0.5 - 0.5 * (d1 + d2) / k, 0., 1.);
  return mix(d1, -d2, h) + k * h * (1. - h);
}

fn opIntersectionBlend(d1: f32, d2: f32, k: f32) -> f32 {
  let h = clamp(0.5 - 0.5 * (d2 - d1) / k, 0., 1.);
  return mix(d2, d1, h) + k * h * (1. - h);
}

fn opDisplace(d1: f32, d2: f32) -> f32 {
  return d1 + d2;
}
//let d = opDisplace(sdfPrimitive3d(p), displacement3d(p));

fn opTwist(p: vec3<f32>, k: f32) -> vec3<f32> {
  let s = sin(k * p.y);
  let c = cos(k * p.y);
  let m = mat2x2<f32>(vec2<f32>(c, s), vec2<f32>(-s, c));
  return vec3<f32>(m * p.xz, p.y);
}
//let d = sdfPrimitive3d(opTwist(p, k));

fn opCheapBend(p: vec3<f32>, k: f32) -> vec3<f32> {
  let s = sin(k * p.x);
  let c = cos(k * p.x);
  let m = mat2x2<f32>(vec2<f32>(c, s), vec2<f32>(-s, c));
  return vec3<f32>(m * p.xy, p.z);
}
//let d = sdfPrimitive3d(opCheapBend(p, k));

fn opTranslate(p: vec3<f32>, t: vec3<f32>) -> vec3<f32> {
  return p - t;
}
//let d = sdfPrimitive3d(opTranslate(p, t));

fn op90RotateX(p: vec3<f32>) -> vec3<f32> {
  return vec3<f32>(p.x, p.z, -p.y);
}

fn op90RotateY(p: vec3<f32>) -> vec3<f32> {
  let r = sqrt(0.5);
  return vec3<f32>(-p.z, p.y, p.x);
}

fn op90RotateZ(p: vec3<f32>) -> vec3<f32> {
  return vec3<f32>(p.y, -p.x, p.z);
}
//let d = sdfPrimitive3d(op90RotateZ(p));

fn opRotateX(p: vec3<f32>, a: f32) -> vec3<f32> {
  let s = sin(a); let c = cos(a);
  return vec3<f32>(p.x, c * p.y + s * p.z, -s * p.y + c * p.z);
}

fn opRotateY(p: vec3<f32>, a: f32) -> vec3<f32> {
  let s = sin(a); let c = cos(a);
  return vec3<f32>(c * p.x - s * p.z, p.y, s * p.x + c * p.z);
}

fn opRotateZ(p: vec3<f32>, a: f32) -> vec3<f32> {
  let s = sin(a); let c = cos(a);
  return vec3<f32>(c * p.x + s * p.y, -s * p.x + c * p.y, p.z);
}
//let d = sdfPrimitive3d(opRotateY(p, a));

fn opRotateE(p: vec3<f32>, e: vec3<f32>, a: f32) -> vec3<f32> {
  let c = cos(a);
  return dot(e, p) * (1. - c) * e - cross(e, p) * sin(a) + c * p;
}
//let d = sdfPrimitive3d(opRotateE(p, normalize(vec3<f32>(1.,0.,.5)), a));

fn opScale(p: vec3<f32>, s: f32) -> vec3<f32> {
  return p / s;
}
//let d = sdfPrimitive3d(opScale(p, s)) * s;

fn opTransform(p: vec3<f32>, transform: mat4x4<f32>) -> vec3<f32> {
  let q = inverse(transform) * vec4<f32>(p, 1.);
}
//let d = sdfPrimitive3d(opTransform(p, transform)) * determinant(transform);

// OR
//let d = sdfPrimitive3d(opScale(opRotateE(opTranslate(p, t), e, a), s)) * s;

fn opSymmetryX(p: vec3<f32>) -> vec3<f32> { return vec3<f32>(abs(p.x), p.y, p.z); }

fn opSymmetryY(p: vec3<f32>) -> vec3<f32> {
  return vec3<f32>(p.x, abs(p.y), p.z);
}

fn opSymmetryZ(p: vec3<f32>) -> vec3<f32> {
  return vec3<f32>(p.x, p.y, abs(p.z));
}
//let d = sdfPrimitive3d(opSymmetryX(p));

fn opInfArray(p: vec3<f32>, c: vec3<f32>) -> vec3<f32> {
  return (p + 0.5 * c % c) - 0.5 * c;
}
//let d = sdfPrimitive3d(opInfArray(p, c));

fn opLimArray(p: vec3<f32>, c: f32, lim: vec3<f32>) -> vec3<f32> {
  return p - c * clamp(round(p / c), -lim, lim);
}
//let d = sdfPrimitive3d(opLimArray(p, c));

fn opElongate(p: vec3<f32>, h: vec3<f32>) -> vec3<f32> {
  return p - clamp(p, -h, h);
}
//let d = sdfPrimitive3d(opElongateFast(p, h));

fn opElongateCorrect(p: vec3<f32>, h: vec3<f32>) -> vec4<f32> {
  let q = abs(p) - h;
  let sgn = 2. * step(vec3<f32>(0.), p) - vec3<f32>(1.);
  return vec4<f32>(sgn * max(q, vec3<f32>(0.)), min(max(q.x, max(q.y, q.z)), 0.));
}
//let p2 = opElongateCorrect(p, h);
//let d = p2.w + sdfPrimitive3d(p2.xyz);

fn opRound(p: vec3<f32>, r: f32) -> f32 {
  return sdfPrimitive3d(p) - r;
}

fn opOnion(d: f32, thickness: f32) -> f32 {
  return abs(d) - thickness;
}
//let d = opOnion(sdfPrimitive3d(p), thickness);

fn opExtrusion(d: f32, z: f32, h: f32) -> f32 {
  let w = vec2<f32>(d, abs(z) - h);
  return min(max(w.x, w.y), 0.) + length(max(w, vec2<f32>(0.)));
}
//let d = opExtrusion(sdfPrimitive2d(p.xy), p.z, h));

fn opRevolution(p: vec3<f32>, o: f32) -> vec2<f32> {
  return vec2<f32>(length(p.xz) - o, p.y);
}
//let d = sdfPrimitive2d(opRevolution(p, h));

fn length4(p: vec3<f32>) -> f32 {
  var q: vec3<f32> = p * p;
  q = q * q;
  return sqrt(sqrt(q.x + q.y + q.z));
}

fn length6(p: vec3<f32>) -> f32 {
  var q: vec3<f32> = p * p * p;
  q = q * q;
  return pow(q.x + q.y + q.z, 1. / 6.);
}

fn length8(p: vec3<f32>) -> f32 {
  var q: vec3<f32> = p * p;
  q = q * q; q = q * q;
  return pow(q.x + q.y + q.z, 1. / 8.);
}