From cc3558de1141eb6afbac821a1673147991162311 Mon Sep 17 00:00:00 2001 From: lukassimon Date: Fri, 23 Oct 2026 23:20:43 +0000 Subject: [PATCH 1/2] Align native dependencies with future-null kinematics --- Cargo.lock | 34 ++++++++++++++-------------- rust/crates/rusticol-core/Cargo.toml | 2 +- rust/loop-reduction/Cargo.toml | 22 +++++++++--------- 3 files changed, 29 insertions(+), 29 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index f2dbbe74f..0400357fd 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -819,7 +819,7 @@ dependencies = [ [[package]] name = "dot-parser" version = "0.5.1" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "either", "pest", @@ -945,7 +945,7 @@ dependencies = [ [[package]] name = "feynkit-amplitude" version = "0.1.0" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "feynkit-graph", "feynkit-model", @@ -959,7 +959,7 @@ dependencies = [ [[package]] name = "feynkit-generator" version = "0.1.0" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "bincode", "bincode-trait-derive", @@ -979,7 +979,7 @@ dependencies = [ [[package]] name = "feynkit-graph" version = "0.1.0" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "feynkit-kinematics", "feynkit-model", @@ -998,7 +998,7 @@ dependencies = [ [[package]] name = "feynkit-kinematics" version = "0.1.0" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "linnet", "numerica", @@ -1011,7 +1011,7 @@ dependencies = [ [[package]] name = "feynkit-model" version = "0.1.0" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "bincode", "bincode-trait-derive", @@ -1099,7 +1099,7 @@ dependencies = [ [[package]] name = "gammaloop-workspace-hack" version = "0.1.0" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "ahash", "arrayvec", @@ -1312,7 +1312,7 @@ dependencies = [ [[package]] name = "idenso" version = "0.3.0" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "ahash", "base64", @@ -1475,7 +1475,7 @@ dependencies = [ [[package]] name = "linnet" version = "0.17.0" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "ahash", "bincode", @@ -1645,7 +1645,7 @@ dependencies = [ [[package]] name = "ninjars" version = "0.1.0" -source = "git+https://github.com/alphal00p/ninjars.git?rev=a069504e3dc884d6e2682aadc75fb12a7f4e7c07#a069504e3dc884d6e2682aadc75fb12a7f4e7c07" +source = "git+https://github.com/alphal00p/ninjars.git?rev=084e466bd464bf4cf847fe5708503d62bc835ca9#084e466bd464bf4cf847fe5708503d62bc835ca9" dependencies = [ "feynkit-kinematics", "symbolica", @@ -1655,7 +1655,7 @@ dependencies = [ [[package]] name = "ninjars-graph" version = "0.1.0" -source = "git+https://github.com/alphal00p/ninjars.git?rev=a069504e3dc884d6e2682aadc75fb12a7f4e7c07#a069504e3dc884d6e2682aadc75fb12a7f4e7c07" +source = "git+https://github.com/alphal00p/ninjars.git?rev=084e466bd464bf4cf847fe5708503d62bc835ca9#084e466bd464bf4cf847fe5708503d62bc835ca9" dependencies = [ "feynkit-amplitude", "feynkit-generator", @@ -1674,7 +1674,7 @@ dependencies = [ [[package]] name = "ninjars-positional" version = "0.1.0" -source = "git+https://github.com/alphal00p/ninjars.git?rev=a069504e3dc884d6e2682aadc75fb12a7f4e7c07#a069504e3dc884d6e2682aadc75fb12a7f4e7c07" +source = "git+https://github.com/alphal00p/ninjars.git?rev=084e466bd464bf4cf847fe5708503d62bc835ca9#084e466bd464bf4cf847fe5708503d62bc835ca9" dependencies = [ "symbolica", "thiserror", @@ -1683,7 +1683,7 @@ dependencies = [ [[package]] name = "ninjars-symbolica" version = "0.1.0" -source = "git+https://github.com/alphal00p/ninjars.git?rev=a069504e3dc884d6e2682aadc75fb12a7f4e7c07#a069504e3dc884d6e2682aadc75fb12a7f4e7c07" +source = "git+https://github.com/alphal00p/ninjars.git?rev=084e466bd464bf4cf847fe5708503d62bc835ca9#084e466bd464bf4cf847fe5708503d62bc835ca9" dependencies = [ "ninjars", "ninjars-positional", @@ -2660,7 +2660,7 @@ checksum = "02ef091d4e8b981dbaa4cc17789d00b53d5d2c001ae689f14e6093f30c529036" [[package]] name = "spenso" version = "0.6.0" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "ahash", "append-only-vec", @@ -2703,7 +2703,7 @@ dependencies = [ [[package]] name = "spenso-hep-lib" version = "0.1.7" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "ahash", "idenso", @@ -2716,7 +2716,7 @@ dependencies = [ [[package]] name = "spenso-macros" version = "0.3.2" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "proc-macro2", "quote", @@ -2796,7 +2796,7 @@ dependencies = [ [[package]] name = "symbolica-utils" version = "0.1.0" -source = "git+https://github.com/alphal00p/gammaloop.git?rev=ee18b8e6386cb36ca88b1917e00ab2ebae93e68e#ee18b8e6386cb36ca88b1917e00ab2ebae93e68e" +source = "git+https://github.com/alphal00p/gammaloop.git?rev=6cbe1b4d296eb78767a55b9ce9cd25ef67feed41#6cbe1b4d296eb78767a55b9ce9cd25ef67feed41" dependencies = [ "bincode", "bincode-trait-derive", diff --git a/rust/crates/rusticol-core/Cargo.toml b/rust/crates/rusticol-core/Cargo.toml index 09c26c288..3782181ee 100644 --- a/rust/crates/rusticol-core/Cargo.toml +++ b/rust/crates/rusticol-core/Cargo.toml @@ -29,7 +29,7 @@ loop-v2 = ["dep:ninjars", "symbolica-runtime"] bincode.workspace = true libloading = { version = "0.9", optional = true } libc.workspace = true -ninjars = { git = "https://github.com/alphal00p/ninjars.git", rev = "a069504e3dc884d6e2682aadc75fb12a7f4e7c07", optional = true } +ninjars = { git = "https://github.com/alphal00p/ninjars.git", rev = "084e466bd464bf4cf847fe5708503d62bc835ca9", optional = true } num-complex = { version = "0.4", optional = true } rayon.workspace = true serde.workspace = true diff --git a/rust/loop-reduction/Cargo.toml b/rust/loop-reduction/Cargo.toml index 1ba05ef78..48022895a 100644 --- a/rust/loop-reduction/Cargo.toml +++ b/rust/loop-reduction/Cargo.toml @@ -9,18 +9,18 @@ license = "0BSD" publish = false [dependencies] -feynkit-amplitude = { git = "https://github.com/alphal00p/gammaloop.git", rev = "ee18b8e6386cb36ca88b1917e00ab2ebae93e68e", default-features = false, features = ["native"] } -feynkit-generator = { git = "https://github.com/alphal00p/gammaloop.git", rev = "ee18b8e6386cb36ca88b1917e00ab2ebae93e68e", default-features = false, features = ["native"] } -feynkit-graph = { git = "https://github.com/alphal00p/gammaloop.git", rev = "ee18b8e6386cb36ca88b1917e00ab2ebae93e68e", default-features = false, features = ["native"] } -idenso = { git = "https://github.com/alphal00p/gammaloop.git", rev = "ee18b8e6386cb36ca88b1917e00ab2ebae93e68e", default-features = false, features = ["native"] } -spenso = { git = "https://github.com/alphal00p/gammaloop.git", rev = "ee18b8e6386cb36ca88b1917e00ab2ebae93e68e", default-features = false, features = ["native", "shadowing"] } -spenso-hep-lib = { git = "https://github.com/alphal00p/gammaloop.git", rev = "ee18b8e6386cb36ca88b1917e00ab2ebae93e68e", default-features = false, features = ["native"] } -feynkit-model = { git = "https://github.com/alphal00p/gammaloop.git", rev = "ee18b8e6386cb36ca88b1917e00ab2ebae93e68e", default-features = false, features = ["native"] } -feynkit-kinematics = { git = "https://github.com/alphal00p/gammaloop.git", rev = "ee18b8e6386cb36ca88b1917e00ab2ebae93e68e", default-features = false, features = ["native"] } +feynkit-amplitude = { git = "https://github.com/alphal00p/gammaloop.git", rev = "6cbe1b4d296eb78767a55b9ce9cd25ef67feed41", default-features = false, features = ["native"] } +feynkit-generator = { git = "https://github.com/alphal00p/gammaloop.git", rev = "6cbe1b4d296eb78767a55b9ce9cd25ef67feed41", default-features = false, features = ["native"] } +feynkit-graph = { git = "https://github.com/alphal00p/gammaloop.git", rev = "6cbe1b4d296eb78767a55b9ce9cd25ef67feed41", default-features = false, features = ["native"] } +idenso = { git = "https://github.com/alphal00p/gammaloop.git", rev = "6cbe1b4d296eb78767a55b9ce9cd25ef67feed41", default-features = false, features = ["native"] } +spenso = { git = "https://github.com/alphal00p/gammaloop.git", rev = "6cbe1b4d296eb78767a55b9ce9cd25ef67feed41", default-features = false, features = ["native", "shadowing"] } +spenso-hep-lib = { git = "https://github.com/alphal00p/gammaloop.git", rev = "6cbe1b4d296eb78767a55b9ce9cd25ef67feed41", default-features = false, features = ["native"] } +feynkit-model = { git = "https://github.com/alphal00p/gammaloop.git", rev = "6cbe1b4d296eb78767a55b9ce9cd25ef67feed41", default-features = false, features = ["native"] } +feynkit-kinematics = { git = "https://github.com/alphal00p/gammaloop.git", rev = "6cbe1b4d296eb78767a55b9ce9cd25ef67feed41", default-features = false, features = ["native"] } rusticol-core = { path = "../crates/rusticol-core", features = ["symbolica-runtime", "python-generation-bridge", "loop-v2"] } -ninjars = { git = "https://github.com/alphal00p/ninjars.git", rev = "a069504e3dc884d6e2682aadc75fb12a7f4e7c07" } -ninjars-graph = { git = "https://github.com/alphal00p/ninjars.git", rev = "a069504e3dc884d6e2682aadc75fb12a7f4e7c07" } -ninjars-symbolica = { git = "https://github.com/alphal00p/ninjars.git", rev = "a069504e3dc884d6e2682aadc75fb12a7f4e7c07", features = ["native-masters"] } +ninjars = { git = "https://github.com/alphal00p/ninjars.git", rev = "084e466bd464bf4cf847fe5708503d62bc835ca9" } +ninjars-graph = { git = "https://github.com/alphal00p/ninjars.git", rev = "084e466bd464bf4cf847fe5708503d62bc835ca9" } +ninjars-symbolica = { git = "https://github.com/alphal00p/ninjars.git", rev = "084e466bd464bf4cf847fe5708503d62bc835ca9", features = ["native-masters"] } oneloop = { git = "https://github.com/alphal00p/oneloopmaster.git", rev = "8ec256fc26cb361b20e10f59c19beb35717ef32c", default-features = false } symbolica = { version = "3.0.1", default-features = false, features = ["integer-gmp", "float-mpfr"] } serde = { version = "1", features = ["derive"] } From 7663e820424ceeda13f1a437af7359480a564e04 Mon Sep 17 00:00:00 2001 From: lukassimon Date: Fri, 23 Oct 2026 23:35:13 +0000 Subject: [PATCH 2/2] Bind constructive null states to prepared tree programs --- rust/loop-reduction/src/lib.rs | 1 + .../src/runtime/coefficients.rs | 13 +- rust/loop-reduction/src/runtime/mod.rs | 22 +- rust/loop-reduction/src/runtime/reduction.rs | 13 +- rust/loop-reduction/src/runtime/states.rs | 32 +- .../src/tree_correlations/runtime.rs | 57 +- rust/loop-reduction/src/tree_null_states.rs | 358 ++++++ .../src/tree_null_states/tests.rs | 1026 +++++++++++++++++ rust/loop-reduction/src/tree_preparation.rs | 54 +- rust/loop-reduction/src/tree_runtime.rs | 64 +- 10 files changed, 1580 insertions(+), 60 deletions(-) create mode 100644 rust/loop-reduction/src/tree_null_states.rs create mode 100644 rust/loop-reduction/src/tree_null_states/tests.rs diff --git a/rust/loop-reduction/src/lib.rs b/rust/loop-reduction/src/lib.rs index 959f9e1c9..db127e609 100644 --- a/rust/loop-reduction/src/lib.rs +++ b/rust/loop-reduction/src/lib.rs @@ -17,6 +17,7 @@ pub mod tree_correlations; mod tree_color; mod tree_components; mod tree_norm; +pub mod tree_null_states; pub mod tree_preparation; mod tree_runtime; pub mod verification; diff --git a/rust/loop-reduction/src/runtime/coefficients.rs b/rust/loop-reduction/src/runtime/coefficients.rs index 1a6707f18..83185d0b3 100644 --- a/rust/loop-reduction/src/runtime/coefficients.rs +++ b/rust/loop-reduction/src/runtime/coefficients.rs @@ -96,18 +96,7 @@ impl CoefficientEvaluator { arithmetic: Arithmetic, bits: u32, ) -> ExecutionResult { - match arithmetic { - Arithmetic::Double if bits != 53 => { - return Err(error("double arithmetic requires exactly 53 bits")); - } - Arithmetic::DoubleDouble if bits != 106 => { - return Err(error("double-double arithmetic requires exactly 106 bits")); - } - Arithmetic::Arbitrary if bits < 2 => { - return Err(error("arbitrary arithmetic requires at least two bits")); - } - _ => {} - } + validate_working_precision(arithmetic, bits)?; let exact = exact_state(state)?; let evaluator = match arithmetic { Arithmetic::Double => Coefficients::Double(batch_evaluator(exact, &0_f64)), diff --git a/rust/loop-reduction/src/runtime/mod.rs b/rust/loop-reduction/src/runtime/mod.rs index 64dc50671..8232dbed2 100644 --- a/rust/loop-reduction/src/runtime/mod.rs +++ b/rust/loop-reduction/src/runtime/mod.rs @@ -29,6 +29,7 @@ pub use master_cache::{MasterCache, MasterStatistics}; pub use ninjars::fast::tensor_layout::monomials as tensor_monomials; pub use reduction::{DecimalMomentum, PolynomialPlan, RawLane, ReductionBatch}; pub use schedule::{CoefficientKernelCache, CoefficientPlan}; +pub(crate) use states::wavefunction_error; pub use states::{ PartialStates, on_shell_wavefunctions_batch, wavefunctions_batch, wavefunctions_partial_batch, }; @@ -45,6 +46,25 @@ pub fn register_scalar_primitives() { } pub type DecimalComplex = [String; 2]; pub type ExecutionResult = Result; + +// One shared working-carrier policy. This is declared attempt precision, not +// a lower bound on the reliable bits of values after numerical cancellation. +pub(crate) fn validate_working_precision(arithmetic: Arithmetic, bits: u32) -> ExecutionResult<()> { + match arithmetic { + Arithmetic::Double if bits != 53 => { + return Err(error("double arithmetic requires exactly 53 bits")); + } + Arithmetic::DoubleDouble if bits != 106 => { + return Err(error("double-double arithmetic requires exactly 106 bits")); + } + Arithmetic::Arbitrary if bits < 2 => { + return Err(error("arbitrary arithmetic requires at least two bits")); + } + _ => {} + } + Ok(()) +} + #[derive(Debug, Clone)] pub struct ExecutionError(rusticol_core::RusticolError); impl ExecutionError { @@ -63,7 +83,7 @@ impl ExecutionError { value.to_string(), )) } - fn context(self, value: impl Display) -> Self { + pub(crate) fn context(self, value: impl Display) -> Self { Self(rusticol_core::RusticolError::with_kind( self.kind(), value.to_string(), diff --git a/rust/loop-reduction/src/runtime/reduction.rs b/rust/loop-reduction/src/runtime/reduction.rs index c5833846b..44eea6567 100644 --- a/rust/loop-reduction/src/runtime/reduction.rs +++ b/rust/loop-reduction/src/runtime/reduction.rs @@ -341,18 +341,7 @@ impl PolynomialPlan { bits: u32, engine: &str, ) -> ExecutionResult { - match arithmetic { - Arithmetic::Double if bits != 53 => { - return Err(error("double arithmetic requires exactly 53 bits")); - } - Arithmetic::DoubleDouble if bits != 106 => { - return Err(error("double-double arithmetic requires exactly 106 bits")); - } - Arithmetic::Arbitrary if bits < 2 => { - return Err(error("arbitrary arithmetic requires at least two bits")); - } - _ => {} - } + validate_working_precision(arithmetic, bits)?; if !matches!(engine, "generic" | "fast" | "fast-tensor") { return Err(error("engine must be generic, fast or fast-tensor")); } diff --git a/rust/loop-reduction/src/runtime/states.rs b/rust/loop-reduction/src/runtime/states.rs index e3406c54a..8f89f161b 100644 --- a/rust/loop-reduction/src/runtime/states.rs +++ b/rust/loop-reduction/src/runtime/states.rs @@ -60,20 +60,24 @@ fn state_components( PhysicalMomentum::from(std::array::from_fn::<_, 4, _>(|i| p[i].re.clone())) .wavefunction(kind, helicity) .map(|wavefunction| wavefunction.into_components()) - .map_err(|failure| match failure { - WavefunctionError::NonFinite => ExecutionError::nonfinite_result(failure), - WavefunctionError::InsufficientPrecision { reason, .. } => { - // Input components may have lower tracked precision after prior - // arithmetic. Preserve the batch's original working precision. - ExecutionError::insufficient_precision( - working_bits, - format!( - "external-state chart needs more than {working_bits} working bits: {reason}" - ), - ) - } - _ => error(failure), - }) + .map_err(|failure| wavefunction_error(failure, working_bits)) +} + +/// Preserve the original attempt precision and typed numerical category for +/// both generic and constructive source-specific FeynKit external states. +pub(crate) fn wavefunction_error(failure: WavefunctionError, working_bits: u32) -> ExecutionError { + match failure { + WavefunctionError::NonFinite => ExecutionError::nonfinite_result(failure), + WavefunctionError::InsufficientPrecision { reason, .. } => { + ExecutionError::insufficient_precision( + working_bits, + format!( + "external-state chart needs more than {working_bits} working bits: {reason}" + ), + ) + } + _ => error(failure), + } } fn evaluate( diff --git a/rust/loop-reduction/src/tree_correlations/runtime.rs b/rust/loop-reduction/src/tree_correlations/runtime.rs index 3cbb88821..1c9260ec5 100644 --- a/rust/loop-reduction/src/tree_correlations/runtime.rs +++ b/rust/loop-reduction/src/tree_correlations/runtime.rs @@ -2,6 +2,8 @@ use std::sync::Arc; +use feynkit_kinematics::KinematicScalar; + use rusticol_core::{ RusticolError, RusticolResult, color::ExactColorCoefficient, preparation::encode_exact_evaluator_state, @@ -16,7 +18,9 @@ use crate::{ Complex, runtime::{Arithmetic, Batch, BatchCarrier, CoefficientEvaluator, Data}, tree_norm::hermitian_norm, + tree_null_states::FutureNullTreeBatch, tree_preparation::{TreeEvaluator, TreeProgram}, + tree_runtime::{RawTreeStateSource, TreeStateSource}, }; fn invalid(message: impl Into) -> RusticolError { @@ -304,15 +308,45 @@ impl TreeColorEvaluator { parameters: Option<&Batch>, helicities: &[i32], ) -> RusticolResult { - let coefficients = match parameters { - Some(p) => self - .tree - .evaluate_coefficients_with_parameters(momenta, p, helicities)?, - None => self.tree.evaluate_coefficients(momenta, helicities)?, - }; + self.fixed_source(&RawTreeStateSource(momenta), parameters, helicities) + } + + fn fixed_source( + &mut self, + source: &dyn TreeStateSource, + parameters: Option<&Batch>, + helicities: &[i32], + ) -> RusticolResult { + let coefficients = self + .tree + .evaluate_coefficients_from_source(source, parameters, helicities)?; self.self_contract(&coefficients) } + /// One physical helicity assignment from a constructive future-null batch + /// bound to this same immutable tree. The result retains directed complex + /// operator values; no source precision or conservation certificate is implied. + pub fn evaluate_future_null_fixed_helicities( + &mut self, + source: &FutureNullTreeBatch, + parameters: Option<&Batch>, + helicities: &[i32], + ) -> RusticolResult { + source.validate_program(&self.program)?; + self.fixed_source(source, parameters, helicities) + } + + /// Stream the existing helicity product for program-bound future-null rows. + /// All external masses were proved identically zero at preparation. + pub fn evaluate_future_null_batch( + &mut self, + source: &FutureNullTreeBatch, + parameters: Option<&Batch>, + ) -> RusticolResult { + source.validate_program(&self.program)?; + self.sum_source(source, parameters) + } + /// Directed `` 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. @@ -366,6 +400,15 @@ impl TreeColorEvaluator { self.sum(momenta, Some(parameters)) } fn sum(&mut self, momenta: &Batch, parameters: Option<&Batch>) -> RusticolResult { + self.sum_source(&RawTreeStateSource(momenta), parameters) + } + + fn sum_source( + &mut self, + source: &dyn TreeStateSource, + parameters: Option<&Batch>, + ) -> RusticolResult { + let momenta = source.momenta(); let helicities: Vec<_> = self .tree .states() @@ -380,7 +423,7 @@ impl TreeColorEvaluator { .zip(&positions) .map(|(h, &p)| h[p]) .collect(); - let value = self.fixed(momenta, parameters, &assignment)?; + let value = self.fixed_source(source, 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 { diff --git a/rust/loop-reduction/src/tree_null_states.rs b/rust/loop-reduction/src/tree_null_states.rs new file mode 100644 index 000000000..9d085cc29 --- /dev/null +++ b/rust/loop-reduction/src/tree_null_states.rs @@ -0,0 +1,358 @@ +// SPDX-License-Identifier: 0BSD + +//! Program-bound constructive massless external states. +//! +//! Admission uses the exact external masses retained by native preparation. +//! The immutable FeynKit owners supply both coefficient momenta and external +//! wavefunctions. This association does not certify event conservation, source +//! precision, gauge completeness or numerical accuracy. Rebuild owners from the +//! original spatial sources when changing working precision. + +use std::{any::TypeId, collections::BTreeSet, sync::Arc}; + +use feynkit_generator::Process; +use feynkit_graph::ExternalState; +use feynkit_kinematics::{FutureNullMomentum, Helicity, KinematicScalar, WavefunctionKind}; +use feynkit_model::{Model, ParticleId}; +use rusticol_core::{RusticolError, RusticolResult}; +use symbolica::{ + atom::AtomCore, + domains::float::{DoubleFloat, Float}, +}; + +use crate::{ + Complex, + runtime::{ + Arithmetic, Batch, BatchCarrier, shape, validate_working_precision, wavefunction_error, + }, + tree_preparation::TreeProgram, + tree_runtime::{ExternalTreeState, TreeEvaluator, TreeStateSource}, +}; + +#[cfg(test)] +mod tests; + +fn error(message: impl std::fmt::Display) -> RusticolError { + RusticolError::invalid_argument(format!("future-null tree states: {message}")) +} + +fn carrier_arithmetic() -> RusticolResult { + let actual = TypeId::of::(); + if actual == TypeId::of::() { + Ok(Arithmetic::Double) + } else if actual == TypeId::of::() { + Ok(Arithmetic::DoubleDouble) + } else if actual == TypeId::of::() { + Ok(Arithmetic::Arbitrary) + } else { + Err(RusticolError::unsupported_precision( + "future-null tree states require a native f64, DoubleFloat or Float batch carrier", + )) + } +} + +/// The complete labelled external domain of one immutable prepared program. +/// +/// Every external mass must specialize identically to exact zero, including +/// scalar spectators. Runtime-dependent or unresolved masses are rejected even +/// when a numerical default is zero. Internal supported masses are unrestricted. +/// Massive external particles require a future mixed-leg owner; the existing +/// raw tree evaluator remains available for its declared domain. +#[derive(Clone)] +pub struct FutureNullTreeDomain { + program: Arc, +} + +impl FutureNullTreeDomain { + pub(crate) fn new(program: Arc) -> RusticolResult { + let legs = &program.external_domain.legs; + if legs.len() != program.momentum_legs().len() + || legs + .iter() + .zip(program.momentum_legs()) + .any(|(leg, label)| leg.label != *label) + || legs + .iter() + .map(|leg| leg.label) + .collect::>() + .len() + != legs.len() + { + return Err(error( + "prepared external occurrence inventory is inconsistent", + )); + } + for leg in legs { + if !leg.resolved_mass.is_zero() { + return Err(RusticolError::unsupported_runtime_capability( + "native-tree-future-null", + format!( + "external occurrence {} has nonzero or unresolved specialized mass {}; all external masses must be identically zero", + leg.label, leg.resolved_mass + ), + )); + } + } + Ok(Self { program }) + } + + /// Native external labels in coefficient-column order, including scalars. + pub fn momentum_legs(&self) -> &[usize] { + self.program.momentum_legs() + } + + /// Prepared model namespace for the retained native particle identities. + pub fn model(&self) -> &Model { + &self.program.external_domain.model + } + + /// Original process supplied to this program's native preparation. + pub fn process(&self) -> &Process { + &self.program.external_domain.process + } + + /// Immutable native occurrence, particle and physical crossing role. + pub fn external_legs(&self) -> impl Iterator + '_ { + self.program + .external_domain + .legs + .iter() + .map(|leg| (leg.label, leg.particle, leg.direction)) + } + + /// Bind complete rows of constructive owners to this particular program. + /// + /// Every row supplies each external label exactly once, in any order. Row + /// identities must be unique. The canonical momentum batch is derived from + /// these same consumed owners, never from an independently supplied vector. + /// Working bits are attempt metadata; tracked Float precision may be lower + /// after arithmetic and is neither padded nor treated as an accuracy proof. + pub fn bind( + &self, + row_ids: Vec, + keyed_legs: Vec)>>, + working_bits: u32, + ) -> RusticolResult> { + validate_working_precision(carrier_arithmetic::()?, working_bits) + .map_err(RusticolError::from)?; + if row_ids.len() != keyed_legs.len() { + return Err(error("row identities and external-state rows differ")); + } + if row_ids.iter().copied().collect::>().len() != row_ids.len() { + return Err(error("row identities must be unique")); + } + let labels = self.momentum_legs(); + let expected = labels.iter().copied().collect::>(); + let columns = labels + .len() + .checked_mul(4) + .ok_or_else(|| error("momentum column dimensions overflow"))?; + let size = shape(row_ids.len(), columns).map_err(RusticolError::from)?; + // Validate every complete label inventory before canonicalizing any row. + for (row, entries) in keyed_legs.iter().enumerate() { + let actual = entries + .iter() + .map(|(label, _)| *label) + .collect::>(); + if entries.len() != labels.len() || actual != expected { + return Err(error(format!( + "lane {} (row {row}) requires a complete bijection of prepared external labels", + row_ids[row] + ))); + } + } + let mut values = Vec::new(); + values + .try_reserve_exact(size) + .map_err(|e| error(format!("momentum allocation failed: {e}")))?; + let mut legs = Vec::new(); + legs.try_reserve_exact(row_ids.len()) + .map_err(|e| error(format!("external-state row allocation failed: {e}")))?; + for entries in keyed_legs { + let mut keyed = entries + .into_iter() + .collect::>(); + let mut row = Vec::new(); + row.try_reserve_exact(labels.len()) + .map_err(|e| error(format!("external-state leg allocation failed: {e}")))?; + for label in labels { + let owner = keyed + .remove(label) + .expect("complete label bijection was checked"); + values.extend( + owner + .as_momentum() + .into_iter() + .map(|value| Complex::new(value.clone(), value.zero())), + ); + row.push(owner); + } + legs.push(row); + } + let momenta = + Batch::typed(values, row_ids, columns, working_bits).map_err(RusticolError::from)?; + Ok(FutureNullTreeBatch { + program: Arc::clone(&self.program), + legs, + momenta, + }) + } + + /// Reuse the prepared coefficient and norm instructions at any supported + /// working precision. No state-specific source compilation is performed. + pub fn evaluator_with_precision_bits( + &self, + arithmetic: Arithmetic, + bits: u32, + ) -> RusticolResult { + Ok(FutureNullTreeEvaluator { + inner: self + .program + .evaluator_with_precision_bits(arithmetic, bits)?, + program: Arc::clone(&self.program), + }) + } +} + +/// Immutable rows whose owners and coefficient momenta share one program. +/// +/// Shared views cannot replace components or retain this association after a +/// scalar conversion. Cloning the raw [`Batch`] returned by [`Self::momenta`] +/// does not create another program-bound source. +/// +/// ```compile_fail +/// use rusticol_loop_reduction::tree_null_states::FutureNullTreeBatch; +/// fn change(batch: &mut FutureNullTreeBatch) { +/// batch.leg(0, 1).unwrap().as_momentum().temporal.value = 7.0; +/// } +/// ``` +pub struct FutureNullTreeBatch { + program: Arc, + legs: Vec>>, + momenta: Batch, +} + +impl FutureNullTreeBatch { + pub fn momenta(&self) -> &Batch { + &self.momenta + } + + /// Shared owner for a row index and native external occurrence label. + pub fn leg(&self, row: usize, label: usize) -> Option<&FutureNullMomentum> { + let column = self + .program + .momentum_legs() + .iter() + .position(|value| *value == label)?; + self.legs.get(row)?.get(column) + } + + /// Identity is deliberately stronger than equal process/shape metadata. + /// Independent preparations cannot exchange admitted sources, including + /// empty batches. Domain clones from the same program are interchangeable. + pub(crate) fn validate_program(&self, program: &Arc) -> RusticolResult<()> { + if !Arc::ptr_eq(&self.program, program) { + return Err(error( + "source belongs to a different prepared program instance", + )); + } + Ok(()) + } +} + +impl TreeStateSource for FutureNullTreeBatch { + fn momenta(&self) -> &Batch { + &self.momenta + } + + fn components(&self, state: &ExternalTreeState, helicity: i32) -> RusticolResult { + // Match all four columns and the immutable prepared role/domain. A raw + // evaluator may use noncanonical tuples; its first column is no proof + // that the remaining columns describe the same physical occurrence. + let (position, leg, prepared) = self + .program + .external_domain + .legs + .iter() + .enumerate() + .filter_map(|(position, leg)| { + let prepared = self.program.external_states().get(leg.state_index?)?; + Some((position, leg, prepared)) + }) + .find(|(position, _, prepared)| { + let columns = std::array::from_fn(|component| position * 4 + component); + state.momentum_columns == columns + && prepared.momentum_columns == columns + && state.kind == prepared.kind + && state.helicities == prepared.helicities + }) + .ok_or_else(|| { + error("external-state tuple or role differs from its prepared occurrence") + })?; + if !prepared.helicities.contains(&helicity) { + return Err(error( + "helicity differs from the prepared external-state domain", + )); + } + let kind = prepared.kind.parse::().map_err(error)?; + let helicity = Helicity::try_from(i8::try_from(helicity).map_err(error)?).map_err(error)?; + let columns = kind.component_count(); + let size = shape(self.momenta.rows(), columns).map_err(RusticolError::from)?; + let mut values = Vec::new(); + values + .try_reserve_exact(size) + .map_err(|e| error(format!("external-component allocation failed: {e}")))?; + for (row, owners) in self.legs.iter().enumerate() { + let components = owners[position] + .wavefunction(kind, helicity) + .map_err(|failure| { + let failure = wavefunction_error(failure, self.momenta.precision_bits()); + let message = format!( + "future-null external state lane {} (row {row}), occurrence {}: {failure}", + self.momenta.ids()[row], + leg.label + ); + RusticolError::from(failure.context(message)) + })?; + values.extend(components.into_components()); + } + Batch::typed( + values, + self.momenta.ids().to_vec(), + columns, + self.momenta.precision_bits(), + ) + .map_err(RusticolError::from) + } +} + +/// Prepared tree evaluator accepting only sources bound to its program. +/// Existing coefficient contraction and helicity-sum algorithms are reused. +pub struct FutureNullTreeEvaluator { + program: Arc, + inner: TreeEvaluator, +} + +impl FutureNullTreeEvaluator { + pub fn evaluate_coefficients( + &mut self, + source: &FutureNullTreeBatch, + parameters: Option<&Batch>, + helicities: &[i32], + ) -> RusticolResult { + source.validate_program(&self.program)?; + self.inner + .evaluate_coefficients_from_source(source, parameters, helicities) + } + + /// Return the unaveraged, unsymmetrized color norm summed over the prepared + /// physical helicities. Result values retain the selected numerical carrier. + pub fn evaluate_batch( + &mut self, + source: &FutureNullTreeBatch, + parameters: Option<&Batch>, + ) -> RusticolResult { + source.validate_program(&self.program)?; + self.inner.evaluate_batch_from_source(source, parameters) + } +} diff --git a/rust/loop-reduction/src/tree_null_states/tests.rs b/rust/loop-reduction/src/tree_null_states/tests.rs new file mode 100644 index 000000000..79051ac50 --- /dev/null +++ b/rust/loop-reduction/src/tree_null_states/tests.rs @@ -0,0 +1,1026 @@ +// SPDX-License-Identifier: 0BSD + +//! Finite component and program-association checks, not accuracy certificates. + +use std::{collections::BTreeMap, sync::Arc}; + +use feynkit_generator::{GenerationOptions, Process}; +use feynkit_graph::ExternalState; +use feynkit_kinematics::{FutureNullMomentum, Helicity, KinematicScalar, ThreeMomentum}; +use feynkit_model::Model; +use rusticol_core::color::{ColorRepresentation, evaluate_nc_terms}; +use rusticol_core::{RusticolErrorKind, RusticolResult}; +use symbolica::{ + atom::{Atom, Symbol}, + domains::{ + float::{DoubleFloat, Float, FloatLike, Real}, + rational::Rational, + }, + symbol, +}; + +use crate::{ + Complex, + runtime::{Arithmetic, Batch, BatchCarrier, Data}, + tree_correlations::{TreeColorInsertion, TreeColorProgram}, + tree_preparation::{RuntimeTreeParameter, TreeFieldClass, TreeParameterDomain, TreeProgram}, + tree_runtime::TreeStateSource, +}; + +fn scalar_definition() -> serde_json::Value { + // The same native scalar-contact fixture class as tree_preparation tests. + // A separate width parameter allows the name ZERO itself to carry mass. + serde_json::json!({ + "name":"owned-scalar-contact", "restriction":null, + "orders":[{"name":"SELF","expansion_order":99,"hierarchy":1}], + "parameters":[ + {"name":"ZERO","lhablock":null,"lhacode":null,"nature":"internal","parameter_type":"real","value":[0.,0.],"expression":null}, + {"name":"g_","lhablock":"SCALAR","lhacode":[1],"nature":"external","parameter_type":"real","value":[1.,0.],"expression":null}, + {"name":"m","lhablock":"SCALAR","lhacode":[2],"nature":"external","parameter_type":"real","value":[4.,0.],"expression":null}, + {"name":"width","lhablock":null,"lhacode":null,"nature":"internal","parameter_type":"real","value":[0.,0.],"expression":null} + ], + "particles":[{"pdg_code":100,"name":"s","antiname":"s","spin":1,"color":1,"mass":"m","width":"width","texname":"s","antitexname":"s","charge":0.,"ghost_number":0,"lepton_number":0,"y_charge":0.,"propagating":true,"goldstoneboson":false,"propagator":"s_prop"}], + "propagators":[{"name":"s_prop","particle":"s","numerator":"1","denominator":"P^2-m^2"}], + "lorentz_structures":[{"name":"SSSS","spins":[1,1,1,1],"structure":"1"}], + "couplings":[{"name":"GC1","expression":"g_","orders":[["SELF",1]],"value":[1.,0.]}], + "vertex_rules":[{"name":"V1","particles":["s","s","s","s"],"color_structures":["1"],"lorentz_structures":["SSSS"],"couplings":[["GC1"]]}], + "functions":[],"form_factors":[] + }) +} + +fn scalar_program( + definition: serde_json::Value, + overrides: BTreeMap, + runtime: &[RuntimeTreeParameter], +) -> RusticolResult> { + TreeProgram::prepare_with_runtime_parameters( + Arc::new(Model::from_json(&definition.to_string()).unwrap()), + &Process::new(["s", "s"], ["s", "s"]), + &GenerationOptions::default().threads(1), + None, + &overrides, + TreeFieldClass::SingletScalars, + runtime, + ) + .map(Arc::new) +} + +fn massless_scalar() -> Arc { + scalar_program( + scalar_definition(), + BTreeMap::from([(symbol!("UFO::m"), Atom::Zero)]), + &[], + ) + .unwrap() +} + +fn qcd(incoming: [&str; 2], outgoing: [&str; 2], runtime_g: bool) -> Arc { + let model = Arc::new(Model::qcd()); + let mut exact = BTreeMap::new(); + if !runtime_g { + exact.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(); + exact.insert(symbol!(&format!("UFO::{}", mass.name)), Atom::Zero); + } + let runtime = if runtime_g { + vec![RuntimeTreeParameter { + symbol: symbol!("UFO::G"), + domain: TreeParameterDomain::Real, + }] + } else { + vec![] + }; + Arc::new( + TreeProgram::prepare_with_runtime_parameters( + model, + &Process::new(incoming, outgoing), + &GenerationOptions::default().threads(1), + None, + &exact, + TreeFieldClass::MasslessQcdAndSingletScalars, + &runtime, + ) + .unwrap(), + ) +} + +fn rational(prototype: &T, numerator: i64, denominator: i64) -> T { + prototype.from_i64(numerator) / prototype.from_i64(denominator) +} + +fn keyed_row( + program: &TreeProgram, + prototype: &T, + boosted: bool, +) -> Vec<(usize, FutureNullMomentum)> { + // Rotation of the existing (5,4,0,3) fixture by cos(phi)=3/5, + // sin(phi)=4/5. The second row has a transverse beta_y=3/5 boost. + // Both retain s=100, t=-20, u=-80, with exact rational spatial sources. + let spatial = if boosted { + [ + [(0, 1), (-15, 4), (5, 1)], + [(0, 1), (-15, 4), (-5, 1)], + [(12, 5), (1, 4), (3, 1)], + [(-12, 5), (-31, 4), (-3, 1)], + ] + } else { + [ + [(0, 1), (0, 1), (5, 1)], + [(0, 1), (0, 1), (-5, 1)], + [(12, 5), (16, 5), (3, 1)], + [(-12, 5), (-16, 5), (-3, 1)], + ] + }; + assert_eq!(program.momentum_legs().len(), spatial.len()); + program + .momentum_legs() + .iter() + .copied() + .zip(spatial) + .map(|(label, coordinates)| { + let [x, y, z] = coordinates.map(|(n, d)| rational(prototype, n, d)); + ( + label, + FutureNullMomentum::from_spatial(ThreeMomentum::new(x, y, z)).unwrap(), + ) + }) + .collect() +} + +fn failure(result: RusticolResult) -> rusticol_core::RusticolError { + result.err().expect("request must fail explicitly") +} + +fn assert_close_values( + actual: &[Complex], + expected: &[Complex], + bits: u32, +) { + assert_eq!(actual.len(), expected.len()); + for (actual, expected) in actual.iter().zip(expected) { + let guard = if bits == 53 { 16 } else { 32 }; + let tolerance = expected.re.from_usize(2).pow(u64::from(bits - guard)).inv(); + let scale = expected.re.norm() + expected.re.one(); + assert!((actual.re.clone() - expected.re.clone()).norm() < tolerance.clone() * &scale); + assert!((actual.im.clone() - expected.im.clone()).norm() < tolerance * scale); + } +} + +fn assert_close_batch(actual: &Batch, expected: &Batch) { + actual.same_rows(expected).unwrap(); + assert_eq!(actual.columns(), expected.columns()); + let bits = actual.precision_bits(); + match (actual.data(), expected.data()) { + (Data::Double(a), Data::Double(b)) => assert_close_values(a, b, bits), + (Data::DoubleDouble(a), Data::DoubleDouble(b)) => assert_close_values(a, b, bits), + (Data::Arbitrary(a), Data::Arbitrary(b)) => assert_close_values(a, b, bits), + _ => panic!("carrier changed"), + } +} + +#[test] +fn external_null_admission_uses_resolved_exact_mass_including_scalars() { + let massive = scalar_program(scalar_definition(), BTreeMap::new(), &[]).unwrap(); + assert!(massive.external_states().is_empty()); + assert_eq!( + failure(massive.future_null_domain()).kind(), + RusticolErrorKind::UnsupportedRuntimeCapability, + ); + // The old raw scalar API still admits a massive spectator inventory. + let raw = Batch::typed(vec![Complex::new(0., 0.); 16], vec![3], 16, 53).unwrap(); + assert!( + massive + .evaluator("double", None) + .unwrap() + .evaluate_batch(&raw) + .is_ok() + ); + assert!(massless_scalar().future_null_domain().is_ok()); + + let mut structural_zero = scalar_definition(); + structural_zero["particles"][0]["mass"] = serde_json::json!("ZERO"); + assert!( + scalar_program(structural_zero.clone(), BTreeMap::new(), &[]) + .unwrap() + .future_null_domain() + .is_ok() + ); + structural_zero["parameters"][0]["value"] = serde_json::json!([3., 0.]); + assert_eq!( + failure( + scalar_program(structural_zero, BTreeMap::new(), &[]) + .unwrap() + .future_null_domain() + ) + .kind(), + RusticolErrorKind::UnsupportedRuntimeCapability + ); + + let mut runtime_zero = scalar_definition(); + runtime_zero["parameters"][2]["value"] = serde_json::json!([0., 0.]); + let runtime = [RuntimeTreeParameter { + symbol: symbol!("UFO::m"), + domain: TreeParameterDomain::Real, + }]; + let dependent = scalar_program(runtime_zero, BTreeMap::new(), &runtime).unwrap(); + assert_eq!( + failure(dependent.future_null_domain()).kind(), + RusticolErrorKind::UnsupportedRuntimeCapability + ); + + // Leave a runtime symbol free; only the exactly vanishing dependent mass + // is admissible, irrespective of that symbol's stored default. + let mut identity = scalar_definition(); + identity["parameters"][2]["nature"] = serde_json::json!("internal"); + identity["parameters"][2]["expression"] = serde_json::json!("g_-g_"); + let runtime = [RuntimeTreeParameter { + symbol: symbol!("UFO::g_"), + domain: TreeParameterDomain::Real, + }]; + assert!( + scalar_program(identity, BTreeMap::new(), &runtime) + .unwrap() + .future_null_domain() + .is_ok() + ); + + let unresolved = scalar_program( + scalar_definition(), + BTreeMap::from([( + symbol!("UFO::m"), + Atom::var(symbol!("null_fixture::unresolved_mass")), + )]), + &[], + ) + .unwrap(); + assert_eq!( + failure(unresolved.future_null_domain()).kind(), + RusticolErrorKind::UnsupportedRuntimeCapability + ); + // Actual quark, antiquark and vector inventories use the same factory. + assert!( + qcd(["u", "u~"], ["g", "g"], false) + .future_null_domain() + .is_ok() + ); +} + +#[test] +fn binding_canonicalizes_complete_labels_and_rejects_row_or_carrier_ambiguity() { + let program = massless_scalar(); + let domain = program.future_null_domain().unwrap(); + let canonical = keyed_row(&program, &0_f64, false); + let mut shuffled = canonical.clone(); + shuffled.rotate_left(1); + shuffled.reverse(); + let batch = domain + .bind(vec![41, 7], vec![shuffled, canonical.clone()], 53) + .unwrap(); + assert_eq!(batch.momenta().ids(), &[41, 7]); + assert_eq!(batch.momenta().columns(), 16); + let Data::Double(values) = batch.momenta().data() else { + panic!("carrier changed") + }; + for row in 0..2 { + for (column, (label, owner)) in canonical.iter().enumerate() { + assert_eq!(batch.leg(row, *label), Some(owner)); + for (component, expected) in owner.as_momentum().components().iter().enumerate() { + assert_eq!( + values[row * 16 + column * 4 + component], + Complex::new(**expected, 0.) + ); + } + } + } + assert!(batch.leg(2, canonical[0].0).is_none()); + assert!(batch.leg(0, usize::MAX).is_none()); + assert!( + domain + .bind( + vec![1], + Vec::)>>::new(), + 53 + ) + .is_err() + ); + assert!( + domain + .bind(vec![1, 1], vec![canonical.clone(); 2], 53) + .is_err() + ); + let mut missing = canonical.clone(); + missing.pop(); + let mut duplicate = canonical.clone(); + duplicate[1].0 = duplicate[0].0; + let mut extra = canonical.clone(); + extra.push((usize::MAX, canonical[0].1.clone())); + let mut foreign = canonical.clone(); + foreign[0].0 = usize::MAX; + for invalid in [missing, duplicate, extra, foreign] { + assert!(domain.bind(vec![41], vec![invalid], 53).is_err()); + } + assert!(domain.bind(vec![41], vec![canonical], 113).is_err()); + assert!(domain.bind::(vec![], vec![], 113).is_err()); + assert!(domain.bind::(vec![], vec![], 106).is_ok()); + assert!(domain.bind::(vec![], vec![], 113).is_err()); + assert!(domain.bind::(vec![], vec![], 1).is_err()); + assert!(domain.bind::(vec![], vec![], 233).is_ok()); +} + +#[test] +fn program_identity_is_checked_even_for_empty_scalar_batches() { + let first = massless_scalar(); + let second = massless_scalar(); + assert!(!Arc::ptr_eq(&first, &second)); + assert_eq!(first.momentum_legs(), second.momentum_legs()); + let first_domain = first.future_null_domain().unwrap(); + let cloned_domain = first_domain.clone(); + let second_domain = second.future_null_domain().unwrap(); + let mut evaluator = first_domain + .evaluator_with_precision_bits(Arithmetic::Double, 53) + .unwrap(); + for empty in [false, true] { + let rows = if empty { + vec![] + } else { + vec![keyed_row(&first, &0_f64, false)] + }; + let ids = if empty { vec![] } else { vec![19] }; + let accepted = cloned_domain.bind(ids.clone(), rows.clone(), 53).unwrap(); + let foreign = second_domain.bind(ids, rows, 53).unwrap(); + let result = evaluator.evaluate_batch(&accepted, None).unwrap(); + assert_eq!(result.rows(), usize::from(!empty)); + assert_eq!(result.columns(), 1); + assert!( + evaluator + .evaluate_coefficients(&accepted, None, &[]) + .is_ok() + ); + assert!(evaluator.evaluate_batch(&foreign, None).is_err()); + assert!( + evaluator + .evaluate_coefficients(&foreign, None, &[]) + .is_err() + ); + assert!( + evaluator + .evaluate_coefficients(&accepted, None, &[1]) + .is_err() + ); + if !empty { + let Data::Double(values) = result.data() else { + panic!("carrier changed") + }; + assert_eq!(values.as_slice(), &[Complex::new(1., 0.)]); + } + } +} + +#[test] +fn owned_state_source_checks_the_entire_prepared_four_column_tuple() { + let program = qcd(["u", "d"], ["u", "d"], false); + let domain = program.future_null_domain().unwrap(); + let source = domain + .bind(vec![71], vec![keyed_row(&program, &0_f64, false)], 53) + .unwrap(); + let empty = domain.bind::(vec![], vec![], 53).unwrap(); + let state = program.external_states()[0].clone(); + assert!(source.components(&state, state.helicities[0]).is_ok()); + assert!(empty.components(&state, state.helicities[0]).is_ok()); + for column in 0..4 { + let mut corrupt = state.clone(); + corrupt.momentum_columns[column] = (corrupt.momentum_columns[column] + 4) % 16; + assert!(source.components(&corrupt, state.helicities[0]).is_err()); + assert!(empty.components(&corrupt, state.helicities[0]).is_err()); + } + // Matching the first column alone must not identify a native occurrence. + let mut swapped = state.clone(); + swapped.momentum_columns.swap(1, 2); + assert!(source.components(&swapped, state.helicities[0]).is_err()); + let mut wrong_role = state.clone(); + wrong_role.kind = "epsilon".into(); + assert!(source.components(&wrong_role, state.helicities[0]).is_err()); + let mut wrong_helicities = state.clone(); + wrong_helicities.helicities = vec![0]; + assert!(source.components(&wrong_helicities, 0).is_err()); +} + +fn analytic_norm( + program: &Arc, + expected: &Rational, + prototype: T, + arithmetic: Arithmetic, + bits: u32, +) { + let domain = program.future_null_domain().unwrap(); + let mut rows = vec![ + keyed_row(program, &prototype, false), + keyed_row(program, &prototype, true), + ]; + rows[1].reverse(); + let source = domain.bind(vec![13, 907], rows, bits).unwrap(); + let mut evaluator = domain + .evaluator_with_precision_bits(arithmetic, bits) + .unwrap(); + let output = evaluator.evaluate_batch(&source, None).unwrap(); + let expected = Batch::typed( + vec![Complex::new(prototype.from_rational(expected), prototype.zero()); 2], + vec![13, 907], + 1, + bits, + ) + .unwrap(); + assert_close_batch(&output, &expected); + assert!(output.nonfinite_rows().is_empty()); + // The coefficient seam consumes components of the same retained owners. + // This is a path-association check; the analytic norm above is independent. + let helicities: Vec<_> = program + .external_states() + .iter() + .enumerate() + .map(|(i, state)| state.helicities[i % state.helicities.len()]) + .collect(); + let components: Vec<_> = program + .external_states() + .iter() + .zip(&helicities) + .map(|(state, &helicity)| source.components(state, helicity).unwrap()) + .collect(); + for ((state, &helicity), actual) in program + .external_states() + .iter() + .zip(&helicities) + .zip(&components) + { + let label = program.momentum_legs()[state.momentum_columns[0] / 4]; + let expected = Batch::typed( + (0..source.momenta().rows()) + .flat_map(|row| { + source + .leg(row, label) + .unwrap() + .wavefunction( + state.kind.parse().unwrap(), + Helicity::try_from(helicity as i8).unwrap(), + ) + .unwrap() + .into_components() + }) + .collect(), + source.momenta().ids().to_vec(), + 4, + bits, + ) + .unwrap(); + assert_close_batch(actual, &expected); + } + let owned = evaluator + .evaluate_coefficients(&source, None, &helicities) + .unwrap(); + let supplied = program + .evaluator_with_precision_bits(arithmetic, bits) + .unwrap() + .evaluate_external_components(source.momenta(), &components.iter().collect::>()) + .unwrap(); + assert_close_batch(&owned, &supplied); +} + +#[test] +fn owned_qcd_crossings_match_analytic_norms_with_fresh_runtime_width_sources() { + // PDG 2024 Cross-Section Formulae, Eqs. 51.5, 51.4, 51.10, 51.12: + // https://pdg.lbl.gov/2024/reviews/rpp2024-rev-cross-section-formulae.pdf + // Spin/color sums restore initial degeneracies; no symmetry or flux factor. + let s = Rational::from(100); + let t = Rational::from(-20); + let u = Rational::from(-80); + let cases = [ + ( + ["u", "d"], + ["u", "d"], + Rational::from(16) * (s.clone().pow(2) + u.clone().pow(2)) / t.clone().pow(2), + ), + ( + ["u", "u~"], + ["d", "d~"], + Rational::from(16) * (t.clone().pow(2) + u.clone().pow(2)) / s.clone().pow(2), + ), + ( + ["u", "g"], + ["u", "g"], + Rational::from(96) + * ((s.clone().pow(2) + u.clone().pow(2)) / t.clone().pow(2) + - Rational::from((4, 9)) * (s.clone().pow(2) + u.clone().pow(2)) + / (s.clone() * &u)), + ), + ( + ["g", "g"], + ["g", "g"], + Rational::from(1152) + * (Rational::from(3) + - u.clone() * &t / s.clone().pow(2) + - s.clone() * &u / t.clone().pow(2) + - s.clone() * &t / u.clone().pow(2)), + ), + ]; + for (incoming, outgoing, expected) in cases { + let program = qcd(incoming, outgoing, false); + let domain = program.future_null_domain().unwrap(); + let legs: Vec<_> = domain.external_legs().collect(); + assert_eq!( + legs.iter().map(|(label, _, _)| *label).collect::>(), + program.momentum_legs() + ); + assert_eq!(domain.process(), &Process::new(incoming, outgoing)); + for ((_, particle, role), (name, expected_role)) in legs.into_iter().zip( + incoming + .into_iter() + .map(|name| (name, ExternalState::Incoming)) + .chain( + outgoing + .into_iter() + .map(|name| (name, ExternalState::Outgoing)), + ), + ) { + assert_eq!(domain.model().particle_by_id(particle).unwrap().name, name); + assert_eq!(role, expected_role); + } + analytic_norm(&program, &expected, 0_f64, Arithmetic::Double, 53); + for bits in [113, 137, 233, 256, 512, 1024] { + analytic_norm( + &program, + &expected, + Float::with_val(bits, 0), + Arithmetic::Arbitrary, + bits, + ); + } + } +} + +fn delicate(prototype: &T, mixed: bool) -> FutureNullMomentum { + let spatial = if mixed { + let huge = prototype.from_usize(2).pow(600); + ThreeMomentum::new(huge.clone().inv(), prototype.zero(), huge) + } else { + let delta = prototype.from_usize(2).pow(30).inv(); + ThreeMomentum::new( + prototype.from_usize(2) * &delta, + prototype.zero(), + -prototype.one() + delta.pow(2), + ) + }; + FutureNullMomentum::from_spatial(spatial).unwrap() +} + +#[test] +fn typed_chart_failures_keep_external_label_row_and_attempt_bits() { + let program = qcd(["u", "d"], ["u", "d"], false); + let domain = program.future_null_domain().unwrap(); + let state = &program.external_states()[0]; + let label = program.momentum_legs()[state.momentum_columns[0] / 4]; + for mixed in [false, true] { + let mut low = keyed_row(&program, &0_f64, false); + low.iter_mut().find(|(id, _)| *id == label).unwrap().1 = delicate(&0_f64, mixed); + let source = domain + .bind( + vec![31, 907], + vec![keyed_row(&program, &0_f64, false), low], + 53, + ) + .unwrap(); + let error = failure(source.components(state, state.helicities[0])); + assert_eq!( + error.kind(), + RusticolErrorKind::InsufficientPrecision { bits: 53 } + ); + assert!( + error + .message() + .contains(&format!("lane 907 (row 1), occurrence {label}:")) + ); + for bits in [113, 137, 233, 256, 512, 1024] { + let prototype = Float::with_val(bits, 0); + let mut row = keyed_row(&program, &prototype, false); + row.iter_mut().find(|(id, _)| *id == label).unwrap().1 = delicate(&prototype, mixed); + let source = domain.bind(vec![907], vec![row], bits).unwrap(); + let output = source.components(state, state.helicities[0]).unwrap(); + assert_eq!(output.ids(), &[907]); + assert_eq!(output.precision_bits(), bits); + assert!(output.nonfinite_rows().is_empty()); + } + } + // Metadata records the attempt, not the lesser tracked precision of an + // intermediate. Do not pad its precision or silently promote the source. + let tracked = Float::with_val(53, 0); + let mut row = keyed_row(&program, &tracked, false); + row.iter_mut().find(|(id, _)| *id == label).unwrap().1 = delicate(&tracked, false); + let source = domain.bind(vec![1021], vec![row], 137).unwrap(); + let error = failure(source.components(state, state.helicities[0])); + assert_eq!( + error.kind(), + RusticolErrorKind::InsufficientPrecision { bits: 137 } + ); + assert!( + error + .message() + .contains(&format!("lane 1021 (row 0), occurrence {label}:")) + ); +} + +#[test] +fn runtime_parameters_preserve_typed_rows_and_sub_binary64_couplings() { + let program = qcd(["u", "d"], ["u", "d"], true); + let domain = program.future_null_domain().unwrap(); + let bits = 137; + let p = Float::with_val(bits, 0); + let source = domain + .bind(vec![29, 11], vec![keyed_row(&program, &p, false); 2], bits) + .unwrap(); + let first = p.one(); + let second = p.one() + Float::with_val(bits, 2).pow(90).inv(); + let parameters = Batch::typed( + vec![ + Complex::new(first.clone(), p.zero()), + Complex::new(second.clone(), p.zero()), + ], + vec![29, 11], + 1, + bits, + ) + .unwrap(); + let mut evaluator = domain + .evaluator_with_precision_bits(Arithmetic::Arbitrary, bits) + .unwrap(); + assert!(evaluator.evaluate_batch(&source, None).is_err()); + let result = evaluator + .evaluate_batch(&source, Some(¶meters)) + .unwrap(); + let expected = Batch::typed( + [first, second] + .into_iter() + .map(|g| Complex::new(Float::with_val(bits, 656) * g.pow(4), p.zero())) + .collect(), + vec![29, 11], + 1, + bits, + ) + .unwrap(); + assert_close_batch(&result, &expected); + let Data::Arbitrary(values) = result.data() else { + panic!("carrier changed") + }; + assert!(values[1].re > values[0].re); + assert!( + evaluator + .evaluate_batch( + &source, + Some(¶meters.take_rows(vec![1, 0], None).unwrap()) + ) + .is_err() + ); + let complex = Batch::typed( + vec![Complex::new(p.one(), p.one()); 2], + vec![29, 11], + 1, + bits, + ) + .unwrap(); + assert!(evaluator.evaluate_batch(&source, Some(&complex)).is_err()); + let mismatch = domain + .evaluator_with_precision_bits(Arithmetic::Arbitrary, 233) + .unwrap() + .evaluate_batch(&source, Some(¶meters)); + assert!(mismatch.is_err()); +} + +fn color_requests(program: &TreeProgram) -> Vec { + let legs = program.color_tensors()[0].legs(); + let [a, b, c] = std::array::from_fn(|i| legs[i].label); + let mut requests = vec![TreeColorInsertion::Identity]; + // Two anchors include fundamental and adjoint Casimirs in the QCD fixture. + for anchor in &legs[..2] { + requests.extend(legs.iter().map(|leg| TreeColorInsertion::Pair { + left: anchor.label, + right: leg.label, + })); + } + requests.extend([ + TreeColorInsertion::OrderedProduct(vec![[a, b], [b, c]]), + TreeColorInsertion::OrderedProduct(vec![[b, c], [a, b]]), + TreeColorInsertion::Anticommutator { + first: [a, b], + second: [b, c], + }, + TreeColorInsertion::StructureConstantTripole { + first: a, + second: b, + third: c, + }, + ]); + requests +} + +fn check_color_conservation( + values: &[Complex], + columns: usize, + legs: &[rusticol_core::color::ColorLeg], + bits: u32, +) { + for row in values.chunks_exact(columns) { + let identity = &row[0]; + assert!(identity.re > identity.re.zero()); + for (anchor_index, anchor) in legs[..2].iter().enumerate() { + let casimir = match anchor.representation { + ColorRepresentation::Adjoint => Rational::from(3), + _ => Rational::from((4, 3)), + }; + let expected = + identity * &Complex::new(identity.re.from_rational(&casimir), identity.re.zero()); + let start = 1 + anchor_index * legs.len(); + assert_close_values( + &row[start + anchor_index..start + anchor_index + 1], + std::slice::from_ref(&expected), + bits, + ); + let sum = row[start..start + legs.len()] + .iter() + .enumerate() + .filter(|(index, _)| *index != anchor_index) + .fold( + Complex::new(identity.re.zero(), identity.re.zero()), + |sum, (_, value)| sum + value, + ); + // Sum of all distinct spectators = -C_i times the identity. + assert_close_values(&[sum], &[-expected], bits); + } + } +} + +fn owned_color_paths( + program: &Arc, + color: &TreeColorProgram, + prototype: T, + arithmetic: Arithmetic, + bits: u32, +) { + let domain = program.future_null_domain().unwrap(); + let exact_row = program + .momentum_legs() + .iter() + .copied() + .zip([[0, 0, 5], [0, 0, -5], [4, 0, 3], [-4, 0, -3]]) + .map(|(label, [x, y, z])| { + ( + label, + FutureNullMomentum::from_spatial(ThreeMomentum::new( + prototype.from_i64(x), + prototype.from_i64(y), + prototype.from_i64(z), + )) + .unwrap(), + ) + }) + .collect(); + let exact = domain.bind(vec![81], vec![exact_row], bits).unwrap(); + let helicities: Vec<_> = program + .external_states() + .iter() + .map(|state| state.helicities[0]) + .collect(); + let mut evaluator = color + .evaluator_with_precision_bits(arithmetic, bits) + .unwrap(); + let owned_fixed = evaluator + .evaluate_future_null_fixed_helicities(&exact, None, &helicities) + .unwrap(); + let raw_fixed = evaluator + .evaluate_fixed_helicities(exact.momenta(), &helicities) + .unwrap(); + assert_close_batch(&owned_fixed, &raw_fixed); + let owned_sum = evaluator.evaluate_future_null_batch(&exact, None).unwrap(); + let raw_sum = evaluator.evaluate_batch(exact.momenta()).unwrap(); + assert_close_batch(&owned_sum, &raw_sum); + + let nonaxis = domain + .bind( + vec![13, 907], + vec![ + keyed_row(program, &prototype, false), + keyed_row(program, &prototype, true), + ], + bits, + ) + .unwrap(); + for output in [ + evaluator + .evaluate_future_null_fixed_helicities(&nonaxis, None, &helicities) + .unwrap(), + evaluator + .evaluate_future_null_batch(&nonaxis, None) + .unwrap(), + ] { + let legs = program.color_tensors()[0].legs(); + match output.data() { + Data::Double(v) => check_color_conservation(v, output.columns(), legs, bits), + Data::DoubleDouble(v) => check_color_conservation(v, output.columns(), legs, bits), + Data::Arbitrary(v) => check_color_conservation(v, output.columns(), legs, bits), + } + } + + // Existing directed supplied-component probe, not a new admitted physical + // spin/source API: phase one ket block by i and retain i times every output. + // The unphased blocks and momenta come from these same bound owners. + let components: Vec<_> = program + .external_states() + .iter() + .zip(&helicities) + .map(|(state, &helicity)| exact.components(state, helicity).unwrap()) + .collect(); + let phase = |batch: &Batch| -> Batch { + fn rotate(values: &[Complex], batch: &Batch) -> Batch { + Batch::typed( + values + .iter() + .map(|z| Complex::new(-z.im.clone(), z.re.clone())) + .collect(), + batch.ids().to_vec(), + batch.columns(), + batch.precision_bits(), + ) + .unwrap() + } + match batch.data() { + Data::Double(v) => rotate(v, batch), + Data::DoubleDouble(v) => rotate(v, batch), + Data::Arbitrary(v) => rotate(v, batch), + } + }; + let mut ket = components.clone(); + ket[0] = phase(&ket[0]); + let directed = evaluator + .evaluate_external_components( + exact.momenta(), + &components.iter().collect::>(), + &ket.iter().collect::>(), + ) + .unwrap(); + assert_close_batch(&directed, &phase(&owned_fixed)); + let ordered_column = 1 + 2 * program.color_tensors()[0].legs().len(); + match directed.data() { + Data::Double(v) => assert_ne!(v[ordered_column].im, 0.), + Data::DoubleDouble(v) => assert_ne!(v[ordered_column].im, v[ordered_column].im.zero()), + Data::Arbitrary(v) => assert_ne!(v[ordered_column].im, v[ordered_column].im.zero()), + } +} + +#[test] +fn owned_color_insertions_reuse_raw_contractions_and_obey_independent_color_identities() { + let tree = qcd(["u", "g"], ["u", "g"], false); + let color = TreeColorProgram::prepare(tree.clone(), color_requests(&tree)).unwrap(); + // Use the existing exact insertion matrix. Verify the chosen ordered probe + // really is non-Hermitian; no duplicate operator or color-network oracle. + let ordered_column = 1 + 2 * tree.color_tensors()[0].legs().len(); + let matrix = color.matrix_nc_terms(ordered_column).unwrap(); + assert!((0..matrix.len()).any(|i| (0..matrix.len()).any(|j| { + evaluate_nc_terms(&matrix[i][j], 3).unwrap() + != evaluate_nc_terms(&matrix[j][i], 3).unwrap().conj() + }))); + owned_color_paths(&tree, &color, 0_f64, Arithmetic::Double, 53); + for bits in [113, 137, 233, 256, 512, 1024] { + owned_color_paths( + &tree, + &color, + Float::with_val(bits, 0), + Arithmetic::Arbitrary, + bits, + ); + } +} + +#[test] +fn owned_color_empty_batches_still_validate_program_parameters_precision_and_helicities() { + let tree = qcd(["u", "d"], ["u", "d"], true); + let other = qcd(["u", "d"], ["u", "d"], true); + let domain = tree.future_null_domain().unwrap(); + let foreign_domain = other.future_null_domain().unwrap(); + let color = + TreeColorProgram::prepare(tree.clone(), vec![TreeColorInsertion::Identity]).unwrap(); + let bits = 137; + let empty = domain.bind::(vec![], vec![], bits).unwrap(); + let foreign = foreign_domain.bind::(vec![], vec![], bits).unwrap(); + let parameters = Batch::typed::(vec![], vec![], 1, bits).unwrap(); + let helicities: Vec<_> = tree + .external_states() + .iter() + .map(|state| state.helicities[0]) + .collect(); + let mut evaluator = color + .evaluator_with_precision_bits(Arithmetic::Arbitrary, bits) + .unwrap(); + let valid = evaluator + .evaluate_future_null_batch(&empty, Some(¶meters)) + .unwrap(); + assert_eq!( + (valid.rows(), valid.columns(), valid.precision_bits()), + (0, 1, bits) + ); + assert!( + evaluator + .evaluate_future_null_fixed_helicities(&empty, Some(¶meters), &helicities) + .is_ok() + ); + assert!( + evaluator + .evaluate_future_null_batch(&foreign, Some(¶meters)) + .is_err() + ); + assert!( + evaluator + .evaluate_future_null_fixed_helicities(&foreign, Some(¶meters), &helicities) + .is_err() + ); + let malformed = [ + None, + Some(Batch::typed::(vec![], vec![], 0, bits).unwrap()), + Some(Batch::typed::(vec![], vec![], 2, bits).unwrap()), + Some(Batch::typed::(vec![], vec![], 1, 233).unwrap()), + // Batch is the raw parameter boundary: an empty wrong carrier with + // matching bit metadata must still be rejected before evaluation. + Some(Batch::typed::(vec![], vec![], 1, bits).unwrap()), + Some( + Batch::typed( + vec![Complex::new( + Float::with_val(bits, 1), + Float::with_val(bits, 0), + )], + vec![7], + 1, + bits, + ) + .unwrap(), + ), + ]; + for parameters in malformed { + assert!( + evaluator + .evaluate_future_null_batch(&empty, parameters.as_ref()) + .is_err() + ); + assert!( + evaluator + .evaluate_future_null_fixed_helicities(&empty, parameters.as_ref(), &helicities) + .is_err() + ); + } + let wrong_source = domain.bind::(vec![], vec![], 233).unwrap(); + assert!( + evaluator + .evaluate_future_null_batch(&wrong_source, Some(¶meters)) + .is_err() + ); + assert!( + evaluator + .evaluate_future_null_fixed_helicities(&wrong_source, Some(¶meters), &helicities) + .is_err() + ); + let wrong_carrier = domain.bind::(vec![], vec![], 53).unwrap(); + assert!( + evaluator + .evaluate_future_null_batch(&wrong_carrier, Some(¶meters)) + .is_err() + ); + assert!( + evaluator + .evaluate_future_null_fixed_helicities(&wrong_carrier, Some(¶meters), &helicities) + .is_err() + ); + let float53 = domain.bind::(vec![], vec![], 53).unwrap(); + let parameters53 = Batch::typed::(vec![], vec![], 1, 53).unwrap(); + let mut double = color + .evaluator_with_precision_bits(Arithmetic::Double, 53) + .unwrap(); + assert!( + double + .evaluate_future_null_batch(&float53, Some(¶meters53)) + .is_err() + ); + assert!( + double + .evaluate_future_null_fixed_helicities(&float53, Some(¶meters53), &helicities) + .is_err() + ); + assert!( + evaluator + .evaluate_future_null_fixed_helicities(&empty, Some(¶meters), &[]) + .is_err() + ); + let mut wrong_helicity = helicities; + wrong_helicity[0] = 0; + assert!( + evaluator + .evaluate_future_null_fixed_helicities(&empty, Some(¶meters), &wrong_helicity) + .is_err() + ); +} diff --git a/rust/loop-reduction/src/tree_preparation.rs b/rust/loop-reduction/src/tree_preparation.rs index aa91809a7..d39c40eaa 100644 --- a/rust/loop-reduction/src/tree_preparation.rs +++ b/rust/loop-reduction/src/tree_preparation.rs @@ -19,7 +19,7 @@ use feynkit_generator::{ FilterScope, GenerationFilter, GenerationOptions, NumeratorGrouping, Process, }; use feynkit_graph::ExternalState; -use feynkit_model::{Model, ParameterCard, ParameterType, Particle}; +use feynkit_model::{Model, ParameterCard, ParameterType, Particle, ParticleId}; use idenso::{representations::Bispinor, tensor::SymbolicTensor}; use rusticol_core::{ RusticolError, RusticolResult, @@ -53,6 +53,7 @@ use crate::{ /// Flux, initial spin/color averages and identical-particle phase-space factors /// belong to the consumer. pub struct TreeProgram { + pub(crate) external_domain: PreparedExternalDomain, field_class: TreeFieldClass, momentum_legs: Vec, diagram_count: usize, @@ -64,6 +65,23 @@ pub struct TreeProgram { color_tensors: Vec, } +// Created only while preparing the model-derived instructions. The shared +// program owns the model namespace and complete external facts, including scalar +// legs; arbitrary saved-instruction evaluators cannot manufacture this domain. +pub(crate) struct PreparedExternalDomain { + pub(crate) model: Arc, + pub(crate) process: Process, + pub(crate) legs: Vec, +} + +pub(crate) struct PreparedExternalLeg { + pub(crate) label: usize, + pub(crate) particle: ParticleId, + pub(crate) direction: ExternalState, + pub(crate) resolved_mass: Atom, + pub(crate) state_index: Option, +} + /// Singlet-scalar entry point delegating to the common amplitude preparation. /// Its existing evaluation methods remain available through [`TreeProgram`]. pub struct ScalarTreeProgram(TreeProgram); @@ -374,6 +392,7 @@ impl TreeProgram { .symmetrize_left_right(false) .symmetrize_external_fermions(false) .numerator_grouping(NumeratorGrouping::None); + let retained_model = Arc::clone(&model); let generated = process.generate_diagrams(model, &options).map_err(error)?; for diagram in &generated.diagrams { for (_, _, edge) in diagram.edges() { @@ -441,6 +460,7 @@ impl TreeProgram { } } let momentum_legs: Vec<_> = amplitude.legs().iter().map(|leg| leg.index).collect(); + let mut prepared_external_legs = Vec::with_capacity(momentum_legs.len()); let source = physical_expression .with_checked_expression(bind_scalars( physical_expression.expression(), @@ -479,6 +499,19 @@ impl TreeProgram { .particle_by_id(leg.particle) .map_err(error)?; let incoming = leg.state == ExternalState::Incoming; + prepared_external_legs.push(PreparedExternalLeg { + label: leg.index, + particle: leg.particle, + direction: leg.state, + // Keep runtime placeholders free: a numerical zero default is + // not a masslessness proof for a parameter that may vary. + resolved_mass: bind_scalars( + &particle.symbolic_mass(amplitude.model()), + &scalar_bindings, + ) + .expand(), + state_index: (particle.spin != 1).then_some(states.len()), + }); let representation = match particle.color { 1 => ColorRepresentation::Singlet, 3 => ColorRepresentation::Fundamental, @@ -606,6 +639,11 @@ impl TreeProgram { .build() .map_err(error)?; Ok(Self { + external_domain: PreparedExternalDomain { + model: retained_model, + process: process.clone(), + legs: prepared_external_legs, + }, field_class, momentum_legs, diagram_count: amplitude.diagrams().len(), @@ -618,10 +656,24 @@ impl TreeProgram { }) } + /// Bind constructive future-null states to this exact immutable program. + /// Every external mass, including scalar masses, must resolve identically + /// to zero with runtime inputs left symbolic. This does not certify event + /// conservation, original exact-source precision, or gauge completeness. + pub fn future_null_domain( + self: &Arc, + ) -> RusticolResult { + crate::tree_null_states::FutureNullTreeDomain::new(Arc::clone(self)) + } + pub fn field_class(&self) -> TreeFieldClass { self.field_class } + pub(crate) fn external_states(&self) -> &[ExternalTreeState] { + &self.states + } + pub fn momentum_legs(&self) -> &[usize] { &self.momentum_legs } diff --git a/rust/loop-reduction/src/tree_runtime.rs b/rust/loop-reduction/src/tree_runtime.rs index ebe224835..d9a170478 100644 --- a/rust/loop-reduction/src/tree_runtime.rs +++ b/rust/loop-reduction/src/tree_runtime.rs @@ -26,6 +26,36 @@ pub struct ExternalTreeState { pub helicities: Vec, } +// Only crate-owned sources enter the shared evaluator. Public raw entry points +// keep their existing semantics; program-bound sources are validated by their +// immutable wrapper before reaching this seam. +pub(crate) trait TreeStateSource { + fn momenta(&self) -> &Batch; + fn components(&self, state: &ExternalTreeState, helicity: i32) -> RusticolResult; +} + +pub(crate) struct RawTreeStateSource<'a>(pub(crate) &'a Batch); + +impl TreeStateSource for RawTreeStateSource<'_> { + fn momenta(&self) -> &Batch { + self.0 + } + + fn components(&self, state: &ExternalTreeState, helicity: i32) -> RusticolResult { + let p = self + .0 + .take_columns(state.momentum_columns.to_vec()) + .map_err(error)?; + wavefunctions_batch( + &p, + vec![helicity; self.0.rows()], + vec![state.kind.clone(); self.0.rows()], + Some(4), + ) + .map_err(RusticolError::from) + } +} + /// Scalar domains assumed while preparing a runtime model parameter. /// Runtime batches must obey these assumptions before numerical evaluation. #[derive(Clone, Copy, Debug, PartialEq, Eq)] @@ -405,6 +435,16 @@ impl TreeEvaluator { parameters: Option<&Batch>, helicities: &[i32], ) -> RusticolResult { + self.evaluate_coefficients_from_source(&RawTreeStateSource(momenta), parameters, helicities) + } + + pub(crate) fn evaluate_coefficients_from_source( + &mut self, + source: &dyn TreeStateSource, + parameters: Option<&Batch>, + helicities: &[i32], + ) -> RusticolResult { + let momenta = source.momenta(); self.validate_momenta(momenta)?; self.validate_parameters(momenta, parameters)?; if helicities.len() != self.states.len() @@ -425,18 +465,7 @@ impl TreeEvaluator { .states .iter() .zip(helicities) - .map(|(state, &helicity)| { - let p = momenta - .take_columns(state.momentum_columns.to_vec()) - .map_err(error)?; - wavefunctions_batch( - &p, - vec![helicity; momenta.rows()], - vec![state.kind.clone(); momenta.rows()], - Some(4), - ) - .map_err(RusticolError::from) - }) + .map(|(state, &helicity)| source.components(state, helicity)) .collect::>>()?; for component in &components { finite_result(component)?; @@ -531,6 +560,15 @@ impl TreeEvaluator { momenta: &Batch, parameters: Option<&Batch>, ) -> RusticolResult { + self.evaluate_batch_from_source(&RawTreeStateSource(momenta), parameters) + } + + pub(crate) fn evaluate_batch_from_source( + &mut self, + source: &dyn TreeStateSource, + parameters: Option<&Batch>, + ) -> RusticolResult { + let momenta = source.momenta(); self.validate_momenta(momenta)?; self.validate_parameters(momenta, parameters)?; if momenta.rows() == 0 { @@ -546,7 +584,7 @@ impl TreeEvaluator { .map(|(state, &position)| state.helicities[position]) .collect::>(); let coefficients = - self.evaluate_coefficients_inner(momenta, parameters, &helicities)?; + self.evaluate_coefficients_from_source(source, parameters, &helicities)?; let norm_inputs = if self.cartesian_norm { cartesian_coefficients(&coefficients)? } else {