Evaluate native tree norms as real Hermitian forms

rust/loop-reduction/README.md

@@ -180,6 +180,16 @@ Nc=3, TR=1/2. The color reducer's formal Nc polynomials do not establish a
 general-Nc model capability. Ordered color connections act on the retained
 tensors before contraction, so higher correlations do not multiply Gram matrices.
 
+Identity norms verify exact Gram Hermiticity and lower the quadratic form to
+real/imaginary coefficient components before numerical evaluation. Thus a norm
+is real by construction, including complex color bases and helicity phases;
+no computed imaginary residue is clipped. The same coefficient executor is
+used, and complete directed complex amplitudes remain available unchanged.
+Real cancellation error and possible negative numerical results are not hidden
+or promoted to an accuracy guarantee. Legacy saved norm programs consuming one
+complex column per coefficient remain accepted; new Hermitian norm programs
+consume adjacent real/imaginary columns instead.
+
 `evaluator_with_precision_bits` constructs the existing typed binary64,
 double-double or arbitrary-precision runtime. Requested bits are initial working
 precision; values retain their propagated precision after cancellation. Quad

rust/loop-reduction/src/lib.rs

@@ -12,11 +12,12 @@
 
 pub mod preparation;
 pub mod runtime;
+pub mod scalar_loop;
 mod tree_color;
 mod tree_components;
-mod tree_runtime;
+mod tree_norm;
 pub mod tree_preparation;
-pub mod scalar_loop;
+mod tree_runtime;
 pub mod verification;
 pub mod virtual_runtime;
 

rust/loop-reduction/src/tree_norm.rs

@@ -0,0 +1,148 @@
+// SPDX-License-Identifier: 0BSD
+
+//! Exact Cartesian lowering of Hermitian color norms. Directed correlations
+//! retain their complex coefficient programs and do not use this lowering.
+
+use rusticol_core::{RusticolError, RusticolResult, color::ExactColorCoefficient};
+use symbolica::{atom::Atom, symbol};
+
+pub(crate) fn hermitian_norm(
+    gram: &[Vec<ExactColorCoefficient>],
+) -> RusticolResult<(Atom, Vec<Atom>)> {
+    if gram.iter().any(|row| row.len() != gram.len()) {
+        return Err(RusticolError::evaluation(
+            "tree norm Gram matrix must be square",
+        ));
+    }
+    for (i, row) in gram.iter().enumerate() {
+        for (j, value) in row.iter().enumerate() {
+            if *value != gram[j][i].conj() {
+                return Err(RusticolError::evaluation(
+                    "tree norm Gram matrix must be exactly Hermitian",
+                ));
+            }
+        }
+    }
+    let inputs: Vec<_> = (0..gram.len())
+        .flat_map(|index| {
+            ["real", "imag"].map(move |part| {
+                Atom::var(symbol!(&format!("rusticol_native_tree::amplitude_{part}_{index}"); Real))
+            })
+        })
+        .collect();
+    let mut norm = Atom::Zero;
+    for i in 0..gram.len() {
+        let (x, y) = (&inputs[2 * i], &inputs[2 * i + 1]);
+        norm += Atom::num(gram[i][i].re.clone()) * (x * x + y * y);
+        for j in i + 1..gram.len() {
+            let (u, v) = (&inputs[2 * j], &inputs[2 * j + 1]);
+            // 2 Re[(x-i y) (R+i I) (u+i v)]. The imaginary
+            // Gram entry contributes +2 I (y u-x v), not its negative.
+            norm += Atom::num(2)
+                * (Atom::num(gram[i][j].re.clone()) * (x * u + y * v)
+                    + Atom::num(gram[i][j].im.clone()) * (y * u - x * v));
+        }
+    }
+    Ok((norm, inputs))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::{
+        Complex,
+        runtime::{Arithmetic, Batch, Data},
+        tree_runtime::TreeEvaluator,
+    };
+    use rusticol_core::preparation::encode_exact_evaluator_state;
+    use symbolica::atom::AtomCore;
+    use symbolica::domains::{
+        float::{Float, FloatLike, Real, SingleFloat},
+        rational::Rational,
+    };
+
+    fn gram() -> Vec<Vec<ExactColorCoefficient>> {
+        vec![
+            vec![
+                Complex::new(2.into(), 0.into()),
+                Complex::new(1.into(), 2.into()),
+            ],
+            vec![
+                Complex::new(1.into(), (-2).into()),
+                Complex::new(3.into(), 0.into()),
+            ],
+        ]
+    }
+
+    #[test]
+    fn complex_hermitian_gram_uses_real_instructions_without_losing_phases() {
+        let matrix = gram();
+        let vector = [
+            Complex::new(Rational::from((1, 3)), Rational::from((2, 7))),
+            Complex::new(Rational::from((-5, 11)), Rational::from((3, 13))),
+        ];
+        // Independent exact directed double sum, without Cartesian lowering.
+        let mut expected = Complex::new(Rational::from(0), Rational::from(0));
+        for i in 0..2 {
+            for j in 0..2 {
+                expected += vector[i].conj() * &matrix[i][j] * &vector[j];
+            }
+        }
+        assert_eq!(expected.im, Rational::from(0));
+        let coefficients = vector
+            .iter()
+            .map(|z| Atom::num(z.re.clone()) + Atom::i() * Atom::num(z.im.clone()))
+            .collect::<Vec<_>>();
+        let exact = Atom::evaluator_multiple(&coefficients, &[] as &[Atom])
+            .cores(1)
+            .build()
+            .unwrap();
+        let (norm, inputs) = hermitian_norm(&matrix).unwrap();
+        let norm = Atom::evaluator_multiple(&[norm], &inputs)
+            .cores(1)
+            .build()
+            .unwrap();
+        let exact = encode_exact_evaluator_state(&exact).unwrap();
+        let norm = encode_exact_evaluator_state(&norm).unwrap();
+        let mut evaluator =
+            TreeEvaluator::new(&exact, &norm, vec![], 0, Arithmetic::Double, 53).unwrap();
+        let p = Batch::typed::<f64>(vec![], vec![9], 0, 53).unwrap();
+        let result = evaluator.evaluate_batch(&p).unwrap();
+        let Data::Double(result) = result.data() else {
+            panic!("carrier changed")
+        };
+        assert_eq!(result[0].im, 0.);
+        assert!((result[0].re - expected.re.to_f64()).abs() < 1e-14);
+        let directed = evaluator.evaluate_coefficients(&p, &[]).unwrap();
+        let Data::Double(directed) = directed.data() else {
+            panic!("carrier changed")
+        };
+        assert!(directed.iter().all(|z| z.im != 0.));
+        for bits in [113, 173, 512] {
+            let mut evaluator =
+                TreeEvaluator::new(&exact, &norm, vec![], 0, Arithmetic::Arbitrary, bits).unwrap();
+            let p = Batch::typed::<Float>(vec![], vec![9], 0, bits).unwrap();
+            let output = evaluator.evaluate_batch(&p).unwrap();
+            let Data::Arbitrary(values) = output.data() else {
+                panic!("carrier changed")
+            };
+            let target = Float::with_val(bits, 0).from_rational(&expected.re);
+            assert_eq!(values[0].im, Float::with_val(bits, 0));
+            assert!(
+                (values[0].re.clone() - target).norm()
+                    < Float::with_val(bits, 2).pow(u64::from(bits - 12)).inv()
+            );
+        }
+    }
+
+    #[test]
+    fn nonhermitian_gram_is_rejected_before_real_lowering() {
+        let mut invalid = gram();
+        invalid[1][0].im = 2.into();
+        assert!(hermitian_norm(&invalid).is_err());
+        let mut invalid = gram();
+        invalid[0][0].im = 1.into();
+        assert!(hermitian_norm(&invalid).is_err());
+        assert!(hermitian_norm(&[vec![]]).is_err());
+    }
+}

