/**
 * chiral_cross_enrich.wgsl — Cross-Enriched Chiral Sidon Filter on GPU
 *
 * Cross-enrichment: each strand can have MULTIPLE chiral types
 * contributing simultaneously (not just one label per strand).
 *
 * With 4 ChiralLabel types × 8 strands = 4^8 = 65,536 configurations.
 * Each config is a 4-bit vector per strand (which types are active).
 *
 * This is why GPU is needed: 65K configs × pairwise checks = millions
 * of operations. The existing dna_braid.wgsl workgroup (256) handles
 * 2^8=256 binary configs; this shader handles 4^8=65536 enriched configs.
 *
 * ChiralLabel types (from BraidStateN.lean):
 *   0 = achiral_stable      → quaternion 1,  Rossby weight 0
 *   1 = chiral_scarred      → quaternion k,  Rossby weight +32768
 *   2 = left_handed_mass    → quaternion i,  Rossby weight +65536
 *   3 = right_handed_vector → quaternion j,  Rossby weight -65536
 *
 * Cross-enrichment: strand s has a 4-bit mask (a,s,l,r) where each bit
 * indicates whether that chiral type is active. The strand's effective
 * quaternion is the SUM of active basis elements, normalized.
 *
 * Rossby drift = sum across all strands of active weights.
 * Helical residue = ⌊step × 25042⌋ mod 28 (golden angle winding).
 */

const WORKGROUP_SIZE: u32 = 256u;
const N_STRANDS: u32 = 8u;
const N_TYPES: u32 = 4u;  // achiral, scarred, left, right
const N_CONFIGS: u32 = 65536u;  // 4^8
const Q16_ONE: u32 = 65536u;
const Q16_HALF: u32 = 32768u;
const GOLDEN_ANGLE: u32 = 25042u;
const EXOTIC_CLASSES: u32 = 28u;
const COUCH_THRESHOLD: u32 = 49152u;  // 0.75 × 65536

// Rossby weights (Q16_16 raw)
const W_ACHIRAL: i32 = 0;
const W_SCARRED: i32 = 32768;   // +0.5
const W_LEFT: i32 = 65536;      // +1
const W_RIGHT: i32 = -65536;    // -1

@group(0) @binding(0) var<storage, read> labels: array<u32, 8>;      // Sidon labels
@group(0) @binding(1) var<storage, read> moduli: array<u32, 9>;      // CRT moduli
@group(0) @binding(2) var<uniform> params: EnrichParams;
@group(0) @binding(3) var<storage, read_write> sidon_results: array<u32>;    // 1=Sidon, 0=collision
@group(0) @binding(4) var<storage, read_write> rossby_drifts: array<i32>;   // Rossby drift per config
@group(0) @binding(5) var<storage, read_write> helical_residues: array<u32>; // Winding mod 28
@group(0) @binding(6) var<storage, read_write> couch_stable: array<u32>;     // COUCH gate result

struct EnrichParams {
    S: u32,               // reflection point
    n_configs: u32,       // 65536
    step: u32,            // braid step (for helical residue)
    _pad: u32,
};

/// Extract chiral type for strand s from config_id.
/// Config_id is a base-4 number: digit s = (config_id >> (2*s)) & 3
fn get_chiral_type(config_id: u32, strand: u32) -> u32 {
    return (config_id >> (strand * 2u)) & 3u;
}

/// Get Rossby weight for a chiral type (returns i32 for signed arithmetic).
fn rossby_weight(chi_type: u32) -> i32 {
    switch chi_type {
        case 0u: { return W_ACHIRAL; }   // achiral_stable
        case 1u: { return W_SCARRED; }   // chiral_scarred
        case 2u: { return W_LEFT; }      // left_handed
        case 3u: { return W_RIGHT; }     // right_handed
        default: { return 0; }
    }
}

/// Get quaternion basis (a,b,c,d) for chiral type.
/// achiral=1=(1,0,0,0), scarred=k=(0,0,0,1), left=i=(0,1,0,0), right=j=(0,0,1,0)
fn quat_basis(chi_type: u32) -> vec4<i32> {
    switch chi_type {
        case 0u: { return vec4<i32>(i32(Q16_ONE), 0, 0, 0); }     // 1
        case 1u: { return vec4<i32>(0, 0, 0, i32(Q16_ONE)); }     // k
        case 2u: { return vec4<i32>(0, i32(Q16_ONE), 0, 0); }     // i
        case 3u: { return vec4<i32>(0, 0, i32(Q16_ONE), 0); }     // j
        default: { return vec4<i32>(0, 0, 0, 0); }
    }
}

