gammalooprs/utils/

mod.rs

use crate::integrands::process::{
    GenericEvaluatorFloat, LmbMultiChannelingSetup, SamplingChannelBridge,
    sampling::maps::SamplingEvaluationError,
};
use crate::model::Model;
use crate::momentum::sample::{
    ExternalFourMomenta, ExternalIndex, ExternalThreeMomenta, LoopMomenta, SubspaceData,
};
use crate::momentum::signature::{ExternalSignature, LoopSignature};
use crate::momentum::{FourMomentum, ThreeMomentum};
use crate::numerator::aind::Aind;
use crate::numerator::ufo::UFO;
use crate::settings::runtime::ParameterizationSettings;
use crate::settings::runtime::SamplingSettings;
use crate::settings::runtime::kinematic::Externals;
use crate::settings::runtime::{ParameterizationMapping, ParameterizationMode};
use crate::utils::hyperdual_utils::new_constant;

use bincode::{Decode, Encode};
use colored::Colorize;
use idenso::representations::initialize;
use itertools::Itertools;
use linnet::half_edge::involution::EdgeIndex;

use rand::Rng;
use ref_ops::{RefAdd, RefDiv, RefMul, RefNeg, RefRem, RefSub};
use rug::float::{Constant, ParseFloatError, Round};
use rug::ops::{CompleteRound, Pow};
use rug::{Assign, Float};
use schemars::JsonSchema;
use serde::{Deserialize, Deserializer, Serialize};
use smallvec::SmallVec;
use spenso::algebra::algebraic_traits::RefOne;
use spenso::algebra::algebraic_traits::RefZero;
use spenso::algebra::complex::Complex;
use spenso::algebra::complex::R;
use spenso::algebra::complex::SymbolicaComplex;
use spenso::algebra::complex::symbolica_traits::ToFloat;
use spenso::algebra::upgrading_arithmetic::TrySmallestUpgrade;
use spenso::network::library::TensorLibraryData;
use spenso::network::library::function_lib::{INBUILTS, Panic, SymbolLib, Wrap};
use spenso::network::library::symbolic::{ExplicitKey, TensorLibrary};
use spenso::network::parsing::ShadowedStructure;
use spenso::structure::concrete_index::ExpandedIndex;
use spenso::tensors::data::StorageTensor;
use spenso::tensors::parametric::ParamTensor;
use spenso::tensors::parametric::to_param::ToAtom;
use spenso_hep_lib::hep_lib_atom;
use symbolica::{
    domains::{
        dual::HyperDual,
        float::{FixedPrecision, Float as SymbolicaFloat, FloatLike as SymFloatLike},
    },
    prelude::*,
};

use std::cell::RefCell;
use std::cmp::{Ord, Ordering};
use std::collections::HashMap;
use std::fmt::{Debug, Display};
use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, Sub, SubAssign};
use std::str::FromStr;
use std::sync::{LazyLock, OnceLock, RwLock};
use std::time::Duration;

use vakint::Vakint;
// use symbolica_community::physics::tensors::library::{
//     gamma5_weyl_data, gamma_data_weyl, proj_m_data_weyl, proj_p_data_weyl, SpensorLibrary,
//     TensorNamespace,
// };
// use symbolica_community::physics::tensors::structure::SpensoStucture;
// use symbolica::domains::Field;
use crate::MAX_LOOP;
use ::tracing::debug;
use typed_index_collections::TiVec;

/// Transient runtime data which contributes no bytes to persisted state.
/// Decoding always starts empty, so numerical bindings are rebuilt at warmup.
#[derive(Clone, Debug, PartialEq)]
pub struct RuntimeCache<T>(Option<T>);

impl<T> Default for RuntimeCache<T> {
    fn default() -> Self {
        Self(None)
    }
}

impl<T> RuntimeCache<T> {
    pub(crate) fn invalidate(&mut self) {
        self.0 = None;
    }

    pub(crate) fn set(&mut self, value: T) {
        self.0 = Some(value);
    }

    pub(crate) fn take(&mut self) -> Option<T> {
        self.0.take()
    }

    pub(crate) fn as_ref(&self) -> Option<&T> {
        self.0.as_ref()
    }

    pub(crate) fn as_mut(&mut self) -> Option<&mut T> {
        self.0.as_mut()
    }
}

impl<T> bincode::Encode for RuntimeCache<T> {
    fn encode<E: bincode::enc::Encoder>(
        &self,
        _encoder: &mut E,
    ) -> std::result::Result<(), bincode::error::EncodeError> {
        Ok(())
    }
}

impl<C, T> bincode::Decode<C> for RuntimeCache<T> {
    fn decode<D: bincode::de::Decoder<Context = C>>(
        _decoder: &mut D,
    ) -> std::result::Result<Self, bincode::error::DecodeError> {
        Ok(Self::default())
    }
}

pub const GIT_VERSION: &str = env!("CARGO_PKG_VERSION");
pub const VERSION: &str = "0.0.1";
pub mod fitting;
pub mod hyperdual_utils;
pub mod progress;
pub mod representations;
pub mod serde_utils;
pub use progress::{long_running_progress_style, long_running_progress_style_with_eta_warmup};
/// can be used instead of commenting out code
#[macro_export]
macro_rules! disable {
    ($($tokens:tt)*) => {};
}

/// Shorthand for tagged debug events.
///
/// `#tag` expands to `tag = true` and the remaining tokens are forwarded to
/// `tracing::debug!` unchanged.
///
/// GammaLoop's logging DSL treats tags as boolean event fields.
///
/// Emitting `#generation` makes directives such as `[{generation}]` eligible to
/// match this event. A negated directive term such as `[{!inspect}]` matches
/// events where the field is absent. If you need value-sensitive filtering, use
/// an explicit field such as `inspect = false` or `mode = "summary"` at the
/// callsite and a directive such as `[{inspect=false}]` or `[{mode=summary}]`.
///
/// Separately, the formatter supports sink-only field prefixes:
/// - `file.<name>` is rendered only in the logfile/json sink
/// - `display.<name>` is rendered only in the stderr/display sink
///
/// For example, `file.integrands = %expr` is useful for attaching large dumps
/// to the logfile without printing them to stderr.
///
/// Pipeline-wide tag set:
///
/// Prefer a small stable vocabulary that cuts across the whole pipeline.
///
/// Primary phase tags:
///
/// | Tag | Meaning |
/// | --- | --- |
/// | `#generation` | Work that builds or prepares runtime objects before Monte Carlo integration starts |
/// | `#integration` | Work performed while evaluating samples or managing the adaptive integration loop |
/// | `#profile` | UV/IR/profile-style diagnostic scans and their analysis |
/// | `#persistence` | Reading or writing saved state, manifests, checkpoints, and exported results |
///
/// Core domain tags:
///
/// | Tag | Meaning |
/// | --- | --- |
/// | `#uv` | Ultraviolet counterterms, UV forests, UV profiles, or UV-specific generation/evaluation logic |
/// | `#ir` | Infrared subtraction, IR profiles, or IR-specific generation/evaluation logic |
/// | `#subtraction` | Threshold, UV, or IR subtraction logic when the main concern is subtraction rather than the specific regime |
/// | `#sampling` | Channel choice, discrete axes, parameterizations, or sample-generation choices |
/// | `#stability` | Precision escalation, retries, instability diagnosis, and numerical safety checks |
/// | `#observables` | Histogramming, event-to-observable projection, and observable snapshot production |
/// | `#selectors` | Event-selection logic and selector decisions |
/// | `#cache` | Cache lookup, reuse, invalidation, or cache-debug instrumentation |
///
/// Common work-unit tags:
///
/// | Tag | Meaning |
/// | --- | --- |
/// | `#graph` | The log is about one graph or graph-local data |
/// | `#group` | The log is about a graph group or another explicitly grouped aggregate |
/// | `#orientation` | The log is about orientation-dependent data or choosing/summing orientations |
/// | `#channel` | The log is about multi-channeling or a particular channel |
/// | `#cut` | The log is about a cut, raised cut, or cut-local computation |
/// | `#event` | The log is about generated/retained event objects or event-group processing |
/// | `#sample` | The log is about one evaluation sample, its coordinates, or per-sample intermediate values |
/// | `#iteration` | The log is about one adaptive integration iteration or iteration-level summaries |
/// | `#term` | The log is about one algebraic term, summand, or term-local contribution |
///
/// Common purpose tags:
///
/// | Tag | Meaning |
/// | --- | --- |
/// | `#solver` | The log is about root-finding, linear solves, fitting, or similar numerical solver state |
/// | `#compile` | The log is about evaluator/code generation or compilation-oriented preparation |
/// | `#inspect` | The log exposes detailed intermediate state for debugging, rather than a high-level milestone |
/// | `#summary` | The log is a compact roll-up rather than a step-by-step trace |
/// | `#dump` | The log emits large or structured payloads such as expressions, tables, or serialized views |
///
/// Tag boundaries:
///
/// - Prefer `#generation` over `#compile` when the message is a broad generation milestone.
/// - Add `#compile` only when the message is specifically about evaluator construction or compilation-like work.
/// - Prefer `#uv` or `#ir` when the regime matters to the user; use `#subtraction` when the subtraction mechanism is the real concern.
/// - Prefer `#inspect` for verbose intermediate values that are mainly useful while debugging internals.
/// - Add `#dump` when the payload is large enough that users may want to suppress it separately from lighter debug logs.
/// - Do not use `#graph`, `#cut`, or `#sample` as substitutes for `#generation` or `#integration`; they refine phase tags rather than replace them.
///
/// Naming guidance:
///
/// - Use phase tags first. They answer "when in the pipeline did this happen?"
/// - Use domain tags second. They answer "what subsystem is this about?"
/// - Use work-unit tags third. They answer "what object is being worked on?"
/// - Use purpose tags last for optional refinement.
/// - Prefer tags that are meaningful across multiple modules.
/// - Avoid tags that merely restate a function name.
/// - Avoid tags tied to one internal representation unless that representation is a stable concept in user-facing debugging.
///
/// `parametric` is usually not a core pipeline tag. It is acceptable as a local refinement when needed, but should not be treated as part of the primary vocabulary unless it becomes a consistently useful cross-cutting concept.
///
/// Examples:
/// - `debug_tags!(#generation, #uv, #graph, #dump; name = %graph.name, "generated graph payload");`
/// - `debug_tags!(#integration, #summary; inspect = false, "iteration summary");`
#[macro_export]
macro_rules! debug_tags {
    (@collect_tags [$($acc:tt)*] # $tag:ident, $($tail:tt)*) => {
        $crate::debug_tags!(@collect_tags [$($acc)* $tag = true,] $($tail)*)
    };
    (@collect_tags [$($acc:tt)*] # $tag:ident; $($rest:tt)*) => {
        $crate::debug_tags!(@collect_fields [$($acc)* $tag = true,] [] $($rest)*)
    };
    (@collect_fields [$($tags:tt)*] [$($fields:tt)*] log.$name:ident = $value:expr, $($tail:tt)*) => {
        $crate::debug_tags!(@collect_fields
            [$($tags)*]
            [$($fields)*
                display.$name = %$crate::LogMessage::log_display(&$value),
                file.$name = %$crate::LogMessage::log_file(&$value),
            ]
            $($tail)*
        )
    };
    (@collect_fields [$($tags:tt)*] [$($fields:tt)*] $prefix:ident.$name:ident = %$value:expr, $($tail:tt)*) => {
        $crate::debug_tags!(@collect_fields [$($tags)*] [$($fields)* $prefix.$name = %$value,] $($tail)*)
    };
    (@collect_fields [$($tags:tt)*] [$($fields:tt)*] $prefix:ident.$name:ident = ?$value:expr, $($tail:tt)*) => {
        $crate::debug_tags!(@collect_fields [$($tags)*] [$($fields)* $prefix.$name = ?$value,] $($tail)*)
    };
    (@collect_fields [$($tags:tt)*] [$($fields:tt)*] $prefix:ident.$name:ident = $value:expr, $($tail:tt)*) => {
        $crate::debug_tags!(@collect_fields [$($tags)*] [$($fields)* $prefix.$name = $value,] $($tail)*)
    };
    (@collect_fields [$($tags:tt)*] [$($fields:tt)*] $name:ident = %$value:expr, $($tail:tt)*) => {
        $crate::debug_tags!(@collect_fields [$($tags)*] [$($fields)* $name = %$value,] $($tail)*)
    };
    (@collect_fields [$($tags:tt)*] [$($fields:tt)*] $name:ident = ?$value:expr, $($tail:tt)*) => {
        $crate::debug_tags!(@collect_fields [$($tags)*] [$($fields)* $name = ?$value,] $($tail)*)
    };
    (@collect_fields [$($tags:tt)*] [$($fields:tt)*] $name:ident = $value:expr, $($tail:tt)*) => {
        $crate::debug_tags!(@collect_fields [$($tags)*] [$($fields)* $name = $value,] $($tail)*)
    };
    (@collect_fields [$($tags:tt)*] [$($fields:tt)*] $name:ident, $($tail:tt)*) => {
        $crate::debug_tags!(@collect_fields [$($tags)*] [$($fields)* $name,] $($tail)*)
    };
    (@collect_fields [$($tags:tt)*] [$($fields:tt)*] $($rest:tt)+) => {
        tracing::debug!($($tags)* $($fields)* $($rest)+)
    };
    ($($input:tt)*) => {
        $crate::debug_tags!(@collect_tags [] $($input)*)
    };
}

#[allow(unused)]
const MAX_DIMENSION: usize = MAX_LOOP * 3;

pub(crate) const ESURFACE_SHIFT_THRESHOLD: f64 = 1.0e-13;
/// Default dimensionless tolerance for E-surface invariant-mass-squared margins.
pub const DEFAULT_ESURFACE_EXISTENCE_THRESHOLD: f64 = 1.0e-7;

pub const LEFT: usize = 0;
pub const RIGHT: usize = 1;

#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
pub enum Side {
    LEFT = 0,
    RIGHT = 1,
}

impl From<Side> for usize {
    fn from(val: Side) -> Self {
        val as usize
    }
}

pub mod linnet_ext;
pub mod symbolica_ext;
pub mod tracing;

#[cfg(test)]
mod tests {
    #[test]
    fn native_hypot_preserves_finite_range() {
        use super::{ArbPrec, F};

        for scale in [1e-200, 1.0, 1e200] {
            for sign in [-1.0, 1.0] {
                let norm = F(sign * 3.0 * scale).hypot(&F(4.0 * scale)).0;
                assert!(norm.is_finite());
                assert!((norm / scale - 5.0).abs() < 1e-14);
            }
        }
        assert_eq!(F(0.0).hypot(&F(-0.0)), F(0.0));
        assert!(F(f64::NAN).hypot(&F(1.0)).0.is_nan());
        assert!(F(1.0).hypot(&F(f64::NAN)).0.is_nan());
        assert!(F(f64::INFINITY).hypot(&F(1.0)).0.is_infinite());
        assert!(F(f64::MAX).hypot(&F(f64::MAX)).0.is_infinite());

        // Native precision must not pass through binary64 for scaling/order.
        let one = F::<ArbPrec>::default().one();
        for exponent in [-2000, 2000] {
            let scale = one.from_i64(2).powi(exponent);
            let norm = (&scale * one.from_i64(3)).hypot(&(&scale * one.from_i64(4)));
            assert!((norm / scale - one.from_i64(5)).abs() < one.from_i64(10).powi(-50));
        }
    }

    #[test]
    fn canonical_arb_materialization_preserves_native_bits_and_checks_range() {
        use super::{ArbPrec, F, FloatLike, QuadFloat, SamplingEvaluationError};
        use rug::Float;

        let low = Float::with_val(1000, 1) >> 400;
        let source = ArbPrec {
            float: Float::with_val(1000, 1) + &low,
        };
        let quad = F::<QuadFloat>::from_arb(&source).unwrap();
        let (lower, upper) = quad.0.mpfr_enclosure(1000);
        assert_eq!(lower, source.float);
        assert_eq!(upper, source.float);
        assert_eq!(
            F::<ArbPrec>::from_arb(&source).unwrap().0.float,
            source.float
        );
        assert_eq!(F::<f64>::from_arb(&source).unwrap().0, 1.0);

        let one = F::<ArbPrec>::default().one();
        let large = one.from_i64(2).powi(2000);
        for value in [large.clone(), -large.clone(), large.inv(), -large.inv()] {
            assert!(value.0.float.is_finite() && value != one.zero());
            for error in [
                F::<f64>::from_arb(&value.0).unwrap_err(),
                F::<QuadFloat>::from_arb(&value.0).unwrap_err(),
            ] {
                assert!(matches!(
                    error.downcast_ref::<SamplingEvaluationError>(),
                    Some(SamplingEvaluationError::Unrepresentable { .. })
                ));
            }
            assert_eq!(F::<ArbPrec>::from_arb(&value.0).unwrap(), value);
        }
        assert_eq!(F::<f64>::from_arb(&one.zero().0).unwrap().0, 0.0);
        let invalid = ArbPrec {
            float: Float::with_val(1000, f64::NAN),
        };
        assert!(matches!(
            F::<ArbPrec>::from_arb(&invalid)
                .unwrap_err()
                .downcast_ref::<SamplingEvaluationError>(),
            Some(SamplingEvaluationError::Unrepresentable { .. })
        ));
    }

    #[test]
    fn canonical_native_embedding_requires_finite_singleton_enclosures() {
        use super::{F, QuadFloat, SamplingEvaluationError};
        use rug::Float;
        use symbolica::domains::float::DoubleFloat;
        for sign in [-1.0, 1.0] {
            let low = sign * 2.0_f64.powi(-400);
            let native = F(QuadFloat::from(DoubleFloat::from_compensated_sum(
                sign, low,
            )));
            let exact = Float::with_val(1000, sign) + Float::with_val(1000, low);
            assert_eq!(native.to_arb_exact().unwrap().0.float, exact);
            // The exact sum spans 1075 bits, although both binary64 limbs are
            // finite. Rounding it to the container would change the source.
            let separated = F(QuadFloat::from(DoubleFloat::from_compensated_sum(
                sign,
                sign * f64::from_bits(1),
            )));
            assert!(matches!(
                separated
                    .to_arb_exact()
                    .unwrap_err()
                    .downcast_ref::<SamplingEvaluationError>(),
                Some(SamplingEvaluationError::Unrepresentable {
                    operation: "canonical source embedding",
                    ..
                })
            ));
        }
        for value in [0.0, -0.0, f64::from_bits(1), f64::MAX] {
            assert_eq!(
                F(value).to_arb_exact().unwrap().0.float,
                Float::with_val(1000, value)
            );
        }
        assert!(F(-0.0).to_arb_exact().unwrap().0.float.is_sign_negative());
        for zero in [0.0, -0.0] {
            let native = F(QuadFloat::from_f64_exact_binary(zero));
            assert_eq!(
                native.to_arb_exact().unwrap().0.float.is_sign_negative(),
                native.0.as_f64().is_sign_negative()
            );
        }
        for value in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
            assert!(F(value).to_arb_exact().is_err());
        }
    }

    #[test]
    fn fixed256_source_scalar_embeds_exactly_without_binary64_range_limits() {
        use super::{ArbPrec, F, FloatLike, SamplingFloat, SamplingPrecision};
        use crate::settings::runtime::Precision;
        use rug::Float;

        assert_eq!(
            SamplingFloat::sampling_precision(),
            SamplingPrecision::Fixed256
        );
        for (physical, sampling) in [
            (Precision::Double, SamplingPrecision::Double),
            (Precision::Quad, SamplingPrecision::Quad),
            (Precision::Arb, SamplingPrecision::Arb),
        ] {
            assert_eq!(SamplingPrecision::from(physical), sampling);
        }
        let one = F::<SamplingFloat>::default().one();
        assert_eq!(one.epsilon(), one.from_i64(2).powi(-255));
        for exponent in [-4096, -255, 0, 255, 4096] {
            for sign in [-1, 1] {
                let native = one.from_i64(sign)
                    * one.from_i64(2).powi(exponent)
                    * (&one + one.from_i64(2).powi(-254));
                let promoted = native.to_arb_exact().unwrap();
                assert_eq!(promoted.0.float, native.0.float);
                assert_eq!(F::<SamplingFloat>::from_arb(&promoted.0).unwrap(), native);
                assert_eq!(promoted.0.float.prec(), 1000);
                assert_eq!(native.0.float.prec(), 256);
            }
        }
        for zero in [0.0, -0.0] {
            let native = F(SamplingFloat::from_f64_exact_binary(zero));
            assert_eq!(
                native.to_arb_exact().unwrap().0.float.is_sign_negative(),
                zero.is_sign_negative()
            );
        }
        for invalid in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
            assert!(
                F(SamplingFloat::from(Float::with_val(256, invalid)))
                    .to_arb_exact()
                    .is_err()
            );
        }
        // Physical inputs remain Arb values; constructing proposal externals at
        // 256 bits does not redefine the original higher-precision kinematics.
        let physical =
            F::<ArbPrec>::default().one() + F::<ArbPrec>::default().from_i64(2).powi(-300);
        assert_eq!(F::<SamplingFloat>::from_arb(&physical.0).unwrap(), one);
        assert_ne!(physical, one.to_arb_exact().unwrap());
    }

    #[test]
    fn fixed_mpfr_rational_coefficients_are_rounded_once() {
        use super::VarFloat;
        use rug::{Float, Integer, Rational};
        use symbolica::domains::rational::Rational as SymbolicaRational;

        fn check<const N: u32>() {
            for (numerator, denominator) in [(0, 7), (1, 7), (-7, 3), (7, -3), (-7, -3)] {
                let exact = Rational::from((numerator, denominator));
                let rational = SymbolicaRational::from(exact.clone());
                let actual = VarFloat::<N>::from(&rational);
                assert_eq!(actual.float, Float::with_val(N, exact));
                assert_eq!(actual.float.prec(), N);
            }
            // Exercise Symbolica's Double and Large integer storage, with
            // either sign. The N-bit boundary also has an analytic rounding
            // result: (2^N+3)/(2^N+1) rounds to 1+2^(1-N), not 1+2^(2-N).
            for bits in [80, N] {
                let power = Integer::from(1) << bits;
                for sign in [-1, 1] {
                    let exact = Rational::from(((power.clone() + 3) * sign, power.clone() + 1));
                    let rational = SymbolicaRational::from(exact.clone());
                    let actual = VarFloat::<N>::from(&rational);
                    assert_eq!(actual.float, Float::with_val(N, exact));
                    if bits == N {
                        let expected =
                            (Float::with_val(N, 1) + (Float::with_val(N, 1) >> (N - 1))) * sign;
                        assert_eq!(actual.float, expected);
                    }
                }
            }
        }
        check::<256>();
        check::<1000>();
    }

    #[test]
    fn fixed_mpfr_symbolica_inputs_use_declared_precision_for_later_arithmetic() {
        use super::{SymbolicaFloat, VarFloat};
        use rug::Float;

        fn check<const N: u32>() {
            for input_precision in [53, N + 31] {
                let input = Float::with_val(input_precision, 1);
                let native = VarFloat::<N>::from(SymbolicaFloat::from(input.clone()));
                assert_eq!(native.float.prec(), N);
                assert_eq!(native.float, input);
                // Expanding a represented 53-bit input does not restore lost
                // information, but subsequent arithmetic must use all N bits.
                let small = VarFloat::<N>::from(Float::with_val(N, 1) >> 200_u32);
                let expected = Float::with_val(N, 1) + &small.float;
                assert_eq!((native + small).float, expected);
            }
            let precise = Float::with_val(N + 31, 1) + (Float::with_val(N + 31, 1) >> (N + 10));
            let native = VarFloat::<N>::from(SymbolicaFloat::from(precise.clone()));
            assert_eq!(native.float, Float::with_val(N, precise));
            assert_eq!(native.float.prec(), N);
            for value in [0.0, -0.0, f64::INFINITY, f64::NEG_INFINITY, f64::NAN] {
                let native = VarFloat::<N>::from(SymbolicaFloat::from(Float::with_val(53, value)));
                assert_eq!(native.float.prec(), N);
                if value.is_nan() {
                    assert!(native.float.is_nan());
                } else {
                    assert_eq!(native.float, value);
                    assert_eq!(native.float.is_sign_negative(), value.is_sign_negative());
                }
            }
        }
        check::<256>();
        check::<1000>();
    }

    #[test]
    fn canonical_scalar_relative_budget_is_directed_and_signed() {
        use super::{ArbPrec, F, FloatLike, QuadFloat, SamplingEvaluationError};
        fn check<T: FloatLike>() {
            let one = F::<ArbPrec>::default().one();
            let source = &one - one.from_i64(2).powi(-54);
            let tolerance = 2.0_f64.powi(-54);
            let loose = f64::from_bits(tolerance.to_bits() + 1);
            for sign in [-1, 1] {
                let native = F::<T>::default().from_i64(sign);
                let source = &source * one.from_i64(sign);
                // |native-source|=2^-54, but tolerance*|source| is strictly
                // smaller by 2^-108. Nearest rounding must not erase that gap.
                let error = native
                    .verify_arb_materialization(&source.0, tolerance)
                    .unwrap_err();
                assert!(matches!(
                    error.downcast_ref::<SamplingEvaluationError>(),
                    Some(SamplingEvaluationError::UncertainGeometry { .. })
                ));
                native.verify_arb_materialization(&source.0, loose).unwrap();
            }
            let zero = F::<T>::default().zero();
            zero.verify_arb_materialization(&one.zero().0, tolerance)
                .unwrap();
            assert!(
                zero.one()
                    .verify_arb_materialization(&one.zero().0, tolerance)
                    .is_err()
            );
            for budget in [0.0, -1.0, f64::INFINITY, f64::NAN] {
                assert!(
                    zero.verify_arb_materialization(&one.zero().0, budget)
                        .is_err()
                );
            }
        }
        check::<f64>();
        check::<QuadFloat>();
        check::<ArbPrec>();
    }

    #[test]
    fn native_mpfr_enclosures_preserve_separated_limbs_and_signs() {
        use super::{ArbPrec, FloatLike, QuadFloat};
        use rug::Float;
        use symbolica::domains::float::DoubleFloat;

        for sign in [-1.0, 1.0] {
            let low = sign * 2.0_f64.powi(-400);
            let quad = QuadFloat(DoubleFloat::from_compensated_sum(sign, low));
            let exact = Float::with_val(512, sign) + Float::with_val(512, low);
            for precision in [24, 106, 512] {
                let (lower, upper) = quad.mpfr_enclosure(precision);
                assert!(lower <= exact && upper >= exact);
                if precision == 512 {
                    assert_eq!(lower, exact);
                    assert_eq!(upper, exact);
                } else {
                    assert!(
                        lower < upper,
                        "a separated limb cannot disappear from a certificate"
                    );
                }
                let arb = ArbPrec {
                    float: exact.clone(),
                };
                let (lower, upper) = arb.mpfr_enclosure(precision);
                assert!(lower <= exact && upper >= exact);
                let native = sign * 1.23456789012345;
                let (lower, upper) = native.mpfr_enclosure(precision);
                let exact_native = Float::with_val(53, native);
                assert!(lower <= exact_native && upper >= exact_native);
            }
        }
    }

    #[test]
    fn normalized_h_functions_integrate_to_one() {
        use super::{F, h, h_dual, hyperdual_utils::new_constant};
        use crate::settings::runtime::{HFunction, HFunctionSettings};
        use symbolica::domains::dual::HyperDual;

        let template = HyperDual::from_values(vec![vec![0], vec![1]], vec![F(1.0), F(0.0)]);
        // With t = sigma e^u, the trapezoidal rule converges exponentially for
        // these profiles, which decay doubly exponentially at both ends.
        let step = 0.01;
        for function in [
            HFunction::Exponential,
            HFunction::PolyExponential,
            HFunction::PolyLeftRightExponential,
        ] {
            // The exponential profile ignores the power.
            for power in [None, Some(1), Some(4), Some(16)] {
                let settings = HFunctionSettings {
                    function: function.clone(),
                    sigma: 1.5,
                    power,
                    ..Default::default()
                };
                let integral = (-4000..=4000)
                    .map(|n| {
                        let t = F(settings.sigma * (n as f64 * step).exp());
                        let value = h(&t, None, None, &settings);
                        let dual = h_dual(&new_constant(&template, &t), None, None, &settings);
                        // Far tails amplify argument rounding but are negligible.
                        let scale = value.0.abs().max(1.0);
                        assert!((dual.values[0].0 - value.0).abs() <= 1e-14 * scale);
                        value.0 * t.0 * step
                    })
                    .sum::<f64>();
                assert!(
                    (integral - 1.0).abs() < 1e-12,
                    "{settings:?} integrates to {integral}"
                );
            }
        }
    }

    #[test]
    fn debug_tags_macro_supports_tags_and_regular_fields() {
        crate::debug_tags!(#integration, #summary;
            inspect = false,
            answer = 42,
            "macro compiles with tag shorthands and regular fields"
        );
    }
}

