diff --git a/python/chiral_cross_enrich.wgsl b/python/chiral_cross_enrich.wgsl new file mode 100644 index 00000000..60947460 --- /dev/null +++ b/python/chiral_cross_enrich.wgsl @@ -0,0 +1,190 @@ +/** + * 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 labels: array; // Sidon labels +@group(0) @binding(1) var moduli: array; // CRT moduli +@group(0) @binding(2) var params: EnrichParams; +@group(0) @binding(3) var sidon_results: array; // 1=Sidon, 0=collision +@group(0) @binding(4) var rossby_drifts: array; // Rossby drift per config +@group(0) @binding(5) var helical_residues: array; // Winding mod 28 +@group(0) @binding(6) var couch_stable: array; // 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 { + switch chi_type { + case 0u: { return vec4(i32(Q16_ONE), 0, 0, 0); } // 1 + case 1u: { return vec4(0, 0, 0, i32(Q16_ONE)); } // k + case 2u: { return vec4(0, i32(Q16_ONE), 0, 0); } // i + case 3u: { return vec4(0, 0, i32(Q16_ONE), 0); } // j + default: { return vec4(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, q2: vec4) -> vec4 { + 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(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: scarred count → self-loop → threshold check. +fn compute_couch(config_id: u32) -> u32 { + 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; + } + } + // self_loop = 57942 * scarred_count / 8 (Q16_16) + let self_loop = (57942u * scarred_count) / N_STRANDS; + // COUCH stable if self_loop < threshold + if (self_loop < COUCH_THRESHOLD) { + return 1u; // stable + } + return 0u; // unstable +} + +/// 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) { + 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, 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); +}