// =====================================================================
// CC BSpline math in WGSL (Phase 6d stage 1).
//
// Foundation for the GPU patch-eval port: the sharpness-aware bicubic
// basis machinery from cc_topology.js's Phase 1b.13a port. Pure
// scalar / 4x4 matrix math, no scene-dependent inputs -- validates
// directly against the bsplineOracle sweep dumped by cc_dump.
//
// CC_GPU_BSPLINE_WGSL defines the helper functions as a string
// constant (for embedding into other shaders via WGSL include-style
// concatenation). Each function mirrors the JS port in
// cc_topology.js, with these subtleties handled:
//
//   - amolite stores Mat44 row-major (m_data[r*4 + c]); WGSL mat4x4
//     is column-major. We use flat array<f32, 16> in row-major order
//     and provide matVecRowMajorTranspose() to match amolite's
//     `MulMatVec(GetTranspose(M), v)` call exactly.
//   - amolite's f32 reassociation matters for byte-equal vs the
//     bsplineOracle; the JS port uses Math.fround() at each step
//     in cc_topology.js. WGSL is f32-native so the equivalent
//     rounding happens automatically.
//
// Public WGSL symbols:
//   matVecRowMajorTranspose(m, v) -> vec4<f32>
//   getMInfinityRM() -> array<f32, 16>
//   evalBSplineBasis(srcT)        -> vec4<f32>
//   evalBSplineBasisDeriv(srcT)   -> vec4<f32>
//   calcMSharpness(rho)           -> array<f32, 16>
//   calcSharpnessMatricesAndEndpoints(sharpness, out_M_inf_hat,
//                                     out_M_s_hat, out_seg0, out_seg1)
//   calcBasisAndDerivForSharpnessRef(t, sharpness, reverseBasis,
//                                    out_basis, out_basisDeriv)
//
// (out_* are pointer<function> parameters since WGSL functions can't
// return multiple values without a wrapping struct.)
'use strict';

