diff --git a/rust/loop-reduction/README.md b/rust/loop-reduction/README.md index 7c03912be..d6e0e4047 100644 --- a/rust/loop-reduction/README.md +++ b/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 diff --git a/rust/loop-reduction/src/lib.rs b/rust/loop-reduction/src/lib.rs index e017935f2..a8e5d263f 100644 --- a/rust/loop-reduction/src/lib.rs +++ b/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; diff --git a/rust/loop-reduction/src/tree_norm.rs b/rust/loop-reduction/src/tree_norm.rs new file mode 100644 index 000000000..85351e1dc --- /dev/null +++ b/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], +) -> RusticolResult<(Atom, Vec)> { + 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![ + 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::>(); + 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::(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::(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()); + } +} diff --git a/rust/loop-reduction/src/tree_preparation.rs b/rust/loop-reduction/src/tree_preparation.rs index 2afe9a0cc..aa91809a7 100644 --- a/rust/loop-reduction/src/tree_preparation.rs +++ b/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 * &litude_inputs[ket]; - } - } + let gram = color_tensors + .iter() + .map(|left| { + color_tensors + .iter() + .map(|right| left.overlap(right)) + .collect::>>() + }) + .collect::>>()?; + 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); diff --git a/rust/loop-reduction/src/tree_runtime.rs b/rust/loop-reduction/src/tree_runtime.rs index c430eeb46..2d5bdec06 100644 --- a/rust/loop-reduction/src/tree_runtime.rs +++ b/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, 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 { Ok(sum) } +fn cartesian_values( + values: &[Complex], + metadata: &Batch, +) -> RusticolResult { + 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 { + 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::>(); 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,