pub trait FloatConvertFrom<U> {
    fn convert_from(x: &U) -> Self;
}

//     fn convert_from(x: &f128::f128) -> f64 {
//         (*x).to_f64().unwrap()
//     }
// }

// impl FloatConvertFrom<f128::f128> for f128::f128 {
//     fn convert_from(x: &f128::f128) -> f128::f128 {
//         *x
//     }
// }

// impl FloatConvertFrom<f64> for f128::f128 {
//     fn convert_from(x: &f64) -> f128::f128 {
//         f128::f128::from_f64(*x).unwrap()
//     }
// }

#[derive(Debug, Clone, PartialEq, PartialOrd, Encode, Decode)]
pub struct VarFloat<const N: u32> {
    #[bincode(with_serde)]
    float: rug::Float,
}

impl<const N: u32> Serialize for VarFloat<N> {
    fn serialize<S>(&self, serializer: S) -> std::result::Result<S::Ok, S::Error>
    where
        S: serde::Serializer,
    {
        let string = self.float.to_string();
        string.serialize(serializer)
    }
}

impl<'de, const N: u32> Deserialize<'de> for VarFloat<N> {
    fn deserialize<D>(deserializer: D) -> std::result::Result<Self, D::Error>
    where
        D: Deserializer<'de>,
    {
        let string: String = serde::Deserialize::deserialize(deserializer)?;
        let val: Self = string
            .parse()
            .unwrap_or_else(|_| panic!("failed to parse arb prec from string: {}", string));

        Ok(val)
    }
}

impl<const N: u32> FromStr for VarFloat<N> {
    type Err = ParseFloatError;
    fn from_str(s: &str) -> std::result::Result<Self, Self::Err> {
        let float = rug::Float::parse(s)?;
        Ok(Self {
            float: rug::Float::with_val(N, float),
        })
    }
}

#[derive(
    Debug,
    Clone,
    Copy,
    Default,
    PartialEq,
    PartialOrd,
    Eq,
    Hash,
    Serialize,
    Deserialize,
    Encode,
    Decode,
)]
#[repr(transparent)]
pub struct QuadFloat(DoubleFloat);

impl std::fmt::Display for QuadFloat {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        std::fmt::Display::fmt(&self.0, f)
    }
}

impl std::fmt::LowerExp for QuadFloat {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        std::fmt::LowerExp::fmt(&self.0, f)
    }
}

impl From<DoubleFloat> for QuadFloat {
    fn from(value: DoubleFloat) -> Self {
        Self(value)
    }
}

impl From<QuadFloat> for DoubleFloat {
    fn from(value: QuadFloat) -> Self {
        value.0
    }
}

impl From<&QuadFloat> for DoubleFloat {
    fn from(value: &QuadFloat) -> Self {
        value.0
    }
}

impl From<SymbolicaFloat> for QuadFloat {
    fn from(value: SymbolicaFloat) -> Self {
        Self(value.to_double_float())
    }
}

impl From<Float> for QuadFloat {
    fn from(value: Float) -> Self {
        Self(SymbolicaFloat::from(value).to_double_float())
    }
}

impl From<&Rational> for QuadFloat {
    fn from(value: &Rational) -> Self {
        Self(value.into())
    }
}

impl<const N: u32> Rem<&VarFloat<N>> for &VarFloat<N> {
    type Output = VarFloat<N>;

    fn rem(self, rhs: &VarFloat<N>) -> Self::Output {
        (&self.float % &rhs.float).complete(N).into()
    }
}

impl<const N: u32> From<Float> for VarFloat<N> {
    fn from(x: Float) -> Self {
        VarFloat {
            float: rug::Float::with_val(N, x),
        }
    }
}

impl<const N: u32> From<&Rational> for VarFloat<N> {
    fn from(x: &Rational) -> Self {
        // Round the exact rational once; separately rounded numerator and
        // denominator can change the result by an ulp.
        Self::from(x.to_multi_prec_float(N))
    }
}

macro_rules! impl_quad_binary_op {
    ($trait:ident, $fn:ident, $op:tt) => {
        impl std::ops::$trait for QuadFloat {
            type Output = Self;

            fn $fn(self, rhs: Self) -> Self::Output {
                Self(self.0 $op rhs.0)
            }
        }

        impl std::ops::$trait<&QuadFloat> for QuadFloat {
            type Output = Self;

            fn $fn(self, rhs: &Self) -> Self::Output {
                Self(self.0 $op rhs.0)
            }
        }

        impl std::ops::$trait<QuadFloat> for &QuadFloat {
            type Output = QuadFloat;

            fn $fn(self, rhs: QuadFloat) -> Self::Output {
                QuadFloat(self.0 $op rhs.0)
            }
        }

        impl std::ops::$trait<&QuadFloat> for &QuadFloat {
            type Output = QuadFloat;

            fn $fn(self, rhs: &QuadFloat) -> Self::Output {
                QuadFloat(self.0 $op rhs.0)
            }
        }
    };
}

macro_rules! impl_quad_assign_op {
    ($trait:ident, $fn:ident, $op:tt) => {
        impl std::ops::$trait for QuadFloat {
            fn $fn(&mut self, rhs: Self) {
                self.0 $op rhs.0;
            }
        }

        impl std::ops::$trait<&QuadFloat> for QuadFloat {
            fn $fn(&mut self, rhs: &Self) {
                self.0 $op rhs.0;
            }
        }
    };
}

impl_quad_binary_op!(Add, add, +);
impl_quad_binary_op!(Sub, sub, -);
impl_quad_binary_op!(Mul, mul, *);
impl_quad_binary_op!(Div, div, /);

impl_quad_assign_op!(AddAssign, add_assign, +=);
impl_quad_assign_op!(SubAssign, sub_assign, -=);
impl_quad_assign_op!(MulAssign, mul_assign, *=);
impl_quad_assign_op!(DivAssign, div_assign, /=);

impl Neg for QuadFloat {
    type Output = Self;

    fn neg(self) -> Self::Output {
        Self(-self.0)
    }
}

impl Neg for &QuadFloat {
    type Output = QuadFloat;

    fn neg(self) -> Self::Output {
        QuadFloat(-self.0)
    }
}

impl Rem<&QuadFloat> for &QuadFloat {
    type Output = QuadFloat;

    fn rem(self, rhs: &QuadFloat) -> Self::Output {
        self.truncating_remainder(rhs)
    }
}

impl<const N: u32> std::ops::Mul for VarFloat<N> {
    type Output = Self;

    fn mul(self, rhs: Self) -> Self::Output {
        (self.float * rhs.float).into()
    }
}

impl<const N: u32> std::ops::Mul<&VarFloat<N>> for VarFloat<N> {
    type Output = Self;

    fn mul(self, rhs: &Self) -> Self::Output {
        (self.float * &rhs.float).into()
    }
}

impl<const N: u32> std::ops::Mul<VarFloat<N>> for &VarFloat<N> {
    type Output = VarFloat<N>;

    fn mul(self, rhs: VarFloat<N>) -> Self::Output {
        (&self.float * &rhs.float).complete(N).into()
    }
}

impl<const N: u32> std::ops::Mul<&VarFloat<N>> for &VarFloat<N> {
    type Output = VarFloat<N>;

    fn mul(self, rhs: &VarFloat<N>) -> Self::Output {
        (&self.float * &rhs.float).complete(N).into()
    }
}

impl<const N: u32> std::ops::MulAssign for VarFloat<N> {
    fn mul_assign(&mut self, rhs: Self) {
        self.float *= rhs.float;
        self.float.set_prec(N);
    }
}

impl<const N: u32> std::ops::MulAssign<&VarFloat<N>> for VarFloat<N> {
    fn mul_assign(&mut self, rhs: &Self) {
        self.float *= &rhs.float;
        self.float.set_prec(N);
    }
}

impl<const N: u32> std::ops::Add for VarFloat<N> {
    type Output = Self;

    fn add(self, rhs: Self) -> Self::Output {
        (self.float + rhs.float).into()
    }
}

impl<const N: u32> std::ops::Add<&VarFloat<N>> for VarFloat<N> {
    type Output = Self;

    fn add(self, rhs: &Self) -> Self::Output {
        (self.float + &rhs.float).into()
    }
}

impl<const N: u32> std::ops::Add<VarFloat<N>> for &VarFloat<N> {
    type Output = VarFloat<N>;

    fn add(self, rhs: VarFloat<N>) -> Self::Output {
        (&self.float + &rhs.float).complete(N).into()
    }
}

impl<const N: u32> std::ops::Add<&VarFloat<N>> for &VarFloat<N> {
    type Output = VarFloat<N>;

    fn add(self, rhs: &VarFloat<N>) -> Self::Output {
        (&self.float + &rhs.float).complete(N).into()
    }
}

impl<const N: u32> std::ops::AddAssign for VarFloat<N> {
    fn add_assign(&mut self, rhs: Self) {
        self.float += rhs.float;
        self.float.set_prec(N);
    }
}

impl<const N: u32> std::ops::AddAssign<&VarFloat<N>> for VarFloat<N> {
    fn add_assign(&mut self, rhs: &Self) {
        self.float += &rhs.float;
        self.float.set_prec(N);
    }
}

impl<const N: u32> std::ops::Sub for VarFloat<N> {
    type Output = Self;

    fn sub(self, rhs: Self) -> Self::Output {
        (self.float - rhs.float).into()
    }
}

impl<const N: u32> std::ops::Sub<&VarFloat<N>> for VarFloat<N> {
    type Output = Self;

    fn sub(self, rhs: &Self) -> Self::Output {
        (self.float - &rhs.float).into()
    }
}

impl<const N: u32> std::ops::Sub<VarFloat<N>> for &VarFloat<N> {
    type Output = VarFloat<N>;

    fn sub(self, rhs: VarFloat<N>) -> Self::Output {
        (&self.float - &rhs.float).complete(N).into()
    }
}

impl<const N: u32> std::ops::Sub<&VarFloat<N>> for &VarFloat<N> {
    type Output = VarFloat<N>;

    fn sub(self, rhs: &VarFloat<N>) -> Self::Output {
        (&self.float - &rhs.float).complete(N).into()
    }
}

impl<const N: u32> std::ops::SubAssign for VarFloat<N> {
    fn sub_assign(&mut self, rhs: Self) {
        self.float -= rhs.float;
        self.float.set_prec(N);
    }
}

impl<const N: u32> std::ops::SubAssign<&VarFloat<N>> for VarFloat<N> {
    fn sub_assign(&mut self, rhs: &Self) {
        self.float -= &rhs.float;
        self.float.set_prec(N);
    }
}

impl<const N: u32> Div for VarFloat<N> {
    type Output = Self;

    fn div(self, rhs: Self) -> Self::Output {
        (self.float / rhs.float).into()
    }
}

impl<const N: u32> Div<&VarFloat<N>> for VarFloat<N> {
    type Output = Self;

    fn div(self, rhs: &Self) -> Self::Output {
        (self.float / &rhs.float).into()
    }
}

impl<const N: u32> Div<VarFloat<N>> for &VarFloat<N> {
    type Output = VarFloat<N>;

    fn div(self, rhs: VarFloat<N>) -> Self::Output {
        (&self.float / &rhs.float).complete(N).into()
    }
}

impl<const N: u32> Div<&VarFloat<N>> for &VarFloat<N> {
    type Output = VarFloat<N>;

    fn div(self, rhs: &VarFloat<N>) -> Self::Output {
        (&self.float / &rhs.float).complete(N).into()
    }
}

impl<const N: u32> std::ops::DivAssign for VarFloat<N> {
    fn div_assign(&mut self, rhs: Self) {
        self.float /= rhs.float;
        self.float.set_prec(N);
    }
}

impl<const N: u32> std::ops::DivAssign<&VarFloat<N>> for VarFloat<N> {
    fn div_assign(&mut self, rhs: &Self) {
        self.float /= &rhs.float;
        self.float.set_prec(N);
    }
}

impl<const N: u32> std::ops::Neg for VarFloat<N> {
    type Output = Self;

    fn neg(self) -> Self::Output {
        (-self.float).into()
    }
}

impl<const N: u32> std::ops::Neg for &VarFloat<N> {
    type Output = VarFloat<N>;

    fn neg(self) -> Self::Output {
        (-&self.float).complete(N).into()
    }
}

impl<const N: u32> std::fmt::Display for VarFloat<N> {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(f, "{}", self.float)
    }
}

impl<const N: u32> std::fmt::LowerExp for VarFloat<N> {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(f, "{:e}", self.float)
    }
}

impl<const N: u32> RefZero for VarFloat<N> {
    fn ref_zero(&self) -> Self {
        Self::new_zero()
    }
}

impl<const N: u32> R for VarFloat<N> {}

impl<const N: u32> RefOne for VarFloat<N> {
    fn ref_one(&self) -> Self {
        self.one()
    }
}

impl<const N: u32> SymFloatLike for VarFloat<N> {
    #[inline]
    fn set_from(&mut self, other: &Self) {
        self.float.assign(&other.float);
    }

    fn mul_add(&self, a: &Self, b: &Self) -> Self {
        (&self.float * &a.float + &b.float).complete(N).into()
    }

    fn is_fully_zero(&self) -> bool {
        self.float.is_zero()
    }

    // fn norm(&self) -> Self {
    //     self.float.clone().abs().into()
    // }

    fn from_i64(&self, a: i64) -> Self {
        VarFloat {
            float: Float::with_val(N, a),
        }
    }

    fn from_usize(&self, a: usize) -> Self {
        VarFloat {
            float: Float::with_val(N, a),
        }
    }

    fn get_precision(&self) -> u32 {
        N
    }

    fn set_precision(&mut self, _precision: u32) {}

    fn zero(&self) -> Self {
        Self::new_zero()
    }

    fn one(&self) -> Self {
        self.from_i64(1)
    }

    fn new_zero() -> Self {
        VarFloat {
            float: Float::new(N),
        }
    }

    fn inv(&self) -> Self {
        self.float.clone().recip().into()
    }

    fn pow(&self, e: u64) -> Self {
        rug::ops::Pow::pow(&self.float, e).complete(N).into()
    }

    fn sample_unit<R: Rng + ?Sized>(&self, rng: &mut R) -> Self {
        let f: f64 = rng.random();
        Float::with_val(N, f).into()
    }

    fn neg(&self) -> Self {
        (-self.float.clone()).into()
    }

    fn get_epsilon(&self) -> f64 {
        2.0f64.powi(1 - (N as i32))
    }

    fn fixed_precision(&self) -> bool {
        true
    }
}

impl<const N: u32> SingleFloat for VarFloat<N> {
    fn is_finite(&self) -> bool {
        self.float.is_finite()
    }

    fn is_one(&self) -> bool {
        self.float == 1.
    }

    fn is_zero(&self) -> bool {
        self.float == 0.
    }

    fn from_rational(&self, rat: &Rational) -> Self {
        rat.into()
    }
}
impl<const N: u32> RealLike for VarFloat<N> {
    fn to_f64(&self) -> f64 {
        self.float.to_f64()
    }

    fn round_to_nearest_integer(&self) -> Integer {
        symbolica::domains::integer::MultiPrecisionInteger::from_raw(
            self.float.clone().round().to_integer().unwrap(),
        )
        .into()
    }

    fn to_usize_clamped(&self) -> usize {
        self.float
            .to_integer()
            .unwrap()
            .to_usize()
            .unwrap_or(usize::MAX)
    }
}

impl<const N: u32> Real for VarFloat<N> {
    fn i(&self) -> Option<Self> {
        None
    }

    fn conj(&self) -> Self {
        self.clone()
    }

    #[inline(always)]
    fn pi(&self) -> Self {
        Float::with_val(N, rug::float::Constant::Pi).into()
    }

    #[inline(always)]
    fn e(&self) -> Self {
        self.one().exp()
    }

    #[inline(always)]
    fn euler(&self) -> Self {
        Float::with_val(N, rug::float::Constant::Euler).into()
    }

    #[inline(always)]
    fn phi(&self) -> Self {
        (self.one() + self.from_i64(5).sqrt()) / Self::from_f64(2.)
    }
    fn atan2(&self, x: &Self) -> Self {
        self.float.clone().atan2(&x.float).into()
    }

    fn powf(&self, e: &Self) -> Self {
        self.float.clone().pow(e.float.clone()).into()
    }

    fn log(&self) -> Self {
        self.float.ln_ref().complete(N).into()
    }
    fn norm(&self) -> Self {
        self.float.clone().abs().into()
    }

    delegate! {
        #[into]
        to self.float.clone(){
            fn sqrt(&self) -> Self;
            fn exp(&self) -> Self;
            fn sin(&self) -> Self;
            fn cos(&self) -> Self;
            fn tan(&self) -> Self;
            fn asin(&self) -> Self;
            fn acos(&self) -> Self;
            fn sinh(&self) -> Self;
            fn cosh(&self) -> Self;
            fn tanh(&self) -> Self;
            fn asinh(&self) -> Self;
            fn acosh(&self) -> Self;
            fn atanh(&self) -> Self;
        }
    }
}

impl<const N: u32> From<VarFloat<N>> for SymbolicaFloat {
    fn from(value: VarFloat<N>) -> Self {
        SymbolicaFloat::from(value.float)
    }
}

impl<const N: u32> From<&VarFloat<N>> for SymbolicaFloat {
    fn from(value: &VarFloat<N>) -> Self {
        SymbolicaFloat::from(value.float.clone())
    }
}

impl<const N: u32> From<SymbolicaFloat> for VarFloat<N> {
    fn from(value: SymbolicaFloat) -> Self {
        // Preserve the represented input at this lane's precision so later
        // owned arithmetic cannot silently retain the input's smaller precision.
        Self::from(value.into_raw())
    }
}

impl From<QuadFloat> for SymbolicaFloat {
    fn from(value: QuadFloat) -> Self {
        value.0.into()
    }
}

impl From<&QuadFloat> for SymbolicaFloat {
    fn from(value: &QuadFloat) -> Self {
        value.0.into()
    }
}

impl RefZero for QuadFloat {
    fn ref_zero(&self) -> Self {
        Self::new_zero()
    }
}

impl R for QuadFloat {}

impl RefOne for QuadFloat {
    fn ref_one(&self) -> Self {
        self.one()
    }
}

impl SymFloatLike for QuadFloat {
    #[inline]
    fn set_from(&mut self, other: &Self) {
        self.0.set_from(&other.0);
    }

    fn mul_add(&self, a: &Self, b: &Self) -> Self {
        Self(self.0.mul_add(&a.0, &b.0))
    }

    fn is_fully_zero(&self) -> bool {
        self.0.is_fully_zero()
    }