const CC_GPU_BSPLINE_WGSL = /* wgsl */ `
// --- scalar helpers ---
fn saturateF(x: f32) -> f32 { return clamp(x, 0.0, 1.0); }
fn safeInvF(x: f32)  -> f32 { return select(0.0, 1.0 / x, x != 0.0); }
fn calcFracF(t: f32) -> f32 { return t - floor(t); }

// --- matrix as row-major flat array ---
// amolite's Mat44 stores m_data[r*4 + c]. WGSL doesn't naturally
// expose row-major matrices, so we keep a flat array<f32, 16> and
// do mat-vec ourselves. matVecRowMajorTranspose(m, v) computes
// (m^T) * v in amolite's row-major semantics, which is what
// CalcBasisAndDerivFromSharpnessMatrices needs.
fn matVecRowMajorTranspose(m: array<f32, 16>, v: vec4<f32>) -> vec4<f32> {
  // (m^T)[r][c] = m[c][r] = m[c*4 + r]
  // result[r] = sum_c (m^T)[r][c] * v[c] = sum_c m[c*4 + r] * v[c]
  return vec4<f32>(
    m[0]*v.x + m[4]*v.y + m[8] *v.z + m[12]*v.w,
    m[1]*v.x + m[5]*v.y + m[9] *v.z + m[13]*v.w,
    m[2]*v.x + m[6]*v.y + m[10]*v.z + m[14]*v.w,
    m[3]*v.x + m[7]*v.y + m[11]*v.z + m[15]*v.w);
}

// Identity helper.
fn mat44IdentityRM() -> array<f32, 16> {
  return array<f32, 16>(
    1.0, 0.0, 0.0, 0.0,
    0.0, 1.0, 0.0, 0.0,
    0.0, 0.0, 1.0, 0.0,
    0.0, 0.0, 0.0, 1.0);
}

// (1-t)*A + t*B for row-major Mat44.
fn mat44ScaledSumRM(a: array<f32, 16>, sa: f32,
                    b: array<f32, 16>, sb: f32) -> array<f32, 16> {
  var ret : array<f32, 16>;
  for (var i: u32 = 0u; i < 16u; i = i + 1u) {
    ret[i] = a[i] * sa + b[i] * sb;
  }
  return ret;
}

// --- BSplineUtil ---

// GetMInfinity (DcciSubdiv.cpp:8798). Identity except At(3, 1) = -1,
// At(3, 2) = 2, At(3, 3) = 0. Row-major.
fn getMInfinityRM() -> array<f32, 16> {
  var m : array<f32, 16> = mat44IdentityRM();
  m[3 * 4u + 1u] = -1.0;
  m[3 * 4u + 2u] =  2.0;
  m[3 * 4u + 3u] =  0.0;
  return m;
}

// EvalBSplineBasis (DcciSubdiv.cpp:8808).
fn evalBSplineBasis(srcT: f32) -> vec4<f32> {
  let t   = saturateF(srcT);
  let tt  = t * t;
  let ttt = t * t * t;
  return vec4<f32>(
    (-ttt + 3.0 * tt - 3.0 * t + 1.0)       / 6.0,
    ( 3.0 * ttt - 6.0 * tt + 4.0)           / 6.0,
    (-3.0 * ttt + 3.0 * tt + 3.0 * t + 1.0) / 6.0,
    ttt                                     / 6.0);
}

fn evalBSplineBasisDeriv(srcT: f32) -> vec4<f32> {
  let t  = saturateF(srcT);
  let tt = t * t;
  return vec4<f32>(
    (-3.0 * tt + 6.0 * t - 3.0) / 6.0,
    ( 9.0 * tt - 12.0 * t)      / 6.0,
    (-9.0 * tt + 6.0 * t + 3.0) / 6.0,
    ( 3.0 * tt)                 / 6.0);
}

// CalcMSharpness (DcciSubdiv.cpp:8848). Row-major.
fn calcMSharpness(rho: f32) -> array<f32, 16> {
  let invOneMinusRho = safeInvF(1.0 - rho);
  let r2 = rho * rho;
  let r1 = rho;
  let scale = (1.0 - rho) / (6.0 * rho);

  var m : array<f32, 16>;
  m[ 0] = scale * (6.0 * r1 * invOneMinusRho);
  m[ 1] = scale * (6.0 * r2 - 5.0 * r1 + 1.0);
  m[ 2] = scale * (-12.0 * r2 + 10.0 * r1 - 2.0);
  m[ 3] = scale * (6.0 * r2 - 5.0 * r1 + 1.0);

  m[ 4] = 0.0;
  m[ 5] = scale * ((2.0 * r2 + 3.0 * r1 + 1.0) * invOneMinusRho);
  m[ 6] = scale * (4.0 * r1 - 2.0);
  m[ 7] = scale * (-2.0 * r1 + 1.0);

  m[ 8] = 0.0;
  m[ 9] = scale * (r1 + 1.0);
  m[10] = scale * ((2.0 * r2 + 6.0 * r1 - 2.0) * invOneMinusRho);
  m[11] = scale * (r1 + 1.0);

  m[12] = 0.0;
  m[13] = scale * (-2.0 * r1 + 1.0);
  m[14] = scale * (4.0 * r1 - 2.0);
  m[15] = scale * ((2.0 * r2 + 3.0 * r1 + 1.0) * invOneMinusRho);

  return m;
}

// CalcSharpnessMatricesAndEndpoints (DcciSubdiv.cpp:8885).
// Output via pointer<function> params since WGSL functions can't
// return multiple values without a struct (which we avoid for inline
// helpers).
struct SharpnessMatrices {
  m_infinity_hat : array<f32, 16>,
  m_s_hat        : array<f32, 16>,
  seg0           : f32,
  seg1           : f32,
};

fn calcSharpnessMatricesAndEndpoints(sharpness: f32) -> SharpnessMatrices {
  let m_infinity = getMInfinityRM();
  let s_floor = floor(sharpness);
  let s_ceil  = ceil(sharpness);

  let rho_floor = pow(2.0, s_floor);
  let rho_ceil  = pow(2.0, s_ceil);

  var m_floor = calcMSharpness(rho_floor);
  let m_ceil  = calcMSharpness(rho_ceil);

  if (sharpness < 1.0) {
    m_floor = mat44IdentityRM();
  }

  let s_frac = calcFracF(sharpness);

  let m_infinity_hat = mat44ScaledSumRM(m_floor,  1.0 - s_frac, m_infinity, s_frac);
  let m_s_hat        = mat44ScaledSumRM(m_floor,  1.0 - s_frac, m_ceil,     s_frac);

  let seg0 = 1.0 - pow(2.0, -s_floor);
  let seg1 = 1.0 - pow(2.0, -s_ceil);

  var out : SharpnessMatrices;
  out.m_infinity_hat = m_infinity_hat;
  out.m_s_hat        = m_s_hat;
  out.seg0           = seg0;
  out.seg1           = seg1;
  return out;
}

// CalcBasisAndDerivFromSharpnessMatrices (DcciSubdiv.cpp:8916).
struct BasisAndDeriv {
  basis      : vec4<f32>,
  basisDeriv : vec4<f32>,
};

fn calcBasisAndDerivFromSharpnessMatrices(
    sm: SharpnessMatrices, t: f32, reverseSharpnessBasis: bool) -> BasisAndDeriv {
  let m_infinity = getMInfinityRM();
  let sharpnessT = select(t, 1.0 - t, reverseSharpnessBasis);

  let basis      = evalBSplineBasis(sharpnessT);
  let basisDeriv = evalBSplineBasisDeriv(sharpnessT);

  var retBasis      : vec4<f32>;
  var retBasisDeriv : vec4<f32>;
  if (sharpnessT < sm.seg0) {
    retBasis      = matVecRowMajorTranspose(m_infinity, basis);
    retBasisDeriv = matVecRowMajorTranspose(m_infinity, basisDeriv);
  } else if (sharpnessT < sm.seg1) {
    retBasis      = matVecRowMajorTranspose(sm.m_infinity_hat, basis);
    retBasisDeriv = matVecRowMajorTranspose(sm.m_infinity_hat, basisDeriv);
  } else {
    retBasis      = matVecRowMajorTranspose(sm.m_s_hat, basis);
    retBasisDeriv = matVecRowMajorTranspose(sm.m_s_hat, basisDeriv);
  }

  if (reverseSharpnessBasis) {
    let tb = retBasis;
    let td = retBasisDeriv;
    retBasis      = vec4<f32>(tb.w, tb.z, tb.y, tb.x);
    retBasisDeriv = vec4<f32>(td.w, td.z, td.y, td.x);
  }

  var out : BasisAndDeriv;
  out.basis      = retBasis;
  out.basisDeriv = retBasisDeriv;
  return out;
}

fn calcBasisAndDerivForSharpnessRef(t: f32, sharpness: f32, reverseSharpnessBasis: bool) -> BasisAndDeriv {
  let sm = calcSharpnessMatricesAndEndpoints(sharpness);
  return calcBasisAndDerivFromSharpnessMatrices(sm, t, reverseSharpnessBasis);
}

// =====================================================================
// EvalBicubicBSplineSharpness (DcciSubdiv.cpp:8998). Vec3 variant.
// Given a 4x4 grid of point indices into ptData + (tx, ty) + sharpness,
// returns (foundP, foundN) where foundN = normalize(cross(dx_tan, dy_tan)).
//
// Sharpness applies on the tx axis (with reverse=true so the sharp
// edge sits along tx=0). ty uses standard B-spline basis.
//
// The shader-side caller must provide:
//   - ptOffsets : array<u32, 16> -- the 16 indices into ptData
//   - ptData    : array<f32>     -- flat xyz Vec3 storage
fn evalBicubicBSplineSharpness_setup(sharpness: f32, tx: f32, ty: f32) -> array<vec4<f32>, 4> {
  let sx       = calcBasisAndDerivForSharpnessRef(tx, sharpness, true);
  let basisY   = evalBSplineBasis(ty);
  let basisDy  = evalBSplineBasisDeriv(ty);
  // Pack as 4 vec4s: basisX, basisDx, basisY, basisDy.
  return array<vec4<f32>, 4>(sx.basis, sx.basisDeriv, basisY, basisDy);
}

// 2-axis sharpness variant. Used by the multi-sharp Fix-3 fallback
// for 2-adjacent (perpendicular creases): sharpnessX > 0 with one
// reverseX, sharpnessY > 0 with one reverseY. The 2-opposite case
// is handled via two separate calls + position averaging at the
// caller (see cc_gpu_patch.js's evalRegular).
fn evalBicubicBSplineSharpness_setup_full(
    sharpnessX: f32, reverseX: bool,
    sharpnessY: f32, reverseY: bool,
    tx: f32, ty: f32) -> array<vec4<f32>, 4> {
  let sx = calcBasisAndDerivForSharpnessRef(tx, sharpnessX, reverseX);
  var basisY  : vec4<f32>;
  var basisDy : vec4<f32>;
  if (sharpnessY > 0.0) {
    let sy = calcBasisAndDerivForSharpnessRef(ty, sharpnessY, reverseY);
    basisY  = sy.basis;
    basisDy = sy.basisDeriv;
  } else {
    basisY  = evalBSplineBasis(ty);
    basisDy = evalBSplineBasisDeriv(ty);
  }
  return array<vec4<f32>, 4>(sx.basis, sx.basisDeriv, basisY, basisDy);
}

struct BicubicEvalOut {
  P : vec3<f32>,
  N : vec3<f32>,
  dx : vec3<f32>,
  dy : vec3<f32>,
};
`;