rust/loop-reduction/src/tree_preparation.rs

@@ -41,6 +41,7 @@ use crate::{
     runtime::{Arithmetic, arithmetic_mode},
     tree_color::{factor_tree_color, index_label},
     tree_components::{TreeComponentVector, contract_tree_state_components},
+    tree_norm::hermitian_norm,
 };
 
 /// Immutable instructions for coherent tree amplitudes and their color norm.
@@ -566,17 +567,16 @@ impl TreeProgram {
             .cores(1)
             .build()
             .map_err(error)?;
-        let amplitude_inputs: Vec<_> = (0..color_tensors.len())
-            .map(|index| Atom::var(symbol!(&format!("rusticol_native_tree::amplitude_{index}"))))
-            .collect();
-        let mut norm = Atom::Zero;
-        for (bra, left) in color_tensors.iter().enumerate() {
-            for (ket, right) in color_tensors.iter().enumerate() {
-                let gram = left.overlap(right)?;
-                let coefficient = Atom::num(gram.re) + Atom::i() * Atom::num(gram.im);
-                norm += amplitude_inputs[bra].conj() * coefficient * &amplitude_inputs[ket];
-            }
-        }
+        let gram = color_tensors
+            .iter()
+            .map(|left| {
+                color_tensors
+                    .iter()
+                    .map(|right| left.overlap(right))
+                    .collect::<RusticolResult<Vec<_>>>()
+            })
+            .collect::<RusticolResult<Vec<_>>>()?;
+        let (norm, amplitude_inputs) = hermitian_norm(&gram)?;
         // Preserve the scalar entry point's direct-norm state by composing the
         // already contracted coefficients. No diagrams or tensors are rebuilt.
         let scalar_norm_state = if field_class == TreeFieldClass::SingletScalars {
@@ -585,7 +585,12 @@ impl TreeProgram {
                     .iter()
                     .position(|input| input.as_view() == term)
                 {
-                    **out = coefficients[index].clone();
+                    let coefficient = &coefficients[index / 2];
+                    **out = if index % 2 == 0 {
+                        (coefficient + coefficient.conj()) * Atom::num((1, 2))
+                    } else {
+                        (coefficient - coefficient.conj()) * Atom::i() * Atom::num((-1, 2))
+                    };
                 }
             });
             let evaluator = Atom::evaluator_multiple(&[expression], &inputs)
@@ -814,7 +819,10 @@ mod tests {
                     (value.re - expected).abs() < 1e-9 * expected,
                     "{incoming:?} -> {outgoing:?}: got {value}, expected {expected}"
                 );
-                assert!(value.im.abs() < 1e-9 * expected);
+                assert_eq!(
+                    value.im, 0.,
+                    "Hermitian norm is constructed as a real quadratic form"
+                );
             }
             check_color_charge_conservation(&program);
             check_gluon_ward_identities(&program, &input, expected);

