export declare const solveWGSL = "\n\n// Shared by every 3D module (avbd3d/gpu/layout.ts)\n\nconst T_NONE = 0;\n\nconst MAX_COLORS = 64u;\nconst NO_COLOR = 255u;\nconst PENDING = 256u;\n\nconst C_PAIRS = 0u;\nconst C_CONTACTS = 1u;\nconst C_PREV_CONTACTS = 2u;\nconst C_OVERFLOW = 3u;\nconst C_CLASHES = 4u;\nconst C_NUM_COLORS = 5u;\n\nconst IA_PAIRS = 0u;\nconst IA_CONTACTS = 3u;\nconst IA_PREV = 6u;\nconst IA_CONSTRAINTS = 9u;\nconst IA_COLOR = 12u;\n\nconst WG = 64u;\nconst COLOR_WG = 256u;\n\nconst BIG = 3e38;\nconst HARD = 1e30;\n\nfn groupsFor(n: u32) -> u32 {\n  return (n + WG - 1u) / WG;\n}\n\n/** lowbias32 integer hash (same as the CPU colouring priorities). */\nfn hash32(v: u32) -> u32 {\n  var x = v;\n  x ^= x >> 16u;\n  x *= 0x7feb352du;\n  x ^= x >> 15u;\n  x *= 0x846ca68bu;\n  x ^= x >> 16u;\n  return x;\n}\n\nconst T_JOINT = 1;\nconst T_SPRING = 2;\n\nconst PENALTY_MIN = 1.0;\nconst PENALTY_MAX = 1e10;\n// Least LDL\u1D40 pivot of the primal solve, relative to its diagonal (wgsl-solve.ts finishBody)\nconst PIVOT_FLOOR = 1e-5;\nconst COLLISION_MARGIN = 0.01;\nconst STICK_THRESH = 0.00001;\n\nconst SHAPE_SPHERE = 1.0;\nconst SHAPE_SAIL = 2.0;\nconst SHAPE_HULL = 3.0;\nconst STICK_BIT = 2147483648u;\nconst C_MANIFOLDS = 6u;\nconst C_PREV_MANIFOLDS = 7u;\n/** Feature key of every contact involving a sphere (one contact per pair). */\nconst SPHERE_FEATURE = 3u << 24u;\n\nconst FLAG_MATCH_NEAREST = 1u;\nconst FLAG_FACE_BIAS = 2u;\nconst FLAG_REUSE_CONTACTS = 4u;\nconst FLAG_START_AT_REST = 8u;\nconst FLAG_MASS_PENALTY = 16u;\nconst REST_SPEED = 0.1;\nconst NEAREST_FRACTION = 0.05;\n\nstruct Body {\n  pos: vec4f,          // xyz, w: friction\n  rot: vec4f,          // orientation quaternion (x, y, z, w)\n  initialPos: vec4f,   // pose at the start of the step (x-)\n  initialRot: vec4f,\n  size: vec4f,         // full widths, w: mass (0 = static)\n  moment: vec4f,       // principal moments, w: bounding radius\n  inertialPos: vec4f,  // inertial target y, w: step it last moved from its reference pose (u32)\n  inertialRot: vec4f,\n  vel: vec4f,          // xyz, w: previous step's velocity along up (adaptive warm start)\n  angVel: vec4f,       // xyz, w: shape (SHAPE_BOX, SHAPE_SPHERE, SHAPE_SAIL, or SHAPE_HULL + hull offset)\n}\n\nstruct Joint {\n  penLin: vec4f,  // xyz, w: linear stiffness (>= HARD: hard); spring: its stiffness\n  penAng: vec4f,  // xyz, w: angular stiffness\n  lamLin: vec4f,  // xyz, w: fracture threshold on |lamAng| (negative: on |lamLin| too)\n  lamAng: vec4f,  // xyz, w: torque arm\n  c0Lin: vec4f,   // C(x-)\n  c0Ang: vec4f,   // xyz: C(x-), w: yield threshold on |lamAng| (BIG: never yields)\n  rA: vec4f,      // xyz: anchor on A (world point when A is the world), w: spring rest length\n  rB: vec4f,      // xyz: anchor on B, w: bend limit in radians (BIG: a plastic joint never breaks on it)\n  rest: vec4f,    // the relative rotation the angle lock holds: rotB = rotA\u00B7rest (quaternion x, y, z, w)\n  restStart: vec4f, // the rest it was made with: its bend is the turn from this to rest\n}\n\nstruct Manifold {\n  ids: vec4u,     // bodyA, bodyB (A > B), first contact, count | generation step << 4\n  geo: vec4f,     // normal (B to A), friction\n}\n\nfn pairCount(m: Manifold) -> u32 {\n  return m.ids.w & 15u;\n}\n\nstruct Contact {\n  rA: vec3f,      // body-local anchors\n  key: u32,       // feature | STICK_BIT\n  rB: vec3f,\n  c0x: f32,       // C(x-) in the pair's basis, including the collision margin\n  pen: vec3f,     // normal, tangent, tangent\n  c0y: f32,\n  lam: vec3f,\n  c0z: f32,\n}\n\nstruct Params {\n  dt: f32,\n  gravity: f32,         // acceleration along up (m/s\u00B2, negative: down)\n  betaLin: f32,\n  betaAng: f32,\n  gamma: f32,\n  alpha: f32,\n  flags: u32,\n  bodyCount: u32,\n  jointCount: u32,\n  colorCap: u32,\n  cellSize: f32,\n  tableMask: u32,\n  maxSmallRadius: f32,\n  largeCount: u32,\n  noCollideCount: u32,\n  pairCapacity: u32,\n  contactCapacity: u32,\n  hashMask: u32,\n  adjFillOffset: u32,\n  adjListOffset: u32,\n  gridCursorOffset: u32,\n  gridSortedOffset: u32,\n  gridCellOffset: u32,\n  stateBOffset: u32,\n  colorHistOffset: u32,\n  colorStartOffset: u32,\n  colorGroups: u32,\n  colorBodiesOffset: u32,\n  rounds: u32,\n  manifoldCapacity: u32,\n  step: u32,           // steps since the solver was created (from 1)\n  primalLanes: u32,    // bytes 0-2: log2 of the colour sizes from which 1, 2, 4 lanes per body suffice\n  wind: vec4f,         // xyz: wind velocity (m/s), w: air pressure coefficient \u00BD\u03C1C_d\n  gust: vec4f,         // x: gust strength (fraction of the wind)\n  up: vec4f,           // xyz: which way is up (unit)\n}\n\nfn qmul(a: vec4f, b: vec4f) -> vec4f {\n  return vec4f(\n    a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,\n    a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,\n    a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w,\n    a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z);\n}\n\nfn qconj(q: vec4f) -> vec4f {\n  return vec4f(-q.xyz, q.w);\n}\n\n/** The demo's quat - quat: 2\u00B7vec(a\u00B7b\u207B\u00B9), a small-angle rotation vector. */\nfn qsub(a: vec4f, b: vec4f) -> vec3f {\n  return qmul(a, qconj(b) / dot(b, b)).xyz * 2.0;\n}\n\n/** The demo's quat + float3: integrate a rotation vector, normalize(a + (v, 0)\u00B7a\u00B7\u00BD). */\nfn qadd(a: vec4f, v: vec3f) -> vec4f {\n  return normalize(a + qmul(vec4f(v, 0.0), a) * 0.5);\n}\n\nfn qrotate(q: vec4f, v: vec3f) -> vec3f {\n  let t = cross(q.xyz, v) * 2.0;\n  return v + t * q.w + cross(q.xyz, t);\n}\n\n/** Contact basis rows (n, t1, t2), as maths.h orthonormal(). */\nfn orthonormal(n: vec3f) -> mat3x3f {\n  var t1 = select(vec3f(0.0, -n.z, n.y), vec3f(-n.y, n.x, 0.0), abs(n.x) > abs(n.z));\n  t1 = normalize(t1);\n  // Columns of the WGSL matrix are the basis rows\n  return mat3x3f(n, t1, cross(n, t1));\n}\n\n\nfn lanesFor(bodies: u32) -> u32 {\n  let rule = params.primalLanes;\n  if (bodies >= (1u << (rule & 255u))) { return 1u; }\n  if (bodies >= (1u << ((rule >> 8u) & 255u))) { return 2u; }\n  if (bodies >= (1u << ((rule >> 16u) & 255u))) { return 4u; }\n  return 8u;\n}\n\nfn colorThreads(bodies: u32) -> u32 {\n  return bodies * lanesFor(bodies);\n}\n\n\n// Per-dispatch constants, selected with a dynamic uniform offset\nstruct PassConstants {\n  color: u32,   // colour solved by this primal dispatch\n  alpha: f32,   // stabilization for this iteration\n  pad0: u32,\n  pad1: u32,\n}\n\n@group(0) @binding(0) var<uniform> params: Params;\n@group(0) @binding(1) var<storage, read_write> bodies: array<Body>;\n@group(0) @binding(2) var<storage, read_write> joints: array<Joint>;\n@group(0) @binding(3) var<storage, read> info: array<vec4i>;  // type, bodyA (-1 = world), bodyB\n@group(0) @binding(4) var<storage, read_write> contacts: array<Contact>;\n@group(0) @binding(8) var<storage, read> manifolds: array<Manifold>;\n@group(0) @binding(5) var<storage, read> adj: array<u32>;\n@group(0) @binding(6) var<storage, read> color: array<u32>;\n@group(0) @binding(7) var<storage, read> counters: array<u32>;\n@group(1) @binding(0) var<uniform> pc: PassConstants;\n\nfn manifoldCount() -> u32 {\n  return min(counters[C_MANIFOLDS], params.manifoldCapacity);\n}\n\n/** outer(u, v): row i, column j = u[i] v[j] (WGSL matrices are column-major). */\nfn outer(u: vec3f, v: vec3f) -> mat3x3f {\n  return mat3x3f(u * v.x, u * v.y, u * v.z);\n}\n\n/**\n * One body's Newton system [lin cross\u1D40; cross ang]\u00B7dx = rhs. cross holds row = angular\n * index, column = linear index.\n */\nstruct Acc {\n  lin: mat3x3f,\n  ang: mat3x3f,\n  cross: mat3x3f,\n  rLin: vec3f,\n  rAng: vec3f,\n}\n\n/** Fold one constraint row in: Jacobian (l, a) w.r.t. the body, stiffness k and force f. */\nfn addRow(acc: ptr<function, Acc>, l: vec3f, a: vec3f, k: f32, f: f32) {\n  (*acc).lin += outer(l, l) * k;\n  (*acc).ang += outer(a, a) * k;\n  (*acc).cross += outer(a, l) * k;\n  (*acc).rLin += l * f;\n  (*acc).rAng += a * f;\n}\n\n// --- Joints and springs -----------------------------------------------------------------------\n\nfn anchorA(k: Joint, a: i32) -> vec3f {\n  if (a < 0) { return k.rA.xyz; }\n  return qrotate(bodies[a].rot, k.rA.xyz) + bodies[a].pos.xyz;\n}\n\nfn rotA(a: i32) -> vec4f {\n  if (a < 0) { return vec4f(0.0, 0.0, 0.0, 1.0); }\n  return bodies[a].rot;\n}\n\nfn jointLinC(k: Joint, a: i32, b: i32) -> vec3f {\n  return anchorA(k, a) - (qrotate(bodies[b].rot, k.rB.xyz) + bodies[b].pos.xyz);\n}\n\n/**\n * The angle lock's error: the turn from rotB to rotA\u00B7rest (the demo's qsub, 2\u00B7vec(q\u00B7rotB\u207B\u00B9)), in\n * length units by the torque arm. The identity rest gives the demo's qsub(rotA, rotB). q and -q\n * are one rotation, so the quaternion between them is taken in the hemisphere with w >= 0: the\n * short way round, and always restoring (rest, or the two bodies' quaternions, may carry either sign).\n */\nfn jointAngC(k: Joint, a: i32, b: i32) -> vec3f {\n  let qb = bodies[b].rot;\n  var d = qmul(qmul(rotA(a), k.rest), qconj(qb) / dot(qb, qb));\n  d = select(d, -d, d.w < 0.0);\n  return d.xyz * 2.0 * k.lamAng.w;\n}\n\n/** Stabilized linear / angular constraints (Eq. 18 on hard rows). */\nfn jointLinRows(k: Joint, a: i32, b: i32, alpha: f32) -> vec3f {\n  return jointLinC(k, a, b) - select(vec3f(0.0), k.c0Lin.xyz * alpha, k.penLin.w >= HARD);\n}\n\nfn jointAngRows(k: Joint, a: i32, b: i32, alpha: f32) -> vec3f {\n  return jointAngC(k, a, b) - select(vec3f(0.0), k.c0Ang.xyz * alpha, k.penAng.w >= HARD);\n}\n\nfn addJoint(acc: ptr<function, Acc>, j: u32, alpha: f32, i: u32) {\n  let t = info[j].x;\n  let a = info[j].y;\n  let b = info[j].z;\n  let k = joints[j];\n  let isA = i32(i) == a;\n  let sg = select(-1.0, 1.0, isA);\n\n  if (t == T_SPRING) {\n    let rAW = qrotate(bodies[a].rot, k.rA.xyz);\n    let rBW = qrotate(bodies[b].rot, k.rB.xyz);\n    let d = (rAW + bodies[a].pos.xyz) - (rBW + bodies[b].pos.xyz);\n    let len = length(d);\n    if (len <= 1.0e-6) { return; }\n    let n = d / len;\n    let stiffness = k.penLin.w;\n    let r = select(rBW, rAW, isA);\n    addRow(acc, n * sg, cross(r, n) * sg, stiffness, stiffness * (len - k.rA.w));\n    return;\n  }\n\n  // Ball-socket rows\n  if (dot(k.penLin.xyz, k.penLin.xyz) > 0.0) {\n    let F = k.penLin.xyz * jointLinRows(k, a, b, alpha) + k.lamLin.xyz;\n    // Jacobian: sg\u00B7I (linear) and skew(-rA) for A, skew(rB) for B (angular); rows below\n    var r: vec3f;\n    if (isA) { r = qrotate(bodies[a].rot, k.rA.xyz); }\n    else { r = -qrotate(bodies[b].rot, k.rB.xyz); }\n    // Row m of skew(-r): (-r) \u00D7 e_m... written out: skew(v) rows (0,-vz,vy), (vz,0,-vx), (-vy,vx,0)\n    let v = -r;\n    addRow(acc, vec3f(sg, 0.0, 0.0), vec3f(0.0, -v.z, v.y), k.penLin.x, F.x);\n    addRow(acc, vec3f(0.0, sg, 0.0), vec3f(v.z, 0.0, -v.x), k.penLin.y, F.y);\n    addRow(acc, vec3f(0.0, 0.0, sg), vec3f(-v.y, v.x, 0.0), k.penLin.z, F.z);\n    // Diagonally lumped geometric stiffness (Sec. 3.5): column norms of -(r\u00B7F) I + r F\u1D40\n    let rf = dot(r, F);\n    let g = vec3f(\n      length(r * F.x - vec3f(rf, 0.0, 0.0)),\n      length(r * F.y - vec3f(0.0, rf, 0.0)),\n      length(r * F.z - vec3f(0.0, 0.0, rf)));\n    (*acc).ang += mat3x3f(vec3f(g.x, 0.0, 0.0), vec3f(0.0, g.y, 0.0), vec3f(0.0, 0.0, g.z));\n  }\n\n  // Angle-lock rows: Jacobian sg\u00B7torqueArm\u00B7I on the angular part\n  if (dot(k.penAng.xyz, k.penAng.xyz) > 0.0) {\n    let F = k.penAng.xyz * jointAngRows(k, a, b, alpha) + k.lamAng.xyz;\n    let s = sg * k.lamAng.w;\n    addRow(acc, vec3f(0.0), vec3f(s, 0.0, 0.0), k.penAng.x, F.x);\n    addRow(acc, vec3f(0.0), vec3f(0.0, s, 0.0), k.penAng.y, F.y);\n    addRow(acc, vec3f(0.0), vec3f(0.0, 0.0, s), k.penAng.z, F.z);\n  }\n}\n\n// --- Contacts ---------------------------------------------------------------------------------\n\n/** What a contact needs of one body, loaded once per pair: rotation and displacement since x-. */\nstruct PairBody {\n  index: u32,\n  rot: vec4f,\n  dLin: vec3f,\n  dAng: vec3f,\n}\n\nfn pairBody(i: u32) -> PairBody {\n  let rot = bodies[i].rot;\n  return PairBody(i, rot, bodies[i].pos.xyz - bodies[i].initialPos.xyz, qsub(rot, bodies[i].initialRot));\n}\n\n/**\n * One contact point evaluated at the current poses: C, the force and the lever arms (callers\n * rebuild their own angular Jacobian rows from these; keeping all six rows live cost\n * registers, and the primal ran 60% slower on mixed piles).\n */\nstruct ContactEval {\n  rAW: vec3f,\n  rBW: vec3f,\n  C: vec3f,\n  F: vec3f,        // cone-clamped force\n  frictionScale: f32,\n  bounds: f32,\n}\n\nfn evalContact(k: Contact, basis: mat3x3f, friction: f32, A: PairBody, B: PairBody, alpha: f32) -> ContactEval {\n  // Lever arms at the current rotation, as the reference does for boxes. Sphere contacts use\n  // the step-start rotation instead (the Taylor point x-): a rolling sphere turns ~0.1 rad per\n  // step, and a lever arm rotated with it puts a false separation into the normal row, so the\n  // sphere sank through the ground while gaining energy.\n  var rotA = A.rot;\n  var rotB = B.rot;\n  if ((k.key & ~STICK_BIT) == SPHERE_FEATURE) {\n    rotA = bodies[A.index].initialRot;\n    rotB = bodies[B.index].initialRot;\n  }\n  var e: ContactEval;\n  e.rAW = qrotate(rotA, k.rA);\n  e.rBW = qrotate(rotB, k.rB);\n\n  // Taylor series approximation of C(x) about x- (Sec. 4), one row per basis vector l, with\n  // angular Jacobians rA \u00D7 l and rB \u00D7 (-l)\n  let n = basis[0];\n  let t1 = basis[1];\n  let t2 = basis[2];\n  let keep = 1.0 - alpha;\n  e.C = vec3f(\n    k.c0x * keep + dot(n, A.dLin) - dot(n, B.dLin) + dot(cross(e.rAW, n), A.dAng) + dot(cross(e.rBW, -n), B.dAng),\n    k.c0y * keep + dot(t1, A.dLin) - dot(t1, B.dLin) + dot(cross(e.rAW, t1), A.dAng) + dot(cross(e.rBW, -t1), B.dAng),\n    k.c0z * keep + dot(t2, A.dLin) - dot(t2, B.dLin) + dot(cross(e.rAW, t2), A.dAng) + dot(cross(e.rBW, -t2), B.dAng));\n  var F = k.pen * e.C + k.lam;\n  // Normal pushes only; friction is clamped to the cone\n  F.x = min(F.x, 0.0);\n  e.bounds = abs(F.x) * friction;\n  e.frictionScale = length(F.yz);\n  if (e.frictionScale > e.bounds && e.frictionScale > 0.0) {\n    F = vec3f(F.x, F.yz * (e.bounds / e.frictionScale));\n  }\n  e.F = F;\n  return e;\n}\n\n// --- Warm start -------------------------------------------------------------------------------\n\n@compute @workgroup_size(64)\nfn warmStartJoints(@builtin(global_invocation_id) gid: vec3u) {\n  let j = gid.x;\n  if (j >= params.jointCount || info[j].x != T_JOINT) { return; }\n  let k = joints[j];\n  let a = info[j].y;\n  let b = info[j].z;\n  // Only the fields that change are stored back (not the anchors or the rests: a joint record is\n  // 160 bytes, and a joint-heavy scene runs this over every one)\n  joints[j].c0Lin = vec4f(jointLinC(k, a, b), 0.0);\n  joints[j].c0Ang = vec4f(jointAngC(k, a, b), k.c0Ang.w);\n  let decay = params.alpha * params.gamma;\n  joints[j].lamLin = vec4f(k.lamLin.xyz * decay, k.lamLin.w);\n  joints[j].lamAng = vec4f(k.lamAng.xyz * decay, k.lamAng.w);\n  // Penalties decay, stay within bounds, and never exceed the material stiffness\n  joints[j].penLin = vec4f(min(clamp(k.penLin.xyz * params.gamma, vec3f(PENALTY_MIN), vec3f(PENALTY_MAX)), vec3f(k.penLin.w)), k.penLin.w);\n  joints[j].penAng = vec4f(min(clamp(k.penAng.xyz * params.gamma, vec3f(PENALTY_MIN), vec3f(PENALTY_MAX)), vec3f(k.penAng.w)), k.penAng.w);\n}\n\n/**\n * The wind on a sail (a thin plate): pressure drag \u00BD\u03C1C_d A (n\u00B7u)|n\u00B7u| along its normal n, u\n * the wind relative to the plate, plus skin friction along the plate (6% of that\n * coefficient), which streams a flag out downwind. The wind gusts in time and across space,\n * sideways too (a few travelling sines), so a flag lying along it still catches it and\n * ripples. The pressure is linearised implicitly in the plate's normal speed (divided by\n * 1 + dt \u2202a/\u2202v): a light plate in real air is stiff, and explicit drag would blow up.\n */\nfn windAccel(i: u32, x: vec3f, rot: vec4f, v: vec3f) -> vec3f {\n  let t = f32(params.step) * params.dt;\n  let w = params.wind.xyz;\n  let speed = length(w);\n  // Gusts travel downwind: phases in the wind's frame (a downwind, c across)\n  let d = select(vec3f(1.0, 0.0, 0.0), w / speed, speed > 1e-6);\n  let side = vec3f(-d.y, d.x, 0.0);\n  let a = dot(x, d);\n  let c = dot(x, side);\n  let along = 0.5 * sin(1.1 * t - 0.35 * a + 0.2 * c) + 0.3 * sin(2.9 * t - 0.8 * a + 0.6 * x.z)\n    + 0.2 * sin(6.1 * t - 1.7 * a + 1.3 * x.z + 0.9 * c);\n  let across = 0.6 * sin(1.7 * t - 0.6 * a + 0.4 * x.z) + 0.4 * sin(4.3 * t - 1.3 * a - 0.7 * x.z);\n  let lift = sin(2.3 * t - 0.9 * a + 1.1 * c);\n  let g = params.gust.x;\n  let u = w * (1.0 + g * along) + (side * (0.25 * across) + vec3f(0.0, 0.0, 0.1 * lift)) * (g * speed) - v;\n  let n = qrotate(rot, vec3f(0.0, 0.0, 1.0));\n  let size = bodies[i].size;\n  let k = params.wind.w * size.x * size.y / size.w; // \u00BD\u03C1C_d A / m\n  let un = dot(u, n);\n  let ut = u - n * un;\n  return n * (k * un * abs(un) / (1.0 + 2.0 * params.dt * k * abs(un))) + ut * (0.06 * k * length(ut));\n}\n\n@compute @workgroup_size(64)\nfn warmStartBodies(@builtin(global_invocation_id) gid: vec3u) {\n  let i = gid.x;\n  if (i >= params.bodyCount) { return; }\n  // Field-wise loads and stores: only the poses change (whole-record copies moved 320 bytes)\n  let pos = bodies[i].pos;\n  let rot = bodies[i].rot;\n  let vel = bodies[i].vel;\n  let angVel = bodies[i].angVel.xyz;\n  let dt = params.dt;\n  let g = params.gravity;\n  let dynamic = bodies[i].size.w > 0.0;\n\n  // Inertial target (Eq. 2), with the wind's push on sails as a second external force\n  var inertialPos = pos.xyz + vel.xyz * dt;\n  if (dynamic) { inertialPos += params.up.xyz * (g * (dt * dt)); }\n  if (dynamic && bodies[i].angVel.w == SHAPE_SAIL) { inertialPos += windAccel(i, pos.xyz, rot, vel.xyz) * (dt * dt); }\n  bodies[i].inertialPos = vec4f(inertialPos, bodies[i].inertialPos.w);\n  bodies[i].inertialRot = qadd(rot, angVel * dt);\n\n  // Adaptive warm start (original VBD paper), along up; vel.w holds last step's velocity along up. With\n  // FLAG_START_AT_REST a body slower than REST_SPEED starts from x- instead: at few\n  // iterations a resting stack never corrects the extrapolated guess, which pumps tall walls\n  // until they buckle (docs/FINDINGS.md, Stage 8m)\n  var w = 0.0;\n  if (abs(g) > 0.0) { w = clamp((dot(vel.xyz, params.up.xyz) - vel.w) / dt * sign(g) / abs(g), 0.0, 1.0); }\n  let slow = length(vel.xyz) + length(angVel) * bodies[i].moment.w < REST_SPEED;\n  let atRest = (params.flags & FLAG_START_AT_REST) != 0u && slow;\n\n  bodies[i].initialPos = pos;\n  bodies[i].initialRot = rot;\n  if (dynamic && !atRest) {\n    bodies[i].pos = vec4f(pos.xyz + vel.xyz * dt + params.up.xyz * (g * (w * dt * dt)), pos.w);\n    bodies[i].rot = qadd(rot, angVel * dt);\n  }\n  // A fixed body with a velocity is kinematic: it moves by it over the step, from where it\n  // started (initialPos), so what it touches sees it slide and is pushed and carried along\n  if (!dynamic && (any(vel.xyz != vec3f(0.0)) || any(angVel != vec3f(0.0)))) {\n    bodies[i].pos = vec4f(pos.xyz + vel.xyz * dt, pos.w);\n    bodies[i].rot = qadd(rot, angVel * dt);\n  }\n}\n\n// --- Primal: one colour -----------------------------------------------------------------------\n\n/** A partial system per body slot handed between lanes (primal). */\nvar<workgroup> partials: array<Acc, 32>;\n\n/**\n * One colour's bodies, lanesFor(colour size) threads per body: each lane accumulates every\n * lanes-th adjacency entry, the partial systems are summed in workgroup memory (halving the\n * lanes each round), and lane 0 solves. Bodies of one colour share no constraint, so writing\n * poses in place never races with a neighbour's read (bodies left clashing by the colouring\n * are the only, counted, exception).\n */\n@compute @workgroup_size(64)\nfn primal(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: vec3u) {\n  let start = color[params.colorStartOffset + pc.color];\n  let count = color[params.colorStartOffset + pc.color + 1u] - start;\n  let lanes = lanesFor(count);\n  let slot = lid.x / lanes;\n  let lane = lid.x % lanes;\n  let k = wid.x * (64u / lanes) + slot;\n  let live = k < count;\n  var i = 0u;\n  var acc: Acc;\n  if (live) {\n    i = color[params.colorBodiesOffset + start + k];\n    acc = accumulate(i, lane, lanes);\n  }\n  // Lanes s..2s-1 hand their sums to lanes 0..s-1 (constant bounds keep the barriers uniform)\n  for (var s = 4u; s > 0u; s >>= 1u) {\n    let merging = s < lanes;\n    if (merging && lane >= s && lane < 2u * s) { partials[slot * s + lane - s] = acc; }\n    workgroupBarrier();\n    if (merging && lane < s) {\n      let other = partials[slot * s + lane];\n      acc.lin += other.lin;\n      acc.ang += other.ang;\n      acc.cross += other.cross;\n      acc.rLin += other.rLin;\n      acc.rAng += other.rAng;\n    }\n    workgroupBarrier();\n  }\n  if (live && lane == 0u) { finishBody(i, acc); }\n}\n\n/**\n * Body i's Newton system from adjacency entries lane, lane + lanes, ... (lane 0 also adds the\n * inertia terms); lanes > 1 split a body's constraints across threads.\n */\nfn accumulate(i: u32, lane: u32, lanes: u32) -> Acc {\n  let pos = bodies[i].pos;\n  let rot = bodies[i].rot;\n  let dt2 = params.dt * params.dt;\n  let m = bodies[i].size.w / dt2;\n  let I = bodies[i].moment.xyz / dt2;\n  var acc: Acc;\n  if (lane == 0u) {\n    acc.lin = mat3x3f(vec3f(m, 0.0, 0.0), vec3f(0.0, m, 0.0), vec3f(0.0, 0.0, m));\n    acc.ang = mat3x3f(vec3f(I.x, 0.0, 0.0), vec3f(0.0, I.y, 0.0), vec3f(0.0, 0.0, I.z));\n    acc.rLin = m * (pos.xyz - bodies[i].inertialPos.xyz);\n    acc.rAng = I * qsub(rot, bodies[i].inertialRot);\n  }\n\n  // One flat loop over joints and contact points: each iteration handles one point, loading\n  // its pair when the previous pair runs out. Only what the point needs stays live across\n  // iterations (partner and own rotation/displacement, the normal): live registers limit\n  // how many threads hide memory latency, and mixed piles run in small, latency-bound\n  // per-colour dispatches.\n  var e = adj[i] + lane;\n  let end = adj[i + 1u];\n  var c = 0u;\n  var cEnd = 0u;\n  var A: PairBody;\n  var B: PairBody;\n  var normal = vec3f(0.0);\n  var friction = 0.0;\n  var isA = false;\n  loop {\n    if (c == cEnd) {\n      if (e >= end) { break; }\n      let id = adj[params.adjListOffset + e];\n      e += lanes;\n      if (id < params.jointCount) {\n        addJoint(&acc, id, pc.alpha, i);\n        continue;\n      }\n      let mf = manifolds[id - params.jointCount];\n      isA = i == mf.ids.x;\n      A = pairBody(mf.ids.x);\n      B = pairBody(mf.ids.y);\n      normal = mf.geo.xyz;\n      friction = mf.geo.w;\n      c = mf.ids.z;\n      cEnd = mf.ids.z + pairCount(mf);\n      if (c == cEnd) { continue; }\n    }\n    let k = contacts[c];\n    c++;\n    let basis = orthonormal(normal);\n    let ev = evalContact(k, basis, friction, A, B, pc.alpha);\n    // This body's rows: linear \u00B1l, angular r \u00D7 (\u00B1l)\n    let sg = select(-1.0, 1.0, isA);\n    let r = select(ev.rBW, ev.rAW, isA);\n    let l0 = basis[0] * sg;\n    let l1 = basis[1] * sg;\n    let l2 = basis[2] * sg;\n    addRow(&acc, l0, cross(r, l0), k.pen.x, ev.F.x);\n    addRow(&acc, l1, cross(r, l1), k.pen.y, ev.F.y);\n    addRow(&acc, l2, cross(r, l2), k.pen.z, ev.F.z);\n  }\n  return acc;\n}\n\n/** Solve body i's 6x6 system and apply the update (Eq. 4). */\nfn finishBody(i: u32, acc: Acc) {\n  let pos = bodies[i].pos;\n  let rot = bodies[i].rot;\n  // LDL\u1D40 solve of the 6x6 SPD system (maths.h solve), lower triangle only. Each pivot is floored\n  // at 1e-5 of its diagonal: in f32 a direction held only by a light body's inertia, under\n  // contacts near PENALTY_MAX (a sliver on its edge under a slab), cancels to noise and can come\n  // out negative, and the step then climbs until the body is Inf, then NaN. The floor stiffens\n  // just that direction, and only when it is that ill-conditioned (pivot.gpu.test.ts).\n  let A11 = acc.lin[0][0];\n  let A21 = acc.lin[0][1]; let A22 = acc.lin[1][1];\n  let A31 = acc.lin[0][2]; let A32 = acc.lin[1][2]; let A33 = acc.lin[2][2];\n  let A41 = acc.cross[0][0]; let A42 = acc.cross[1][0]; let A43 = acc.cross[2][0]; let A44 = acc.ang[0][0];\n  let A51 = acc.cross[0][1]; let A52 = acc.cross[1][1]; let A53 = acc.cross[2][1]; let A54 = acc.ang[0][1]; let A55 = acc.ang[1][1];\n  let A61 = acc.cross[0][2]; let A62 = acc.cross[1][2]; let A63 = acc.cross[2][2]; let A64 = acc.ang[0][2]; let A65 = acc.ang[1][2]; let A66 = acc.ang[2][2];\n\n  let D1 = A11;\n  let L21 = A21 / D1;\n  let L31 = A31 / D1;\n  let L41 = A41 / D1;\n  let L51 = A51 / D1;\n  let L61 = A61 / D1;\n  let D2 = max(A22 - L21 * L21 * D1, A22 * PIVOT_FLOOR);\n  let L32 = (A32 - L21 * L31 * D1) / D2;\n  let L42 = (A42 - L21 * L41 * D1) / D2;\n  let L52 = (A52 - L21 * L51 * D1) / D2;\n  let L62 = (A62 - L21 * L61 * D1) / D2;\n  let D3 = max(A33 - (L31 * L31 * D1 + L32 * L32 * D2), A33 * PIVOT_FLOOR);\n  let L43 = (A43 - L31 * L41 * D1 - L32 * L42 * D2) / D3;\n  let L53 = (A53 - L31 * L51 * D1 - L32 * L52 * D2) / D3;\n  let L63 = (A63 - L31 * L61 * D1 - L32 * L62 * D2) / D3;\n  let D4 = max(A44 - (L41 * L41 * D1 + L42 * L42 * D2 + L43 * L43 * D3), A44 * PIVOT_FLOOR);\n  let L54 = (A54 - L41 * L51 * D1 - L42 * L52 * D2 - L43 * L53 * D3) / D4;\n  let L64 = (A64 - L41 * L61 * D1 - L42 * L62 * D2 - L43 * L63 * D3) / D4;\n  let D5 = max(A55 - (L51 * L51 * D1 + L52 * L52 * D2 + L53 * L53 * D3 + L54 * L54 * D4), A55 * PIVOT_FLOOR);\n  let L65 = (A65 - L51 * L61 * D1 - L52 * L62 * D2 - L53 * L63 * D3 - L54 * L64 * D4) / D5;\n  let D6 = max(A66 - (L61 * L61 * D1 + L62 * L62 * D2 + L63 * L63 * D3 + L64 * L64 * D4 + L65 * L65 * D5), A66 * PIVOT_FLOOR);\n\n  let y1 = acc.rLin.x;\n  let y2 = acc.rLin.y - L21 * y1;\n  let y3 = acc.rLin.z - L31 * y1 - L32 * y2;\n  let y4 = acc.rAng.x - L41 * y1 - L42 * y2 - L43 * y3;\n  let y5 = acc.rAng.y - L51 * y1 - L52 * y2 - L53 * y3 - L54 * y4;\n  let y6 = acc.rAng.z - L61 * y1 - L62 * y2 - L63 * y3 - L64 * y4 - L65 * y5;\n\n  let w3 = y6 / D6;\n  let w2 = y5 / D5 - L65 * w3;\n  let w1 = y4 / D4 - L54 * w2 - L64 * w3;\n  let v3 = y3 / D3 - L43 * w1 - L53 * w2 - L63 * w3;\n  let v2 = y2 / D2 - L32 * v3 - L42 * w1 - L52 * w2 - L62 * w3;\n  let v1 = y1 / D1 - L21 * v2 - L31 * v3 - L41 * w1 - L51 * w2 - L61 * w3;\n\n  // dx = -A\u207B\u00B9 rhs (Eq. 4)\n  bodies[i].pos = vec4f(pos.xyz - vec3f(v1, v2, v3), pos.w);\n  bodies[i].rot = qadd(rot, -vec3f(w1, w2, w3));\n}\n\n// --- Dual -------------------------------------------------------------------------------------\n\nfn dualJoint(j: u32) {\n  let t = info[j].x;\n  if (t != T_JOINT) { return; }\n  let a = info[j].y;\n  let b = info[j].z;\n  let k = joints[j];\n  let alpha = pc.alpha;\n  var penLin = k.penLin;\n  var penAng = k.penAng;\n  var lamLin = k.lamLin;\n  var lamAng = k.lamAng;\n  if (dot(penLin.xyz, penLin.xyz) > 0.0) {\n    let C = jointLinRows(k, a, b, alpha);\n    if (penLin.w >= HARD) { lamLin = vec4f(penLin.xyz * C + lamLin.xyz, lamLin.w); }\n    penLin = vec4f(min(penLin.xyz + abs(C) * params.betaLin, vec3f(min(penLin.w, PENALTY_MAX))), penLin.w);\n  }\n  let frac = lamLin.w;\n  let limit = abs(frac);\n  var yielded = false;\n  if (dot(penAng.xyz, penAng.xyz) > 0.0) {\n    let C = jointAngRows(k, a, b, alpha);\n    if (penAng.w >= HARD) {\n      lamAng = vec4f(penAng.xyz * C + lamAng.xyz, lamAng.w);\n      // Yield (a plastic joint, see layout.ts J_YIELD): the force asked of a rigid angle lock\n      // is past yield but short of fracture, so the bend becomes the new rest (the error and\n      // the step's reference error with it: else stabilization, Eq. 18, would pull 99% of the\n      // bend back) and the force carried is cut back to yield. Perfectly plastic: the joint\n      // holds the bend with nothing more than yield, and lets the rest go on bending.\n      let yieldForce = k.c0Ang.w;\n      if (yieldForce < BIG) {\n        let force = length(lamAng.xyz);\n        if (force > yieldForce && !(limit < BIG && force > limit)) {\n          lamAng = vec4f(lamAng.xyz * (yieldForce / force), lamAng.w);\n          yielded = true;\n        }\n      }\n    }\n    // A yielding row is at its bound, as a clamped force: its penalty is not ramped\n    if (!yielded) {\n      penAng = vec4f(min(penAng.xyz + abs(C) * params.betaAng, vec3f(min(penAng.w, PENALTY_MAX))), penAng.w);\n    }\n  }\n  // A joint that gave moves its rest to where the bodies are. If it has then bent past its limit\n  // (layout.ts J_BREAK_BEND: the turn from the rest it started with to this one) it breaks, which\n  // is how a sustained overload tears a plastic joint: yield cuts its force back before the\n  // fracture test below sees more than a sudden jump\n  var rest = k.rest;\n  var bentOff = false;\n  if (yielded) {\n    rest = normalize(qmul(qconj(rotA(a)), bodies[b].rot));\n    if (k.rB.w < BIG) {\n      let d = qmul(qconj(k.restStart), rest);\n      bentOff = 2.0 * atan2(length(d.xyz), abs(d.w)) > k.rB.w;\n    }\n  }\n  // Fracture: the joint stops acting for good (the CPU deletes it). A negative threshold (not in\n  // the paper: GpuSolver3D.appendJoints) breaks on the linear force too, at the same limit\n  let linear = frac < 0.0 && dot(lamLin.xyz, lamLin.xyz) > limit * limit;\n  if ((limit < BIG && (dot(lamAng.xyz, lamAng.xyz) > limit * limit || linear)) || bentOff) {\n    penLin = vec4f(0.0);\n    penAng = vec4f(0.0);\n    lamLin = vec4f(0.0, 0.0, 0.0, frac);\n    lamAng = vec4f(0.0, 0.0, 0.0, lamAng.w);\n  }\n  joints[j].penLin = penLin;\n  joints[j].penAng = penAng;\n  joints[j].lamLin = lamLin;\n  joints[j].lamAng = lamAng;\n  if (yielded) {\n    joints[j].rest = rest;\n    joints[j].c0Ang = vec4f(0.0, 0.0, 0.0, k.c0Ang.w);\n  }\n}\n\nfn dualManifold(m: u32) {\n  let mf = manifolds[m];\n  let A = pairBody(mf.ids.x);\n  let B = pairBody(mf.ids.y);\n  let basis = orthonormal(mf.geo.xyz);\n  let end = mf.ids.z + pairCount(mf);\n  for (var c = mf.ids.z; c < end; c++) {\n    let k = contacts[c];\n    let e = evalContact(k, basis, mf.geo.w, A, B, pc.alpha);\n    // Write back only what changes (lambda, penalty, stick): the dual is bandwidth-bound\n    contacts[c].lam = e.F;\n    // Ramp the penalty where the force is within its bounds (Eq. 16)\n    var pen = k.pen;\n    if (e.F.x < 0.0) { pen.x = min(pen.x + params.betaLin * abs(e.C.x), PENALTY_MAX); }\n    if (e.frictionScale <= e.bounds) {\n      pen.y = min(pen.y + params.betaLin * abs(e.C.y), PENALTY_MAX);\n      pen.z = min(pen.z + params.betaLin * abs(e.C.z), PENALTY_MAX);\n      let stick = length(e.C.yz) < STICK_THRESH;\n      contacts[c].key = (k.key & ~STICK_BIT) | select(0u, STICK_BIT, stick);\n    }\n    contacts[c].pen = pen;\n  }\n}\n\n/** Dual update of every constraint: joints first, then contact pairs (one dispatch). */\n@compute @workgroup_size(64)\nfn dual(@builtin(global_invocation_id) gid: vec3u) {\n  let id = gid.x;\n  if (id < params.jointCount) { dualJoint(id); }\n  else if (id - params.jointCount < manifoldCount()) { dualManifold(id - params.jointCount); }\n}\n\n@compute @workgroup_size(64)\nfn updateVelocities(@builtin(global_invocation_id) gid: vec3u) {\n  let i = gid.x;\n  if (i >= params.bodyCount) { return; }\n  // This step's starting velocity along up, for the next step's adaptive warm start\n  let prevUp = dot(bodies[i].vel.xyz, params.up.xyz);\n  if (bodies[i].size.w > 0.0) {\n    bodies[i].vel = vec4f((bodies[i].pos.xyz - bodies[i].initialPos.xyz) / params.dt, prevUp);\n    bodies[i].angVel = vec4f(qsub(bodies[i].rot, bodies[i].initialRot) / params.dt, bodies[i].angVel.w);\n  } else {\n    bodies[i].vel.w = prevUp;\n  }\n}\n";
//# sourceMappingURL=wgsl-solve.d.ts.map