// =====================================================================
// CC_GPU_BICUBIC_TEST_WGSL  (Phase 6d stage 2)
//
// Test shader: for each test case, read 16 control points + (tx, ty,
// sharpness) and compute (P, N, dx_tan, dy_tan). The bicubic eval
// lives inline since WGSL function args don't accept variable-size
// arrays without storage-buffer indirection.
//
// Layout:
//   inputs[i].sharpness, tx, ty, _pad
//   inputs[i].ptOffsets[16] (u32 indices into ptData)
//   ptData : array<f32> flat xyz Vec3
//   outputs[i].P, N, dx, dy (vec3 each, with padding)
// =====================================================================
const CC_GPU_BICUBIC_TEST_WGSL = /* wgsl */ `
${CC_GPU_BSPLINE_WGSL}

struct BicubicCaseIn {
  sharpness  : f32,
  tx         : f32,
  ty         : f32,
  _pad0      : f32,
  ptOffsets  : array<u32, 16>,
};

struct BicubicCaseOut {
  P  : vec4<f32>,
  N  : vec4<f32>,
  dx : vec4<f32>,
  dy : vec4<f32>,
};

@group(0) @binding(0) var<storage, read>       inputs  : array<BicubicCaseIn>;
@group(0) @binding(1) var<storage, read>       ptData  : array<f32>;
@group(0) @binding(2) var<storage, read_write> outputs : array<BicubicCaseOut>;
@group(0) @binding(3) var<uniform>             NUM     : vec4<u32>;

@compute @workgroup_size(64)
fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
  let i = gid.x;
  if (i >= NUM.x) { return; }
  let tc = inputs[i];

  let basis = evalBicubicBSplineSharpness_setup(tc.sharpness, tc.tx, tc.ty);
  let basisX  = basis[0];
  let basisDx = basis[1];
  let basisY  = basis[2];
  let basisDy = basis[3];

  var sumP  = vec3<f32>(0.0);
  var sumDx = vec3<f32>(0.0);
  var sumDy = vec3<f32>(0.0);

  for (var r: u32 = 0u; r < 4u; r = r + 1u) {
    let wY  = basisY[r];
    let wDy = basisDy[r];
    for (var c: u32 = 0u; c < 4u; c = c + 1u) {
      let wX  = basisX[c];
      let wDx = basisDx[c];
      let off = tc.ptOffsets[r * 4u + c];
      let P   = vec3<f32>(ptData[off*3u + 0u], ptData[off*3u + 1u], ptData[off*3u + 2u]);
      let wxy  = wX  * wY;
      let wdxy = wDx * wY;
      let wxdy = wX  * wDy;
      sumP  = sumP  + P * wxy;
      sumDx = sumDx + P * wdxy;
      sumDy = sumDy + P * wxdy;
    }
  }

  let nRaw = cross(sumDx, sumDy);
  let nLen = length(nRaw);
  let N    = select(vec3<f32>(0.0), nRaw / nLen, nLen > 0.0);

  outputs[i].P  = vec4<f32>(sumP,  0.0);
  outputs[i].N  = vec4<f32>(N,     0.0);
  outputs[i].dx = vec4<f32>(sumDx, 0.0);
  outputs[i].dy = vec4<f32>(sumDy, 0.0);
}
`;