    fn from_i64(&self, a: i64) -> Self {
        Self(self.0.from_i64(a))
    }

    fn from_usize(&self, a: usize) -> Self {
        Self(self.0.from_usize(a))
    }

    fn get_precision(&self) -> u32 {
        self.0.get_precision()
    }

    fn set_precision(&mut self, precision: u32) {
        self.0.set_precision(precision);
    }

    fn zero(&self) -> Self {
        Self(self.0.zero())
    }

    fn one(&self) -> Self {
        Self(self.0.one())
    }

    fn new_zero() -> Self {
        Self(DoubleFloat::default().zero())
    }

    fn inv(&self) -> Self {
        Self(self.0.inv())
    }

    fn pow(&self, e: u64) -> Self {
        // Old direct wrapper around Symbolica's DoubleFloat::pow -> Df64::powi.
        // Keep this inactive until the upstream bug is fixed in Symbolica.
        // Self(self.0.pow(e))

        // Temporary local workaround kept active in GammaLoop for now.
        // Equivalent explicit implementation:
        // let mut base = *self;
        // let mut acc = Self::one();
        // let mut exponent = e;
        // while exponent > 0 {
        //     if exponent & 1 == 1 {
        //         acc *= base;
        //     }
        //     exponent >>= 1;
        //     if exponent > 0 {
        //         base *= base;
        //     }
        // }
        // acc
        self.pow_u64_via_mul(e)
    }

    fn sample_unit<R: Rng + ?Sized>(&self, rng: &mut R) -> Self {
        Self(self.0.sample_unit(rng))
    }

    fn neg(&self) -> Self {
        Self(self.0.neg())
    }

    fn get_epsilon(&self) -> f64 {
        self.0.get_epsilon()
    }

    fn fixed_precision(&self) -> bool {
        self.0.fixed_precision()
    }
}

impl SingleFloat for QuadFloat {
    fn is_finite(&self) -> bool {
        self.0.is_finite()
    }

    fn is_one(&self) -> bool {
        self.0.is_one()
    }

    fn is_zero(&self) -> bool {
        self.0.is_zero()
    }

    fn from_rational(&self, rat: &Rational) -> Self {
        rat.into()
    }
}

impl RealLike for QuadFloat {
    fn to_f64(&self) -> f64 {
        self.0.to_f64()
    }

    fn round_to_nearest_integer(&self) -> Integer {
        self.0.round_to_nearest_integer()
    }

    fn to_usize_clamped(&self) -> usize {
        self.0.to_usize_clamped()
    }
}

impl Real for QuadFloat {
    fn i(&self) -> Option<Self> {
        None
    }

    fn conj(&self) -> Self {
        *self
    }

    #[inline(always)]
    fn pi(&self) -> Self {
        Self(self.0.pi())
    }

    #[inline(always)]
    fn e(&self) -> Self {
        Self(self.0.e())
    }

    #[inline(always)]
    fn euler(&self) -> Self {
        Self(self.0.euler())
    }

    #[inline(always)]
    fn phi(&self) -> Self {
        Self(self.0.phi())
    }

    fn atan2(&self, x: &Self) -> Self {
        Self(self.0.atan2(&x.0))
    }

    fn powf(&self, e: &Self) -> Self {
        Self(self.0.powf(&e.0))
    }

    fn log(&self) -> Self {
        Self(self.0.log())
    }

    fn norm(&self) -> Self {
        Self(self.0.norm())
    }

    fn sqrt(&self) -> Self {
        Self(self.0.sqrt())
    }

    fn exp(&self) -> Self {
        Self(self.0.exp())
    }

    fn sin(&self) -> Self {
        Self(self.0.sin())
    }

    fn cos(&self) -> Self {
        Self(self.0.cos())
    }

    fn tan(&self) -> Self {
        Self(self.0.tan())
    }

    fn asin(&self) -> Self {
        Self(self.0.asin())
    }

    fn acos(&self) -> Self {
        Self(self.0.acos())
    }

    fn sinh(&self) -> Self {
        Self(self.0.sinh())
    }

    fn cosh(&self) -> Self {
        Self(self.0.cosh())
    }

    fn tanh(&self) -> Self {
        Self(self.0.tanh())
    }

    fn asinh(&self) -> Self {
        Self(self.0.asinh())
    }

    fn acosh(&self) -> Self {
        Self(self.0.acosh())
    }

    fn atanh(&self) -> Self {
        Self(self.0.atanh())
    }
}

impl FloatLike for f128 {
    fn mpfr_enclosure(&self, precision: u32) -> (Float, Float) {
        let limbs = self.0.into_inner();
        let high = Float::with_val(53, limbs.hi());
        let low = Float::with_val(53, limbs.lo());
        (
            Float::with_val_round(precision, &high + &low, Round::Down).0,
            Float::with_val_round(precision, &high + &low, Round::Up).0,
        )
    }

    fn E(&self) -> Self {
        Self::E()
    }

    fn PIHALF(&self) -> Self {
        Self::PIHALF()
    }

    fn SQRT_2(&self) -> Self {
        Self::from_f64(2.0).sqrt()
    }

    fn SQRT_2_HALF(&self) -> Self {
        Self::from_f64(2.0).sqrt() / Self::from_f64(2.0)
    }

    fn rem_euclid(&self, rhs: &Self) -> Self {
        let r = self.ref_rem(rhs);
        if r < r.zero() { r + rhs } else { r }
    }

    fn FRAC_1_PI(&self) -> Self {
        Self::FRAC_1_PI()
    }

    fn PI(&self) -> Self {
        Self::PI()
    }

    fn TAU(&self) -> Self {
        Self::TAU()
    }

    fn from_f64_exact_binary(x: f64) -> Self {
        Self::from_f64_exact_binary(x)
    }

    fn from_arb(value: &ArbPrec) -> Self {
        Self::from(value.float.clone())
    }

    fn sampling_precision() -> SamplingPrecision {
        SamplingPrecision::Quad
    }

    fn sampling_bridge_cache(
        setup: &LmbMultiChannelingSetup,
    ) -> &RuntimeCache<std::result::Result<SamplingChannelBridge<Self>, SamplingEvaluationError>>
    {
        &setup.sampling_bridge_quad
    }

    fn sampling_bridge_cache_mut(
        setup: &mut LmbMultiChannelingSetup,
    ) -> &mut RuntimeCache<std::result::Result<SamplingChannelBridge<Self>, SamplingEvaluationError>>
    {
        &mut setup.sampling_bridge_quad
    }

    fn from_f64(x: f64) -> Self {
        // There are two reasonable f64 -> higher-precision policies:
        // preserve the exact binary64 value or reinterpret the visible decimal
        // spelling of the f64. GammaLoop currently chooses the decimal route here
        // because these upcasts are overwhelmingly user-authored settings, and
        // values like 0.1 are less surprising when they retain their decimal
        // semantics instead of exposing the hidden binary64 tail. The exact-binary
        // helper is kept alongside this for callers that need a faithful embedding
        // of an already-computed f64.
        Self::from_f64_decimal(x)
    }

    fn into_f64(&self) -> f64 {
        self.as_f64()
    }

    fn is_nan(&self) -> bool {
        self.is_nan_value()
    }

    fn is_infinite(&self) -> bool {
        self.is_infinite_value()
    }

    fn floor(&self) -> Self {
        self.floor_value()
    }

    fn try_extract_externals_from_cache(
        externals: &Externals,
    ) -> Option<&TiVec<ExternalIndex, FourMomentum<F<Self>>>> {
        match externals {
            Externals::Constant { f_128_cache, .. } => f_128_cache.as_ref(),
        }
    }

    fn epsilon(&self) -> Self {
        self.machine_epsilon()
    }
}

// Only the runtime cache hooks differ between the two fixed MPFR lanes.
macro_rules! impl_mpfr_float_like {
    ($scalar:ty, $precision:ident, $bridge:ident, $externals:ident => $external_cache:expr) => {
        impl FloatLike for $scalar {
            fn mpfr_enclosure(&self, precision: u32) -> (Float, Float) {
                (
                    Float::with_val_round(precision, &self.float, Round::Down).0,
                    Float::with_val_round(precision, &self.float, Round::Up).0,
                )
            }

            fn E(&self) -> Self {
                Self::E()
            }

            fn PIHALF(&self) -> Self {
                Self::PIHALF()
            }

            fn SQRT_2(&self) -> Self {
                self.from_i64(2).sqrt()
            }

            fn SQRT_2_HALF(&self) -> Self {
                self.from_i64(2).sqrt() / self.from_i64(2)
            }

            fn rem_euclid(&self, rhs: &Self) -> Self {
                let r = self.ref_rem(rhs);
                if r < r.zero() { r + rhs } else { r }
            }

            fn FRAC_1_PI(&self) -> Self {
                Self::FRAC_1_PI()
            }

            fn PI(&self) -> Self {
                Self::PI()
            }

            fn TAU(&self) -> Self {
                Self::TAU()
            }

            fn from_f64_exact_binary(x: f64) -> Self {
                Self::from_f64_exact_binary(x)
            }

            fn from_arb(value: &ArbPrec) -> Self {
                Self::from(value.float.clone())
            }

            fn sampling_precision() -> SamplingPrecision {
                SamplingPrecision::$precision
            }

            fn sampling_bridge_cache(
                setup: &LmbMultiChannelingSetup,
            ) -> &RuntimeCache<
                std::result::Result<SamplingChannelBridge<Self>, SamplingEvaluationError>,
            > {
                &setup.$bridge
            }

            fn sampling_bridge_cache_mut(
                setup: &mut LmbMultiChannelingSetup,
            ) -> &mut RuntimeCache<
                std::result::Result<SamplingChannelBridge<Self>, SamplingEvaluationError>,
            > {
                &mut setup.$bridge
            }

            fn from_f64(x: f64) -> Self {
                // There are two reasonable f64 -> higher-precision policies:
                // preserve the exact binary64 value or reinterpret the visible decimal
                // spelling of the f64. GammaLoop currently chooses the decimal route here
                // because these upcasts are overwhelmingly user-authored settings, and
                // values like 0.1 are less surprising when they retain their decimal
                // semantics instead of exposing the hidden binary64 tail. The exact-binary
                // helper is kept alongside this for callers that need a faithful embedding
                // of an already-computed f64.
                VarFloat::from_f64(x)
            }

            fn into_f64(&self) -> f64 {
                self.to_f64()
            }

            fn is_nan(&self) -> bool {
                self.float.is_nan()
            }

            fn is_infinite(&self) -> bool {
                self.float.is_infinite()
            }

            fn floor(&self) -> Self {
                self.float.clone().floor().into()
            }

            fn try_extract_externals_from_cache(
                $externals: &Externals,
            ) -> Option<&TiVec<ExternalIndex, FourMomentum<F<Self>>>> {
                $external_cache
            }

            fn epsilon(&self) -> Self {
                self.machine_epsilon()
            }
        }

        impl PrecisionUpgradable for $scalar {
            type Higher = ArbPrec;
            type Lower = f128;

            fn higher(&self) -> Self::Higher {
                ArbPrec::from(self.float.clone())
            }

            fn lower(&self) -> Self::Lower {
                f128::from(self.float.clone())
            }
        }
    };
}

impl_mpfr_float_like!(ArbPrec, Arb, sampling_bridge_arb, externals => match externals {
    Externals::Constant { arb_cache, .. } => arb_cache.as_ref(),
});
impl_mpfr_float_like!(SamplingFloat, Fixed256, sampling_bridge_fixed256, _externals => None);

impl<const N: u32> VarFloat<N> {
    fn machine_epsilon(&self) -> Self {
        self.from_i64(2).pow((N - 1) as u64).inv()
    }

    fn one() -> Self {
        VarFloat {
            float: rug::Float::with_val(N, 1.0),
        }
    }
    #[allow(non_snake_case)]
    fn E() -> Self {
        Self::one().exp()
    }

    #[allow(non_snake_case)]
    fn PIHALF() -> Self {
        Self::PI() / Self::from_f64(2.0)
    }

    #[allow(non_snake_case)]
    fn PI() -> Self {
        VarFloat {
            float: rug::Float::with_val(N, Constant::Pi),
        }
    }

    #[allow(non_snake_case)]
    fn TAU() -> Self {
        let mut tau = Self::PI() + Self::PI();
        tau.float.set_prec(N);
        tau
    }

    #[allow(non_snake_case)]
    fn FRAC_1_PI() -> Self {
        Self::PI().inv()
    }

    #[allow(dead_code)]
    pub(crate) fn from_f64_exact_binary(x: f64) -> Self {
        VarFloat {
            float: rug::Float::with_val(N, x),
        }
    }

    pub(crate) fn from_f64_decimal(x: f64) -> Self {
        if !x.is_finite() {
            return Self::from_f64_exact_binary(x);
        }

        // Decimal parsing can acquire a process-wide locale lock. Cache repeated
        // settings/model values per worker, including the target precision in the
        // key: a thread-local static is shared by all const-generic instantiations.
        thread_local! {
            static CACHE: RefCell<HashMap<(u32, u64), Float>> = RefCell::default();
        }
        CACHE.with_borrow_mut(|cache| {
            let key = (N, x.to_bits());
            if let Some(float) = cache.get(&key) {
                return Self {
                    float: float.clone(),
                };
            }
            let valid = Float::parse(format!("{}", x)).unwrap();
            let float = Float::with_val(N, valid);
            // Dynamic inputs must not cause unbounded retention. Eviction only
            // repeats the identical conversion; it never changes its semantics.
            if cache.len() == 1024 {
                cache.clear();
            }
            cache.insert(key, float.clone());
            Self { float }
        })
    }

    pub(crate) fn from_f64(x: f64) -> Self {
        // There are two reasonable f64 -> higher-precision policies:
        // preserve the exact binary64 value or reinterpret the visible decimal
        // spelling of the f64. GammaLoop currently chooses the decimal route here
        // because these upcasts are overwhelmingly user-authored settings, and
        // values like 0.1 are less surprising when they retain their decimal
        // semantics instead of exposing the hidden binary64 tail. The exact-binary
        // helper is kept alongside this for callers that need a faithful embedding
        // of an already-computed f64.
        Self::from_f64_decimal(x)
    }

    pub(crate) fn to_f64(&self) -> f64 {
        self.float.to_f64()
    }

    #[allow(dead_code)]
    fn as_f64(&self) -> f64 {
        self.float.to_f64()
    }

    #[allow(dead_code)]
    fn is_nan_value(&self) -> bool {
        self.float.is_nan()
    }

    #[allow(dead_code)]
    fn is_infinite_value(&self) -> bool {
        self.float.is_infinite()
    }

    #[allow(dead_code)]
    fn floor_value(&self) -> Self {
        self.float.clone().floor().into()
    }
}

impl QuadFloat {
    fn pow_u64_via_mul(&self, mut exponent: u64) -> Self {
        let mut base = *self;
        let mut acc = Self::one();

        while exponent > 0 {
            if exponent & 1 == 1 {
                acc *= base;
            }
            exponent >>= 1;
            if exponent > 0 {
                base *= base;
            }
        }

        acc
    }

    fn machine_epsilon(&self) -> Self {
        self.from_i64(2)
            .pow((self.get_precision() - 1) as u64)
            .inv()
    }

    #[allow(dead_code)]
    fn one() -> Self {
        Self(DoubleFloat::from(1.0))
    }

    #[allow(non_snake_case)]
    fn E() -> Self {
        Self(DoubleFloat::default().e())
    }

    #[allow(non_snake_case)]
    fn PIHALF() -> Self {
        Self::PI() / Self::from_f64(2.0)
    }

    #[allow(non_snake_case)]
    fn PI() -> Self {
        Self(DoubleFloat::default().pi())
    }

    #[allow(non_snake_case)]
    fn TAU() -> Self {
        Self::PI() + Self::PI()
    }

    #[allow(non_snake_case)]
    fn FRAC_1_PI() -> Self {
        Self::PI().inv()
    }

    #[allow(dead_code)]
    pub(crate) fn from_f64_exact_binary(x: f64) -> Self {
        Self(DoubleFloat::from(x))
    }

    pub(crate) fn from_f64_decimal(x: f64) -> Self {
        if !x.is_finite() {
            return Self::from_f64_exact_binary(x);
        }

        // Keep decimal promotion (not exact-binary embedding), but avoid the
        // parser's locale lock for repeated settings/model values on each worker.
        // Bit keys distinguish signed zero and need no model/settings invalidation.
        thread_local! {
            static CACHE: RefCell<HashMap<u64, QuadFloat>> = RefCell::default();
        }
        CACHE.with_borrow_mut(|cache| {
            let key = x.to_bits();
            if let Some(value) = cache.get(&key) {
                return *value;
            }
            let value = SymbolicaFloat::parse(
                &format!("{}", x),
                Some(DoubleFloat::default().get_precision()),
            )
            .unwrap()
            .into();
            if cache.len() == 1024 {
                cache.clear();
            }
            cache.insert(key, value);
            value
        })
    }

    fn from_f64(x: f64) -> Self {
        Self::from_f64_decimal(x)
    }

    fn floor_via_symbolica(&self) -> Self {
        SymbolicaFloat::from(self).into_raw().floor().into()
    }

    fn trunc_toward_zero(&self) -> Self {
        if *self < Self::new_zero() {
            -(-self).floor_via_symbolica()
        } else {
            self.floor_via_symbolica()
        }
    }

    fn truncating_remainder(&self, rhs: &Self) -> Self {
        if rhs.is_zero() {
            return Self::from_f64_exact_binary(f64::NAN);
        }

        *self - ((*self / *rhs).trunc_toward_zero() * *rhs)
    }

    fn as_f64(&self) -> f64 {
        self.0.to_f64()
    }

    fn is_nan_value(&self) -> bool {
        self.0.to_f64().is_nan()
    }

    fn is_infinite_value(&self) -> bool {
        self.0.to_f64().is_infinite()
    }

    fn floor_value(&self) -> Self {
        self.floor_via_symbolica()
    }
}

impl<const N: u32> Default for VarFloat<N> {
    fn default() -> Self {
        VarFloat {
            float: rug::Float::with_val(N, 0.0),
        }
    }
}

impl PrecisionUpgradable for f128 {
    type Higher = ArbPrec;
    type Lower = f64;

    fn higher(&self) -> Self::Higher {
        ArbPrec::from(SymbolicaFloat::from(self))
    }

    fn lower(&self) -> Self::Lower {
        self.as_f64()
    }
}

impl<T: Real + PrecisionUpgradable, H: Real, L: Real> PrecisionUpgradable for Complex<T>
where
    T: PrecisionUpgradable<Higher = H, Lower = L>,
{
    type Higher = Complex<H>;
    type Lower = Complex<L>;

    fn higher(&self) -> Self::Higher {
        Complex::new(self.re.higher(), self.im.higher())
    }

    fn lower(&self) -> Self::Lower {
        Complex::new(self.re.lower(), self.im.lower())
    }
}

// #[allow(non_camel_case_types)]
// pub type f256 = VarFloat<243>;

pub trait PrecisionUpgradable {
    type Higher;
    type Lower;

    fn higher(&self) -> Self::Higher;
    fn lower(&self) -> Self::Lower;
}

pub trait FloatLike:
    Real
    + Send
    + Sync
    + 'static
    +R
    +Default
    + Clone
    + PartialOrd
    + RealLike
    + for<'a> RefAdd<&'a Self, Output = Self>
    // + for<'a> RefMutAdd<&'a Self, Output = Self>
    + RefAdd<Self, Output = Self>
    // + RefMutAdd<Self, Output = Self>
    + for<'a> RefMul<&'a Self, Output = Self>
    // + for<'a> RefMutMul<&'a Self, Output = Self>
    + RefMul<Self, Output = Self>
    // + RefMutMul<Self, Output = Self>
    + for<'a> RefSub<&'a Self, Output = Self>
    // + for<'a> RefMutSub<&'a Self, Output = Self>
    + RefSub<Self, Output = Self>
    // + RefMutSub<Self, Output = Self>
    + for<'a> RefDiv<&'a Self, Output = Self>
    // + for<'a> RefMutDiv<&'a Self, Output = Self> f64 doesn't have RefMutDiv
    + RefDiv<Self, Output = Self>
    + for<'a> RefRem<&'a Self, Output = Self>
    // + RefMutDiv<Self, Output = Self>
    + RefNeg<Output = Self>
    + RefZero
    + RefOne
    // + RefMutNeg<Output = Self> f64 doesn't have RefMutNeg
    + PrecisionUpgradable
    + Serialize
    + Display
    + GenericEvaluatorFloat
    + Into<symbolica::domains::float::Float>
{

    #[allow(non_snake_case)]
    fn PI(&self) -> Self;
    #[allow(non_snake_case)]
    fn E(&self) -> Self;
    #[allow(non_snake_case)]
    fn TAU(&self) -> Self;
    #[allow(non_snake_case)]
    fn SQRT_2(&self) -> Self;
    #[allow(non_snake_case)]
    fn SQRT_2_HALF(&self) -> Self;
    #[allow(non_snake_case)]
    fn PIHALF(&self) -> Self;
    #[allow(non_snake_case)]
    fn FRAC_1_PI(&self) -> Self;

    fn from_f64(x: f64) -> Self;

    /// Preserve the exact stored binary64 value of an original Monte Carlo
    /// coordinate. User-authored settings retain the existing decimal policy.
    fn from_f64_exact_binary(x: f64) -> Self;

    /// Round directly from the canonical Arb draw, never through another lane.
    /// Materialization uses the checked F::from_arb boundary below.
    fn from_arb(value: &ArbPrec) -> Self;

    /// Enclose the exact represented native value at the requested MPFR
    /// precision. Directed rounding includes both separated Quad limbs; this
    /// is a certificate boundary, unlike ordinary nearest-rounded conversion.
    fn mpfr_enclosure(&self, precision: u32) -> (Float, Float);

    /// Native sampling lane, distinct from the allowed physical stability levels.
    fn sampling_precision() -> SamplingPrecision;
    fn sampling_bridge_cache(setup: &LmbMultiChannelingSetup) -> &RuntimeCache<std::result::Result<SamplingChannelBridge<Self>, SamplingEvaluationError>>;
    fn sampling_bridge_cache_mut(setup: &mut LmbMultiChannelingSetup) -> &mut RuntimeCache<std::result::Result<SamplingChannelBridge<Self>, SamplingEvaluationError>>;


    #[allow(clippy::wrong_self_convention)]
    fn into_f64(&self) -> f64; // for inverse gamma in tropical sampling

    fn is_nan(&self) -> bool;

    fn is_infinite(&self) -> bool;

    fn floor(&self) -> Self;

    fn square(&self) -> Self {
        self.pow(2)
    }

    fn powi(&self, n: i32) -> Self {
        let absn = n.unsigned_abs() as u64;
        if n.is_negative() {
            self.pow(absn).inv()
        } else {
            self.pow(absn)
        }
    }

    fn epsilon(&self) -> Self;

    fn less_than_epsilon(&self) -> bool {
        self < &self.epsilon()
    }

    fn positive(&self) -> bool {
        self > &self.zero()
    }

    fn max_value(&self) -> Self {
        Self::from_f64(f64::MAX)
    }

    fn min_value(&self) -> Self {
        Self::from_f64(f64::MIN)
    }

    fn ln(&self) -> Self {
       self.log()
    }

    fn rem_euclid(&self, rhs: &Self) -> Self;

    fn try_extract_externals_from_cache(externals: &Externals) -> Option<&TiVec<ExternalIndex, FourMomentum<F<Self>>>>;
}

#[derive(
    Debug,
    Clone,
    PartialEq,
    PartialOrd,
    Copy,
    Default,
    Serialize,
    Deserialize,
    Encode,
    Decode,
    Hash,
    JsonSchema,
)]
pub struct F<T: FloatLike>(pub T);