rust/loop-reduction/src/tree_runtime.rs

@@ -41,7 +41,11 @@ pub enum TreeParameterDomain {
 /// all momentum columns, four components per external state, then the runtime
 /// parameter columns in their declared order. Its outputs are the color-basis
 /// coefficients; the norm program consumes those outputs, including their
-/// complex phases, and returns one scalar column. Evaluation methods without
+/// complex phases, and returns one scalar column. Its input layout may contain
+/// one complex column per coefficient (legacy) or two real columns per
+/// coefficient, ordered real then imaginary (exact Hermitian norm lowering).
+/// The latter uses real arithmetic expressions before any numerical summation;
+/// no imaginary part of a computed contraction is discarded. Methods without
 /// a parameter batch require an empty runtime parameter domain.
 pub struct TreeEvaluator {
     amplitude: CoefficientEvaluator,
@@ -50,6 +54,7 @@ pub struct TreeEvaluator {
     parameter_domains: Vec<TreeParameterDomain>,
     momentum_columns: usize,
     coefficient_columns: usize,
+    cartesian_norm: bool,
     arithmetic: Arithmetic,
     bits: u32,
 }
@@ -153,6 +158,39 @@ fn add_norms(left: &Batch, right: &Batch) -> RusticolResult<Batch> {
     Ok(sum)
 }
 
+fn cartesian_values<T: BatchCarrier>(
+    values: &[Complex<T>],
+    metadata: &Batch,
+) -> RusticolResult<Batch> {
+    let values = values
+        .iter()
+        .flat_map(|value| {
+            [
+                Complex::new(value.re.clone(), value.re.zero()),
+                Complex::new(value.im.clone(), value.im.zero()),
+            ]
+        })
+        .collect();
+    Batch::typed(
+        values,
+        metadata.ids().to_vec(),
+        metadata
+            .columns()
+            .checked_mul(2)
+            .ok_or_else(|| error("Cartesian norm dimensions overflow"))?,
+        metadata.precision_bits(),
+    )
+    .map_err(error)
+}
+
+fn cartesian_coefficients(coefficients: &Batch) -> RusticolResult<Batch> {
+    match coefficients.data() {
+        Data::Double(values) => cartesian_values(values, coefficients),
+        Data::DoubleDouble(values) => cartesian_values(values, coefficients),
+        Data::Arbitrary(values) => cartesian_values(values, coefficients),
+    }
+}
+
 impl TreeEvaluator {
     pub fn new(
         coefficient_state: &[u8],
@@ -238,9 +276,12 @@ impl TreeEvaluator {
         let coefficient = exact_state(coefficient_state).map_err(error)?;
         let norm = exact_state(norm_state).map_err(error)?;
         let coefficient_columns = coefficient.get_output_len();
+        let cartesian_norm = coefficient_columns
+            .checked_mul(2)
+            .is_some_and(|columns| columns == norm.get_input_len());
         if coefficient.get_input_len() != inputs
             || coefficient_columns == 0
-            || norm.get_input_len() != coefficient_columns
+            || (!cartesian_norm && norm.get_input_len() != coefficient_columns)
             || norm.get_output_len() != 1
         {
             return Err(error(
@@ -260,6 +301,7 @@ impl TreeEvaluator {
             parameter_domains,
             momentum_columns,
             coefficient_columns,
+            cartesian_norm,
             arithmetic,
             bits,
         })
@@ -496,7 +538,12 @@ impl TreeEvaluator {
                 .collect::<Vec<_>>();
             let coefficients =
                 self.evaluate_coefficients_inner(momenta, parameters, &helicities)?;
-            let norm = self.norm.evaluate_batch(&coefficients).map_err(error)?;
+            let norm_inputs = if self.cartesian_norm {
+                cartesian_coefficients(&coefficients)?
+            } else {
+                coefficients
+            };
+            let norm = self.norm.evaluate_batch(&norm_inputs).map_err(error)?;
             finite(&norm)?;
             accumulated = Some(match accumulated {
                 None => norm,