async function gpuComputeCcBicubic(device, cases, ptData16Flat) {
  // cases: [{sharpness, tx, ty, ptOffsets: [16 ints]}]
  // ptData16Flat: Float32Array of (numPts * 3) flat-xyz; offsets index into this.
  const numCases = cases.length;

  // Per-case: 4 floats (sharpness, tx, ty, pad) + 16 u32 = 20 words = 80 bytes.
  const caseBytes = 80;
  const inBytes = numCases * caseBytes;
  const inBuf = new ArrayBuffer(inBytes);
  const inF32 = new Float32Array(inBuf);
  const inU32 = new Uint32Array(inBuf);
  for (let i = 0; i < numCases; i++) {
    const base = i * 20;
    inF32[base + 0] = cases[i].sharpness;
    inF32[base + 1] = cases[i].tx;
    inF32[base + 2] = cases[i].ty;
    inF32[base + 3] = 0;
    for (let k = 0; k < 16; k++) inU32[base + 4 + k] = cases[i].ptOffsets[k];
  }

  const inputsBuf = device.createBuffer({
    size: Math.max(inBytes, 16),
    usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
    label: 'cc_gpu_bicubic/inputs',
  });
  device.queue.writeBuffer(inputsBuf, 0, inBuf);

  const ptBuf = device.createBuffer({
    size: Math.max(ptData16Flat.byteLength, 16),
    usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
    label: 'cc_gpu_bicubic/ptData',
  });
  device.queue.writeBuffer(ptBuf, 0, ptData16Flat.buffer, ptData16Flat.byteOffset, ptData16Flat.byteLength);

  // Outputs: 4 vec4 per case = 64 bytes.
  const outBytes = numCases * 64;
  const outputsBuf = device.createBuffer({
    size: Math.max(outBytes, 16),
    usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC,
    label: 'cc_gpu_bicubic/outputs',
  });

  const numBuf = device.createBuffer({
    size: 16, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
  });
  device.queue.writeBuffer(numBuf, 0, new Uint32Array([numCases, 0, 0, 0]).buffer);

  const layout = device.createBindGroupLayout({
    entries: [
      { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'read-only-storage' } },
      { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'read-only-storage' } },
      { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage'           } },
      { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'uniform'           } },
    ],
  });
  const pipeline = device.createComputePipeline({
    layout:  device.createPipelineLayout({ bindGroupLayouts: [layout] }),
    compute: { module: device.createShaderModule({ code: CC_GPU_BICUBIC_TEST_WGSL }), entryPoint: 'main' },
  });
  const bindGroup = device.createBindGroup({
    layout,
    entries: [
      { binding: 0, resource: { buffer: inputsBuf  } },
      { binding: 1, resource: { buffer: ptBuf      } },
      { binding: 2, resource: { buffer: outputsBuf } },
      { binding: 3, resource: { buffer: numBuf     } },
    ],
  });

  const enc  = device.createCommandEncoder();
  const pass = enc.beginComputePass();
  pass.setPipeline(pipeline);
  pass.setBindGroup(0, bindGroup);
  pass.dispatchWorkgroups(Math.ceil(numCases / 64));
  pass.end();
  const rb = device.createBuffer({
    size: outBytes, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST,
  });
  enc.copyBufferToBuffer(outputsBuf, 0, rb, 0, outBytes);
  device.queue.submit([enc.finish()]);
  await rb.mapAsync(GPUMapMode.READ);
  const flat = new Float32Array(rb.getMappedRange().slice(0, outBytes));
  rb.unmap();

  inputsBuf.destroy(); ptBuf.destroy(); outputsBuf.destroy();
  numBuf.destroy(); rb.destroy();

  const out = new Array(numCases);
  for (let i = 0; i < numCases; i++) {
    const base = i * 16; // 4 vec4 per case, 4 floats per vec4
    out[i] = {
      P:  [flat[base+0],  flat[base+1],  flat[base+2]],
      N:  [flat[base+4],  flat[base+5],  flat[base+6]],
      dx: [flat[base+8],  flat[base+9],  flat[base+10]],
      dy: [flat[base+12], flat[base+13], flat[base+14]],
    };
  }
  return out;
}