/// Hamilton product of two Q16_16 quaternions.
/// q1 = (a1,b1,c1,d1), q2 = (a2,b2,c2,d2)
/// Result scaled by Q16_ONE (divide at end).
fn quat_mul(q1: vec4<i32>, q2: vec4<i32>) -> vec4<i32> {
    let a1 = q1.x; let b1 = q1.y; let c1 = q1.z; let d1 = q1.w;
    let a2 = q2.x; let b2 = q2.y; let c2 = q2.z; let d2 = q2.w;
    let a = (a1*a2 - b1*b2 - c1*c2 - d1*d2) / i32(Q16_ONE);
    let b = (a1*b2 + b1*a2 + c1*d2 - d1*c2) / i32(Q16_ONE);
    let c = (a1*c2 - b1*d2 + c1*a2 + d1*b2) / i32(Q16_ONE);
    let d = (a1*d2 + b1*c2 - c1*b2 + d1*a2) / i32(Q16_ONE);
    return vec4<i32>(a, b, c, d);
}

/// Compute Rossby drift for a config (sum of weights across strands).
fn compute_rossby(config_id: u32) -> i32 {
    var total: i32 = 0;
    for (var s: u32 = 0u; s < N_STRANDS; s = s + 1u) {
        let chi = get_chiral_type(config_id, s);
        total = total + rossby_weight(chi);
    }
    return total;
}

/// Compute helical residue (winding number mod 28).
fn compute_helical(step: u32) -> u32 {
    return ((step * GOLDEN_ANGLE) / Q16_ONE) % EXOTIC_CLASSES;
}

/// COUCH stability: two-stage check.
/// Stage A: Rossby/Kelvin regime (drift=0 → Kelvin → boundary-trapped → FAIL)
/// Stage B: scarred contention (self-loop < threshold)
fn compute_couch(config_id: u32) -> u32 {
    // Stage A: Rossby drift must be non-zero (Kelvin regime fails)
    let drift = compute_rossby(config_id);
    if (drift == 0) {
        // Kelvin regime: boundary-trapped, no mixing, zero dissipation
        // (proven: rossby_energy_dissipation_rate requires isActive)
        return 0u;  // COUCH FAILS — stuck like AVX-512 (self_loop=0.885)
    }

    // Stage B: scarred contention check
    var scarred_count: u32 = 0u;
    for (var s: u32 = 0u; s < N_STRANDS; s = s + 1u) {
        if (get_chiral_type(config_id, s) == 1u) {  // chiral_scarred
            scarred_count = scarred_count + 1u;
        }
    }
    let self_loop = (57942u * scarred_count) / N_STRANDS;
    if (self_loop < COUCH_THRESHOLD) {
        return 1u;  // COUCH PASSES — Rossby active + low contention
    }
    return 0u;  // COUCH FAILS — too much scarred contention
}

/// Main kernel: each thread processes one cross-enriched configuration.
/// 65K configs → 256 threads/workgroup × 256 workgroups = 65536 threads.
@compute @workgroup_size(WORKGROUP_SIZE)
fn cross_enriched_filter(@builtin(global_invocation_id) gid: vec3<u32>) {
    let config_id = gid.x;
    if (config_id >= params.n_configs) {
        return;
    }

    // 1. Compute Rossby drift
    let drift = compute_rossby(config_id);
    rossby_drifts[config_id] = drift;

    // 2. Compute helical residue
    let residue = compute_helical(params.step + config_id);
    helical_residues[config_id] = residue;

    // 3. COUCH stability check
    let couch = compute_couch(config_id);
    couch_stable[config_id] = couch;

    // Skip Sidon check if COUCH fails (cheap filter first)
    if (couch == 0u) {
        sidon_results[config_id] = 0u;
        return;
    }

    // 4. Quaternion Sidon check (expensive, only for COUCH-passing configs)
    // Get quaternion for each strand
    var quats: array<vec4<i32>, 8>;
    for (var s: u32 = 0u; s < N_STRANDS; s = s + 1u) {
        let chi = get_chiral_type(config_id, s);
        quats[s] = quat_basis(chi);
    }

    // Check all pairwise quaternion products for collisions
    var collisions: u32 = 0u;
    for (var i: u32 = 0u; i < N_STRANDS; i = i + 1u) {
        for (var j: u32 = i; j < N_STRANDS; j = j + 1u) {
            let prod_ij = quat_mul(quats[i], quats[j]);

            // Check against all other pairs
            for (var k: u32 = 0u; k < N_STRANDS; k = k + 1u) {
                for (var l: u32 = k; l < N_STRANDS; l = l + 1u) {
                    if (i * 8u + j < k * 8u + l) {
                        let prod_kl = quat_mul(quats[k], quats[l]);
                        // Collision: products are equal
                        if (prod_ij == prod_kl) {
                            collisions = collisions + 1u;
                        }
                    }
                }
            }
        }
    }

    sidon_results[config_id] = select(1u, 0u, collisions > 0u);
}