impl<T: FloatLike> ToCoefficient for F<T> {
    fn to_coefficient(self) -> Coefficient {
        let z = self.zero();
        Coefficient::Float(symbolica::domains::float::Complex {
            re: self.into(),
            im: z.into(),
        })
    }
}

use symbolica::evaluate::{EvaluationDomain, ExportNumber};

#[cfg(test)]
mod evaluation_domain_tests;

impl ExportNumber for QuadFloat {
    fn export(&self) -> String {
        self.to_string()
    }

    fn is_real(&self) -> bool {
        true
    }

    fn to_complex_double(&self) -> symbolica::domains::float::Complex<f64> {
        symbolica::domains::float::Complex::new(self.0.to_f64(), 0.0)
    }
}

impl<const N: u32> ExportNumber for VarFloat<N> {
    fn export(&self) -> String {
        self.to_string()
    }

    fn is_real(&self) -> bool {
        true
    }

    fn to_complex_double(&self) -> symbolica::domains::float::Complex<f64> {
        symbolica::domains::float::Complex::new(self.to_f64(), 0.0)
    }
}

impl FixedPrecision for QuadFloat {
    const BINARY_PRECISION: usize = <DoubleFloat as FixedPrecision>::BINARY_PRECISION;
}

impl<const N: u32> FixedPrecision for VarFloat<N> {
    const BINARY_PRECISION: usize = N as usize;
}

impl EvaluationDomain for QuadFloat {
    const FIXED_PRECISION: Option<u32> = <DoubleFloat as EvaluationDomain>::FIXED_PRECISION;

    fn resolve_function(
        tags: &[AtomView],
        info: &EvaluationInfo,
    ) -> Option<Box<dyn ExternalFunction<Self>>> {
        if let Some(f) = info.get_evaluator::<Self>(tags) {
            return Some(f);
        }

        let f = DoubleFloat::resolve_function(tags, info)?;
        Some(Box::new(move |args: &[Self]| {
            let args: SmallVec<[DoubleFloat; 4]> = args.iter().map(|x| x.0).collect();
            Self(f(&args))
        }))
    }

    fn try_from_complex_float(
        f: symbolica::domains::float::Complex<symbolica::domains::float::Float>,
    ) -> Result<Self, String> {
        if f.is_real() {
            Ok(Self(f.re.to_double_float()))
        } else {
            Err(format!(
                "Cannot convert from Complex<Float> to {} because the result is not real",
                std::any::type_name::<Self>()
            ))
        }
    }
}

impl<const N: u32> EvaluationDomain for VarFloat<N> {
    const FIXED_PRECISION: Option<u32> = Some(N);

    fn resolve_function(
        tags: &[AtomView],
        info: &EvaluationInfo,
    ) -> Option<Box<dyn ExternalFunction<Self>>> {
        if let Some(f) = info.get_evaluator::<Self>(tags) {
            return Some(f);
        }

        let f = SymbolicaFloat::resolve_function(tags, info)?;
        Some(Box::new(move |args: &[Self]| {
            let args: SmallVec<[SymbolicaFloat; 4]> = args.iter().map(Into::into).collect();
            Self::from(f(&args))
        }))
    }

    fn try_from_complex_float(
        f: symbolica::domains::float::Complex<symbolica::domains::float::Float>,
    ) -> Result<Self, String> {
        if f.is_real() {
            Ok(Self::from(f.re))
        } else {
            Err(format!(
                "Cannot convert from Complex<Float> to {} because the result is not real",
                std::any::type_name::<Self>()
            ))
        }
    }
}

impl<T: FloatLike + ExportNumber> ExportNumber for F<T> {
    fn export(&self) -> String {
        self.0.to_string()
    }

    fn is_real(&self) -> bool {
        self.0.is_real()
    }

    fn to_complex_double(&self) -> symbolica::domains::float::Complex<f64> {
        self.0.to_complex_double()
    }
}

impl<T: FloatLike + FixedPrecision> FixedPrecision for F<T> {
    const BINARY_PRECISION: usize = T::BINARY_PRECISION;
}

impl<T: FloatLike + EvaluationDomain> EvaluationDomain for F<T> {
    const FIXED_PRECISION: Option<u32> = T::FIXED_PRECISION;

    fn resolve_function(
        tags: &[AtomView],
        info: &EvaluationInfo,
    ) -> Option<Box<dyn ExternalFunction<Self>>> {
        if let Some(f) = info.get_evaluator::<Self>(tags) {
            return Some(f);
        }

        let f = T::resolve_function(tags, info)?;
        Some(Box::new(move |args: &[Self]| {
            let args: SmallVec<[T; 4]> = args.iter().map(|x| x.0.clone()).collect();
            Self(f(&args))
        }))
    }

    fn try_from_complex_float(
        f: symbolica::domains::float::Complex<symbolica::domains::float::Float>,
    ) -> Result<Self, String> {
        if f.is_real() {
            T::try_from_complex_float(f).map(F)
        } else {
            Err(format!(
                "Cannot convert from Complex<Float> to {} because the result is not real",
                std::any::type_name::<Self>()
            ))
        }
    }
}

impl<T: TensorLibraryData + FloatLike> TensorLibraryData for F<T> {
    fn one() -> Self {
        F(<T as TensorLibraryData>::one())
    }

    fn minus_one() -> Self {
        F(<T as TensorLibraryData>::minus_one())
    }

    fn zero() -> Self {
        F(<T as TensorLibraryData>::zero())
    }
}

impl From<F<f64>> for Coefficient {
    fn from(value: F<f64>) -> Self {
        Coefficient::from(value.0)
    }
}

impl<'a> From<&'a F<f64>> for Coefficient {
    fn from(x: &'a F<f64>) -> Self {
        x.0.into()
    }
}

impl ToFloat for F<f64> {
    fn to_float(&self) -> symbolica::domains::float::Float {
        symbolica::domains::float::Float::with_val(53, self.0)
    }
}

impl ToAtom for F<f64> {
    fn to_atom(self) -> Atom {
        Atom::num(self.0)
    }
}

pub trait ToCoefficient {
    fn to_coefficient(self) -> Coefficient;
}

impl<T: FloatLike> From<F<T>> for symbolica::domains::float::Float {
    fn from(value: F<T>) -> Self {
        value.0.into()
    }
}

impl<T: FloatLike> ToCoefficient for Complex<F<T>> {
    fn to_coefficient(self) -> Coefficient {
        Coefficient::Float(symbolica::domains::float::Complex {
            re: self.re.0.into(),
            im: self.im.0.into(),
        })
    }
}

impl TrySmallestUpgrade<F<f64>> for Atom {
    type LCM = Atom;

    fn try_upgrade(&'_ self) -> Option<std::borrow::Cow<'_, Self::LCM>> {
        Some(std::borrow::Cow::Borrowed(self))
    }
}

impl TrySmallestUpgrade<Atom> for F<f64> {
    type LCM = Atom;

    fn try_upgrade(&'_ self) -> Option<std::borrow::Cow<'_, Self::LCM>> {
        <f64 as TrySmallestUpgrade<Atom>>::try_upgrade(&self.0)
    }
}

impl<T: FloatLike> R for F<T> {}

impl<T: FloatLike> Rem<&F<T>> for &F<T> {
    type Output = F<T>;

    fn rem(self, rhs: &F<T>) -> Self::Output {
        F(self.0.ref_rem(&rhs.0))
    }
}

impl<T: FloatLike> RefZero<F<T>> for &F<T> {
    fn ref_zero(&self) -> F<T> {
        F(self.0.ref_zero())
    }
}

impl<T: FloatLike> RealLike for F<T> {
    delegate! {
        to self.0{
            fn to_usize_clamped(&self)->usize;
            fn to_f64(&self)->f64;
            fn round_to_nearest_integer(&self)->Integer;
        }
    }
}

impl<T: FloatLike> SingleFloat for F<T> {
    delegate! {
        to self.0{
            fn is_zero(&self)->bool;
            fn is_one(&self)->bool;
            fn is_finite(&self)->bool;
        }
    }
    fn from_rational(&self, rat: &Rational) -> Self {
        F(self.0.from_rational(rat))
    }
}

impl<T: FloatLike> PrecisionUpgradable for F<T>
where
    T::Higher: FloatLike,
    T::Lower: FloatLike,
{
    type Higher = F<T::Higher>;
    type Lower = F<T::Lower>;

    fn higher(&self) -> Self::Higher {
        F(self.0.higher())
    }

    fn lower(&self) -> Self::Lower {
        F(self.0.lower())
    }
}

impl<T: FloatLike> RefZero for F<T> {
    fn ref_zero(&self) -> Self {
        F(self.0.zero())
    }
}

impl<T: FloatLike> RefOne for F<T> {
    fn ref_one(&self) -> Self {
        F(self.0.one())
    }
}

impl<T: FloatLike> TrySmallestUpgrade<F<T>> for F<T> {
    type LCM = F<T>;
    fn try_upgrade(&'_ self) -> Option<std::borrow::Cow<'_, Self::LCM>> {
        Some(std::borrow::Cow::Borrowed(self))
    }
}

impl<T: FloatLike> TrySmallestUpgrade<F<T>> for Complex<F<T>> {
    type LCM = Complex<F<T>>;
    fn try_upgrade(&'_ self) -> Option<std::borrow::Cow<'_, Self::LCM>> {
        Some(std::borrow::Cow::Borrowed(self))
    }
}

impl<T: FloatLike> TrySmallestUpgrade<Complex<F<T>>> for F<T> {
    type LCM = Complex<F<T>>;
    fn try_upgrade(&'_ self) -> Option<std::borrow::Cow<'_, Self::LCM>> {
        let z = self.ref_zero();

        Some(std::borrow::Cow::Owned(Complex::new(self.clone(), z)))
    }
}

// impl<T:FloatLike> TrySmallestUpgrade<Complex<F<T>>> for F<T> {
//     type LCM = Complex<F<T>>;
//     fn try_upgrade(&self) -> Option<std::borrow::Cow<Self::LCM>> {
//         Some(std::borrow::Cow::Borrowed(self))
//     }
// }

impl<'a, T: FloatLike> From<&'a Rational> for F<T> {
    fn from(x: &'a Rational) -> Self {
        F(T::new_zero().from_rational(x))
    }
}

impl<T: FloatLike> From<T> for F<T> {
    fn from(x: T) -> Self {
        F(x)
    }
}

impl<T: FloatLike> std::fmt::Display for F<T> {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(f, "{}", self.0)
    }
}

impl<T: FloatLike> std::fmt::LowerExp for F<T> {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        let rendered = match f.precision() {
            Some(precision) => format!("{:.*e}", precision, self.0),
            None => format!("{:e}", self.0),
        };

        let (is_nonnegative, body) = match rendered.strip_prefix('-') {
            Some(body) => (false, body),
            None => (true, rendered.as_str()),
        };

        f.pad_integral(is_nonnegative, "", body)
    }
}

impl<T: FloatLike> SymFloatLike for F<T> {
    #[inline]
    fn set_from(&mut self, other: &Self) {
        self.0.set_from(&other.0);
    }

    fn mul_add(&self, a: &Self, b: &Self) -> Self {
        F(self.0.mul_add(&a.0, &b.0))
    }

    fn is_fully_zero(&self) -> bool {
        self.0.is_fully_zero()
    }

    fn new_zero() -> Self {
        F(T::new_zero())
    }
    fn sample_unit<R: Rng + ?Sized>(&self, rng: &mut R) -> Self {
        F(self.0.sample_unit(rng))
    }
    fn neg(&self) -> Self {
        F(self.0.ref_neg())
    }

    fn zero(&self) -> Self {
        F(self.0.zero())
    }
    fn one(&self) -> Self {
        F(self.0.one())
    }
    // fn norm(&self) -> Self;
    fn from_usize(&self, x: usize) -> Self {
        F(self.0.from_usize(x))
    }
    fn from_i64(&self, x: i64) -> Self {
        F(self.0.from_i64(x))
    }
    fn pow(&self, n: u64) -> Self {
        F(self.0.pow(n))
    }
    fn inv(&self) -> Self {
        F(self.0.inv())
    }
    fn get_precision(&self) -> u32 {
        self.0.get_precision()
    }

    fn set_precision(&mut self, precision: u32) {
        self.0.set_precision(precision);
    }
    fn get_epsilon(&self) -> f64 {
        self.0.get_epsilon()
    }
    fn fixed_precision(&self) -> bool {
        self.0.fixed_precision()
    }
}

impl<T: FloatLike + Constructible> Constructible for F<T> {
    fn new_from_i64(a: i64) -> Self {
        F(T::new_from_i64(a))
    }

    fn new_from_usize(a: usize) -> Self {
        F(T::new_from_usize(a))
    }

    fn new_one() -> Self {
        F(T::new_one())
    }

    fn new_sample_unit<R: Rng + ?Sized>(rng: &mut R) -> Self {
        F(T::new_sample_unit(rng))
    }
}

impl<T: FloatLike> Real for F<T> {
    fn atan2(&self, x: &Self) -> Self {
        F(self.0.atan2(&x.0))
    }

    fn conj(&self) -> Self {
        F(self.0.conj())
    }

    fn i(&self) -> Option<Self> {
        None
    }

    fn powf(&self, e: &Self) -> Self {
        F(self.0.powf(&e.0))
    }

    fn norm(&self) -> Self {
        F(self.0.norm())
    }

    fn e(&self) -> Self {
        F(self.0.e())
    }
    fn phi(&self) -> Self {
        F(self.0.phi())
    }
    fn euler(&self) -> Self {
        F(self.0.euler())
    }
    fn pi(&self) -> Self {
        F(self.0.pi())
    }
    fn sqrt(&self) -> Self {
        F(self.0.sqrt())
    }
    fn log(&self) -> Self {
        F(self.0.log())
    }
    fn exp(&self) -> Self {
        F(self.0.exp())
    }
    fn sin(&self) -> Self {
        F(self.0.sin())
    }
    fn cos(&self) -> Self {
        F(self.0.cos())
    }
    fn tan(&self) -> Self {
        F(self.0.tan())
    }
    fn asin(&self) -> Self {
        F(self.0.asin())
    }
    fn acos(&self) -> Self {
        F(self.0.acos())
    }
    fn sinh(&self) -> Self {
        F(self.0.sinh())
    }
    fn cosh(&self) -> Self {
        F(self.0.cosh())
    }
    fn tanh(&self) -> Self {
        F(self.0.tanh())
    }
    fn asinh(&self) -> Self {
        F(self.0.asinh())
    }
    fn acosh(&self) -> Self {
        F(self.0.acosh())
    }
    fn atanh(&self) -> Self {
        F(self.0.atanh())
    }
}

use delegate::delegate;

impl<T: FloatLike> F<T> {
    pub(crate) fn max(self, other: F<T>) -> F<T> {
        if self < other { other } else { self }
    }

    pub fn negate(&mut self) {
        self.0 = -self.0.clone();
    }

    pub(crate) fn from_ff64(x: F<f64>) -> Self {
        F(T::from_f64(x.0))
    }

    pub(crate) fn zero(&self) -> Self {
        F(self.0.zero())
    }

    pub(crate) fn one(&self) -> Self {
        F(self.0.one())
    }

    #[allow(clippy::wrong_self_convention)]
    pub(crate) fn from_usize(&self, x: usize) -> Self {
        F(self.0.from_usize(x))
    }

    #[allow(clippy::wrong_self_convention)]
    pub(crate) fn from_i64(&self, x: i64) -> Self {
        F(self.0.from_i64(x))
    }

    pub(crate) fn higher(&self) -> F<T::Higher>
    where
        T::Higher: FloatLike,
    {
        F(self.0.higher())
    }

    pub(crate) fn lower(&self) -> F<T::Lower>
    where
        T::Lower: FloatLike,
    {
        F(self.0.lower())
    }

    pub fn from_f64(x: f64) -> Self {
        F(T::from_f64(x))
    }

    /// Embed the exact represented native value in the canonical Arb container.
    /// In particular, two separated Quad limbs need not fit in 1000 bits.
    pub(crate) fn to_arb_exact(&self) -> eyre::Result<F<ArbPrec>> {
        let (lower, upper) = self.0.mpfr_enclosure(1000);
        if !lower.is_finite() || !upper.is_finite() || lower != upper {
            return Err(SamplingEvaluationError::Unrepresentable {
                operation: "canonical source embedding",
                detail: format!(
                    "{:?} scalar {self} is not exactly representable at 1000 bits",
                    T::sampling_precision()
                ),
            }
            .into());
        }
        // Opposite signed zeros compare equal in the directed enclosure. At
        // exact zero only, the binary64 sign extraction is itself lossless and
        // preserves the native signed-zero convention; nonzero values never narrow.
        let exact = if lower == 0 {
            Float::with_val(1000, self.0.into_f64())
        } else {
            lower
        };
        Ok(F(ArbPrec::from(exact)))
    }

    /// Materialize one canonical scalar without accepting lost range as zero.
    /// Relative factor accuracy and physical positivity belong to their owners.
    pub(crate) fn from_arb(value: &ArbPrec) -> eyre::Result<Self> {
        if !value.float.is_finite() {
            return Err(SamplingEvaluationError::Unrepresentable {
                operation: "canonical draw materialization",
                detail: "nonfinite canonical scalar".into(),
            }
            .into());
        }
        let native = T::from_arb(value);
        if !native.is_finite() || (native == native.zero() && value.float != 0) {
            return Err(SamplingEvaluationError::Unrepresentable {
                operation: "canonical draw materialization",
                detail: format!(
                    "canonical scalar {value} is outside {:?} range",
                    T::sampling_precision()
                ),
            }
            .into());
        }
        Ok(F(native))
    }

    /// Bound materialization error against the original represented scalar.
    /// Sign/positivity requirements of a weight or root remain with its owner.
    pub(crate) fn verify_arb_materialization(
        &self,
        source: &ArbPrec,
        tolerance: f64,
    ) -> eyre::Result<()> {
        const PRECISION: u32 = 2048;
        if !tolerance.is_finite() || tolerance <= 0.0 {
            return Err(eyre::eyre!(
                "materialization requires a finite positive relative budget"
            ));
        }
        if !self.0.is_finite() || !source.float.is_finite() {
            return Err(SamplingEvaluationError::Unrepresentable {
                operation: "canonical draw materialization",
                detail: "nonfinite source or native scalar in relative check".into(),
            }
            .into());
        }
        let (source_lo, source_hi) = source.mpfr_enclosure(PRECISION);
        let (native_lo, native_hi) = self.0.mpfr_enclosure(PRECISION);
        let error = Float::with_val_round(PRECISION, &native_hi - &source_lo, Round::Up)
            .0
            .max(&Float::with_val_round(PRECISION, &source_hi - &native_lo, Round::Up).0);
        let magnitude = if source_lo >= 0 {
            source_lo
        } else if source_hi <= 0 {
            -source_hi
        } else {
            Float::with_val(PRECISION, 0)
        };
        // Borrowed operands keep arithmetic unevaluated until directed rounding.
        let allowed = Float::with_val_round(PRECISION, &magnitude * tolerance, Round::Down).0;
        if !error.is_finite() || !allowed.is_finite() || error > allowed {
            return Err(SamplingEvaluationError::UncertainGeometry {
                detail: format!("canonical scalar {source} materialized as {self} exceeds relative budget {tolerance} in {:?}", T::sampling_precision()),
            }.into());
        }
        Ok(())
    }

    #[allow(clippy::wrong_self_convention)]
    pub(crate) fn into_ff64(&self) -> F<f64> {
        F(self.0.into_f64())
    }

    pub(crate) fn abs(&self) -> Self {
        F(self.0.norm())
    }

    /// Euclidean magnitude without squaring unscaled native components.
    pub(crate) fn hypot(&self, other: &Self) -> Self {
        if !self.0.is_finite() || !other.0.is_finite() {
            return self.abs() + other.abs();
        }
        let scale = self.abs().max(other.abs());
        if scale == scale.zero() {
            return scale;
        }
        ((self / &scale).square() + (other / &scale).square()).sqrt() * scale
    }

    pub(crate) fn sqrt(&self) -> Self {
        F(self.0.sqrt())
    }

    pub(crate) fn powf(&self, e: &Self) -> Self {
        F(self.0.powf(&e.0))
    }

    pub(crate) fn log10(&self) -> Self {
        let ten = self.from_i64(10);
        self.ln() / ten.ln()
    }

    pub(crate) fn complex_sqrt(&self) -> Complex<Self> {
        if self.positive() {
            Complex::new(self.sqrt(), self.zero())
        } else {
            Complex::new(self.zero(), (-self).sqrt())
        }
    }

    pub(crate) fn rem_euclid(&self, rhs: &Self) -> Self {
        F(self.0.rem_euclid(&rhs.0))
    }