// =====================================================================
// CC_GPU_BSPLINE_ORACLE_WGSL
//
// Test shader: takes a per-case input {t, sharpness, reverse} array
// and writes back (basis, basisDeriv) for each. Used to validate the
// WGSL BSpline math against the cc_dump-emitted bsplineOracle sweep.
// =====================================================================
const CC_GPU_BSPLINE_ORACLE_WGSL = /* wgsl */ `
${CC_GPU_BSPLINE_WGSL}

struct TestCaseIn {
  t         : f32,
  sharpness : f32,
  reverse   : u32,
  _pad      : u32,
};

struct TestCaseOut {
  basis      : vec4<f32>,
  basisDeriv : vec4<f32>,
};

@group(0) @binding(0) var<storage, read>          inputs  : array<TestCaseIn>;
@group(0) @binding(1) var<storage, read_write>    outputs : array<TestCaseOut>;
@group(0) @binding(2) var<uniform>                NUM     : vec4<u32>;  // x = numCases

@compute @workgroup_size(64)
fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
  let i = gid.x;
  if (i >= NUM.x) { return; }
  let tc = inputs[i];
  let r  = calcBasisAndDerivForSharpnessRef(tc.t, tc.sharpness, tc.reverse != 0u);
  outputs[i].basis      = r.basis;
  outputs[i].basisDeriv = r.basisDeriv;
}
`;

