Improve resident trimmed ICP moment precision

CHANGELOG.md

@@ -21,6 +21,14 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`).
 
 ### Changed
 
+- Improve resident trimmed ICP moment precision with a default DoubleSingle
+  route: integer product residuals and preserved f32 rounding carry two words
+  through coordinate normalization, centroid division and covariance products.
+  Explicit `GpuIcpTrimPrecision::CompensatedF32` retains the earlier arithmetic.
+  Fraction one keeps canonical untrimmed results; matching, all-gated fitness,
+  CPU Horn solving and fixed 176-byte readbacks remain. Large gated registration
+  regressions now enforce the unchanged CPU pose-agreement tolerance.
+
 - Add explicit resident GPU ICP trimming under `register-icp-gpu`. Stable
   distance/source-index selection retains floor(trim_fraction*gated_count)
   pairs for moments while all gated pairs still define fitness/stopping.

crates/spatialrust-gpu/src/kernels/icp_reduction.rs

@@ -165,12 +165,19 @@ impl GpuIcpInput {
         transform: Isometry3<f32>,
         max_distance: f32,
     ) -> SpatialResult<GpuIcpMoments> {
-        self.moments_impl(runtime, transform, max_distance, None)
+        self.moments_impl(
+            runtime,
+            transform,
+            max_distance,
+            None,
+            trim::GpuIcpTrimPrecision::DoubleSingle,
+        )
     }
 
     /// Submits stable GPU selection of the lowest-distance floor(fraction*N)
     /// gated pairs, breaking distance ties by original source index. Moments
-    /// use retained pairs; fitness uses ALL gated pairs. Binary64 fraction
+    /// use retained pairs with DoubleSingle normalized arithmetic; fitness
+    /// uses ALL gated pairs. Binary64 fraction
     /// multiplication/floor is emulated with integer limbs on the device.
     /// No match/point payload readback or host wait occurs. A fraction of one
     /// delegates to the canonical untrimmed kernels; no CPU fallback occurs.
@@ -180,15 +187,42 @@ impl GpuIcpInput {
         transform: Isometry3<f32>,
         max_distance: f32,
         fraction: f64,
+    ) -> SpatialResult<GpuIcpTrimmedMoments> {
+        self.trimmed_moments_resident_with_precision(
+            runtime,
+            transform,
+            max_distance,
+            fraction,
+            trim::GpuIcpTrimPrecision::DoubleSingle,
+        )
+    }
+
+    /// Explicit choice of retained-moment arithmetic. DoubleSingle carries
+    /// two f32 words through normalized offsets, division and products; f32
+    /// transforms/search stay unchanged. CompensatedF32 retains the earlier
+    /// coordinate/product precision. Fraction one uses canonical ordinary
+    /// moments for either choice. Transfers/ownership/status rules are identical.
+    pub fn trimmed_moments_resident_with_precision(
+        &self,
+        runtime: &WgpuRuntime,
+        transform: Isometry3<f32>,
+        max_distance: f32,
+        fraction: f64,
+        precision: trim::GpuIcpTrimPrecision,
     ) -> SpatialResult<GpuIcpTrimmedMoments> {
         if !fraction.is_finite() || fraction <= 0.0 || fraction > 1.0 {
             return Err(SpatialError::InvalidArgument(
                 "trim_fraction must be finite and in (0, 1]".to_owned(),
             ));
         }
         let trimmed = fraction < 1.0;
-        let inner =
-            self.moments_impl(runtime, transform, max_distance, trimmed.then_some(fraction))?;
+        let inner = self.moments_impl(
+            runtime,
+            transform,
+            max_distance,
+            trimmed.then_some(fraction),
+            precision,
+        )?;
         Ok(GpuIcpTrimmedMoments { inner, trimmed })
     }
 
@@ -198,6 +232,7 @@ impl GpuIcpInput {
         transform: Isometry3<f32>,
         max_distance: f32,
         fraction: Option<f64>,
+        precision: trim::GpuIcpTrimPrecision,
     ) -> SpatialResult<GpuIcpMoments> {
         let q = transform.rotation();
         let t = transform.translation();
@@ -266,10 +301,10 @@ impl GpuIcpInput {
         let partials = scratch.acquire("icp-partials", u64::from(groups) * 144);
         let selection = fraction.map(|_| scratch.acquire("icp-trim-selection", 32));
         let cache = if fraction.is_some() {
-            &runtime
+            runtime
                 .icp_trim_pipelines
                 .get_or_init(|| trim::TrimPipelines::new(runtime.device()))
-                .reduction
+                .reduction(runtime.device(), precision)
         } else {
             runtime
                 .icp_reduction_pipelines

crates/spatialrust-gpu/src/kernels/icp_trim.rs

@@ -3,6 +3,19 @@ use super::*;
 
 pub(super) const SHADER: &str = include_str!("../shaders/icp_trim.wgsl");
 
+/// Arithmetic used for retained ICP moments; transform/search remain f32.
+#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
+pub enum GpuIcpTrimPrecision {
+    /// Canonical f32 normalized coordinates/products with compensated sums.
+    /// Lower arithmetic cost, but coordinate/product rounding can bias moments.
+    CompensatedF32,
+    /// Two f32 words through offsets, normalization, mean division and products.
+    /// Carries quotient/product residuals to reduce moment bias; hardware f32
+    /// transforms/search and their gate/tie precision remain unchanged.
+    #[default]
+    DoubleSingle,
+}
+
 #[repr(C)]
 #[derive(Clone, Copy, Pod, Zeroable)]
 pub(super) struct TrimParams {
@@ -48,7 +61,8 @@ impl TrimParams {
 }
 
 pub(crate) struct TrimPipelines {
-    pub(super) reduction: IcpReductionPipelines,
+    reduction: std::sync::OnceLock<IcpReductionPipelines>,
+    precise: std::sync::OnceLock<IcpReductionPipelines>,
     layout: wgpu::BindGroupLayout,
     histogram: wgpu::ComputePipeline,
     select: wgpu::ComputePipeline,
@@ -96,13 +110,29 @@ impl TrimPipelines {
             })
         });
         Self {
-            reduction: IcpReductionPipelines::with_shader(device, &trimmed_shader(), true),
+            reduction: std::sync::OnceLock::new(),
+            precise: std::sync::OnceLock::new(),
             layout,
             histogram,
             select,
         }
     }
 
+    pub(super) fn reduction(
+        &self,
+        device: &wgpu::Device,
+        precision: GpuIcpTrimPrecision,
+    ) -> &IcpReductionPipelines {
+        match precision {
+            GpuIcpTrimPrecision::CompensatedF32 => self.reduction.get_or_init(|| {
+                IcpReductionPipelines::with_shader(device, &trimmed_shader(), true)
+            }),
+            GpuIcpTrimPrecision::DoubleSingle => self.precise.get_or_init(|| {
+                IcpReductionPipelines::with_shader(device, &precise_shader(), true)
+            }),
+        }
+    }
+
     pub(super) fn encode(
         &self,
         runtime: &WgpuRuntime,
@@ -192,6 +222,64 @@ fn trimmed_shader() -> String {
     shader
 }
 
+fn precise_shader() -> String {
+    let mut shader = trimmed_shader();
+    shader.push_str(include_str!("../shaders/icp_trim_moments.wgsl"));
+    replace_once(&mut shader, "var value=0.0;", "var value=vec2<f32>(0.0);");
+    replace_once(
+        &mut shader,
+        "value=((p-anchor)/params.scale)[term];",
+        "value=ds_offset(p[term],anchor[term],params.scale);",
+    );
+    replace_once(
+        &mut shader,
+        "value=((q-params.target_anchor.xyz)/params.scale)[term-3u];",
+        "value=ds_offset(q[term-3u],params.target_anchor[term-3u],params.scale);",
+    );
+    replace_once(&mut shader, "value=dot3(d,d);", "value=vec2<f32>(dot3(d,d),0.0);");
+    replace_once(&mut shader, "if !finite(value) {", "if !finite(value.x) || !finite(value.y) {");
+    replace_once(
+        &mut shader,
+        "total=add(total,value);",
+        "if term==15u { total=add(total,value.x); } else { total=precise_merge(total,value); }",
+    );
+    replace_once(&mut shader, "result.values[term]=merge(vec2<f32>(total.x/f32(count),0.0),vec2<f32>(total.y/f32(count),0.0));", "result.values[term]=ds_divide(total,f32(count));");
+    replace_once(
+        &mut shader,
+        "let mean_s=result.values[row].x+result.values[row].y;",
+        "let mean_s=result.values[row];",
+    );
+    replace_once(
+        &mut shader,
+        "let mean_t=result.values[col+3u].x+result.values[col+3u].y;",
+        "let mean_t=result.values[col+3u];",
+    );
+    replace_once(
+        &mut shader,
+        "let s=(queries[i].xyz-result.source_anchor.xyz)/params.scale;",
+        "let s=ds_offset(queries[i][row],result.source_anchor[row],params.scale);",
+    );
+    replace_once(&mut shader, "let t=(target_points[matches[i].index].xyz-result.target_anchor.xyz)/params.scale;", "let t=ds_offset(target_points[matches[i].index][col],result.target_anchor[col],params.scale);");
+    replace_once(
+        &mut shader,
+        "total=add(total,(s[row]-mean_s)*(t[col]-mean_t));",
+        "total=precise_merge(total,ds_product(precise_merge(s,-mean_s),precise_merge(t,-mean_t)));",
+    );
+    replace_once(&mut shader, "if lane<stride { sums[lane]=merge(sums[lane],sums[lane+stride]); counts[lane]+=counts[lane+stride]; }", "if lane<stride { if term==15u { sums[lane]=merge(sums[lane],sums[lane+stride]); } else { sums[lane]=precise_merge(sums[lane],sums[lane+stride]); } counts[lane]+=counts[lane+stride]; }");
+    let group_merge =
+        "for (var g=0u;g<params.groups;g+=1u) { total=merge(total,partials[g].values[term]); }";
+    assert_eq!(shader.matches(group_merge).count(), 2);
+    shader=shader.replacen(group_merge,"for (var g=0u;g<params.groups;g+=1u) { if term==15u { total=merge(total,partials[g].values[term]); } else { total=precise_merge(total,partials[g].values[term]); } }",1);
+    replace_once(
+        &mut shader,
+        "if lane<stride { sums[lane]=merge(sums[lane],sums[lane+stride]); }",
+        "if lane<stride { sums[lane]=precise_merge(sums[lane],sums[lane+stride]); }",
+    );
+    // The remaining group merge belongs to covariance finalization.
+    replace_once(&mut shader, "for (var g=0u;g<params.groups;g+=1u) { total=merge(total,partials[g].values[term]); }", "for (var g=0u;g<params.groups;g+=1u) { total=precise_merge(total,partials[g].values[term]); }");
+    shader
+}
+
 #[cfg(all(test, not(target_arch = "wasm32")))]
 #[path = "icp_trim_tests.rs"]
 mod tests;

crates/spatialrust-gpu/src/kernels/icp_trim_math_tests.rs

@@ -0,0 +1,100 @@
+use super::*;
+
+#[test]
+fn double_single_arithmetic_matches_independent_binary64_values() {
+    let runtime = WgpuRuntime::new_headless().unwrap();
+    let mut inputs: Vec<[f32; 4]> = vec![
+        [1.0, 0.0, 3.0, 0.0],
+        [1.0, 2.0f32.powi(-24), 1.0, -2.0f32.powi(-24)],
+        [16777216.0, 0.0, 1.0, 0.0],
+        [8.2388, 1e-7, 209642.0, 0.0],
+        [-2548.44, 0.00007, 209642.0, 0.0],
+        [-64.0, 0.0, 128.7, 0.0],
+        [-0.78, 1e-8, 131.25, 1e-7],
+    ];
+    let mut seed = 20261013u32;
+    for i in 0..512 {
+        seed = seed.wrapping_mul(1664525).wrapping_add(1013904223);
+        let a = f32::from_bits(
+            ((100 + (i % 50)) as u32) << 23 | (seed & 0x7fffff) | ((i as u32 & 1) << 31),
+        );
+        seed = seed.wrapping_mul(1664525).wrapping_add(1013904223);
+        let b = f32::from_bits(
+            ((100 + ((i * 7) % 50)) as u32) << 23
+                | (seed & 0x7fffff)
+                | (((i as u32 >> 1) & 1) << 31),
+        );
+        inputs.push([a, a * 2.0f32.powi(-27), b, -b * 2.0f32.powi(-27)]);
+    }
+    let source = include_str!("../shaders/icp_reduction.wgsl");
+    let start = source.find("fn add(").unwrap();
+    let end = source.find("fn valid_match(").unwrap();
+    let shader=format!("struct Params {{ unused:u32,count:u32,pad0:u32,pad1:u32 }}\n@group(0) @binding(0) var<uniform> params:Params;\n@group(0) @binding(1) var<storage,read> inputs:array<vec4<f32>>;\n@group(0) @binding(2) var<storage,read_write> output:array<vec2<f32>>;\n{}\n{}\n@compute @workgroup_size(64) fn test_math(@builtin(global_invocation_id) id:vec3<u32>) {{ if id.x>=params.count {{ return; }} let v=inputs[id.x]; output[4u*id.x]=precise_add(v.xy,v.z); output[4u*id.x+1u]=two_product(v.x,v.z); output[4u*id.x+2u]=ds_divide(v.xy,v.z); output[4u*id.x+3u]=ds_product(v.xy,v.zw); }}", &source[start..end],include_str!("../shaders/icp_trim_moments.wgsl"));
+    let module = runtime.device().create_shader_module(wgpu::ShaderModuleDescriptor {
+        label: None,
+        source: wgpu::ShaderSource::Wgsl(shader.into()),
+    });
+    let pipeline = runtime.device().create_compute_pipeline(&wgpu::ComputePipelineDescriptor {
+        label: None,
+        layout: None,
+        module: &module,
+        entry_point: Some("test_math"),
+        compilation_options: Default::default(),
+        cache: None,
+    });
+    let uniform = runtime.device().create_buffer_init(&wgpu::util::BufferInitDescriptor {
+        label: None,
+        contents: bytemuck::cast_slice(&[0u32, inputs.len() as u32, 0, 0]),
+        usage: wgpu::BufferUsages::UNIFORM,
+    });
+    let input = runtime.device().create_buffer_init(&wgpu::util::BufferInitDescriptor {
+        label: None,
+        contents: bytemuck::cast_slice(&inputs),
+        usage: wgpu::BufferUsages::STORAGE,
+    });
+    let output = runtime.device().create_buffer(&wgpu::BufferDescriptor {
+        label: None,
+        size: inputs.len() as u64 * 32,
+        usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC,
+        mapped_at_creation: false,
+    });
+    let entries: Vec<_> = [&uniform, &input, &output]
+        .iter()
+        .enumerate()
+        .map(|(binding, b)| wgpu::BindGroupEntry {
+            binding: binding as u32,
+            resource: b.as_entire_binding(),
+        })
+        .collect();
+    let group = runtime.device().create_bind_group(&wgpu::BindGroupDescriptor {
+        label: None,
+        layout: &pipeline.get_bind_group_layout(0),
+        entries: &entries,
+    });
+    let mut encoder = runtime.device().create_command_encoder(&Default::default());
+    dispatch(&mut encoder, &pipeline, &group, (inputs.len() as u32).div_ceil(64), 1);
+    runtime.queue().submit(Some(encoder.finish()));
+    let words = read(&runtime, &output, output.size());
+    let mut maximum_relative_error = 0.0_f64;
+    for (i, v) in inputs.iter().enumerate() {
+        let a = f64::from(v[0]) + f64::from(v[1]);
+        let b = f64::from(v[2]) + f64::from(v[3]);
+        let expected =
+            [a + f64::from(v[2]), f64::from(v[0]) * f64::from(v[2]), a / f64::from(v[2]), a * b];
+        for (op, expected) in expected.into_iter().enumerate() {
+            let hi = f32::from_bits(words[8 * i + 2 * op]);
+            let lo = f32::from_bits(words[8 * i + 2 * op + 1]);
+            let actual = f64::from(hi) + f64::from(lo);
+            maximum_relative_error =
+                maximum_relative_error.max((actual - expected).abs() / expected.abs().max(1e-12));
+            assert!(
+                (actual - expected).abs() <= 1e-11 * expected.abs().max(1e-12),
+                "input={i} operation={op}"
+            );
+        }
+    }
+    eprintln!(
+        "GPU_ICP_PRECISE_MATH comparisons={} maximum_relative_error={maximum_relative_error}",
+        inputs.len() * 4
+    );
+}