    #[allow(non_snake_case)]
    pub(crate) fn PI(&self) -> Self {
        F(self.0.PI())
    }
    #[allow(non_snake_case)]
    pub fn E(&self) -> Self {
        F(self.0.E())
    }
    #[allow(non_snake_case)]
    pub fn TAU(&self) -> Self {
        F(self.0.TAU())
    }
    #[allow(non_snake_case)]
    pub fn PIHALF(&self) -> Self {
        F(self.0.PIHALF())
    }
    #[allow(non_snake_case)]
    pub fn SQRT_2(&self) -> Self {
        F(self.0.SQRT_2())
    }
    #[allow(non_snake_case)]
    pub fn SQRT_2_HALF(&self) -> Self {
        F(self.0.SQRT_2_HALF())
    }
    #[allow(non_snake_case)]
    pub(crate) fn FRAC_1_PI(&self) -> Self {
        F(self.0.FRAC_1_PI())
    }
    #[allow(clippy::wrong_self_convention)]
    pub(crate) fn into_f64(&self) -> f64 {
        self.0.into_f64()
    }
    pub(crate) fn square(&self) -> Self {
        F(self.0.square())
    }
    pub(crate) fn powi(&self, n: i32) -> Self {
        F(self.0.powi(n))
    }
    pub(crate) fn epsilon(&self) -> Self {
        F(self.0.epsilon())
    }
    pub(crate) fn less_than_epsilon(&self) -> bool {
        self.0.less_than_epsilon()
    }
    pub(crate) fn positive(&self) -> bool {
        self.0.positive()
    }
    pub(crate) fn max_value(&self) -> Self {
        F(self.0.max_value())
    }
    pub(crate) fn min_value(&self) -> Self {
        F(self.0.min_value())
    }
    pub(crate) fn ln(&self) -> Self {
        F(self.0.ln())
    }
    pub(crate) fn inv(&self) -> Self {
        F(self.0.inv())
    }
    pub(crate) fn is_nan(&self) -> bool {
        self.0.is_nan()
    }
    pub(crate) fn is_infinite(&self) -> bool {
        self.0.is_infinite()
    }
    pub fn floor(&self) -> Self {
        F(self.0.floor())
    }
}

// impl CompiledEvaluatorFloat for F<f64> {
//     fn evaluate(
//         eval: &mut symbolica::evaluate::CompiledEvaluator,
//         args: &[Self],
//         out: &mut [Self],
//     ) {
//         // cast to f64
//         let args_f64: Vec<f64> = args.iter().map(|x| x.0).collect_vec();
//         let mut out_f64 = out.iter().map(|x| x.0).collect_vec();

//         eval.evaluate_double(&args_f64, &mut out_f64);

//         // write the result to out
//         out.iter_mut()
//             .zip(out_f64)
//             .for_each(|(out_ff64, out_f64)| *out_ff64 = F(out_f64));
//     }
// }

impl<T: FloatLike> Add<F<T>> for F<T> {
    type Output = F<T>;
    fn add(self, rhs: F<T>) -> Self::Output {
        F(self.0 + rhs.0)
    }
}

impl<T: FloatLike> Add<&F<T>> for F<T> {
    type Output = F<T>;
    fn add(self, rhs: &F<T>) -> Self::Output {
        F(self.0 + &rhs.0)
    }
}

impl<T: FloatLike> Add<&F<T>> for &F<T> {
    type Output = F<T>;
    fn add(self, rhs: &F<T>) -> Self::Output {
        F(self.0.ref_add(&rhs.0))
    }
}

impl<T: FloatLike> Add<F<T>> for &F<T> {
    type Output = F<T>;
    fn add(self, rhs: F<T>) -> Self::Output {
        F(self.0.ref_add(rhs.0))
    }
}

impl<T: FloatLike> AddAssign<&F<T>> for F<T> {
    fn add_assign(&mut self, rhs: &F<T>) {
        self.0 += &rhs.0;
    }
}

impl<T: FloatLike> AddAssign<F<T>> for F<T> {
    fn add_assign(&mut self, rhs: F<T>) {
        self.0 += rhs.0;
    }
}

impl<T: FloatLike> Sub<F<T>> for F<T> {
    type Output = F<T>;
    fn sub(self, rhs: F<T>) -> Self::Output {
        F(self.0 - rhs.0)
    }
}

impl<T: FloatLike> Sub<&F<T>> for F<T> {
    type Output = F<T>;
    fn sub(self, rhs: &F<T>) -> Self::Output {
        F(self.0 - &rhs.0)
    }
}

impl<T: FloatLike> Sub<&F<T>> for &F<T> {
    type Output = F<T>;
    fn sub(self, rhs: &F<T>) -> Self::Output {
        F(self.0.ref_sub(&rhs.0))
    }
}

impl<T: FloatLike> Sub<F<T>> for &F<T> {
    type Output = F<T>;
    fn sub(self, rhs: F<T>) -> Self::Output {
        F(self.0.ref_sub(rhs.0))
    }
}

impl<T: FloatLike> SubAssign<&F<T>> for F<T> {
    fn sub_assign(&mut self, rhs: &F<T>) {
        self.0 -= &rhs.0;
    }
}

impl<T: FloatLike> SubAssign<F<T>> for F<T> {
    fn sub_assign(&mut self, rhs: F<T>) {
        self.0 -= rhs.0;
    }
}

impl<T: FloatLike> Mul<F<T>> for F<T> {
    type Output = F<T>;
    fn mul(self, rhs: F<T>) -> Self::Output {
        F(self.0 * rhs.0)
    }
}

impl<T: FloatLike> Mul<&F<T>> for F<T> {
    type Output = F<T>;
    fn mul(self, rhs: &F<T>) -> Self::Output {
        F(self.0 * &rhs.0)
    }
}

impl<T: FloatLike> Mul<&F<T>> for &F<T> {
    type Output = F<T>;
    fn mul(self, rhs: &F<T>) -> Self::Output {
        F(self.0.ref_mul(&rhs.0))
    }
}

impl<T: FloatLike> Mul<F<T>> for &F<T> {
    type Output = F<T>;
    fn mul(self, rhs: F<T>) -> Self::Output {
        F(self.0.ref_mul(rhs.0))
    }
}

impl<T: FloatLike> MulAssign<&F<T>> for F<T> {
    fn mul_assign(&mut self, rhs: &F<T>) {
        self.0 *= &rhs.0;
    }
}

impl<T: FloatLike> MulAssign<F<T>> for F<T> {
    fn mul_assign(&mut self, rhs: F<T>) {
        self.0 *= rhs.0;
    }
}

impl<T: FloatLike> Div<F<T>> for F<T> {
    type Output = F<T>;
    fn div(self, rhs: F<T>) -> Self::Output {
        F(self.0 / rhs.0)
    }
}

impl<T: FloatLike> Div<&F<T>> for F<T> {
    type Output = F<T>;
    fn div(self, rhs: &F<T>) -> Self::Output {
        F(self.0 / &rhs.0)
    }
}

impl<T: FloatLike> Div<&F<T>> for &F<T> {
    type Output = F<T>;
    fn div(self, rhs: &F<T>) -> Self::Output {
        F(self.0.ref_div(&rhs.0))
    }
}

impl<T: FloatLike> Div<F<T>> for &F<T> {
    type Output = F<T>;
    fn div(self, rhs: F<T>) -> Self::Output {
        F(self.0.ref_div(rhs.0))
    }
}

impl<T: FloatLike> DivAssign<&F<T>> for F<T> {
    fn div_assign(&mut self, rhs: &F<T>) {
        self.0 /= &rhs.0;
    }
}

impl<T: FloatLike> DivAssign<F<T>> for F<T> {
    fn div_assign(&mut self, rhs: F<T>) {
        self.0 /= rhs.0;
    }
}

impl<T: FloatLike> Neg for F<T> {
    type Output = F<T>;
    fn neg(self) -> Self::Output {
        F(-self.0)
    }
}

impl<T: FloatLike> Neg for &F<T> {
    type Output = F<T>;
    fn neg(self) -> Self::Output {
        F(self.0.ref_neg())
    }
}

pub trait RefDefault {
    fn default(&self) -> Self;
}

impl<T: FloatLike> RefDefault for T {
    fn default(&self) -> Self {
        self.zero()
    }
}

impl<T: FloatLike> RefDefault for F<T> {
    fn default(&self) -> Self {
        F(self.0.default())
    }
}
impl PrecisionUpgradable for f64 {
    type Higher = f128;
    type Lower = f64;

    fn higher(&self) -> Self::Higher {
        f128::from_f64(*self)
    }

    fn lower(&self) -> Self::Lower {
        *self
    }
}

impl FloatLike for f64 {
    fn mpfr_enclosure(&self, precision: u32) -> (Float, Float) {
        (
            Float::with_val_round(precision, *self, Round::Down).0,
            Float::with_val_round(precision, *self, Round::Up).0,
        )
    }

    fn PI(&self) -> Self {
        std::f64::consts::PI
    }

    fn SQRT_2(&self) -> Self {
        std::f64::consts::SQRT_2
    }

    fn SQRT_2_HALF(&self) -> Self {
        std::f64::consts::SQRT_2 / 2.0
    }

    fn PIHALF(&self) -> Self {
        std::f64::consts::PI / 2.0
    }

    fn E(&self) -> Self {
        std::f64::consts::E
    }

    fn TAU(&self) -> Self {
        std::f64::consts::TAU
    }

    fn FRAC_1_PI(&self) -> Self {
        std::f64::consts::FRAC_1_PI
    }

    fn from_f64_exact_binary(x: f64) -> Self {
        x
    }

    fn from_arb(value: &ArbPrec) -> Self {
        value.float.to_f64()
    }

    fn sampling_precision() -> SamplingPrecision {
        SamplingPrecision::Double
    }

    fn sampling_bridge_cache(
        setup: &LmbMultiChannelingSetup,
    ) -> &RuntimeCache<std::result::Result<SamplingChannelBridge<Self>, SamplingEvaluationError>>
    {
        &setup.sampling_bridge
    }

    fn sampling_bridge_cache_mut(
        setup: &mut LmbMultiChannelingSetup,
    ) -> &mut RuntimeCache<std::result::Result<SamplingChannelBridge<Self>, SamplingEvaluationError>>
    {
        &mut setup.sampling_bridge
    }

    fn from_f64(x: f64) -> Self {
        x
    }

    fn into_f64(&self) -> f64 {
        *self
    }

    fn is_nan(&self) -> bool {
        f64::is_nan(*self)
    }

    fn is_infinite(&self) -> bool {
        f64::is_infinite(*self)
    }

    fn floor(&self) -> Self {
        f64::floor(*self)
    }

    fn rem_euclid(&self, rhs: &Self) -> Self {
        f64::rem_euclid(*self, *rhs)
    }

    fn try_extract_externals_from_cache(
        externals: &Externals,
    ) -> Option<&TiVec<ExternalIndex, FourMomentum<F<Self>>>> {
        match externals {
            Externals::Constant { f_64_cache, .. } => f_64_cache.as_ref(),
        }
    }

    fn epsilon(&self) -> Self {
        f64::EPSILON
    }
}
impl From<F<f64>> for f64 {
    fn from(value: F<f64>) -> Self {
        value.0
    }
}

impl From<F<f64>> for Rational {
    fn from(value: F<f64>) -> Self {
        value.0.try_into().unwrap()
    }
}

#[allow(non_camel_case_types)]
// Keep this one-line fallback close to the active alias: `VarFloat<113>` is a
// slower but safe binary128-like replacement when debugging `DoubleFloat`
// issues. The effective epsilon follows the scalar, but the explicit native
// evaluator/parameter implementations must also be reviewed when flipping it.
// pub type f128 = VarFloat<113>;
pub type f128 = QuadFloat;
pub type ArbPrec = VarFloat<1000>;
/// Fixed source arithmetic; completed values embed exactly in the Arb container.
pub type SamplingFloat = VarFloat<256>;

/// Numerical sampling lanes include a source-only precision which never becomes
/// a physical stability level or a persisted runtime setting.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum SamplingPrecision {
    Double,
    Quad,
    Fixed256,
    Arb,
}

impl From<crate::settings::runtime::Precision> for SamplingPrecision {
    fn from(precision: crate::settings::runtime::Precision) -> Self {
        use crate::settings::runtime::Precision;
        match precision {
            Precision::Double => Self::Double,
            Precision::Quad => Self::Quad,
            Precision::Arb => Self::Arb,
        }
    }
}

impl Display for SamplingPrecision {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.write_str(match self {
            Self::Double => "f64",
            Self::Quad => "f128",
            Self::Fixed256 => "fixed256",
            Self::Arb => "arb1000",
        })
    }
}

/// An iterator which iterates two other iterators simultaneously
#[derive(Clone, Debug)]
#[must_use = "iterator adaptors are lazy and do nothing unless consumed"]
pub struct ZipEq<I, J> {
    a: I,
    b: J,
}

/// An iterator which iterates two other iterators simultaneously and checks
/// if the sizes are equal in debug mode.
#[allow(unused)]
pub(crate) fn zip_eq<I, J>(i: I, j: J) -> ZipEq<I::IntoIter, J::IntoIter>
where
    I: IntoIterator,
    J: IntoIterator,
{
    ZipEq {
        a: i.into_iter(),
        b: j.into_iter(),
    }
}

impl<I, J> Iterator for ZipEq<I, J>
where
    I: Iterator,
    J: Iterator,
{
    type Item = (I::Item, J::Item);

    fn next(&mut self) -> Option<Self::Item> {
        match (self.a.next(), self.b.next()) {
            (None, None) => None,
            (Some(a), Some(b)) => Some((a, b)),
            (None, Some(_)) => {
                #[cfg(debug_assertions)]
                panic!("Unequal length of iterators; first iterator finished first");
                #[cfg(not(debug_assertions))]
                None
            }
            (Some(_), None) => {
                #[cfg(debug_assertions)]
                panic!("Unequal length of iterators; second iterator finished first");
                #[cfg(not(debug_assertions))]
                None
            }
        }
    }

    fn size_hint(&self) -> (usize, Option<usize>) {
        let sa = self.a.size_hint();
        let sb = self.b.size_hint();
        (sa.0.min(sb.0), sa.1.zip(sb.1).map(|(ua, ub)| ua.min(ub)))
    }
}

impl<I, J> ExactSizeIterator for ZipEq<I, J>
where
    I: ExactSizeIterator,
    J: ExactSizeIterator,
{
}

pub(crate) fn parse_python_expression(expression: &str) -> Atom {
    initialize();
    let _ = UFO.metric;
    let processed_string = String::from(expression)
        .replace("**", "^")
        .replace("cmath.sqrt", "sqrt")
        .replace("cmath.pi", "pi")
        .replace("math.sqrt", "sqrt")
        .replace("math.pi", "pi");

    parse!(processed_string)
}

#[test]
fn python_expression_imaginary_unit_is_numeric() {
    // The pinned Symbolica parser accepts both UFO export spellings directly.
    for expression in ["𝑖", "1𝑖", "(𝑖)^5", "1+𝑖-1"] {
        assert_eq!(parse_python_expression(expression), Atom::i());
    }
    assert_eq!(parse_python_expression("-𝑖/3"), -Atom::i() / 3);
}

/// Format a mean ± sdev as mean(sdev) with the correct number of digits.
/// Based on the Python package gvar.
pub fn format_uncertainty(mean: F<f64>, sdev: F<f64>) -> String {
    let mean = mean.0;
    let sdev = sdev.0;

    fn ndec(x: f64, offset: usize) -> i32 {
        let mut ans = (offset as f64 - x.log10()) as i32;
        if ans > 0 && x * 10.0.powi(ans) >= [0.5, 9.5, 99.5][offset] {
            ans -= 1;
        }
        if ans < 0 { 0 } else { ans }
    }
    let v = mean;
    let dv = sdev.abs();

    // special cases
    if v.is_nan() || dv.is_nan() {
        format!("{:e} ± {:e}", v, dv)
    } else if dv.is_infinite() {
        format!("{:e} ± inf", v)
    } else if v.is_zero() && !(1e-4..1e5).contains(&dv) {
        if dv.is_zero() {
            "0(0)".to_owned()
        } else {
            let e = format!("{:.1e}", dv);
            let mut ans = e.split('e');
            let e1 = ans.next().unwrap();
            let e2 = ans.next().unwrap();
            "0.0(".to_owned() + e1 + ")e" + e2
        }
    } else if v.is_zero() {
        if dv >= 9.95 {
            format!("0({:.0})", dv)
        } else if dv >= 0.995 {
            format!("0.0({:.1})", dv)
        } else {
            let ndecimal = ndec(dv, 2);
            format!(
                "{:.*}({:.0})",
                ndecimal as usize,
                v,
                dv * (10.).powi(ndecimal)
            )
        }
    } else if dv.is_zero() {
        let e = format!("{:e}", v);
        let mut ans = e.split('e');
        let e1 = ans.next().unwrap();
        let e2 = ans.next().unwrap();
        if e2 != "0" {
            e1.to_owned() + "(0)e" + e2
        } else {
            e1.to_owned() + "(0)"
        }
    } else if dv > 1e4 * v.abs() {
        format!("{:.1e} ± {:.2e}", v, dv)
    } else if v.abs() >= 1e6 || v.abs() < 1e-5 {
        // exponential notation for large |self.mean|
        let exponent = v.abs().log10().floor();
        let fac = (10.0).powf(&exponent);
        let mantissa = format_uncertainty(F(v / fac), F(dv / fac));
        let e = format!("{:.0e}", fac);
        let mut ee = e.split('e');
        mantissa + "e" + ee.nth(1).unwrap()
    }
    // normal cases
    else if dv >= 9.95 {
        if v.abs() >= 9.5 {
            format!("{:.0}({:.0})", v, dv)
        } else {
            let ndecimal = ndec(v.abs(), 1);
            format!("{:.*}({:.*})", ndecimal as usize, v, ndecimal as usize, dv)
        }
    } else if dv >= 0.995 {
        if v.abs() >= 0.95 {
            format!("{:.1}({:.1})", v, dv)
        } else {
            let ndecimal = ndec(v.abs(), 1);
            format!("{:.*}({:.*})", ndecimal as usize, v, ndecimal as usize, dv)
        }
    } else {
        let ndecimal = ndec(v.abs(), 1).max(ndec(dv, 2));
        format!(
            "{:.*}({:.0})",
            ndecimal as usize,
            v,
            dv * (10.).powi(ndecimal)
        )
    }
}

/// Compare two slices, selecting on length first
#[allow(unused)]
pub(crate) fn compare_slice<T: Ord>(slice1: &[T], slice2: &[T]) -> Ordering {
    match slice1.len().cmp(&slice2.len()) {
        Ordering::Equal => (),
        non_eq => return non_eq,
    }

    let l = slice1.len();
    // Slice to the loop iteration range to enable bound check
    // elimination in the compiler
    let lhs = &slice1[..l];
    let rhs = &slice2[..l];

    for i in 0..l {
        match lhs[i].cmp(&rhs[i]) {
            Ordering::Equal => (),
            non_eq => return non_eq,
        }
    }

    Ordering::Equal
}

pub trait Signum {
    fn multiply_sign(&self, sign: i8) -> Self;
}

// impl Signum for f128::f128 {
//     #[inline]
//     fn multiply_sign(&self, sign: i8) -> f128::f128 {
//         match sign {
//             1 => *self,
//             0 => f128::f128::zero(),
//             -1 => self.neg(),
//             _ => unreachable!("Sign should be -1,0,1"),
//         }
//     }
// }

impl Signum for f64 {
    #[inline]
    fn multiply_sign(&self, sign: i8) -> f64 {
        match sign {
            1 => *self,
            0 => self.zero(),
            -1 => -self,
            _ => unreachable!("Sign should be -1,0,1"),
        }
    }
}

impl Signum for f32 {
    #[inline]
    fn multiply_sign(&self, sign: i8) -> Self {
        match sign {
            1 => *self,
            0 => 0.0,
            -1 => self.neg(),
            _ => unreachable!("Sign should be -1,0,1"),
        }
    }
}

impl<T: FloatLike> Signum for Complex<F<T>> {
    #[inline]
    fn multiply_sign(&self, sign: i8) -> Complex<F<T>> {
        match sign {
            1 => self.clone(),
            0 => self.ref_zero(),
            -1 => -self.clone(),
            _ => unreachable!("Sign should be -1,0,1"),
        }
    }
}

impl<T: FloatLike> Signum for FourMomentum<F<T>> {
    #[inline]
    fn multiply_sign(&self, sign: i8) -> FourMomentum<F<T>> {
        match sign {
            1 => self.clone(),
            0 => self.ref_zero(),
            -1 => -self,
            _ => unreachable!("Sign should be -1,0,1"),
        }
    }
}

#[allow(unused)]
#[inline]
/// Invert with better precision
pub(crate) fn finv<T: FloatLike>(c: Complex<F<T>>) -> Complex<F<T>> {
    let norm = c.norm_squared();
    c.conj() / norm
}

#[allow(unused)]
#[inline]
pub(crate) fn powi<T: FloatLike>(c: Complex<F<T>>, n: i32) -> Complex<F<T>> {
    if n.is_negative() {
        let u = -n as u64;
        finv(c.pow(u))
    } else {
        let u = n as u64;
        c.pow(u)
    }
}

#[allow(unused)]
pub(crate) fn evaluate_signature<T>(
    signature: &[i8],
    momenta: &[FourMomentum<F<T>>],
) -> FourMomentum<F<T>>
where
    T: FloatLike,
{
    let mut momentum = momenta[0].zero();
    for (&sign, mom) in zip_eq(signature, momenta) {
        match sign {
            0 => {}
            1 => momentum += mom,
            -1 => momentum -= mom,
            _ => {
                #[cfg(debug_assertions)]
                panic!("Sign should be -1,0,1")
            }
        }
    }

    momentum
}

#[allow(unused)]
#[inline]
pub(crate) fn pinch_dampening_function<T: FloatLike>(
    dampening_arg: F<T>,
    delta_t: F<T>,
    powers: (u64, u64),
    multiplier: f64,
) -> F<T> {
    // Make sure the function is even in t-tstar
    assert!(powers.1.is_multiple_of(2));
    let a = dampening_arg.pow(powers.0);
    &a / (&a + F::<T>::from_f64(multiplier) * delta_t.pow(powers.1))
}

pub(crate) fn h_dual<T: FloatLike>(
    t: &HyperDual<F<T>>,
    tstar: Option<HyperDual<F<T>>>,
    sigma: Option<F<T>>,
    h_function_settings: &crate::settings::runtime::HFunctionSettings,
) -> HyperDual<F<T>> {
    let sigma = sigma.unwrap_or_else(|| F::<T>::from_f64(h_function_settings.sigma));
    let sig = new_constant(t, &sigma);
    let two = new_constant(t, &sigma.from_usize(2));
    let power = h_function_settings.power;
    let prefactor = || {
        let normalisation = new_constant(t, &h_function_settings.normalization(&sigma));
        match power {
            None | Some(0) => normalisation.inv(),
            Some(p) => (t.clone() / &sig).inv().pow(p as u64) / normalisation,
        }
    };
    match h_function_settings.function {
        crate::settings::runtime::HFunction::Exponential => {
            // This profile has no power.
            (-(t.clone() * t) / (sig.clone() * &sig)).exp()
                / new_constant(t, &h_function_settings.normalization(&sigma))
        }
        crate::settings::runtime::HFunction::PolyExponential => {
            // Normalized by the integral of (t/sigma)^{-p} exp(2-t^2/sigma^2-sigma^2/t^2).
            prefactor()
                * (two
                    - (t.clone() * t) / (sig.clone() * &sig)
                    - (sig.clone() * &sig) / (t.clone() * t))
                    .exp()
        }
        crate::settings::runtime::HFunction::PolyLeftRightExponential => {
            // Normalized by the integral of (t/sigma)^{-p} exp(-((t^2/sigma^2+1)/(t/sigma)-2)).
            prefactor()
                * (two - ((t.clone() * t) / (sig.clone() * &sig) + t.one()) / (t.clone() / &sig))
                    .exp()
        }
        crate::settings::runtime::HFunction::ExponentialCT => {
            let delta_t_sq = (tstar.clone().unwrap() - t) * (tstar.clone().unwrap() - t);
            let tstar_sq = tstar.clone().unwrap() * tstar.unwrap();
            // info!("dampener: {}", dampener);
            // info!("delta_t_sq: {}", delta_t_sq);
            // info!("tstar_sq: {}", tstar_sq);
            // info!(
            //     "Exp arg: {}",
            //     -sig.inv() * (delta_t_sq / tstar_sq + sig * sig * (dampener * dampener))
            // );
            // info!(
            //     "result: {}",
            //     (-sig.inv() * (delta_t_sq / tstar_sq + sig * sig * (dampener * dampener))).exp()
            // );
            if h_function_settings.enabled_dampening {
                let dampener = delta_t_sq.clone() / (delta_t_sq.clone() - &tstar_sq);
                (-sig.inv()
                    * (delta_t_sq.clone() / tstar_sq
                        + sig.clone() * sig * (dampener.clone() * dampener)))
                    .exp()
            } else {
                (-sig.inv() * (delta_t_sq / tstar_sq)).exp()
            }
        }
    }
}

pub(crate) fn h<T: FloatLike>(
    t: &F<T>,
    tstar: Option<F<T>>,
    sigma: Option<F<T>>,
    h_function_settings: &crate::settings::runtime::HFunctionSettings,
) -> F<T> {
    let sig = sigma.unwrap_or_else(|| F::<T>::from_f64(h_function_settings.sigma));
    let two = sig.from_usize(2);
    let power = h_function_settings.power;
    let prefactor = || {
        let normalisation = h_function_settings.normalization(&sig);
        match power {
            None | Some(0) => normalisation.inv(),
            Some(p) => (t / &sig).powi(-(p as i32)) / normalisation,
        }
    };
    match h_function_settings.function {
        crate::settings::runtime::HFunction::Exponential => {
            // This profile has no power.
            (-(t.square()) / (sig.square())).exp() / h_function_settings.normalization(&sig)
        }
        crate::settings::runtime::HFunction::PolyExponential => {
            // Normalized by the integral of (t/sigma)^{-p} exp(2-t^2/sigma^2-sigma^2/t^2).
            prefactor()
                * (two - (t.square()) / (sig.square()) - (sig.square()) / (t.square())).exp()
        }
        crate::settings::runtime::HFunction::PolyLeftRightExponential => {
            // Normalized by the integral of (t/sigma)^{-p} exp(-((t^2/sigma^2+1)/(t/sigma)-2)).
            prefactor() * (two - ((t.square()) / (sig.square()) + t.one()) / (t / &sig)).exp()
        }
        crate::settings::runtime::HFunction::ExponentialCT => {
            let delta_t_sq = (tstar.clone().unwrap() - t).square();
            let tstar_sq = tstar.unwrap().square();
            // info!("dampener: {}", dampener);
            // info!("delta_t_sq: {}", delta_t_sq);
            // info!("tstar_sq: {}", tstar_sq);
            // info!(
            //     "Exp arg: {}",
            //     -sig.inv() * (delta_t_sq / tstar_sq + sig * sig * (dampener * dampener))
            // );
            // info!(
            //     "result: {}",
            //     (-sig.inv() * (delta_t_sq / tstar_sq + sig * sig * (dampener * dampener))).exp()
            // );
            if h_function_settings.enabled_dampening {
                let dampener = delta_t_sq.clone() / (delta_t_sq.clone() - &tstar_sq);
                (-sig.inv() * (delta_t_sq.clone() / tstar_sq + sig.square() * (dampener.square())))
                    .exp()
            } else {
                (-sig.inv() * (delta_t_sq / tstar_sq)).exp()
            }
        }
    }
}

/// Calculate the determinant of any complex-valued input matrix using LU-decomposition.
/// Original C-code by W. Gong and D.E. Soper.
#[allow(unused)]
pub(crate) fn determinant<T: FloatLike>(bb: &[Complex<F<T>>], dimension: usize) -> Complex<F<T>> {
    let one = bb[0].re.one();
    let zero = one.zero();
    // Define matrix related variables.
    let mut determinant = Complex::new(one.clone(), zero.clone());
    let mut indx = [0; MAX_DIMENSION];
    let mut d = 1; // initialize parity parameter

    // Inintialize the matrix to be decomposed with the transferred matrix b.
    let mut aa = bb.to_vec();

    // Define parameters used in decomposition.
    let mut imax = 0;
    let mut flag = 1;
    let mut dumc;
    let mut sum;

    let mut aamax;
    let mut dumr;
    let mut vv = vec![zero.clone(); MAX_DIMENSION];

    // Get the implicit scaling information.
    for i in 0..dimension {
        aamax = zero.clone();
        for j in 0..dimension {
            let r = aa[i * dimension + j].norm_squared();
            if r > aamax {
                aamax = r;
            }
        }
        // Set a flag to check if the determinant is zero.
        if aamax.is_zero() {
            flag = 0;
        }
        // Save the scaling.
        vv[i] = aamax.inv();
    }
    if flag == 1 {
        for j in 0..dimension {
            for i in 0..j {
                sum = aa[i * dimension + j].clone();
                for k in 0..i {
                    sum -= &aa[i * dimension + k] * &aa[k * dimension + j];
                }
                aa[i * dimension + j] = sum;
            }
            //Initialize for the search for largest pivot element.
            aamax = zero.clone();
            for i in j..dimension {
                sum = aa[i * dimension + j].clone();
                for k in 0..j {
                    sum -= &aa[i * dimension + k] * &aa[k * dimension + j];
                }
                aa[i * dimension + j] = sum.clone();
                // Figure of merit for the pivot.
                dumr = &vv[i] * sum.norm_squared();
                // Is it better than the best so far?
                if dumr >= aamax {
                    imax = i;
                    aamax = dumr;
                }
            }
            // See if we need to interchange rows.
            if j != imax {
                for k in 0..dimension {
                    dumc = aa[imax * dimension + k].clone();
                    aa[imax * dimension + k] = aa[j * dimension + k].clone();
                    aa[j * dimension + k] = dumc.clone();
                }
                // Change the parity of d.
                d = -d;
                // Interchange the scale factor.
                vv[imax] = vv[j].clone();
            }
            indx[j] = imax;
            if j + 1 != dimension {
                dumc = aa[j * dimension + j].inv();
                for i in j + 1..dimension {
                    aa[i * dimension + j] *= &dumc;
                }
            }
        }
    }
    // Calculate the determinant using the decomposed matrix.
    if flag == 0 {
        determinant = Complex::new(zero.clone(), zero.clone());
    } else {
        // Multiply the diagonal elements.
        for diagonal in 0..dimension {
            determinant *= &aa[diagonal * dimension + diagonal];
        }
        determinant *= &zero.from_i64(d);
    }
    determinant
}

#[allow(unused)]
pub(crate) fn next_combination_with_replacement(state: &mut [usize], max_entry: usize) -> bool {
    for i in (0..state.len()).rev() {
        if state[i] < max_entry {
            state[i] += 1;
            for j in i + 1..state.len() {
                state[j] = state[i]
            }
            return true;
        }
    }
    false
}

pub(crate) fn compute_loop_part<T: FloatLike>(
    loop_signature: &LoopSignature,
    loop_moms: &LoopMomenta<F<T>>,
) -> ThreeMomentum<F<T>> {
    loop_signature.apply_typed(loop_moms)
}

pub(crate) fn compute_loop_part_subspace<T: FloatLike>(
    loop_signature: &LoopSignature,
    loop_moms: &LoopMomenta<F<T>>,
    subspace: &SubspaceData,
) -> ThreeMomentum<F<T>> {
    let projected: LoopSignature = subspace.project_loop_signature(loop_signature).collect();
    projected.apply_typed(loop_moms)
}

pub(crate) fn compute_shift_part<T: FloatLike>(
    external_signature: &ExternalSignature,
    external_moms: &ExternalFourMomenta<F<T>>,
) -> FourMomentum<F<T>> {
    external_signature
        .apply_typed::<FourMomentum<F<T>>, ExternalIndex, ExternalFourMomenta<F<T>>>(external_moms)
}

pub(crate) fn compute_shift_part_subspace<T: FloatLike>(
    loop_signature: &LoopSignature,
    external_signature: &ExternalSignature,
    loop_moms: &LoopMomenta<F<T>>,
    external_moms: &ExternalThreeMomenta<F<T>>,
    subspace: &SubspaceData,
) -> ThreeMomentum<F<T>> {
    let projected_loop_complement: LoopSignature = subspace
        .project_complement_loop_signature(loop_signature)
        .collect();

    let loop_part = projected_loop_complement.apply_typed(loop_moms);
    let external_part: ThreeMomentum<F<T>> = external_signature
        .apply_typed::<ThreeMomentum<F<T>>, ExternalIndex, ExternalThreeMomenta<F<T>>>(
            external_moms,
        );
    loop_part + external_part
}

pub(crate) fn compute_t_part_of_shift_part<T: FloatLike>(
    external_signature: &ExternalSignature,
    external_moms: &ExternalFourMomenta<F<T>>,
) -> F<T> {
    // external_signature.panic_validate_basis(external_moms);
    external_signature
        .apply_iter(external_moms.iter().map(|m| m.temporal.value.clone()))
        .unwrap_or(external_moms[ExternalIndex(0)].temporal.value.zero())
}

// Bilinear form for E-surface defined as sqrt[(k+p1)^2+m1sq] + sqrt[(k+p2)^2+m2sq] + e_shift
// The Bilinear system then reads 4 k.a.k + 4 k.n + C = 0
#[allow(unused, clippy::type_complexity)]
pub(crate) fn one_loop_e_surface_bilinear_form<T: FloatLike>(
    p1: &[F<T>; 3],
    p2: &[F<T>; 3],
    m1_sq: F<T>,
    m2_sq: F<T>,
    e_shift: F<T>,
) -> ([[F<T>; 3]; 3], [F<T>; 3], F<T>) {
    let zero = e_shift.zero();
    let two = zero.from_i64(2);
    let e_shift_sq = e_shift.square();
    let p1_sq = p1[0].square() + p1[1].square() + p1[2].square();
    let p2_sq = p2[0].square() + p2[1].square() + p2[2].square();

    let zeros = [zero.clone(), zero.clone(), zero.clone()];
    let mut a = [zeros.clone(), zeros.clone(), zeros.clone()];
    a[0][0] = (&p1[0] - &p2[0] - &e_shift) * (&p2[0] - &p1[0] - &e_shift);
    a[0][1] = (&p1[0] - &p2[0]) * (&p2[1] - &p1[1]);
    a[1][0] = a[0][1].clone();
    a[0][2] = (&p1[0] - &p2[0]) * (&p2[2] - &p1[2]);
    a[2][0] = a[0][2].clone();
    a[1][1] = (&p1[1] - &p2[1] - &e_shift) * (&p2[1] - &p1[1] - &e_shift);
    a[1][2] = (&p1[1] - &p2[1]) * (&p2[2] - &p1[2]);
    a[2][1] = a[1][2].clone();
    a[2][2] = (&p1[2] - &p2[2] - &e_shift) * (&p2[2] - &p1[2] - &e_shift);

    let mut b = zeros.clone();
    b[0] =
        (&p2[0] - &p1[0]) * (&m1_sq - &m2_sq + &p1_sq - &p2_sq) + &e_shift_sq * (&p1[0] + &p2[0]);
    b[1] =
        (&p2[1] - &p1[1]) * (&m1_sq - &m2_sq + &p1_sq - &p2_sq) + &e_shift_sq * (&p1[1] + &p2[1]);
    b[2] =
        (&p2[2] - &p1[2]) * (&m1_sq - &m2_sq + &p1_sq - &p2_sq) + &e_shift_sq * (&p1[2] + &p2[2]);

    let c = -&e_shift_sq * &e_shift_sq + &two * &e_shift_sq * (&m1_sq + &m2_sq + &p1_sq + &p2_sq)
        - (&m1_sq - &m2_sq + &p1_sq - &p2_sq) * (&m1_sq - &m2_sq + &p1_sq - &p2_sq);

    (a, b, c)
}

use color_eyre::Result;
use eyre::eyre;
#[allow(unused)]
use std::fmt::LowerExp;

pub trait ApproxEq<U: LowerExp, T: LowerExp>: LowerExp {
    fn approx_eq(&self, other: &U, tolerance: &T) -> bool;

