diff --git a/rust/loop-reduction/src/lib.rs b/rust/loop-reduction/src/lib.rs
index f119e9b56..959f9e1c9 100644
--- a/rust/loop-reduction/src/lib.rs
+++ b/rust/loop-reduction/src/lib.rs
@@ -13,6 +13,7 @@
 pub mod preparation;
 pub mod runtime;
 pub mod scalar_loop;
+pub mod tree_correlations;
 mod tree_color;
 mod tree_components;
 mod tree_norm;
diff --git a/rust/loop-reduction/src/tree_correlations.rs b/rust/loop-reduction/src/tree_correlations.rs
new file mode 100644
index 000000000..b8658c713
--- /dev/null
+++ b/rust/loop-reduction/src/tree_correlations.rs
@@ -0,0 +1,117 @@
+// SPDX-License-Identifier: 0BSD
+
+//! Exact color insertions on the retained coherent tree basis.
+//!
+//! Matrices are directed overlaps of complete operator-transformed tensors,
+//! never products of Gram matrices. Numerical amplitudes use the prepared
+//! model's SU(3), `TR=1/2` domain; formal Nc matrices are diagnostic algebra.
+//! Results include no initial averages, flux or final-state symmetry factors.
+//! This module supplies color bilinears, not an open-index spin tensor API.
+
+use std::sync::Arc;
+
+use rusticol_core::{RusticolError, RusticolResult, color::ColorNcTerms};
+
+use crate::{runtime::Arithmetic, tree_preparation::TreeProgram};
+
+mod matrices;
+mod runtime;
+#[cfg(test)]
+mod tests;
+
+pub use runtime::TreeColorEvaluator;
+
+/// All labels refer to original, typed all-outgoing color legs in the tree.
+/// Crossing is already encoded by the prepared tree; no extra sign is added.
+#[derive(Clone, Debug, PartialEq, Eq)]
+#[non_exhaustive]
+pub enum TreeColorInsertion {
+    Identity,
+    /// `T_left.T_right`; coincident labels give the quadratic Casimir.
+    Pair {
+        left: i64,
+        right: i64,
+    },
+    /// Mathematical left-to-right order, applied rightmost first on the ket.
+    /// Empty products are rejected; request `Identity` explicitly instead.
+    OrderedProduct(Vec<[i64; 2]>),
+    /// `AB+BA`, with no factor one half and no repeated-factor deduplication.
+    Anticommutator {
+        first: [i64; 2],
+        second: [i64; 2],
+    },
+    /// `f_abc T_first^a T_second^b T_third^c` for three distinct colored legs.
+    /// The exact action is `i[P_first,second,P_second,third]`, without an
+    /// additional factor of i in the definition of the returned operator.
+    StructureConstantTripole {
+        first: i64,
+        second: i64,
+        third: i64,
+    },
+}
+
+/// Immutable basis identity, exact insertion matrices and bilinear programs.
+/// Owning the tree prevents coefficients from an unrelated equal-size basis
+/// from being used accidentally. Request order and multiplicity are retained.
+pub struct TreeColorProgram {
+    tree: Arc<TreeProgram>,
+    requests: Vec<TreeColorInsertion>,
+    matrices: Vec<Vec<Vec<ColorNcTerms>>>,
+    self_state: Vec<u8>,
+    directed_state: Vec<u8>,
+}
+
+impl TreeColorProgram {
+    pub fn prepare(
+        tree: Arc<TreeProgram>,
+        requests: Vec<TreeColorInsertion>,
+    ) -> RusticolResult<Self> {
+        let matrices = matrices::prepare_matrices(tree.color_tensors(), &requests)?;
+        let numerical = matrices
+            .iter()
+            .map(|matrix| matrices::at_nc3(matrix))
+            .collect::<RusticolResult<Vec<_>>>()?;
+        let (self_state, directed_state) = runtime::compile_states(&numerical)?;
+        Ok(Self {
+            tree,
+            requests,
+            matrices,
+            self_state,
+            directed_state,
+        })
+    }
+
+    pub fn tree(&self) -> &TreeProgram {
+        &self.tree
+    }
+    pub fn requests(&self) -> &[TreeColorInsertion] {
+        &self.requests
+    }
+    pub fn basis_len(&self) -> usize {
+        self.tree.color_tensors().len()
+    }
+
+    /// Directed `<C_b|O|C_c>` entries as exact rational complex Laurent terms.
+    /// Formal Nc dependence does not extend the numerical model beyond SU(3).
+    pub fn matrix_nc_terms(&self, request: usize) -> RusticolResult<&[Vec<ColorNcTerms>]> {
+        self.matrices
+            .get(request)
+            .map(Vec::as_slice)
+            .ok_or_else(|| {
+                RusticolError::invalid_argument(
+                    "native tree color matrix request index is out of range",
+                )
+            })
+    }
+
+    /// Rebuild numerical coefficients from exact instructions at the requested
+    /// working bits. `Arbitrary,113` is the multiprecision 113-bit tier;
+    /// `DoubleDouble,106` remains the distinct compatibility carrier.
+    pub fn evaluator_with_precision_bits(
+        &self,
+        arithmetic: Arithmetic,
+        bits: u32,
+    ) -> RusticolResult<TreeColorEvaluator> {
+        TreeColorEvaluator::new(self, arithmetic, bits)
+    }
+}
diff --git a/rust/loop-reduction/src/tree_correlations/matrices.rs b/rust/loop-reduction/src/tree_correlations/matrices.rs
new file mode 100644
index 000000000..fee761907
--- /dev/null
+++ b/rust/loop-reduction/src/tree_correlations/matrices.rs
@@ -0,0 +1,552 @@
+// SPDX-License-Identifier: 0BSD
+
+//! Directed exact insertion matrices, built from full ket tensor actions.
+//!
+//! Rows are bras and columns are kets in the caller's retained basis order.
+//! The basis need not be orthogonal or linearly independent. In particular,
+//! overlap matrices are never multiplied or inverted to represent products.
+
+use rusticol_core::{
+    RusticolError, RusticolResult,
+    color::{
+        ColorDipole, ColorDipoleProduct, ColorNcTerms, ConnectedColorTensor, ExactColorCoefficient,
+        evaluate_nc_terms,
+    },
+};
+use symbolica::domains::{float::Complex, rational::Rational};
+
+use super::TreeColorInsertion;
+
+/// Preserve request order and multiplicity, including repeated factors.
+pub(super) fn prepare_matrices(
+    basis: &[ConnectedColorTensor],
+    requests: &[TreeColorInsertion],
+) -> RusticolResult<Vec<Vec<Vec<ColorNcTerms>>>> {
+    let Some(first) = basis.first() else {
+        return Err(RusticolError::invalid_argument(
+            "tree color insertion matrices require a nonempty basis",
+        ));
+    };
+    if basis.iter().any(|tensor| tensor.legs() != first.legs()) {
+        return Err(RusticolError::invalid_argument(
+            "tree color basis tensors must have identical typed legs",
+        ));
+    }
+    if requests.is_empty() {
+        return Err(RusticolError::invalid_argument(
+            "tree color insertion requests must not be empty",
+        ));
+    }
+    let legs = first.legs();
+    requests
+        .iter()
+        .map(|request| match request {
+            TreeColorInsertion::Identity => basis
+                .iter()
+                .map(|bra| basis.iter().map(|ket| bra.overlap_nc_terms(ket)).collect())
+                .collect(),
+            TreeColorInsertion::Pair { left, right } => {
+                let action = ColorDipole::new(legs.to_vec(), *left, *right)?;
+                from_action(basis, |ket| action.apply_connected(ket))
+            }
+            TreeColorInsertion::OrderedProduct(pairs) => {
+                if pairs.is_empty() {
+                    return Err(RusticolError::invalid_argument(
+                        "an ordered tree color product must contain at least one pair; use Identity for the empty product",
+                    ));
+                }
+                let action = ColorDipoleProduct::new(legs.to_vec(), pairs.clone())?;
+                from_action(basis, |ket| action.apply_connected(ket))
+            }
+            TreeColorInsertion::Anticommutator { first, second } => {
+                let ab = ColorDipoleProduct::new(legs.to_vec(), vec![*first, *second])?;
+                let ba = ColorDipoleProduct::new(legs.to_vec(), vec![*second, *first])?;
+                let mut matrix = from_action(basis, |ket| ab.apply_connected(ket))?;
+                let reverse = from_action(basis, |ket| ba.apply_connected(ket))?;
+                combine(&mut matrix, reverse, false);
+                Ok(matrix)
+            }
+            TreeColorInsertion::StructureConstantTripole { first, second, third } => {
+                if first == second || first == third || second == third {
+                    return Err(RusticolError::invalid_argument(
+                        "a structure-constant tree color tripole requires three distinct colored legs",
+                    ));
+                }
+                // The core validates that every label is an original colored
+                // leg and applies mathematical products right-to-left on the
+                // ket. Its independently checked physical convention gives
+                // [P_ij,P_jk] = -i f_abc T_i^a T_j^b T_k^c.
+                let ab = ColorDipoleProduct::new(
+                    legs.to_vec(),
+                    vec![[*first, *second], [*second, *third]],
+                )?;
+                let ba = ColorDipoleProduct::new(
+                    legs.to_vec(),
+                    vec![[*second, *third], [*first, *second]],
+                )?;
+                let mut matrix = from_action(basis, |ket| ab.apply_connected(ket))?;
+                let reverse = from_action(basis, |ket| ba.apply_connected(ket))?;
+                combine(&mut matrix, reverse, true);
+                for terms in matrix.iter_mut().flatten() {
+                    for coefficient in terms.values_mut() {
+                        *coefficient =
+                            Complex::new(-coefficient.im.clone(), coefficient.re.clone());
+                    }
+                }
+                Ok(matrix)
+            }
+        })
+        .collect()
+}
+
+/// Apply the complete action once per ket, then overlap every bra with it.
+fn from_action(
+    basis: &[ConnectedColorTensor],
+    mut action: impl FnMut(&ConnectedColorTensor) -> RusticolResult<ConnectedColorTensor>,
+) -> RusticolResult<Vec<Vec<ColorNcTerms>>> {
+    let mut matrix = vec![vec![ColorNcTerms::new(); basis.len()]; basis.len()];
+    for (column, ket) in basis.iter().enumerate() {
+        let image = action(ket)?;
+        for (row, bra) in basis.iter().enumerate() {
+            matrix[row][column] = bra.overlap_nc_terms(&image)?;
+        }
+    }
+    Ok(matrix)
+}
+
+/// Entrywise exact addition/subtraction; both matrices were constructed from
+/// the same basis. This is not multiplication of operator overlap matrices.
+fn combine(left: &mut [Vec<ColorNcTerms>], right: Vec<Vec<ColorNcTerms>>, subtract: bool) {
+    for (left_row, right_row) in left.iter_mut().zip(right) {
+        for (left_terms, right_terms) in left_row.iter_mut().zip(right_row) {
+            for (power, coefficient) in right_terms {
+                let entry = left_terms
+                    .entry(power)
+                    .or_insert_with(|| Complex::new(Rational::from(0), Rational::from(0)));
+                if subtract {
+                    entry.re -= coefficient.re;
+                    entry.im -= coefficient.im;
+                } else {
+                    entry.re += coefficient.re;
+                    entry.im += coefficient.im;
+                }
+            }
+            left_terms.retain(|_, coefficient| {
+                coefficient.re != Rational::from(0) || coefficient.im != Rational::from(0)
+            });
+        }
+    }
+}
+
+/// The retained TreeProgram amplitude domain is Nc=3. Keep the conversion
+/// exact, including complex phases and negative Laurent powers.
+pub(super) fn at_nc3(
+    matrix: &[Vec<ColorNcTerms>],
+) -> RusticolResult<Vec<Vec<ExactColorCoefficient>>> {
+    if matrix.is_empty() || matrix.iter().any(|row| row.len() != matrix.len()) {
+        return Err(RusticolError::invalid_argument(
+            "a tree color insertion matrix must be nonempty and square",
+        ));
+    }
+    matrix
+        .iter()
+        .map(|row| {
+            row.iter()
+                .map(|terms| evaluate_nc_terms(terms, 3))
+                .collect()
+        })
+        .collect()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use ColorRepresentation::{Antifundamental, Fundamental, Singlet};
+    use TreeColorInsertion::{
+        Anticommutator, Identity, OrderedProduct, Pair, StructureConstantTripole,
+    };
+    use rusticol_core::{
+        RusticolErrorKind,
+        color::{ColorLeg, ColorRepresentation, ColorTensor, OpenColorString, TensorTerm},
+    };
+
+    fn value(re: i64, im: i64) -> ExactColorCoefficient {
+        Complex::new(Rational::from(re), Rational::from(im))
+    }
+
+    fn polynomial(terms: &[(i64, i64, i64, i64)]) -> ColorNcTerms {
+        terms
+            .iter()
+            .map(|&(power, re, im, denominator)| {
+                (
+                    power,
+                    Complex::new(
+                        Rational::from((re, denominator)),
+                        Rational::from((im, denominator)),
+                    ),
+                )
+            })
+            .collect()
+    }
+
+    fn tensor(legs: &[ColorLeg], pairs: &[[i64; 2]]) -> ConnectedColorTensor {
+        ConnectedColorTensor::new(
+            legs.to_vec(),
+            vec![TensorTerm {
+                coefficient: value(1, 0),
+                tensor: ColorTensor::new(
+                    pairs
+                        .iter()
+                        .map(|&[fundamental, antifundamental]| OpenColorString {
+                            fundamental,
+                            adjoints: vec![],
+                            antifundamental,
+                        })
+                        .collect(),
+                    vec![],
+                ),
+                nc_power: 0,
+            }],
+            0,
+        )
+        .unwrap()
+    }
+
+    fn singlet(q: i64, qb: i64) -> ConnectedColorTensor {
+        tensor(
+            &[
+                ColorLeg {
+                    label: q,
+                    representation: Fundamental,
+                },
+                ColorLeg {
+                    label: qb,
+                    representation: Antifundamental,
+                },
+            ],
+            &[[q, qb]],
+        )
+    }
+
+    fn four_quarks() -> Vec<ConnectedColorTensor> {
+        let legs = [
+            ColorLeg {
+                label: 1,
+                representation: Fundamental,
+            },
+            ColorLeg {
+                label: 2,
+                representation: Antifundamental,
+            },
+            ColorLeg {
+                label: 3,
+                representation: Fundamental,
+            },
+            ColorLeg {
+                label: 4,
+                representation: Antifundamental,
+            },
+        ];
+        vec![
+            tensor(&legs, &[[1, 2], [3, 4]]),
+            tensor(&legs, &[[1, 4], [3, 2]]),
+        ]
+    }
+
+    #[test]
+    fn singlet_casimirs_and_request_multiplicity_are_exact() {
+        let basis = [singlet(-13, 29)];
+        let pair = Pair {
+            left: -13,
+            right: 29,
+        };
+        let matrices = prepare_matrices(
+            &basis,
+            &[
+                Identity,
+                Pair {
+                    left: -13,
+                    right: -13,
+                },
+                pair.clone(),
+                pair,
+            ],
+        )
+        .unwrap();
+        assert_eq!(matrices[0][0][0], polynomial(&[(1, 1, 0, 1)]));
+        // <delta|T_q^2|delta> = Nc*CF = (Nc^2-1)/2.
+        assert_eq!(
+            matrices[1][0][0],
+            polynomial(&[(2, 1, 0, 2), (0, -1, 0, 2)])
+        );
+        assert_eq!(
+            matrices[2][0][0],
+            polynomial(&[(2, -1, 0, 2), (0, 1, 0, 2)])
+        );
+        assert_eq!(matrices[2], matrices[3]);
+        assert_eq!(at_nc3(&matrices[1]).unwrap(), vec![vec![value(4, 0)]]);
+        assert_eq!(at_nc3(&matrices[2]).unwrap(), vec![vec![value(-4, 0)]]);
+    }
+
+    #[test]
+    fn four_quark_products_and_tripole_retain_directed_complex_signs() {
+        let basis = four_quarks();
+        let matrices = prepare_matrices(
+            &basis,
+            &[
+                OrderedProduct(vec![[1, 2], [2, 3]]),
+                OrderedProduct(vec![[2, 3], [1, 2]]),
+                StructureConstantTripole {
+                    first: 1,
+                    second: 2,
+                    third: 3,
+                },
+                StructureConstantTripole {
+                    first: 2,
+                    second: 3,
+                    third: 1,
+                },
+                StructureConstantTripole {
+                    first: 2,
+                    second: 1,
+                    third: 3,
+                },
+            ],
+        )
+        .unwrap();
+        // Independent direct/crossed delta-basis formulas, with tr(t^a t^b)
+        // = delta_ab/2. AB is not BA; the product acts on the ket.
+        let ab = vec![
+            vec![
+                ColorNcTerms::new(),
+                polynomial(&[(3, 1, 0, 4), (1, -2, 0, 4), (-1, 1, 0, 4)]),
+            ],
+            vec![
+                polynomial(&[(1, -1, 0, 4), (-1, 1, 0, 4)]),
+                ColorNcTerms::new(),
+            ],
+        ];
+        assert_eq!(matrices[0], ab);
+        for row in 0..2 {
+            for column in 0..2 {
+                assert_eq!(matrices[1][row][column], ab[column][row]);
+            }
+        }
+        // F_123 has off-diagonal values +/-i*Nc*(Nc^2-1)/4.
+        let tripole = vec![
+            vec![
+                ColorNcTerms::new(),
+                polynomial(&[(3, 0, 1, 4), (1, 0, -1, 4)]),
+            ],
+            vec![
+                polynomial(&[(3, 0, -1, 4), (1, 0, 1, 4)]),
+                ColorNcTerms::new(),
+            ],
+        ];
+        assert_eq!(matrices[2], tripole);
+        assert_eq!(matrices[3], tripole);
+        for row in 0..2 {
+            for column in 0..2 {
+                let reversed: ColorNcTerms = tripole[row][column]
+                    .iter()
+                    .map(|(&power, coefficient)| (power, -coefficient.clone()))
+                    .collect();
+                assert_eq!(matrices[4][row][column], reversed);
+            }
+        }
+        assert_eq!(
+            at_nc3(&matrices[2]).unwrap(),
+            vec![
+                vec![value(0, 0), value(0, 6)],
+                vec![value(0, -6), value(0, 0)]
+            ],
+        );
+    }
+
+    #[test]
+    fn repeated_factors_and_anticommutators_have_no_implicit_half() {
+        let basis = [singlet(1, 2)];
+        let matrices = prepare_matrices(
+            &basis,
+            &[
+                OrderedProduct(vec![[1, 2], [1, 2]]),
+                Anticommutator {
+                    first: [1, 2],
+                    second: [1, 2],
+                },
+                OrderedProduct(vec![[1, 2]; 3]),
+            ],
+        )
+        .unwrap();
+        assert_eq!(
+            matrices[0][0][0],
+            polynomial(&[(3, 1, 0, 4), (1, -2, 0, 4), (-1, 1, 0, 4)])
+        );
+        assert_eq!(
+            matrices[1][0][0],
+            polynomial(&[(3, 1, 0, 2), (1, -2, 0, 2), (-1, 1, 0, 2)])
+        );
+        assert_eq!(
+            matrices[2][0][0],
+            polynomial(&[(4, -1, 0, 8), (2, 3, 0, 8), (0, -3, 0, 8), (-2, 1, 0, 8)])
+        );
+
+        let matrices = prepare_matrices(
+            &four_quarks(),
+            &[Anticommutator {
+                first: [1, 2],
+                second: [2, 3],
+            }],
+        )
+        .unwrap();
+        let off_diagonal = polynomial(&[(3, 1, 0, 4), (1, -3, 0, 4), (-1, 2, 0, 4)]);
+        assert_eq!(
+            matrices[0],
+            vec![
+                vec![ColorNcTerms::new(), off_diagonal.clone()],
+                vec![off_diagonal, ColorNcTerms::new()],
+            ]
+        );
+    }
+
+    #[test]
+    fn redundant_complex_basis_keeps_phases_and_exact_normalizations() {
+        let original = four_quarks();
+        let mut terms = vec![];
+        for (source, coefficient) in original.iter().zip([value(2, 1), value(-1, 2)]) {
+            let mut term = source.terms()[0].clone();
+            term.coefficient = coefficient;
+            term.nc_power = -1;
+            terms.push(term);
+        }
+        // The third vector equals [(2+i) C0+(-1+2i) C1]/(2 Nc).
+        // This basis is rank deficient, with nontrivial complex overlaps.
+        let redundant = ConnectedColorTensor::new(original[0].legs().to_vec(), terms, 2).unwrap();
+        let mut basis = original.clone();
+        basis.push(redundant);
+        let requests = [
+            Identity,
+            Pair { left: 1, right: 3 },
+            OrderedProduct(vec![[1, 2], [2, 3]]),
+            StructureConstantTripole {
+                first: 1,
+                second: 2,
+                third: 3,
+            },
+        ];
+        let base = prepare_matrices(&original, &requests).unwrap();
+        let actual = prepare_matrices(&basis, &requests).unwrap();
+        for nc in [2, 3, 5] {
+            let denominator = Rational::from(2 * i64::from(nc));
+            let weights = [
+                [
+                    value(1, 0),
+                    value(0, 0),
+                    Complex::new(
+                        Rational::from(2) / &denominator,
+                        Rational::from(1) / &denominator,
+                    ),
+                ],
+                [
+                    value(0, 0),
+                    value(1, 0),
+                    Complex::new(
+                        Rational::from(-1) / &denominator,
+                        Rational::from(2) / &denominator,
+                    ),
+                ],
+            ];
+            for (small, large) in base.iter().zip(&actual) {
+                for row in 0..3 {
+                    for column in 0..3 {
+                        let mut expected = value(0, 0);
+                        for a in 0..2 {
+                            for b in 0..2 {
+                                expected += weights[a][row].conj()
+                                    * evaluate_nc_terms(&small[a][b], nc).unwrap()
+                                    * &weights[b][column];
+                            }
+                        }
+                        assert_eq!(
+                            evaluate_nc_terms(&large[row][column], nc).unwrap(),
+                            expected
+                        );
+                    }
+                }
+            }
+        }
+        assert_ne!(at_nc3(&actual[0]).unwrap()[0][2].im, Rational::from(0));
+    }
+
+    #[test]
+    fn zero_overlap_does_not_remove_a_nonzero_insertion() {
+        let basis = four_quarks();
+        // C1-C0/Nc is orthogonal to C0, while T1.T3 couples them.
+        let mut terms = basis[1].terms().to_vec();
+        let mut subtraction = basis[0].terms()[0].clone();
+        subtraction.coefficient = value(-1, 0);
+        subtraction.nc_power = -1;
+        terms.push(subtraction);
+        let orthogonal = ConnectedColorTensor::new(basis[0].legs().to_vec(), terms, 0).unwrap();
+        let matrices = prepare_matrices(
+            &[basis[0].clone(), orthogonal],
+            &[Identity, Pair { left: 1, right: 3 }],
+        )
+        .unwrap();
+        assert!(matrices[0][0][1].is_empty());
+        assert!(matrices[0][1][0].is_empty());
+        let off_diagonal = polynomial(&[(2, 1, 0, 2), (0, -1, 0, 2)]);
+        assert_eq!(matrices[1][0][1], off_diagonal);
+        assert_eq!(matrices[1][1][0], off_diagonal);
+    }
+
+    #[test]
+    fn malformed_basis_and_unsupported_insertions_fail_explicitly() {
+        let basis = [singlet(1, 2)];
+        for requests in [
+            vec![],
+            vec![OrderedProduct(vec![])],
+            vec![Pair { left: 1, right: 9 }],
+            vec![StructureConstantTripole {
+                first: 1,
+                second: 1,
+                third: 2,
+            }],
+            vec![StructureConstantTripole {
+                first: 1,
+                second: 2,
+                third: 9,
+            }],
+        ] {
+            assert_eq!(
+                prepare_matrices(&basis, &requests).unwrap_err().kind(),
+                RusticolErrorKind::InvalidArgument
+            );
+        }
+        assert!(prepare_matrices(&[], &[Identity]).is_err());
+        assert!(prepare_matrices(&[basis[0].clone(), singlet(3, 4)], &[Identity]).is_err());
+        let mut legs = basis[0].legs().to_vec();
+        legs.push(ColorLeg {
+            label: 3,
+            representation: Singlet,
+        });
+        let with_singlet = ConnectedColorTensor::new(legs, basis[0].terms().to_vec(), 0).unwrap();
+        assert!(prepare_matrices(&[with_singlet.clone()], &[Pair { left: 1, right: 3 }]).is_err());
+        assert!(
+            prepare_matrices(
+                &[with_singlet],
+                &[StructureConstantTripole {
+                    first: 1,
+                    second: 2,
+                    third: 3
+                }]
+            )
+            .is_err()
+        );
+        let odd = ConnectedColorTensor::new(basis[0].legs().to_vec(), basis[0].terms().to_vec(), 1)
+            .unwrap();
+        assert!(prepare_matrices(&[basis[0].clone(), odd], &[Identity]).is_err());
+        assert!(at_nc3(&[]).is_err());
+        assert!(at_nc3(&[vec![]]).is_err());
+    }
+}
diff --git a/rust/loop-reduction/src/tree_correlations/runtime.rs b/rust/loop-reduction/src/tree_correlations/runtime.rs
new file mode 100644
index 000000000..3cbb88821
--- /dev/null
+++ b/rust/loop-reduction/src/tree_correlations/runtime.rs
@@ -0,0 +1,407 @@
+// SPDX-License-Identifier: 0BSD
+
+use std::sync::Arc;
+
+use rusticol_core::{
+    RusticolError, RusticolResult, color::ExactColorCoefficient,
+    preparation::encode_exact_evaluator_state,
+};
+use symbolica::{
+    atom::{Atom, AtomCore},
+    symbol,
+};
+
+use super::{TreeColorInsertion, TreeColorProgram};
+use crate::{
+    Complex,
+    runtime::{Arithmetic, Batch, BatchCarrier, CoefficientEvaluator, Data},
+    tree_norm::hermitian_norm,
+    tree_preparation::{TreeEvaluator, TreeProgram},
+};
+
+fn invalid(message: impl Into<String>) -> RusticolError {
+    RusticolError::invalid_argument(message)
+}
+
+fn coordinate_symbols(prefix: &str, count: usize) -> Vec<Atom> {
+    (0..count)
+        .flat_map(|i| {
+            ["real", "imag"]
+                .map(move |part| Atom::var(symbol!(&format!("{prefix}_{part}_{i}"); Real)))
+        })
+        .collect()
+}
+
+/// Exact Cartesian expansion of conj(B) K A, including complex off-diagonals.
+fn directed(matrix: &[Vec<ExactColorCoefficient>], bra: &[Atom], ket: &[Atom]) -> (Atom, Atom) {
+    let (mut real, mut imag) = (Atom::Zero, Atom::Zero);
+    for (b, row) in matrix.iter().enumerate() {
+        let (x, y) = (&bra[2 * b], &bra[2 * b + 1]);
+        for (c, k) in row.iter().enumerate() {
+            let (u, v) = (&ket[2 * c], &ket[2 * c + 1]);
+            let p = x * u + y * v;
+            let q = x * v - y * u;
+            real += Atom::num(k.re.clone()) * &p - Atom::num(k.im.clone()) * &q;
+            imag += Atom::num(k.im.clone()) * p + Atom::num(k.re.clone()) * q;
+        }
+    }
+    (real, imag)
+}
+
+pub(super) fn compile_states(
+    matrices: &[Vec<Vec<ExactColorCoefficient>>],
+) -> RusticolResult<(Vec<u8>, Vec<u8>)> {
+    let n = matrices
+        .first()
+        .map(Vec::len)
+        .filter(|&n| n > 0)
+        .ok_or_else(|| invalid("color contraction requires nonempty matrices"))?;
+    n.checked_mul(4)
+        .ok_or_else(|| invalid("color Cartesian dimensions overflow"))?;
+    matrices
+        .len()
+        .checked_mul(2)
+        .ok_or_else(|| invalid("color output dimensions overflow"))?;
+    if matrices
+        .iter()
+        .any(|m| m.len() != n || m.iter().any(|r| r.len() != n))
+    {
+        return Err(invalid(
+            "color insertion matrices must have one common square dimension",
+        ));
+    }
+    let same = coordinate_symbols("rusticol_native_tree::amplitude", n);
+    let bra = coordinate_symbols("rusticol_native_color::bra", n);
+    let ket = coordinate_symbols("rusticol_native_color::ket", n);
+    let (mut self_outputs, mut directed_outputs) = (Vec::new(), Vec::new());
+    for matrix in matrices {
+        let is_hermitian = (0..n).all(|b| (0..n).all(|c| matrix[b][c] == matrix[c][b].conj()));
+        let (real, imag) = if is_hermitian {
+            // Establish reality algebraically before creating instructions.
+            // The same kernel is used by the ordinary prepared tree norm.
+            let (value, inputs) = hermitian_norm(matrix)?;
+            debug_assert_eq!(inputs, same);
+            (value, Atom::Zero)
+        } else {
+            directed(matrix, &same, &same)
+        };
+        self_outputs.extend([real, imag]);
+        let (real, imag) = directed(matrix, &bra, &ket);
+        directed_outputs.extend([real, imag]);
+    }
+    let encode = |outputs: &[Atom], inputs: &[Atom]| {
+        let evaluator = Atom::evaluator_multiple(outputs, inputs)
+            .cores(1)
+            .build()
+            .map_err(|e| {
+                RusticolError::evaluation(format!("native color bilinear preparation: {e}"))
+            })?;
+        encode_exact_evaluator_state(&evaluator)
+    };
+    Ok((
+        encode(&self_outputs, &same)?,
+        encode(&directed_outputs, &[bra, ket].concat())?,
+    ))
+}
+
+fn cartesian_values<T: BatchCarrier>(
+    values: &[Complex<T>],
+    batch: &Batch,
+) -> RusticolResult<Batch> {
+    let values = values
+        .iter()
+        .flat_map(|z| {
+            [
+                Complex::new(z.re.clone(), z.re.zero()),
+                Complex::new(z.im.clone(), z.im.zero()),
+            ]
+        })
+        .collect();
+    let columns = batch
+        .columns()
+        .checked_mul(2)
+        .ok_or_else(|| invalid("Cartesian dimensions overflow"))?;
+    Batch::typed(
+        values,
+        batch.ids().to_vec(),
+        columns,
+        batch.precision_bits(),
+    )
+    .map_err(Into::into)
+}
+
+fn cartesian(batch: &Batch) -> RusticolResult<Batch> {
+    match batch.data() {
+        Data::Double(v) => cartesian_values(v, batch),
+        Data::DoubleDouble(v) => cartesian_values(v, batch),
+        Data::Arbitrary(v) => cartesian_values(v, batch),
+    }
+}
+
+fn pack_values<T: BatchCarrier>(
+    values: &[Complex<T>],
+    batch: &Batch,
+    requests: usize,
+) -> RusticolResult<Batch> {
+    let mut output = Vec::with_capacity(values.len() / 2);
+    for (index, parts) in values.chunks_exact(2).enumerate() {
+        let id = batch.ids()[index / requests];
+        let request = index % requests;
+        if parts.iter().any(|z| !z.re.is_finite() || !z.im.is_finite()) {
+            return Err(RusticolError::nonfinite_result(format!(
+                "native color bilinear produced a nonfinite value in row {id}, request {request}"
+            )));
+        }
+        if parts.iter().any(|z| !z.im.is_fully_zero()) {
+            return Err(RusticolError::integrity(format!(
+                "Cartesian color expression is not real in row {id}, request {request}"
+            )));
+        }
+        // These are separate exact real expressions, not a projection of the
+        // complex contraction. Preserve each Float and its tracked precision.
+        output.push(Complex::new(parts[0].re.clone(), parts[1].re.clone()));
+    }
+    Batch::typed(
+        output,
+        batch.ids().to_vec(),
+        requests,
+        batch.precision_bits(),
+    )
+    .map_err(Into::into)
+}
+
+fn pack(batch: Batch, requests: usize) -> RusticolResult<Batch> {
+    if requests.checked_mul(2) != Some(batch.columns()) {
+        return Err(RusticolError::integrity(
+            "native color Cartesian output shape changed",
+        ));
+    }
+    match batch.data() {
+        Data::Double(v) => pack_values(v, &batch, requests),
+        Data::DoubleDouble(v) => pack_values(v, &batch, requests),
+        Data::Arbitrary(v) => pack_values(v, &batch, requests),
+    }
+}
+
+fn add_values<T: BatchCarrier>(
+    a: &[Complex<T>],
+    b: &[Complex<T>],
+    batch: &Batch,
+) -> RusticolResult<Batch> {
+    let values = a
+        .iter()
+        .zip(b)
+        .enumerate()
+        .map(|(i, (a, b))| {
+            let value = a + b;
+            if !value.re.is_finite() || !value.im.is_finite() {
+                Err(RusticolError::nonfinite_result(format!(
+                    "native color helicity sum is nonfinite in row {}, request {}",
+                    batch.ids()[i / batch.columns()],
+                    i % batch.columns()
+                )))
+            } else {
+                Ok(value)
+            }
+        })
+        .collect::<RusticolResult<Vec<_>>>()?;
+    Batch::typed(
+        values,
+        batch.ids().to_vec(),
+        batch.columns(),
+        batch.precision_bits(),
+    )
+    .map_err(Into::into)
+}
+
+fn add(a: &Batch, b: &Batch) -> RusticolResult<Batch> {
+    a.same_rows(b)?;
+    if a.columns() != b.columns() {
+        return Err(RusticolError::integrity("color sum columns changed"));
+    }
+    match (a.data(), b.data()) {
+        (Data::Double(a_values), Data::Double(b_values)) => add_values(a_values, b_values, a),
+        (Data::DoubleDouble(a_values), Data::DoubleDouble(b_values)) => {
+            add_values(a_values, b_values, a)
+        }
+        (Data::Arbitrary(a_values), Data::Arbitrary(b_values)) => add_values(a_values, b_values, a),
+        _ => unreachable!("same_rows checked the carrier"),
+    }
+}
+
+/// Evaluates color insertions with the same owned tree and numerical carrier.
+/// Each output row has one complex column per requested insertion, in order.
+pub struct TreeColorEvaluator {
+    program: Arc<TreeProgram>,
+    requests: Vec<TreeColorInsertion>,
+    tree: TreeEvaluator,
+    same: CoefficientEvaluator,
+    directed: CoefficientEvaluator,
+}
+
+impl TreeColorEvaluator {
+    pub(super) fn new(
+        program: &TreeColorProgram,
+        arithmetic: Arithmetic,
+        bits: u32,
+    ) -> RusticolResult<Self> {
+        let tree = program
+            .tree
+            .evaluator_with_precision_bits(arithmetic, bits)?;
+        if tree.coefficient_count() != program.basis_len() {
+            return Err(RusticolError::integrity(
+                "tree coefficient and color basis dimensions differ",
+            ));
+        }
+        Ok(Self {
+            program: Arc::clone(&program.tree),
+            requests: program.requests.clone(),
+            tree,
+            same: CoefficientEvaluator::with_precision_bits(&program.self_state, arithmetic, bits)?,
+            directed: CoefficientEvaluator::with_precision_bits(
+                &program.directed_state,
+                arithmetic,
+                bits,
+            )?,
+        })
+    }
+
+    pub fn tree(&self) -> &TreeProgram {
+        &self.program
+    }
+    pub fn requests(&self) -> &[TreeColorInsertion] {
+        &self.requests
+    }
+    pub fn precision_bits(&self) -> u32 {
+        self.same.precision_bits()
+    }
+
+    fn self_contract(&mut self, coefficients: &Batch) -> RusticolResult<Batch> {
+        let inputs = cartesian(coefficients)?;
+        let output = self.same.evaluate_batch(&inputs)?;
+        pack(output, self.requests.len())
+    }
+
+    /// One physical helicity assignment, with the same amplitude on both sides.
+    pub fn evaluate_fixed_helicities(
+        &mut self,
+        momenta: &Batch,
+        helicities: &[i32],
+    ) -> RusticolResult<Batch> {
+        self.fixed(momenta, None, helicities)
+    }
+    pub fn evaluate_fixed_helicities_with_parameters(
+        &mut self,
+        momenta: &Batch,
+        parameters: &Batch,
+        helicities: &[i32],
+    ) -> RusticolResult<Batch> {
+        self.fixed(momenta, Some(parameters), helicities)
+    }
+    fn fixed(
+        &mut self,
+        momenta: &Batch,
+        parameters: Option<&Batch>,
+        helicities: &[i32],
+    ) -> RusticolResult<Batch> {
+        let coefficients = match parameters {
+            Some(p) => self
+                .tree
+                .evaluate_coefficients_with_parameters(momenta, p, helicities)?,
+            None => self.tree.evaluate_coefficients(momenta, helicities)?,
+        };
+        self.self_contract(&coefficients)
+    }
+
+    /// Directed `<bra|O|ket>` for independent supplied external component blocks.
+    /// Each block uses the prepared tree's four-component/adjoint convention.
+    /// Even for Hermitian O, the directed result is kept complex.
+    pub fn evaluate_external_components(
+        &mut self,
+        momenta: &Batch,
+        bra: &[&Batch],
+        ket: &[&Batch],
+    ) -> RusticolResult<Batch> {
+        self.components(momenta, None, bra, ket)
+    }
+    pub fn evaluate_external_components_with_parameters(
+        &mut self,
+        momenta: &Batch,
+        parameters: &Batch,
+        bra: &[&Batch],
+        ket: &[&Batch],
+    ) -> RusticolResult<Batch> {
+        self.components(momenta, Some(parameters), bra, ket)
+    }
+    fn components(
+        &mut self,
+        momenta: &Batch,
+        parameters: Option<&Batch>,
+        bra: &[&Batch],
+        ket: &[&Batch],
+    ) -> RusticolResult<Batch> {
+        let mut evaluate = |states: &[&Batch]| match parameters {
+            Some(p) => self
+                .tree
+                .evaluate_external_components_with_parameters(momenta, p, states),
+            None => self.tree.evaluate_external_components(momenta, states),
+        };
+        let bra = evaluate(bra)?;
+        let ket = evaluate(ket)?;
+        bra.same_rows(&ket)?;
+        let inputs = Batch::concat_columns(&[&cartesian(&bra)?, &cartesian(&ket)?])?;
+        let output = self.directed.evaluate_batch(&inputs)?;
+        pack(output, self.requests.len())
+    }
+
+    /// Stream the full declared physical helicity product, without averages.
+    pub fn evaluate_batch(&mut self, momenta: &Batch) -> RusticolResult<Batch> {
+        self.sum(momenta, None)
+    }
+    pub fn evaluate_batch_with_parameters(
+        &mut self,
+        momenta: &Batch,
+        parameters: &Batch,
+    ) -> RusticolResult<Batch> {
+        self.sum(momenta, Some(parameters))
+    }
+    fn sum(&mut self, momenta: &Batch, parameters: Option<&Batch>) -> RusticolResult<Batch> {
+        let helicities: Vec<_> = self
+            .tree
+            .states()
+            .iter()
+            .map(|s| s.helicities.clone())
+            .collect();
+        let mut positions = vec![0; helicities.len()];
+        let mut accumulated: Option<Batch> = None;
+        loop {
+            let assignment: Vec<_> = helicities
+                .iter()
+                .zip(&positions)
+                .map(|(h, &p)| h[p])
+                .collect();
+            let value = self.fixed(momenta, parameters, &assignment)?;
+            // Validate all tree metadata even for zero rows, then return a
+            // correctly shaped typed empty result without enumerating spins.
+            if momenta.rows() == 0 {
+                return Ok(value);
+            }
+            accumulated = Some(match accumulated {
+                None => value,
+                Some(a) => add(&a, &value)?,
+            });
+            let mut p = positions.len();
+            while p > 0 {
+                p -= 1;
+                positions[p] += 1;
+                if positions[p] < helicities[p].len() {
+                    break;
+                }
+                positions[p] = 0;
+            }
+            if p == 0 && positions.first().is_none_or(|&x| x == 0) {
+                return Ok(accumulated.expect("one assignment in every declared helicity product"));
+            }
+        }
+    }
+}
diff --git a/rust/loop-reduction/src/tree_correlations/tests.rs b/rust/loop-reduction/src/tree_correlations/tests.rs
new file mode 100644
index 000000000..ef7b03260
--- /dev/null
+++ b/rust/loop-reduction/src/tree_correlations/tests.rs
@@ -0,0 +1,551 @@
+// SPDX-License-Identifier: 0BSD
+
+use std::{collections::BTreeMap, sync::Arc};
+
+use feynkit_generator::{GenerationOptions, Process};
+use feynkit_model::Model;
+use rusticol_core::{RusticolErrorKind, color::ExactColorCoefficient};
+use symbolica::{
+    atom::{Atom, Symbol},
+    domains::{
+        float::{DoubleFloat, Float, FloatLike, Real, SingleFloat},
+        rational::Rational,
+    },
+    symbol,
+};
+
+use super::*;
+use crate::{
+    Complex,
+    runtime::{Arithmetic, Batch, BatchCarrier, CoefficientEvaluator, Data, wavefunctions_batch},
+    tree_preparation::{TreeFieldClass, TreeProgram},
+};
+
+fn z(re: i64, im: i64) -> ExactColorCoefficient {
+    Complex::new(re.into(), im.into())
+}
+
+fn exact_bilinear(
+    matrix: &[Vec<ExactColorCoefficient>],
+    bra: &[ExactColorCoefficient],
+    ket: &[ExactColorCoefficient],
+) -> ExactColorCoefficient {
+    let mut result = z(0, 0);
+    for b in 0..bra.len() {
+        for c in 0..ket.len() {
+            result += bra[b].conj() * &matrix[b][c] * &ket[c];
+        }
+    }
+    result
+}
+
+fn coordinates<T: BatchCarrier>(
+    bra: &[ExactColorCoefficient],
+    ket: Option<&[ExactColorCoefficient]>,
+    prototype: T,
+    bits: u32,
+) -> Batch {
+    let values: Vec<_> = bra
+        .iter()
+        .chain(ket.into_iter().flatten())
+        .flat_map(|z| [&z.re, &z.im])
+        .map(|r| Complex::new(prototype.from_rational(r), prototype.zero()))
+        .collect();
+    let columns = values.len();
+    Batch::typed(values, vec![91], columns, bits).unwrap()
+}
+
+#[test]
+fn cartesian_bilinears_match_independent_exact_complex_double_sum_all_widths() {
+    let matrices = vec![
+        vec![vec![z(2, 0), z(1, 2)], vec![z(1, -2), z(3, 0)]],
+        vec![vec![z(1, 1), z(-2, 3)], vec![z(4, -1), z(0, 2)]],
+    ];
+    let bra = vec![
+        Complex::new((1, 3).into(), (2, 7).into()),
+        Complex::new((-5, 11).into(), (3, 13).into()),
+    ];
+    let ket = vec![
+        Complex::new((4, 9).into(), (-2, 5).into()),
+        Complex::new((7, 17).into(), (-3, 19).into()),
+    ];
+    let (same, directed) = runtime::compile_states(&matrices).unwrap();
+    for (state, other) in [(&same, None), (&directed, Some(ket.as_slice()))] {
+        let expected: Vec<_> = matrices
+            .iter()
+            .map(|m| exact_bilinear(m, &bra, other.unwrap_or(&bra)))
+            .flat_map(|z| [z.re, z.im])
+            .collect();
+        let out = CoefficientEvaluator::with_precision_bits(state, Arithmetic::Double, 53)
+            .unwrap()
+            .evaluate_batch(&coordinates(&bra, other, 0_f64, 53))
+            .unwrap();
+        let Data::Double(v) = out.data() else {
+            panic!("carrier changed")
+        };
+        for (actual, expected) in v.iter().zip(&expected) {
+            assert_eq!(actual.im, 0.);
+            assert!((actual.re - expected.to_f64()).abs() < 1e-13);
+        }
+        let out = CoefficientEvaluator::with_precision_bits(state, Arithmetic::DoubleDouble, 106)
+            .unwrap()
+            .evaluate_batch(&coordinates(&bra, other, DoubleFloat::from(0.), 106))
+            .unwrap();
+        let Data::DoubleDouble(values) = out.data() else {
+            panic!("carrier changed")
+        };
+        let prototype = DoubleFloat::from(0.);
+        for (actual, expected) in values.iter().zip(&expected) {
+            assert!(actual.im.is_fully_zero());
+            assert!(
+                (actual.re.clone() - prototype.from_rational(expected)).norm()
+                    < DoubleFloat::from(2.).pow(90).inv()
+            );
+        }
+        for bits in [113, 137, 256, 512, 1024] {
+            let prototype = Float::with_val(bits, 0);
+            let out = CoefficientEvaluator::with_precision_bits(state, Arithmetic::Arbitrary, bits)
+                .unwrap()
+                .evaluate_batch(&coordinates(&bra, other, prototype.clone(), bits))
+                .unwrap();
+            assert_eq!(out.precision_bits(), bits);
+            let Data::Arbitrary(v) = out.data() else {
+                panic!("carrier changed")
+            };
+            for (actual, expected) in v.iter().zip(&expected) {
+                assert!(actual.im.is_fully_zero());
+                assert!(
+                    (actual.re.clone() - prototype.from_rational(expected)).norm()
+                        < Float::with_val(bits, 2).pow(u64::from(bits - 20)).inv()
+                );
+            }
+        }
+    }
+    // A Hermitian operator on independent vectors is not constrained to real.
+    assert_ne!(
+        exact_bilinear(&matrices[0], &bra, &ket).im,
+        Rational::from(0)
+    );
+    // A non-Hermitian self product retains its physical imaginary component.
+    assert_ne!(
+        exact_bilinear(&matrices[1], &bra, &bra).im,
+        Rational::from(0)
+    );
+}
+
+#[test]
+fn exact_hermitian_lowering_preserves_sub_binary64_phase_difference() {
+    let matrix = vec![vec![z(0, 0), z(0, 1)], vec![z(0, -1), z(0, 0)]];
+    let epsilon = Rational::from((1, 2)).pow(80);
+    let vector = vec![
+        z(1, 1),
+        Complex::new(1.into(), Rational::from(1) + &epsilon),
+    ];
+    let expected = exact_bilinear(&matrix, &vector, &vector);
+    assert_eq!(
+        expected,
+        Complex::new(-Rational::from(2) * epsilon, Rational::from(0))
+    );
+    let (same, _) = runtime::compile_states(&[matrix]).unwrap();
+    let low = CoefficientEvaluator::with_precision_bits(&same, Arithmetic::Double, 53)
+        .unwrap()
+        .evaluate_batch(&coordinates(&vector, None, 0_f64, 53))
+        .unwrap();
+    let Data::Double(low) = low.data() else {
+        panic!()
+    };
+    assert_eq!(low[0].re, 0.);
+    for bits in [113, 137, 256, 512, 1024] {
+        let zero = Float::with_val(bits, 0);
+        let output = CoefficientEvaluator::with_precision_bits(&same, Arithmetic::Arbitrary, bits)
+            .unwrap()
+            .evaluate_batch(&coordinates(&vector, None, zero.clone(), bits))
+            .unwrap();
+        let Data::Arbitrary(v) = output.data() else {
+            panic!()
+        };
+        assert_ne!(v[0].re, zero);
+        assert!(
+            (v[0].re.clone() - zero.from_rational(&expected.re)).norm()
+                < Float::with_val(bits, 2).pow(u64::from(bits - 10)).inv()
+        );
+        assert!(v[1].re.is_fully_zero() && v[0].im.is_fully_zero() && v[1].im.is_fully_zero());
+    }
+}
+
+fn qcd(incoming: [&str; 2], outgoing: [&str; 2], runtime_g: bool) -> Arc<TreeProgram> {
+    let model = Arc::new(Model::qcd());
+    let mut parameters: BTreeMap<Symbol, Atom> = BTreeMap::new();
+    if !runtime_g {
+        parameters.insert(symbol!("UFO::G"), Atom::one());
+    }
+    for particle in model.particles().iter().filter(|p| p.spin == 2) {
+        let mass = model.parameter_by_id(particle.mass).unwrap();
+        parameters.insert(symbol!(&format!("UFO::{}", mass.name)), Atom::Zero);
+    }
+    let runtime = if runtime_g {
+        vec![crate::tree_preparation::RuntimeTreeParameter {
+            symbol: symbol!("UFO::G"),
+            domain: crate::tree_preparation::TreeParameterDomain::Real,
+        }]
+    } else {
+        vec![]
+    };
+    Arc::new(
+        TreeProgram::prepare_with_runtime_parameters(
+            model,
+            &Process::new(incoming, outgoing),
+            &GenerationOptions::default().threads(1),
+            None,
+            &parameters,
+            TreeFieldClass::MasslessQcdAndSingletScalars,
+            &runtime,
+        )
+        .unwrap(),
+    )
+}
+
+fn momenta<T: BatchCarrier>(prototype: T, bits: u32) -> Batch {
+    let values = [5, 0, 0, 5, 5, 0, 0, -5, 5, 4, 0, 3, 5, -4, 0, -3]
+        .into_iter()
+        .map(|v| {
+            Complex::new(
+                prototype.from_rational(&Rational::from(v)),
+                prototype.zero(),
+            )
+        })
+        .collect();
+    Batch::typed(values, vec![17], 16, bits).unwrap()
+}
+
+#[test]
+fn generated_color_identities_crossings_and_runtime_parameters() {
+    for (incoming, outgoing) in [
+        (["u", "d"], ["u", "d"]),
+        (["u", "u~"], ["d", "d~"]),
+        (["u", "g"], ["u", "g"]),
+        (["g", "g"], ["g", "g"]),
+    ] {
+        let tree = qcd(incoming, outgoing, true);
+        let legs = tree.color_tensors()[0].legs();
+        let first = legs[0].label;
+        let mut requests = vec![TreeColorInsertion::Identity];
+        requests.extend(legs.iter().map(|leg| TreeColorInsertion::Pair {
+            left: first,
+            right: leg.label,
+        }));
+        requests.push(TreeColorInsertion::Identity);
+        let color = TreeColorProgram::prepare(tree.clone(), requests.clone()).unwrap();
+        assert!(std::ptr::eq(color.tree(), tree.as_ref()));
+        assert_eq!(color.requests(), requests);
+        assert!(color.matrix_nc_terms(requests.len()).is_err());
+        let p = momenta(0_f64, 53);
+        let params = Batch::typed(vec![Complex::new(1., 0.)], vec![17], 1, 53).unwrap();
+        let norm = tree
+            .evaluator_with_precision_bits(Arithmetic::Double, 53)
+            .unwrap()
+            .evaluate_batch_with_parameters(&p, &params)
+            .unwrap();
+        let mut evaluator = color
+            .evaluator_with_precision_bits(Arithmetic::Double, 53)
+            .unwrap();
+        let result = evaluator
+            .evaluate_batch_with_parameters(&p, &params)
+            .unwrap();
+        let (Data::Double(v), Data::Double(norm)) = (result.data(), norm.data()) else {
+            panic!()
+        };
+        assert_eq!(result.ids(), [17]);
+        // A multi-output evaluator can share/reassociate exact subexpressions
+        // differently from the one-output norm. Only numerical rounding may differ.
+        assert!((v[0].re - norm[0].re).abs() <= 8. * f64::EPSILON * norm[0].re.abs());
+        assert_eq!(v[0].im, norm[0].im);
+        assert_eq!(v[0], v[v.len() - 1]);
+        assert!(v.iter().all(|z| z.im == 0.));
+        let casimir = match legs[0].representation {
+            rusticol_core::color::ColorRepresentation::Adjoint => 3.,
+            _ => 4. / 3.,
+        };
+        assert!((v[1].re - casimir * v[0].re).abs() < 1e-10 * v[0].re.abs());
+        assert!(v[1..=legs.len()].iter().map(|z| z.re).sum::<f64>().abs() < 1e-10 * v[0].re.abs());
+        let params2 = Batch::typed(vec![Complex::new(2., 0.)], vec![17], 1, 53).unwrap();
+        let result2 = evaluator
+            .evaluate_batch_with_parameters(&p, &params2)
+            .unwrap();
+        let Data::Double(v2) = result2.data() else {
+            panic!()
+        };
+        for (a, b) in v.iter().zip(v2.iter()) {
+            assert!((b.re - 16. * a.re).abs() < 1e-9 * v[0].re.abs());
+        }
+    }
+}
+
+#[test]
+fn independent_external_components_keep_directed_complex_phase() {
+    let tree = qcd(["u", "d"], ["u", "d"], false);
+    let color =
+        TreeColorProgram::prepare(tree.clone(), vec![TreeColorInsertion::Identity]).unwrap();
+    let p = momenta(0_f64, 53);
+    let mut native = tree
+        .evaluator_with_precision_bits(Arithmetic::Double, 53)
+        .unwrap();
+    let components: Vec<_> = native
+        .states()
+        .iter()
+        .map(|state| {
+            wavefunctions_batch(
+                &p.take_columns(state.momentum_columns.to_vec()).unwrap(),
+                vec![state.helicities[0]],
+                vec![state.kind.clone()],
+                Some(4),
+            )
+            .unwrap()
+        })
+        .collect();
+    let helicities: Vec<_> = native.states().iter().map(|s| s.helicities[0]).collect();
+    let mut evaluator = color
+        .evaluator_with_precision_bits(Arithmetic::Double, 53)
+        .unwrap();
+    let norm = evaluator
+        .evaluate_fixed_helicities(&p, &helicities)
+        .unwrap();
+    let Data::Double(norm) = norm.data() else {
+        panic!()
+    };
+    assert!(norm[0].re > 0.);
+    let mut bra = components.clone();
+    let mut ket = components.clone();
+    let phase = |batch: &Batch, z: Complex<f64>| {
+        let Data::Double(v) = batch.data() else {
+            panic!()
+        };
+        Batch::typed(
+            v.iter().map(|v| v * &z).collect(),
+            batch.ids().to_vec(),
+            4,
+            53,
+        )
+        .unwrap()
+    };
+    bra[0] = phase(&bra[0], Complex::new(0., 2.));
+    ket[0] = phase(&ket[0], Complex::new(3., -4.));
+    let output = evaluator
+        .evaluate_external_components(
+            &p,
+            &bra.iter().collect::<Vec<_>>(),
+            &ket.iter().collect::<Vec<_>>(),
+        )
+        .unwrap();
+    let Data::Double(value) = output.data() else {
+        panic!()
+    };
+    assert!((value[0].re + 8. * norm[0].re).abs() < 1e-10 * norm[0].re);
+    assert!((value[0].im + 6. * norm[0].re).abs() < 1e-10 * norm[0].re);
+    let reverse = evaluator
+        .evaluate_external_components(
+            &p,
+            &ket.iter().collect::<Vec<_>>(),
+            &bra.iter().collect::<Vec<_>>(),
+        )
+        .unwrap();
+    let Data::Double(reverse) = reverse.data() else {
+        panic!()
+    };
+    assert!((reverse[0].re - value[0].re).abs() < 1e-10 * norm[0].re);
+    assert!((reverse[0].im + value[0].im).abs() < 1e-10 * norm[0].re);
+    // Metadata checks include independent component blocks on empty rows.
+    let empty = Batch::typed::<f64>(vec![], vec![], 16, 53).unwrap();
+    let valid: Vec<_> = (0..components.len())
+        .map(|_| Batch::typed::<f64>(vec![], vec![], 4, 53).unwrap())
+        .collect();
+    let invalid = Batch::typed::<f64>(vec![], vec![], 3, 53).unwrap();
+    let empty_result = evaluator
+        .evaluate_external_components(
+            &empty,
+            &valid.iter().collect::<Vec<_>>(),
+            &valid.iter().collect::<Vec<_>>(),
+        )
+        .unwrap();
+    assert_eq!((empty_result.rows(), empty_result.columns()), (0, 1));
+    let mut refs: Vec<_> = valid.iter().collect();
+    refs[0] = &invalid;
+    assert!(
+        evaluator
+            .evaluate_external_components(&empty, &valid.iter().collect::<Vec<_>>(), &refs)
+            .is_err()
+    );
+    // The native coefficient method remains available and is not replaced by a density.
+    assert_eq!(
+        native
+            .evaluate_coefficients(&p, &helicities)
+            .unwrap()
+            .columns(),
+        color.basis_len()
+    );
+}
+
+#[test]
+fn arbitrary_directed_facade_retains_a_subbinary64_external_phase() {
+    let bits = 137;
+    let tree = qcd(["u", "d"], ["u", "d"], false);
+    let color =
+        TreeColorProgram::prepare(tree.clone(), vec![TreeColorInsertion::Identity]).unwrap();
+    let p = momenta(Float::with_val(bits, 0), bits);
+    let native = tree
+        .evaluator_with_precision_bits(Arithmetic::Arbitrary, bits)
+        .unwrap();
+    let states = native.states();
+    let helicities: Vec<_> = states.iter().map(|s| s.helicities[0]).collect();
+    let components: Vec<_> = states
+        .iter()
+        .map(|state| {
+            wavefunctions_batch(
+                &p.take_columns(state.momentum_columns.to_vec()).unwrap(),
+                vec![state.helicities[0]],
+                vec![state.kind.clone()],
+                Some(4),
+            )
+            .unwrap()
+        })
+        .collect();
+    let mut evaluator = color
+        .evaluator_with_precision_bits(Arithmetic::Arbitrary, bits)
+        .unwrap();
+    let norm = evaluator
+        .evaluate_fixed_helicities(&p, &helicities)
+        .unwrap();
+    let Data::Arbitrary(norm) = norm.data() else {
+        panic!("carrier changed")
+    };
+    let epsilon = Float::with_val(bits, 2).pow(80).inv();
+    let mut bra = components.clone();
+    let mut ket = components.clone();
+    let phase = |batch: &Batch, z: Complex<Float>| {
+        let Data::Arbitrary(v) = batch.data() else {
+            panic!("carrier changed")
+        };
+        Batch::typed(
+            v.iter().map(|v| v * &z).collect(),
+            batch.ids().to_vec(),
+            4,
+            bits,
+        )
+        .unwrap()
+    };
+    bra[0] = phase(
+        &bra[0],
+        Complex::new(Float::with_val(bits, 0), Float::with_val(bits, 2)),
+    );
+    ket[0] = phase(
+        &ket[0],
+        Complex::new(
+            Float::with_val(bits, 3),
+            -(Float::with_val(bits, 4) + &epsilon),
+        ),
+    );
+    let value = evaluator
+        .evaluate_external_components(
+            &p,
+            &bra.iter().collect::<Vec<_>>(),
+            &ket.iter().collect::<Vec<_>>(),
+        )
+        .unwrap();
+    assert_eq!(
+        (value.precision_bits(), value.columns(), value.ids()),
+        (bits, 1, &[17][..])
+    );
+    let Data::Arbitrary(values) = value.data() else {
+        panic!("carrier changed")
+    };
+    let expected_re =
+        -(Float::with_val(bits, 8) + Float::with_val(bits, 2) * &epsilon) * &norm[0].re;
+    let expected_im = -Float::with_val(bits, 6) * &norm[0].re;
+    let tolerance = Float::with_val(bits, 2).pow(115).inv() * &norm[0].re;
+    assert!((values[0].re.clone() - expected_re).norm() < tolerance);
+    assert!((values[0].im.clone() - expected_im).norm() < tolerance);
+    // Promoting binary64 components would erase this exact phase perturbation.
+    assert!(
+        (values[0].re.clone() + Float::with_val(bits, 8) * &norm[0].re).norm()
+            > epsilon * &norm[0].re
+    );
+}
+
+#[test]
+fn facade_precision_empty_metadata_and_retryable_overflow() {
+    let tree = qcd(["u", "d"], ["u", "d"], true);
+    let program = TreeColorProgram::prepare(tree, vec![TreeColorInsertion::Identity]).unwrap();
+    let mut double = program
+        .evaluator_with_precision_bits(Arithmetic::Double, 53)
+        .unwrap();
+    let p = momenta(0_f64, 53);
+    let params = |g| Batch::typed(vec![Complex::new(g, 0.)], vec![17], 1, 53).unwrap();
+    let error = double
+        .evaluate_batch_with_parameters(&p, &params(1e100))
+        .err()
+        .unwrap();
+    assert_eq!(error.kind(), RusticolErrorKind::NonFiniteResult);
+    assert_ne!(
+        double
+            .evaluate_batch_with_parameters(&p, &params(f64::INFINITY))
+            .err()
+            .unwrap()
+            .kind(),
+        RusticolErrorKind::NonFiniteResult
+    );
+    let wrong_ids = Batch::typed(vec![Complex::new(1., 0.)], vec![18], 1, 53).unwrap();
+    assert!(
+        double
+            .evaluate_batch_with_parameters(&p, &wrong_ids)
+            .is_err()
+    );
+    for bits in [113, 137, 256, 512, 1024] {
+        let mut evaluator = program
+            .evaluator_with_precision_bits(Arithmetic::Arbitrary, bits)
+            .unwrap();
+        let p = momenta(Float::with_val(bits, 0), bits);
+        let parameters = Batch::typed(
+            vec![Complex::new(
+                Float::parse("1e100", Some(bits)).unwrap(),
+                Float::with_val(bits, 0),
+            )],
+            vec![17],
+            1,
+            bits,
+        )
+        .unwrap();
+        let value = evaluator
+            .evaluate_batch_with_parameters(&p, &parameters)
+            .unwrap();
+        assert!(value.nonfinite_rows().is_empty());
+        assert_eq!(value.precision_bits(), bits);
+        let Data::Arbitrary(v) = value.data() else {
+            panic!("carrier changed")
+        };
+        assert!(v[0].im.is_fully_zero());
+        let empty = Batch::typed::<Float>(vec![], vec![], 16, bits).unwrap();
+        let ep = Batch::typed::<Float>(vec![], vec![], 1, bits).unwrap();
+        let result = evaluator
+            .evaluate_batch_with_parameters(&empty, &ep)
+            .unwrap();
+        assert_eq!((result.rows(), result.columns()), (0, 1));
+        assert!(evaluator.evaluate_batch(&empty).is_err());
+        let wrong = Batch::typed::<Float>(vec![], vec![], 16, bits + 1).unwrap();
+        assert!(
+            evaluator
+                .evaluate_batch_with_parameters(&wrong, &ep)
+                .is_err()
+        );
+        let wrong = Batch::typed::<f64>(vec![], vec![], 16, bits).unwrap();
+        assert!(
+            evaluator
+                .evaluate_batch_with_parameters(&wrong, &ep)
+                .is_err()
+        );
+    }
+    assert!(
+        program
+            .evaluator_with_precision_bits(Arithmetic::Double, 113)
+            .is_err()
+    );
+}