// =====================================================================
// gpuComputeCcBSplineOracle
//
// Run the bsplineOracle sweep on GPU. Returns array of {basis,
// basisDeriv} matching the structure of refJson.subdiv.bsplineOracle.
// basis sweep (NOT the sharpness sweep -- that one tests the matrices
// directly; would need a different shader to dump out M_infinity_hat
// + M_s_hat + seg0/seg1).
// =====================================================================
async function gpuComputeCcBSplineOracle(device, basisCases) {
  // basisCases: array of {t, sharpness, reverse} (numbers 0/1 for reverse).
  const numCases = basisCases.length;

  // Pack inputs as a buffer of {t: f32, sharpness: f32, reverse: u32, pad: u32}.
  const inBytes = numCases * 16;
  const inBuf = new ArrayBuffer(inBytes);
  const inF32 = new Float32Array(inBuf);
  const inU32 = new Uint32Array(inBuf);
  for (let i = 0; i < numCases; i++) {
    inF32[i * 4 + 0] = basisCases[i].t;
    inF32[i * 4 + 1] = basisCases[i].sharpness;
    inU32[i * 4 + 2] = basisCases[i].reverse ? 1 : 0;
    inU32[i * 4 + 3] = 0;
  }

  const inputsBuf = device.createBuffer({
    size: Math.max(inBytes, 16),
    usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
    label: 'cc_gpu_bspline/inputs',
  });
  device.queue.writeBuffer(inputsBuf, 0, inBuf);

  // Outputs: 8 floats per case (basis vec4 + basisDeriv vec4).
  const outBytes = numCases * 32;
  const outputsBuf = device.createBuffer({
    size: Math.max(outBytes, 16),
    usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC,
    label: 'cc_gpu_bspline/outputs',
  });

  const numBuf = device.createBuffer({
    size: 16, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST, label: 'cc_gpu_bspline/numCases',
  });
  device.queue.writeBuffer(numBuf, 0, new Uint32Array([numCases, 0, 0, 0]).buffer);

  const layout = device.createBindGroupLayout({
    entries: [
      { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'read-only-storage' } },
      { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage'           } },
      { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'uniform'           } },
    ],
  });
  const pipeline = device.createComputePipeline({
    layout:  device.createPipelineLayout({ bindGroupLayouts: [layout] }),
    compute: { module: device.createShaderModule({ code: CC_GPU_BSPLINE_ORACLE_WGSL }), entryPoint: 'main' },
  });
  const bindGroup = device.createBindGroup({
    layout,
    entries: [
      { binding: 0, resource: { buffer: inputsBuf  } },
      { binding: 1, resource: { buffer: outputsBuf } },
      { binding: 2, resource: { buffer: numBuf     } },
    ],
  });

  const enc  = device.createCommandEncoder();
  const pass = enc.beginComputePass();
  pass.setPipeline(pipeline);
  pass.setBindGroup(0, bindGroup);
  pass.dispatchWorkgroups(Math.ceil(numCases / 64));
  pass.end();

  const rb = device.createBuffer({
    size: outBytes, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST,
  });
  enc.copyBufferToBuffer(outputsBuf, 0, rb, 0, outBytes);
  device.queue.submit([enc.finish()]);
  await rb.mapAsync(GPUMapMode.READ);
  const flat = new Float32Array(rb.getMappedRange().slice(0, outBytes));
  rb.unmap();

  inputsBuf.destroy(); outputsBuf.destroy(); numBuf.destroy(); rb.destroy();

  // Unpack into [{basis: [b0,b1,b2,b3], basisDeriv: [...] }]
  const out = new Array(numCases);
  for (let i = 0; i < numCases; i++) {
    out[i] = {
      basis:      [flat[i*8+0], flat[i*8+1], flat[i*8+2], flat[i*8+3]],
      basisDeriv: [flat[i*8+4], flat[i*8+5], flat[i*8+6], flat[i*8+7]],
    };
  }
  return out;
}

if (typeof module !== 'undefined' && module.exports) {
  module.exports = {
    CC_GPU_BSPLINE_WGSL,
    CC_GPU_BSPLINE_ORACLE_WGSL,
    CC_GPU_BICUBIC_TEST_WGSL,
    gpuComputeCcBSplineOracle,
    gpuComputeCcBicubic,
  };
}