    fn approx_eq_slice(lhs: &[Self], rhs: &[U], tolerance: &T) -> bool
    where
        Self: Sized,
    {
        lhs.iter()
            .zip_eq(rhs)
            .all(|(l, r)| l.approx_eq(r, tolerance))
    }

    fn approx_eq_iterator<'a, I, J>(lhs: I, rhs: J, tolerance: &'a T) -> bool
    where
        Self: Sized + 'a,
        U: 'a,
        I: IntoIterator<Item = &'a Self>,
        J: IntoIterator<Item = &'a U>,
    {
        lhs.into_iter()
            .zip_eq(rhs)
            .all(|(l, r)| l.approx_eq(r, tolerance))
    }

    fn assert_approx_eq(&self, other: &U, tolerance: &T) {
        assert!(
            self.approx_eq(other, tolerance),
            "assert_approx_eq failed: \n{:+e} != \n{:+e} with tolerance {:+e}",
            self,
            other,
            tolerance
        )
    }
    fn approx_eq_res(&self, other: &U, tolerance: &T) -> Result<()> {
        if self.approx_eq(other, tolerance) {
            Ok(())
        } else {
            Err(eyre!(
                "assert_approx_eq failed: \n{:+e} != \n{:+e} with tolerance {:+e}",
                self,
                other,
                tolerance
            ))
        }
    }
}

// pub trait ApproxEqable: Real+PartialOrd+ for<'a> RefSub<&'a Self,Output = Self>+IsZero{}

// impl<T: Real+PartialOrd+ for<'a> RefSub<&'a T,Output = T>+IsZero> ApproxEqable for T{}

impl<T: FloatLike> ApproxEq<F<T>, F<T>> for F<T> {
    fn approx_eq(&self, other: &F<T>, tolerance: &F<T>) -> bool {
        if other.is_zero() {
            self.norm() < tolerance.clone()
        } else {
            ((self.ref_sub(other)) / other).norm() < tolerance.clone()
        }
    }
}

impl<T: FloatLike> ApproxEq<Complex<F<T>>, F<T>> for Complex<F<T>> {
    fn approx_eq(&self, other: &Complex<F<T>>, tolerance: &F<T>) -> bool {
        if !self.norm().re.approx_eq(&other.norm().re, tolerance) {
            return false;
        } else if self.norm().is_zero() || other.norm().is_zero() {
            return true;
        }
        let two_pi = self.re.PI() + self.re.PI();
        let arg_self = self.arg().rem_euclid(&two_pi);
        let arg_other = other.arg().rem_euclid(&two_pi);
        if !arg_self.approx_eq(&arg_other, tolerance) {
            return false;
        }
        true
    }
    fn approx_eq_res(&self, other: &Complex<F<T>>, tolerance: &F<T>) -> Result<()> {
        if !self.norm().re.approx_eq(&other.norm().re, tolerance) {
            return Err(eyre!(
                "Norms are not approximately equal: \n{:+e} != \n{:+e} with tolerance {:+e}",
                &self.norm().re,
                other.norm().re,
                tolerance
            ));
        } else if self.norm().is_zero() || other.norm().is_zero() {
            return Ok(());
        }

        let two_pi = self.re.PI() + self.re.PI();
        let arg_self = self.arg().rem_euclid(&two_pi);
        let arg_other = other.arg().rem_euclid(&two_pi);
        // let arg_diff = (&self.arg() - &other.arg()).rem_euclid(&two_pi);
        // let arg_zero = self.re.zero();
        if !arg_self.approx_eq(&arg_other, tolerance) {
            return Err(eyre!(
                "Phases are not approximately equal: \n{:+e} - \n{:+e}= \n{:+e}!=0 with tolerance {:+e}",
                arg_self,
                arg_other,
                &arg_self - &arg_other,
                tolerance
            ));
        }
        Ok(())
    }
}

impl<T: FloatLike> ApproxEq<F<T>, F<T>> for Complex<F<T>> {
    fn approx_eq(&self, other: &F<T>, tolerance: &F<T>) -> bool {
        self.re.approx_eq(other, tolerance) && self.im.approx_eq(tolerance, tolerance)
    }
    fn approx_eq_res(&self, other: &F<T>, tolerance: &F<T>) -> Result<()> {
        if self.im.approx_eq(tolerance, tolerance) {
            return Err(eyre!(
                "Non-zero imaginary part: \n{:+e} with tolerance {:+e}",
                &self.im,
                tolerance
            ));
        }
        if !self.re.approx_eq(other, tolerance) {
            return Err(eyre!(
                "Real parts are not approximately equal: \n{:+e} != \n{:+e} with tolerance {:+e}",
                &self.re,
                other,
                tolerance
            ));
        }
        Ok(())
    }
}

impl<T: FloatLike> ApproxEq<Complex<F<T>>, F<T>> for F<T> {
    fn approx_eq(&self, other: &Complex<F<T>>, tolerance: &F<T>) -> bool {
        other.re.approx_eq(self, tolerance) && other.im.approx_eq(tolerance, tolerance)
    }

    fn approx_eq_res(&self, other: &Complex<F<T>>, tolerance: &F<T>) -> Result<()> {
        if other.im.approx_eq(tolerance, tolerance) {
            return Err(eyre!(
                "Non-zero imaginary part: \n{:+e} with tolerance {:+e}",
                &other.im,
                tolerance
            ));
        }
        if !other.re.approx_eq(self, tolerance) {
            return Err(eyre!(
                "Real parts are not approximately equal: \n{:+e} != \n{:+e} with tolerance {:+e}",
                &other.re,
                self,
                tolerance
            ));
        }
        Ok(())
    }
}

#[allow(unused)]
pub(crate) fn one_loop_eval_e_surf<T: FloatLike>(
    k: &[F<T>; 3],
    p1: &[F<T>; 3],
    p2: &[F<T>; 3],
    m1_sq: F<T>,
    m2_sq: F<T>,
    e_shift: F<T>,
) -> F<T> {
    ((&k[0] + &p1[0]) * (&k[0] + &p1[0])
        + (&k[1] + &p1[1]) * (&k[1] + &p1[1])
        + (&k[2] + &p1[2]) * (&k[2] + &p1[2])
        + m1_sq)
        .sqrt()
        + ((&k[0] + &p2[0]) * (&k[0] + &p2[0])
            + (&k[1] + &p2[1]) * (&k[1] + &p2[1])
            + (&k[2] + &p2[2]) * (&k[2] + &p2[2])
            + m2_sq)
            .sqrt()
        + e_shift
}

#[allow(unused)]
pub(crate) fn one_loop_eval_e_surf_k_derivative<T: FloatLike>(
    k: &[F<T>; 3],
    p1: &[F<T>; 3],
    p2: &[F<T>; 3],
    m1_sq: F<T>,
    m2_sq: F<T>,
) -> [F<T>; 3] {
    let e1 = ((&k[0] + &p1[0]) * (&k[0] + &p1[0])
        + (&k[1] + &p1[1]) * (&k[1] + &p1[1])
        + (&k[2] + &p1[2]) * (&k[2] + &p1[2])
        + m1_sq)
        .sqrt();
    let e2 = ((&k[0] + &p2[0]) * (&k[0] + &p2[0])
        + (&k[1] + &p2[1]) * (&k[1] + &p2[1])
        + (&k[2] + &p2[2]) * (&k[2] + &p2[2])
        + m2_sq)
        .sqrt();
    [
        (&k[0] + &p1[0]) / &e1 + (&k[0] + &p2[0]) / &e2,
        (&k[1] + &p1[1]) / &e1 + (&k[1] + &p2[1]) / &e2,
        (&k[2] + &p1[2]) / &e1 + (&k[2] + &p2[2]) / &e2,
    ]
}

#[allow(unused)]
pub(crate) fn one_loop_get_e_surf_t_scaling<T: FloatLike>(
    k: &[F<T>; 3],
    p1: &[F<T>; 3],
    p2: &[F<T>; 3],
    m1_sq: F<T>,
    m2_sq: F<T>,
    e_shift: F<T>,
) -> [F<T>; 2] {
    let zero = e_shift.zero();
    let one = zero.one();
    let (a, b, c_coef) = one_loop_e_surface_bilinear_form(p1, p2, m1_sq, m2_sq, e_shift);
    let mut a_coef = zero.clone();
    for i in 0..=2 {
        for j in 0..=2 {
            a_coef += &k[i] * &a[i][j] * &k[j];
        }
    }
    a_coef *= zero.from_i64(4);
    let mut b_coef = zero.clone();
    for i in 0..=2 {
        b_coef += &k[i] * &b[i];
    }
    b_coef *= zero.from_i64(4);
    let discr = b_coef.square() - zero.from_i64(4) * &a_coef * &c_coef;
    if discr < zero {
        [zero.clone(), zero.clone()]
    } else {
        [
            (-&b_coef + discr.sqrt()) / (zero.from_i64(2) * &a_coef),
            (-&b_coef - discr.sqrt()) / (zero.from_i64(2) * &a_coef),
        ]
    }
}

pub(crate) fn box_muller<T: FloatLike>(x1: F<T>, x2: F<T>) -> (F<T>, F<T>) {
    let r = (-x1.from_i64(2) * x1.log()).sqrt();
    let theta = r.from_i64(2) * r.PI() * x2;
    (r.clone() * theta.cos(), r * theta.sin())
}

pub(crate) fn compute_surface_and_volume<T: FloatLike>(n_dim: usize, radius: F<T>) -> (F<T>, F<T>) {
    let mut surface = radius.from_i64(2);
    let one = radius.one();
    let mut volume = one.clone();
    for i in 1..n_dim + 1 {
        (surface, volume) = (
            one.from_i64(2) * one.PI() * volume,
            surface / one.from_i64(i as i64),
        );
    }
    (
        surface * radius.pow(n_dim as u64),
        volume * radius.pow(n_dim as u64),
    )
}

pub(crate) fn get_n_dim_for_n_loop_momenta(
    settings: &SamplingSettings,
    n_loop_momenta: usize,
    n_edges: Option<usize>, // for tropical parameterization, we need to know the number of edges
) -> usize {
    if let Some(parameterization_settings) = settings.get_parameterization_settings() {
        match parameterization_settings.mode {
            ParameterizationMode::HyperSphericalFlat => {
                // Because we use Box-Muller, we need to have an even number of angular dimensions
                let mut n_dim = 3 * n_loop_momenta;
                if n_dim % 2 == 1 {
                    n_dim += 1;
                }
                n_dim + 1
            }
            ParameterizationMode::HyperSpherical
            | ParameterizationMode::Cartesian
            | ParameterizationMode::Spherical
            | ParameterizationMode::RelativeSpherical
            | ParameterizationMode::SphericalCommonRadial
            | ParameterizationMode::SphericalProductCommonRadial
            | ParameterizationMode::MomentumSpace => 3 * n_loop_momenta,
        }
    } else {
        let tropical_part = 2 * n_edges.expect("No tropical subgraph table generated, please run without tropical sampling or regenerate with tables") - 1;
        let d_l = 3 * n_loop_momenta;
        if d_l % 2 == 1 {
            tropical_part + d_l + 1
        } else {
            tropical_part + d_l
        }
    }
}

pub(crate) fn global_parameterize<T: FloatLike>(
    x: &[F<T>],
    e_cm: F<T>,
    settings: &ParameterizationSettings,
) -> (Vec<[F<T>; 3]>, F<T>) {
    debug!("b: {}", settings.b);
    let zero = e_cm.zero();
    let one = zero.one();
    match settings.mode {
        ParameterizationMode::HyperSpherical | ParameterizationMode::HyperSphericalFlat => {
            let mut jac = one.clone();

            let radius: F<T> = match settings.mapping {
                ParameterizationMapping::Log => {
                    // r = e_cm * ln(1 + b*x/(1-x))
                    let b = F::<T>::from_f64(settings.b);
                    let radius = &e_cm * (&one + &b * &x[0] / (&one - &x[0])).log();
                    jac *= &e_cm * &b / (&one - &x[0]) / (&one + &x[0] * (&b - &one));
                    radius
                }
                ParameterizationMapping::Power => {
                    // r = e_cm * b * (x/(1-x))^p; p=1 reproduces the linear map.
                    let b = F::<T>::from_f64(settings.b);
                    let power = F::<T>::from_f64(settings.power);
                    let odds = &x[0] / (&one - &x[0]);
                    let power_minus_one = &power - &one;
                    let radius = &e_cm * &b * odds.powf(&power);
                    jac *=
                        &e_cm * &b * &power * odds.powf(&power_minus_one) / (&one - &x[0]).square();
                    radius
                }
                ParameterizationMapping::Linear => {
                    // r = e_cm * b * x/(1-x)
                    let b = F::<T>::from_f64(settings.b);
                    let radius = &e_cm * &b * &x[0] / (&one - &x[0]);
                    jac *= (&e_cm * &b + &radius).powi(2) / &e_cm / &b;
                    radius
                }
            };
            match settings.mode {
                ParameterizationMode::HyperSpherical => {
                    let phi = zero.from_i64(2) * zero.PI() * &x[1];
                    jac *= zero.from_i64(2) * zero.PI();

                    let mut cos_thetas = Vec::with_capacity(x.len() - 2);
                    let mut sin_thetas = Vec::with_capacity(x.len() - 2);

                    for (i, xi) in x[2..].iter().enumerate() {
                        let cos_theta = -&one + zero.from_i64(2) * xi;
                        jac *= zero.from_i64(2);
                        let sin_theta = (&one - cos_theta.square()).sqrt();
                        // Uniform cos(theta_i) leaves the angular factor
                        // sin(theta_i)^(D-3-i) for this Cartesian ordering.
                        let angular_power = x.len() - 3 - i;
                        if angular_power > 0 {
                            jac *= sin_theta.pow(angular_power as u64);
                        }
                        cos_thetas.push(cos_theta);
                        sin_thetas.push(sin_theta);
                    }

                    let mut concatenated_vecs = Vec::with_capacity(x.len() / 3);
                    let mut base = radius.clone();
                    for (cos_theta, sin_theta) in cos_thetas.iter().zip(sin_thetas.iter()) {
                        concatenated_vecs.push(&base * cos_theta);
                        base *= sin_theta;
                    }
                    concatenated_vecs.push(&base * phi.cos());
                    concatenated_vecs.push(&base * phi.sin());

                    jac *= radius.pow((x.len() - 1) as u64); // hyperspherical coords

                    (
                        concatenated_vecs
                            .chunks(3)
                            .map(|v| [v[0].clone(), v[1].clone(), v[2].clone()])
                            .collect(),
                        jac,
                    )
                }
                ParameterizationMode::HyperSphericalFlat => {
                    // As we will use Box Muller we expect an even number of random variables
                    assert!(x[1..].len().is_multiple_of(2));
                    let mut normal_distributed_xs = vec![];
                    for x_pair in x[1..].chunks(2) {
                        let (z1, z2) = box_muller(x_pair[0].clone(), x_pair[1].clone());
                        normal_distributed_xs.push(z1);
                        normal_distributed_xs.push(z2);
                    }
                    // ignore the last variable generated if we had to pad to get the 3*n_loop_momenta
                    if normal_distributed_xs.len() % 3 != 0 {
                        normal_distributed_xs.pop();
                    }
                    let curr_norm = normal_distributed_xs[..]
                        .iter()
                        .map(|x| x.square())
                        .reduce(|acc, e| acc + &e)
                        .unwrap_or(zero.clone())
                        .sqrt();
                    let surface =
                        compute_surface_and_volume(normal_distributed_xs.len() - 1, radius.clone())
                            .0;
                    jac *= surface;
                    let rescaling_factor = &radius / &curr_norm;
                    (
                        normal_distributed_xs
                            .chunks(3)
                            .map(|v| {
                                [
                                    &v[0] * &rescaling_factor,
                                    &v[1] * &rescaling_factor,
                                    &v[2] * &rescaling_factor,
                                ]
                            })
                            .collect(),
                        jac,
                    )
                }
                _ => unreachable!(),
            }
        }
        ParameterizationMode::Cartesian | ParameterizationMode::Spherical => {
            let mut jac = one.clone();
            let mut vecs = Vec::with_capacity(x.len() / 3);
            for xi in x.chunks(3) {
                let (vec_i, jac_i) = parameterize3d(xi, e_cm.clone(), settings);
                vecs.push(vec_i);
                jac *= jac_i;
            }
            (vecs, jac)
        }
        ParameterizationMode::RelativeSpherical => {
            if x.len() != 6 {
                panic!(
                    "relative_spherical parameterization currently requires exactly two loop three-momenta"
                );
            }
            let spherical_settings = ParameterizationSettings {
                mode: ParameterizationMode::Spherical,
                mapping: settings.mapping.clone(),
                b: settings.b,
                power: settings.power,
                lmb_basis_ids: Default::default(),
                sampling_channels: Default::default(),
            };
            let (common, common_jac) = parameterize3d(&x[0..3], e_cm.clone(), &spherical_settings);
            let (relative, relative_jac) =
                parameterize3d(&x[3..6], e_cm.clone(), &spherical_settings);
            let half = &one / one.from_i64(2);
            let first = [
                &common[0] + &half * &relative[0],
                &common[1] + &half * &relative[1],
                &common[2] + &half * &relative[2],
            ];
            let second = [
                &common[0] - &half * &relative[0],
                &common[1] - &half * &relative[1],
                &common[2] - &half * &relative[2],
            ];
            (vec![first, second], common_jac * relative_jac)
        }
        ParameterizationMode::SphericalCommonRadial => {
            if x.is_empty() || !x.len().is_multiple_of(3) {
                panic!(
                    "spherical_common_radial parameterization requires complete loop three-momenta"
                );
            }
            let n_loop_momenta = x.len() / 3;

            let mut jac = one.clone();
            let radius = match settings.mapping {
                ParameterizationMapping::Log => {
                    let b = F::<T>::from_f64(settings.b);
                    let radius = &e_cm * (&one + &b * &x[0] / (&one - &x[0])).log();
                    jac *= &e_cm * &b / (&one - &x[0]) / (&one + &x[0] * (&b - &one));
                    radius
                }
                ParameterizationMapping::Power => {
                    let b = F::<T>::from_f64(settings.b);
                    let power = F::<T>::from_f64(settings.power);
                    let odds = &x[0] / (&one - &x[0]);
                    let power_minus_one = &power - &one;
                    let radius = &e_cm * &b * odds.powf(&power);
                    jac *=
                        &e_cm * &b * &power * odds.powf(&power_minus_one) / (&one - &x[0]).square();
                    radius
                }
                ParameterizationMapping::Linear => {
                    let b = F::<T>::from_f64(settings.b);
                    let radius = &e_cm * &b * &x[0] / (&one - &x[0]);
                    jac *= (&e_cm * &b + &radius).powi(2) / &e_cm / &b;
                    radius
                }
            };

            let mut remaining_fraction = one.clone();
            let mut stick_jac = one.clone();
            let mut radii = Vec::with_capacity(n_loop_momenta);
            for i in 0..n_loop_momenta {
                let fraction = if i + 1 == n_loop_momenta {
                    remaining_fraction.clone()
                } else {
                    stick_jac *= &remaining_fraction;
                    let fraction = &remaining_fraction * &x[1 + i];
                    remaining_fraction *= &one - &x[1 + i];
                    fraction
                };
                radii.push(&radius * fraction);
            }

            jac *= radius.pow((n_loop_momenta - 1) as u64) * stick_jac;
            for radius_i in &radii {
                jac *= radius_i.square();
            }

            let angle_offset = n_loop_momenta;
            let mut momenta = Vec::with_capacity(n_loop_momenta);
            for (i, radius_i) in radii.iter().enumerate() {
                let phi = zero.from_i64(2) * zero.PI() * &x[angle_offset + 2 * i];
                jac *= zero.from_i64(2) * zero.PI();
                let cos_theta = -&one + zero.from_i64(2) * &x[angle_offset + 2 * i + 1];
                jac *= zero.from_i64(2);
                let sin_theta = (&one - cos_theta.square()).sqrt();
                momenta.push([
                    radius_i * &sin_theta * phi.cos(),
                    radius_i * &sin_theta * phi.sin(),
                    radius_i * cos_theta,
                ]);
            }

            (momenta, jac)
        }
        ParameterizationMode::SphericalProductCommonRadial => {
            if x.len() != 6 {
                panic!(
                    "spherical_product_common_radial parameterization currently requires exactly two loop three-momenta"
                );
            }
            let mut branch_x = x.to_vec();
            let branch_settings = if x[0] < (&one / one.from_i64(2)) {
                branch_x[0] = &x[0] * zero.from_i64(2);
                ParameterizationSettings {
                    mode: ParameterizationMode::Spherical,
                    mapping: settings.mapping.clone(),
                    b: settings.b,
                    power: settings.power,
                    lmb_basis_ids: Default::default(),
                    sampling_channels: Default::default(),
                }
            } else {
                branch_x[0] = (&x[0] - (&one / one.from_i64(2))) * zero.from_i64(2);
                ParameterizationSettings {
                    mode: ParameterizationMode::SphericalCommonRadial,
                    mapping: settings.mapping.clone(),
                    b: settings.b,
                    power: settings.power,
                    lmb_basis_ids: Default::default(),
                    sampling_channels: Default::default(),
                }
            };
            let (momenta, jac) = global_parameterize(&branch_x, e_cm, &branch_settings);
            (momenta, jac * zero.from_i64(2))
        }
        ParameterizationMode::MomentumSpace => (x.as_chunks::<3>().0.to_vec(), one),
    }
}

#[allow(unused)]
pub(crate) fn global_inv_parameterize<T: FloatLike>(
    moms: &[ThreeMomentum<F<T>>],
    e_cm: F<T>,
    settings: &ParameterizationSettings,
) -> (Vec<F<T>>, F<T>) {
    let one = e_cm.one();
    let zero = one.zero();

    match settings.mode {
        ParameterizationMode::HyperSpherical => {
            let mut inv_jac = one.clone();
            let mut xs = Vec::with_capacity(moms.len() * 3);

            let cartesian_xs = moms
                .iter()
                .flat_map(|lv| lv.clone().into_iter())
                .collect::<Vec<F<T>>>();

            let mut k_r_sq = cartesian_xs
                .iter()
                .map(|xi| xi.square())
                .reduce(|acc, e| acc + &e)
                .unwrap_or(zero.clone());
            // cover the degenerate case
            if k_r_sq.is_zero() {
                return (vec![zero.clone(); cartesian_xs.len()], zero);
            }
            let k_r = k_r_sq.sqrt();
            match settings.mapping {
                ParameterizationMapping::Log => {
                    let b = F::<T>::from_f64(settings.b);
                    let x1 = &one - &b / (-&one + &b + (&k_r / &e_cm).exp());
                    inv_jac /= e_cm * &b / (&one - &x1) / (&one + &x1 * (&b - &one));
                    xs.push(x1);
                }
                ParameterizationMapping::Power => {
                    let b = F::<T>::from_f64(settings.b);
                    let power = F::<T>::from_f64(settings.power);
                    let inv_power = &one / &power;
                    let odds = (&k_r / (&e_cm * &b)).powf(&inv_power);
                    let x1 = &odds / (&one + &odds);
                    let power_minus_one = &power - &one;
                    inv_jac /=
                        &e_cm * &b * &power * odds.powf(&power_minus_one) / (&one - &x1).square();
                    xs.push(x1);
                }
                ParameterizationMapping::Linear => {
                    let b = F::<T>::from_f64(settings.b);
                    inv_jac /= (&e_cm * &b + &k_r).powi(2) / &e_cm / &b;
                    xs.push(&k_r / (&e_cm * &b + &k_r));
                }
            };

            let y = cartesian_xs[cartesian_xs.len() - 2].clone();
            let x = cartesian_xs[cartesian_xs.len() - 1].clone();
            let xphi = if x < zero {
                &one + (&one / one.from_i64(2)) * zero.FRAC_1_PI() * x.atan2(&y)
            } else {
                (&one / one.from_i64(2)) * zero.FRAC_1_PI() * x.atan2(&y)
            };
            xs.push(xphi);
            inv_jac /= zero.from_i64(2) * zero.PI();

            for (i, x) in cartesian_xs[..cartesian_xs.len() - 2].iter().enumerate() {
                xs.push((&one / one.from_i64(2)) * (&one + x / k_r_sq.sqrt()));
                inv_jac /= zero.from_i64(2);
                let angular_power = cartesian_xs.len() - 3 - i;
                if angular_power > 0 {
                    inv_jac /= (&one - (x * x / &k_r_sq)).sqrt().pow(angular_power as u64);
                }
                k_r_sq -= x * x;
            }

            inv_jac /= k_r.powi((cartesian_xs.len() - 1) as i32);

            (xs, inv_jac)
        }
        ParameterizationMode::HyperSphericalFlat => {
            panic!(
                "Inverse of flat hyperspherical sampling is not available since it is not bijective."
            );
        }
        ParameterizationMode::Cartesian | ParameterizationMode::Spherical => {
            let mut inv_jac = one;
            let mut xs = Vec::with_capacity(moms.len() * 3);
            for (i, mom) in moms.iter().enumerate() {
                let (xs_i, inv_jac_i) = inv_parametrize3d(mom, e_cm.clone(), settings);

                xs.extend(xs_i);
                inv_jac *= inv_jac_i;
            }
            (xs, inv_jac)
        }
        ParameterizationMode::RelativeSpherical => {
            if moms.len() != 2 {
                panic!(
                    "relative_spherical inverse parameterization currently requires exactly two loop three-momenta"
                );
            }
            let half = &one / one.from_i64(2);
            let common = ThreeMomentum::new(
                (&moms[0].px + &moms[1].px) * &half,
                (&moms[0].py + &moms[1].py) * &half,
                (&moms[0].pz + &moms[1].pz) * &half,
            );
            let relative = ThreeMomentum::new(
                &moms[0].px - &moms[1].px,
                &moms[0].py - &moms[1].py,
                &moms[0].pz - &moms[1].pz,
            );
            let spherical_settings = ParameterizationSettings {
                mode: ParameterizationMode::Spherical,
                mapping: settings.mapping.clone(),
                b: settings.b,
                power: settings.power,
                lmb_basis_ids: Default::default(),
                sampling_channels: Default::default(),
            };
            let (common_xs, common_inv_jac) =
                inv_parametrize3d(&common, e_cm.clone(), &spherical_settings);
            let (relative_xs, relative_inv_jac) =
                inv_parametrize3d(&relative, e_cm.clone(), &spherical_settings);
            (
                common_xs.into_iter().chain(relative_xs).collect(),
                common_inv_jac * relative_inv_jac,
            )
        }
        ParameterizationMode::SphericalCommonRadial => {
            if moms.is_empty() {
                panic!(
                    "spherical_common_radial inverse parameterization requires at least one loop three-momentum"
                );
            }

            let phi_x = |x: &F<T>, y: &F<T>| {
                if y < &zero {
                    &one + (&one / one.from_i64(2)) * zero.FRAC_1_PI() * y.atan2(x)
                } else {
                    (&one / one.from_i64(2)) * zero.FRAC_1_PI() * y.atan2(x)
                }
            };

            let radii_sq = moms
                .iter()
                .map(|mom| mom.px.square() + mom.py.square() + mom.pz.square())
                .collect::<Vec<_>>();
            let radii = radii_sq
                .iter()
                .map(|radius_sq| radius_sq.sqrt())
                .collect::<Vec<_>>();
            let radius = radii
                .iter()
                .fold(zero.clone(), |acc, radius_i| acc + radius_i);

            let mut xs = Vec::with_capacity(3 * moms.len());
            if radius.is_zero() || radii.iter().any(F::is_zero) {
                xs.resize(3 * moms.len(), zero.clone());
                return (xs, zero);
            }

            let mut inv_jac = one.clone();
            let radial_x = match settings.mapping {
                ParameterizationMapping::Log => {
                    let b = F::<T>::from_f64(settings.b);
                    let x1 = &one - &b / (-&one + &b + (&radius / &e_cm).exp());
                    inv_jac /= &e_cm * &b / (&one - &x1) / (&one + &x1 * (&b - &one));
                    x1
                }
                ParameterizationMapping::Power => {
                    let b = F::<T>::from_f64(settings.b);
                    let power = F::<T>::from_f64(settings.power);
                    let inv_power = &one / &power;
                    let odds = (&radius / (&e_cm * &b)).powf(&inv_power);
                    let x1 = &odds / (&one + &odds);
                    let power_minus_one = &power - &one;
                    inv_jac /=
                        &e_cm * &b * &power * odds.powf(&power_minus_one) / (&one - &x1).square();
                    x1
                }
                ParameterizationMapping::Linear => {
                    let b = F::<T>::from_f64(settings.b);
                    inv_jac /= (&e_cm * &b + &radius).powi(2) / &e_cm / &b;
                    &radius / (&e_cm * &b + &radius)
                }
            };

            xs.push(radial_x);
            let mut remaining_fraction = one.clone();
            let mut stick_jac = one.clone();
            for (i, radius_i) in radii.iter().enumerate().take(moms.len() - 1) {
                let fraction = radius_i / &radius;
                xs.push(&fraction / &remaining_fraction);
                stick_jac *= &remaining_fraction;
                remaining_fraction -= fraction;
                if remaining_fraction.is_zero() && i + 2 < moms.len() {
                    return (vec![zero.clone(); 3 * moms.len()], zero);
                }
            }

            inv_jac /= radius.pow((moms.len() - 1) as u64) * stick_jac;
            for radius_sq in &radii_sq {
                inv_jac /= radius_sq;
            }

            for (mom, radius_i) in moms.iter().zip(radii.iter()) {
                xs.push(phi_x(&mom.px, &mom.py));
                inv_jac /= zero.from_i64(2) * zero.PI();
                xs.push((&one / one.from_i64(2)) * (&one + &mom.pz / radius_i));
                inv_jac /= zero.from_i64(2);
            }

            (xs, inv_jac)
        }
        ParameterizationMode::SphericalProductCommonRadial => {
            if moms.len() != 2 {
                panic!(
                    "spherical_product_common_radial inverse parameterization currently requires exactly two loop three-momenta"
                );
            }
            let spherical_settings = ParameterizationSettings {
                mode: ParameterizationMode::Spherical,
                mapping: settings.mapping.clone(),
                b: settings.b,
                power: settings.power,
                lmb_basis_ids: Default::default(),
                sampling_channels: Default::default(),
            };
            let common_radial_settings = ParameterizationSettings {
                mode: ParameterizationMode::SphericalCommonRadial,
                mapping: settings.mapping.clone(),
                b: settings.b,
                power: settings.power,
                lmb_basis_ids: Default::default(),
                sampling_channels: Default::default(),
            };
            let (mut xs, product_inv_jac) =
                global_inv_parameterize(moms, e_cm.clone(), &spherical_settings);
            let (_, common_inv_jac) = global_inv_parameterize(moms, e_cm, &common_radial_settings);
            if let Some(first_x) = xs.first_mut() {
                *first_x = &*first_x / zero.from_i64(2);
            }
            (xs, product_inv_jac + common_inv_jac)
        }
        ParameterizationMode::MomentumSpace => (
            moms.iter()
                .flat_map(|mom| [mom.px.clone(), mom.py.clone(), mom.pz.clone()])
                .collect(),
            one,
        ),
    }
}

/// Map a vector in the unit hypercube to the infinite hypercube.
/// Also compute the Jacobian.
pub(crate) fn parameterize3d<T: FloatLike>(
    x: &[F<T>],
    e_cm: F<T>,
    settings: &ParameterizationSettings,
) -> ([F<T>; 3], F<T>) {
    let zero = e_cm.zero();
    let one = zero.one();
    let mut l_space = [zero.clone(), zero.clone(), zero.clone()];
    let mut jac = one.clone();

    match settings.mode {
        ParameterizationMode::Cartesian => match settings.mapping {
            ParameterizationMapping::Log => {
                for i in 0..3 {
                    let x = x[i].clone();
                    l_space[i] = &e_cm * (&x / (&one - &x)).log();
                    jac *= &e_cm / (&x - &x * &x);
                }
            }
            ParameterizationMapping::Linear => {
                for i in 0..3 {
                    let x = x[i].clone();
                    l_space[i] = &e_cm * (&one / (&one - &x) - &one / &x);
                    jac *= &e_cm * (&one / (&x * &x) + &one / ((&one - &x) * (&one - &x)));
                }
            }
            ParameterizationMapping::Power => {
                panic!("Power radial mapping is only supported for spherical coordinates");
            }
        },
        ParameterizationMode::Spherical => {
            let radius = match settings.mapping {
                ParameterizationMapping::Log => {
                    // r = &e_cm * ln(1 + b*&x/(1-&x))
                    let x = x[0].clone();
                    let b = F::<T>::from_f64(settings.b);
                    let radius = &e_cm * (&one + &b * &x / (&one - &x)).log();
                    jac *= &e_cm * &b / (&one - &x) / (&one + &x * (&b - &one));

                    radius
                }
                ParameterizationMapping::Power => {
                    // r = e_cm * b * (x/(1-x))^p; p=1 reproduces the linear map.
                    let x = x[0].clone();
                    let b = F::<T>::from_f64(settings.b);
                    let power = F::<T>::from_f64(settings.power);
                    let odds = &x / (&one - &x);
                    let power_minus_one = &power - &one;
                    let radius = &e_cm * &b * odds.powf(&power);
                    jac *= &e_cm * &b * &power * odds.powf(&power_minus_one) / (&one - &x).square();

                    radius
                }
                ParameterizationMapping::Linear => {
                    // r = &e_cm * b * x/(1-x)
                    let b = F::<T>::from_f64(settings.b);
                    let radius = &e_cm * &b * &x[0] / (&one - &x[0]);
                    jac *= (&e_cm * &b + &radius).powi(2) / &e_cm / &b;
                    radius
                }
            };
            let phi = zero.from_i64(2) * zero.PI() * &x[1];
            jac *= zero.from_i64(2) * zero.PI();

            let cos_theta = -&one + zero.from_i64(2) * &x[2]; // out of range
            jac *= zero.from_i64(2);
            let sin_theta = (&one - cos_theta.square()).sqrt();

            l_space[0] = &radius * &sin_theta * phi.cos();
            l_space[1] = &radius * &sin_theta * phi.sin();
            l_space[2] = &radius * &cos_theta;

            jac *= radius.square(); // spherical coord
        }
        _ => {
            panic!(
                "Inappropriate parameterization mapping specified for parameterize: {:?}.",
                settings.mode.clone()
            );
        }
    }

    (l_space, jac)
}

pub(crate) fn inv_parametrize3d<T: FloatLike>(
    mom: &ThreeMomentum<F<T>>,
    e_cm: F<T>,
    settings: &ParameterizationSettings,
) -> ([F<T>; 3], F<T>) {
    let one = e_cm.one();
    let zero = one.zero();
    if settings.mode == ParameterizationMode::Cartesian {
        let two = one.from_i64(2);
        let mut coordinates = [zero.clone(), zero.clone(), zero.clone()];
        let mut jac = one.clone();
        for (coordinate, momentum) in coordinates.iter_mut().zip([&mom.px, &mom.py, &mom.pz]) {
            let scaled = momentum / &e_cm;
            let magnitude = scaled.abs();
            // Invert the negative half first to avoid overflow or cancellation in the tails.
            let lower = match settings.mapping {
                ParameterizationMapping::Log => {
                    let exponential = (-magnitude).exp();
                    &exponential / (&one + &exponential)
                }
                ParameterizationMapping::Linear => {
                    if magnitude > one {
                        let reciprocal = &one / &magnitude;
                        let twice_reciprocal = &two * reciprocal;
                        &twice_reciprocal
                            / ((&one + twice_reciprocal.square()).sqrt() + &one + &twice_reciprocal)
                    } else {
                        &two / ((magnitude.square() + two.square()).sqrt() + magnitude + &two)
                    }
                }
                ParameterizationMapping::Power => {
                    panic!("Power radial mapping is only supported for spherical coordinates");
                }
            };
            let upper = &one - &lower;
            jac *= match settings.mapping {
                ParameterizationMapping::Log => &lower * &upper / &e_cm,
                ParameterizationMapping::Linear => {
                    lower.square() * upper.square() / (&e_cm * (lower.square() + upper.square()))
                }
                ParameterizationMapping::Power => unreachable!(),
            };
            *coordinate = if scaled < zero { lower } else { upper };
        }
        return (coordinates, jac);
    }
    if settings.mode != ParameterizationMode::Spherical {
        panic!("Inverse mapping requires cartesian or spherical coordinates");
    }

    let mut jac = one.clone();

    let x = &mom.px;
    let y = &mom.py;
    let z = &mom.pz;

    let k_r_sq = x.square() + y.square() + z.square();
    let k_r = k_r_sq.sqrt();

    let x2 = if y < &zero {
        &one + (&one / one.from_i64(2)) * zero.FRAC_1_PI() * y.atan2(x)
    } else {
        (&one / one.from_i64(2)) * zero.FRAC_1_PI() * y.atan2(x)
    };

    // cover the degenerate case
    if k_r_sq.is_zero() {
        return ([zero.clone(), x2.clone(), zero.clone()], zero.clone());
    }

    let x1 = match settings.mapping {
        ParameterizationMapping::Log => {
            let b = F::<T>::from_f64(settings.b);
            let x1 = &one - &b / (-&one + &b + (&k_r / &e_cm).exp());
            jac /= &e_cm * &b / (&one - &x1) / (&one + &x1 * (&b - &one));
            x1
        }
        ParameterizationMapping::Power => {
            let b = F::<T>::from_f64(settings.b);
            let power = F::<T>::from_f64(settings.power);
            let inv_power = &one / &power;
            let odds = (&k_r / (&e_cm * &b)).powf(&inv_power);
            let x1 = &odds / (&one + &odds);
            let power_minus_one = &power - &one;
            jac /= &e_cm * &b * &power * odds.powf(&power_minus_one) / (&one - &x1).square();
            x1
        }
        ParameterizationMapping::Linear => {
            let b = F::<T>::from_f64(settings.b);
            jac /= (&e_cm * &b + &k_r).powi(2) / &e_cm / &b;
            &k_r / (&e_cm * &b + &k_r)
        }
    };

    let x3 = (&one / one.from_i64(2)) * (&one + z / &k_r);
    jac /= zero.from_i64(2) * zero.PI();
    jac /= zero.from_i64(2);
    jac /= k_r.square();

    let x = [x1, x2, x3];

    (x, jac)
}

pub const MINUTE: usize = 60;
pub const HOUR: usize = 3_600;
pub const DAY: usize = 86_400;
pub const WEEK: usize = 604_800;
pub(crate) fn format_wdhms(seconds: usize) -> String {
    let mut compound_duration = vec![];
    if seconds == 0 {
        compound_duration.push("0s".to_string());
        return compound_duration.join(" ");
    }

    let mut sec = seconds % WEEK;
    // weeks
    let ws = seconds / WEEK;
    if ws != 0 {
        compound_duration.push(format!("{ws}w"));
    }

    // days
    let ds = sec / DAY;
    sec %= DAY;
    if ds != 0 {
        compound_duration.push(format!("{ds}d"));
    }

    // hours
    let hs = sec / HOUR;
    sec %= HOUR;
    if hs != 0 {
        compound_duration.push(format!("{hs}h"));
    }

    // minutes
    let ms = sec / MINUTE;
    sec %= MINUTE;
    if ms != 0 {
        compound_duration.push(format!("{ms}m"));
    }

    // seconds
    if sec != 0 {
        compound_duration.push(format!("{sec}s"));
    }

    compound_duration.join(" ")
}

pub(crate) fn format_wdhms_from_duration(duration: Duration) -> String {
    format_wdhms(duration.as_secs() as usize)
}

#[allow(unused)]
pub(crate) fn inv_3x3_sig_matrix(mat: [[isize; 3]; 3]) -> [[isize; 3]; 3] {
    let denom = -mat[0][2] * mat[1][1] * mat[2][0]
        + mat[0][1] * mat[1][2] * mat[2][0]
        + mat[0][2] * mat[1][0] * mat[2][1]
        - mat[0][0] * mat[1][2] * mat[2][1]
        - mat[0][1] * mat[1][0] * mat[2][2]
        + mat[0][0] * mat[1][1] * mat[2][2];
    if denom != 1 && denom != -1 {
        panic!("Non invertible signature matrix.");
    }
    let mut inv_mat = [[0; 3]; 3];
    inv_mat[0][0] = (-mat[1][2] * mat[2][1] + mat[1][1] * mat[2][2]) * denom;
    inv_mat[0][1] = (mat[0][2] * mat[2][1] - mat[0][1] * mat[2][2]) * denom;
    inv_mat[0][2] = (-mat[0][2] * mat[1][1] + mat[0][1] * mat[1][2]) * denom;
    inv_mat[1][0] = (mat[1][2] * mat[2][0] - mat[1][0] * mat[2][2]) * denom;
    inv_mat[1][1] = (-mat[0][2] * mat[2][0] + mat[0][0] * mat[2][2]) * denom;
    inv_mat[1][2] = (mat[0][2] * mat[1][0] - mat[0][0] * mat[1][2]) * denom;
    inv_mat[2][0] = (-mat[1][1] * mat[2][0] + mat[1][0] * mat[2][1]) * denom;
    inv_mat[2][1] = (mat[0][1] * mat[2][0] - mat[0][0] * mat[2][1]) * denom;
    inv_mat[2][2] = (-mat[0][1] * mat[1][0] + mat[0][0] * mat[1][1]) * denom;

    inv_mat
}

pub mod index_vec;

#[allow(unused)]
pub(crate) fn format_for_compare_digits(x: F<f64>, y: F<f64>) -> (String, String) {
    let mut string_x = format!("{:.16e}", x);
    let mut string_y = format!("{:.16e}", y);

    #[allow(clippy::comparison_chain)]
    if string_x.len() > string_y.len() {
        for _ in 0..(string_x.len() - string_y.len()) {
            string_y.push(' ');
        }
    } else if string_y.len() > string_x.len() {
        for _ in 0..(string_y.len() - string_x.len()) {
            string_x.push(' ');
        }
    }

    let string_vec = string_x
        .chars()
        .zip(string_y.chars())
        .map(|(char_x, char_y)| {
            if char_x == char_y {
                (char_x.to_string().green(), char_y.to_string().green())
            } else {
                (char_x.to_string().red(), char_y.to_string().red())
            }
        })
        .collect_vec();

    let string_x = string_vec.iter().map(|(x, _)| x).join("");
    let string_y = string_vec.iter().map(|(_, y)| y).join("");

    (string_x, string_y)
}

#[allow(unused)]
pub(crate) fn format_evaluation_time(time: Duration) -> String {
    let seconds = time.as_secs_f64();
    let (value, unit) = if seconds >= 1.0 {
        (seconds, "s")
    } else if seconds >= 1.0e-3 {
        (seconds * 1.0e3, "ms")
    } else {
        (seconds * 1.0e6, "µs")
    };

    let precision = if value >= 100.0 {
        0
    } else if value >= 10.0 {
        1
    } else {
        2
    };

    format!("{value:.precision$} {unit}")
}

pub(crate) fn duration_from_secs_f64_saturating(time: f64) -> Duration {
    if time.is_nan() || time <= 0.0 {
        return Duration::ZERO;
    }

    if time.is_infinite() {
        return Duration::MAX;
    }

    let max_seconds = Duration::MAX.as_secs();
    if time >= max_seconds as f64 {
        return Duration::MAX;
    }

    let seconds = time.trunc() as u64;
    let fractional = (time - seconds as f64).clamp(0.0, 0.999_999_999_999);
    let nanos = (fractional * 1_000_000_000.0) as u32;
    Duration::new(seconds, nanos)
}

pub(crate) fn format_evaluation_time_from_f64(time: f64) -> String {
    format_evaluation_time(duration_from_secs_f64_saturating(time))
}

pub(crate) fn normalize_tabled_separator_rows(rendered: &str) -> String {
    let original_lines = rendered.lines().collect_vec();
    let visible_lines = original_lines
        .iter()
        .map(|line| strip_ansi_escape_codes(line).chars().collect_vec())
        .collect_vec();

    original_lines
        .iter()
        .enumerate()
        .map(|(row_index, line)| {
            let visible_line = &visible_lines[row_index];
            if is_top_border_row(visible_line) {
                return rebuild_top_border_row(&visible_lines, row_index);
            }

            if is_internal_separator_row(visible_line) {
                return rebuild_internal_separator_row(&visible_lines, row_index);
            }

            if is_bottom_border_row(visible_line) {
                return rebuild_bottom_border_row(&visible_lines, row_index);
            }

            (*line).to_string()
        })
        .join("\n")
}

fn rebuild_top_border_row(visible_lines: &[Vec<char>], row_index: usize) -> String {
    rebuild_border_row(visible_lines, row_index, '╭', '╮', |_, _, has_below| {
        if has_below { '┬' } else { '─' }
    })
}

fn rebuild_internal_separator_row(visible_lines: &[Vec<char>], row_index: usize) -> String {
    rebuild_border_row(
        visible_lines,
        row_index,
        '├',
        '┤',
        |row_index, column, has_below| match (
            row_index > 0 && has_vertical_border(&visible_lines[row_index - 1], column),
            has_below,
        ) {
            (false, false) => '─',
            (false, true) => '┬',
            (true, false) => '┴',
            (true, true) => '┼',
        },
    )
}

fn rebuild_bottom_border_row(visible_lines: &[Vec<char>], row_index: usize) -> String {
    rebuild_border_row(
        visible_lines,
        row_index,
        '╰',
        '╯',
        |row_index, column, _| {
            if row_index > 0 && has_vertical_border(&visible_lines[row_index - 1], column) {
                '┴'
            } else {
                '─'
            }
        },
    )
}

fn rebuild_border_row<F>(
    visible_lines: &[Vec<char>],
    row_index: usize,
    left: char,
    right: char,
    mut intersection_for_column: F,
) -> String
where
    F: FnMut(usize, usize, bool) -> char,
{
    let border = &visible_lines[row_index];
    if border.len() < 2 {
        return border.iter().collect();
    }

    let last = border.len() - 1;
    let mut rebuilt = String::with_capacity(border.len());
    rebuilt.push(left);

    for column in 1..last {
        let has_vertical_below = row_index + 1 < visible_lines.len()
            && has_vertical_border(&visible_lines[row_index + 1], column);
        rebuilt.push(intersection_for_column(
            row_index,
            column,
            has_vertical_below,
        ));
    }

    rebuilt.push(right);
    rebuilt
}

fn has_vertical_border(line: &[char], column: usize) -> bool {
    line.get(column)
        .copied()
        .is_some_and(|ch| matches!(ch, '│' | '├' | '┤' | '┬' | '┴' | '┼'))
}

fn is_internal_separator_row(line: &[char]) -> bool {
    matches!(line.first(), Some('├')) && matches!(line.last(), Some('┤'))
}

fn is_top_border_row(line: &[char]) -> bool {
    matches!(line.first(), Some('╭')) && matches!(line.last(), Some('╮'))
}

fn is_bottom_border_row(line: &[char]) -> bool {
    matches!(line.first(), Some('╰')) && matches!(line.last(), Some('╯'))
}

fn strip_ansi_escape_codes(line: &str) -> String {
    let mut stripped = String::with_capacity(line.len());
    let mut chars = line.chars().peekable();

    while let Some(ch) = chars.next() {
        if ch == '\u{1b}' && chars.peek() == Some(&'[') {
            chars.next();
            for code in chars.by_ref() {
                if ('@'..='~').contains(&code) {
                    break;
                }
            }
            continue;
        }

        stripped.push(ch);
    }

    stripped
}

pub(crate) fn into_complex_ff64<T: FloatLike>(c: &Complex<F<T>>) -> Complex<F<f64>> {
    Complex::new(c.re.into_ff64(), c.im.into_ff64())
}

#[test]
fn exact_binary_cube_promotion_preserves_original_value_and_settings_policy() {
    fn check<T: FloatLike>() {
        let one = F::<T>::from_f64(1.0);
        let exact = F(T::from_f64_exact_binary(0.1));
        let expected = one.from_i64(3602879701896397) / one.from_i64(2).powi(55);
        assert_eq!(exact, expected);
        assert_eq!(exact.into_f64().to_bits(), 0.1_f64.to_bits());
        if one.epsilon() < F::<T>::from_f64(f64::EPSILON) {
            let settings_value = F::<T>::from_f64(0.1);
            assert!(exact > settings_value);
            assert!((settings_value - &one / one.from_i64(10)).abs() <= one.epsilon());
        }
    }
    check::<f64>();
    check::<QuadFloat>();
    check::<ArbPrec>();
}

#[test]
fn decimal_promotion_cache_preserves_parsing_and_precision() {
    // Include non-binary decimals, signed zero, subnormals and overflow-range
    // exponents. Repeat after more unique inputs than the cache can retain.
    let mut inputs = vec![
        0.0,
        -0.0,
        0.1,
        -0.3,
        1.0,
        0.5,
        f64::MIN_POSITIVE,
        f64::from_bits(1),
        f64::MAX,
        -f64::MAX,
        f64::INFINITY,
        f64::NEG_INFINITY,
        f64::NAN,
    ];
    inputs.extend((0..1200).map(|i| f64::from_bits(0x3fb999999999999a + i)));
    for _ in 0..2 {
        for &input in &inputs {
            let expected: QuadFloat = if input.is_finite() {
                SymbolicaFloat::parse(
                    &input.to_string(),
                    Some(DoubleFloat::default().get_precision()),
                )
                .unwrap()
                .into()
            } else {
                QuadFloat::from_f64_exact_binary(input)
            };
            // Hit the cache as well as populating it. Conversion back to f64
            // additionally checks the sign of zero and non-finite handling.
            for _ in 0..2 {
                let actual = QuadFloat::from_f64_decimal(input);
                if input.is_nan() {
                    assert!(F(actual).is_nan());
                } else {
                    assert_eq!(actual, expected);
                    assert_eq!(
                        F(actual).into_f64().to_bits(),
                        F(expected).into_f64().to_bits()
                    );
                }
            }
            // Interleave precisions on one thread; the numeric spelling alone
            // must never let a lower-precision result satisfy a higher request.
            for precision in [256, 1000, 256, 1000] {
                let actual = match precision {
                    256 => VarFloat::<256>::from_f64_decimal(input).float,
                    1000 => VarFloat::<1000>::from_f64_decimal(input).float,
                    _ => unreachable!(),
                };
                let expected = if input.is_finite() {
                    Float::with_val(precision, Float::parse(input.to_string()).unwrap())
                } else {
                    Float::with_val(precision, input)
                };
                assert_eq!(actual.prec(), precision);
                if input.is_nan() {
                    assert!(actual.is_nan());
                } else {
                    assert_eq!(actual, expected);
                    assert_eq!(actual.is_sign_negative(), expected.is_sign_negative());
                }
            }
        }
    }
}

#[test]
fn complex_compare() {
    let ltd = Complex::new(0.11773583919739394, -0.22157450463964778).map(F);

    let cff = Complex::new(0.11773583589023458, -0.22157450446824836).map(F);

    ltd.approx_eq_res(&cff, &F(0.00000001)).unwrap();
}

#[cfg(test)]
mod formatting_tests {
    use super::{ArbPrec, F, FloatLike, f128};

    fn assert_has_scientific_exponent(rendered: &str) {
        let (mantissa, exponent) = rendered
            .split_once('e')
            .unwrap_or_else(|| panic!("missing exponent marker in: {rendered}"));
        assert!(
            !mantissa.is_empty(),
            "missing mantissa in scientific output: {rendered}"
        );
        assert!(
            exponent.parse::<i32>().is_ok(),
            "invalid exponent in scientific output: {rendered}"
        );
    }

    fn assert_fractional_digits(rendered: &str, expected_digits: usize) {
        let (mantissa, _) = rendered
            .split_once('e')
            .unwrap_or_else(|| panic!("missing exponent marker in: {rendered}"));
        let (_, fractional) = mantissa
            .trim_start_matches(['+', '-'])
            .split_once('.')
            .unwrap_or_else(|| panic!("missing decimal point in mantissa: {rendered}"));
        assert_eq!(
            fractional.len(),
            expected_digits,
            "unexpected number of fractional digits in: {rendered}"
        );
    }

    fn assert_scientific_formatting_for_backend<T: FloatLike>(x: F<T>) {
        let plain = format!("{:e}", x.clone());
        assert_has_scientific_exponent(&plain);

        let plus = format!("{:+e}", x.clone());
        assert_has_scientific_exponent(&plus);
        assert!(plus.starts_with('+') || plus.starts_with('-'));

        let precision = format!("{:+.16e}", x.clone());
        assert_has_scientific_exponent(&precision);
        assert!(
            precision.contains('.'),
            "missing decimal point in scientific output: {precision}"
        );

        let padded = format!("{:>24.6e}", x);
        assert!(
            padded.len() >= 24,
            "formatted width should be at least requested width: {padded}"
        );
        assert_has_scientific_exponent(padded.trim_start());
    }

    #[test]
    fn display_matches_underlying_f64() {
        let value = F(12.5f64);
        assert_eq!(format!("{}", value), "12.5");
    }

    #[test]
    fn lower_exp_keeps_exponent_with_precision() {
        let rendered = format!("{:+.16e}", F(1.23456789e-8f64));
        assert_has_scientific_exponent(&rendered);
        assert_fractional_digits(&rendered, 16);
    }

    #[test]
    fn lower_exp_width_and_alignment_behave_as_expected() {
        let rendered = format!("{:>20.4e}", F(12.345f64));
        assert_eq!(rendered.len(), 20);
        assert_has_scientific_exponent(rendered.trim_start());

        let left_aligned = format!("{:<20.4e}", F(12.345f64));
        assert_eq!(left_aligned.len(), 20);
        assert_has_scientific_exponent(left_aligned.trim_end());
    }

    #[test]
    fn lower_exp_sign_and_zero_pad_behave_as_expected() {
        let positive = format!("{:+020.4e}", F(12.345f64));
        assert_eq!(positive.len(), 20);
        assert!(positive.starts_with('+'));
        assert_has_scientific_exponent(&positive);

        let negative = format!("{:+020.4e}", F(-12.345f64));
        assert_eq!(negative.len(), 20);
        assert!(negative.starts_with('-'));
        assert_has_scientific_exponent(&negative);
    }

    #[test]
    fn lower_exp_precision_controls_fractional_digits() {
        let p2 = format!("{:.2e}", F(1.0f64 / 3.0f64));
        assert_fractional_digits(&p2, 2);

        let p9 = format!("{:.9e}", F(1.0f64 / 3.0f64));
        assert_fractional_digits(&p9, 9);
    }

    #[test]
    fn lower_exp_works_across_float_backends() {
        assert_scientific_formatting_for_backend(F::<f64>::from_f64(1.23456789e-8));
        assert_scientific_formatting_for_backend(F::<f128>::from_f64(1.23456789e-8));
        assert_scientific_formatting_for_backend(F::<ArbPrec>::from_f64(1.23456789e-8));
    }

    #[test]
    fn quad_square_matches_multiplication() {
        let zero = F::<f128>::from_f64(0.0);
        let zero_sq = zero.square();
        assert!(!zero_sq.is_nan());
        assert!(!zero_sq.is_infinite());
        assert_eq!(zero_sq, zero);

        let x = F::<f128>::from_f64(500.0);
        let x_sq = x.square();
        assert!(!x_sq.is_nan());
        assert!(!x_sq.is_infinite());
        assert_eq!(x_sq, x * x);
        assert_eq!(x_sq.sqrt(), x);
    }
}

pub mod symbols;
pub use symbols::{GS, W_};

type TensorLibStore = RwLock<TensorLibrary<ParamTensor<ExplicitKey<Aind>>, Aind>>;
type FunLibStore = SymbolLib<ParamTensor<ShadowedStructure<Aind>>, Wrap>;

pub static TENSORLIB: LazyLock<TensorLibStore> = LazyLock::new(|| RwLock::new(hep_lib_atom()));

pub static FUN_LIB: LazyLock<FunLibStore> = LazyLock::new(|| {
    let mut lib = Wrap::new_lib();
    lib.insert(INBUILTS.conj, |a: ParamTensor<ShadowedStructure<Aind>>| {
        a.map_data_self(|x| x.conj())
    });
    lib
});

pub static PARAM_FUN_LIB: LazyLock<SymbolLib<ParamTensor<ShadowedStructure<Aind>>, Panic>> =
    LazyLock::new(|| {
        let mut lib = Panic::new_lib();
        lib.insert(INBUILTS.conj, |a: ParamTensor<ShadowedStructure<Aind>>| {
            a.map_data_self(|x| x.conj())
        });
        lib
    });

pub static VAKINT: OnceLock<Result<Vakint>> = OnceLock::new();

pub fn init_vakint() -> Result<()> {
    let r: &Result<Vakint, eyre::Report> =
        VAKINT.get_or_init(|| Vakint::new().map_err(|e| eyre!(e.to_string())));

    match r {
        Ok(_) => Ok(()),
        Err(e) => Err(eyre!(e.to_string())),
    }
}

pub fn vakint() -> Result<&'static Vakint> {
    let r = VAKINT
        .get()
        .ok_or_else(|| eyre!("Vakint not initialised"))?;
    match r {
        Ok(v) => Ok(v),
        Err(e) => Err(eyre!(e.to_string())),
    }
}

impl<T: FloatLike> momtrop::float::MomTropFloat for F<T> {
    #[inline]
    fn PI(&self) -> Self {
        self.PI()
    }

    #[inline]
    fn abs(&self) -> Self {
        self.abs()
    }

    #[inline]
    fn cos(&self) -> Self {
        <F<T> as Real>::cos(self)
    }

    #[inline]
    fn exp(&self) -> Self {
        <F<T> as Real>::exp(self)
    }

    #[inline]
    fn one(&self) -> Self {
        <F<T> as SymFloatLike>::one(self)
    }

    #[inline]
    fn from_f64(&self, value: f64) -> Self {
        F::from_f64(value)
    }

    #[inline]
    fn from_isize(&self, value: isize) -> Self {
        Self(self.0.from_i64(value as i64))
    }

    #[inline]
    fn inv(&self) -> Self {
        <F<T> as SymFloatLike>::inv(self)
    }

    #[inline]
    fn ln(&self) -> Self {
        <F<T> as Real>::log(self)
    }

    #[inline]
    fn powf(&self, power: &Self) -> Self {
        <F<T> as Real>::powf(self, power)
    }

    #[inline]
    fn sin(&self) -> Self {
        <F<T> as Real>::sin(self)
    }

    #[inline]
    fn sqrt(&self) -> Self {
        <F<T> as Real>::sqrt(self)
    }

    #[inline]
    fn zero(&self) -> Self {
        <F<T> as SymFloatLike>::zero(self)
    }

    #[inline]
    fn to_f64(&self) -> f64 {
        self.into_f64()
    }
}

pub trait Length {
    fn len(&self) -> usize;
    fn is_empty(&self) -> bool {
        self.len() == 0
    }
}

impl<I, T> Length for TiVec<I, T> {
    fn len(&self) -> usize {
        self.len()
    }

    fn is_empty(&self) -> bool {
        self.is_empty()
    }
}

impl<T> Length for Vec<T> {
    fn len(&self) -> usize {
        self.len()
    }
}

pub(crate) fn ose_atom_from_index(index: EdgeIndex) -> Atom {
    function!(
        GS.ose,
        usize::from(index) as i64 // Atom::from(FlatIndex::from(0))
    )
}

pub(crate) fn cut_energy(index: EdgeIndex) -> Atom {
    function!(GS.energy, usize::from(index) as i64)
}

pub(crate) fn external_energy_atom_from_index(index: EdgeIndex) -> Atom {
    GS.emr_mom(index, Atom::from(ExpandedIndex::from_iter([0])))
}

pub mod newton_solver;
use include_dir::{Dir, include_dir};
static BUILTIN_MODELS: Dir = include_dir!("$CARGO_MANIFEST_DIR/../../assets/models/json");

pub fn load_generic_model(name: &str) -> Model {
    if let Some(file) = BUILTIN_MODELS.get_file(format!("{}/{}.json", name, name)) {
        Model::from_str(file.contents_utf8().unwrap().into(), "json").unwrap()
    } else {
        panic!("Model {} not found in built-in models.", name);
    }
}

#[cfg(test)]
pub mod test_utils;

#[cfg(feature = "python")]
impl<'py> pyo3::IntoPyObject<'py> for F<f64> {
    type Target = pyo3::types::PyFloat;
    type Output = pyo3::Bound<'py, Self::Target>;
    type Error = std::convert::Infallible;

    #[inline]
    fn into_pyobject(self, py: pyo3::Python<'py>) -> Result<Self::Output, Self::Error> {
        Ok(pyo3::types::PyFloat::new(py, self.0))
    }
}

#[cfg(feature = "python")]
impl<'py> pyo3::IntoPyObject<'py> for &F<f64> {
    type Target = pyo3::types::PyFloat;
    type Output = pyo3::Bound<'py, Self::Target>;
    type Error = std::convert::Infallible;

    #[inline]
    fn into_pyobject(self, py: pyo3::Python<'py>) -> Result<Self::Output, Self::Error> {
        (*self).into_pyobject(py)
    }
}

#[cfg(feature = "python")]
impl<'a, 'py> pyo3::FromPyObject<'a, 'py> for F<f64> {
    type Error = pyo3::PyErr;

    fn extract(obj: pyo3::Borrowed<'a, 'py, pyo3::types::PyAny>) -> Result<Self, Self::Error> {
        obj.extract::<f64>().map(F)
    }
}