gammalooprs/utils/

mod.rs

use crate::integrands::process::GenericEvaluatorFloat;
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};
use rug::ops::{CompleteRound, Pow};
use rug::{Assign, Float};
use schemars::JsonSchema;
use serde::{Deserialize, Deserializer, Serialize};
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, PanicMissingConcrete, SymbolLib};
use spenso::network::library::symbolic::{ExplicitKey, TensorLibrary};
use spenso::network::parsing::ShadowedStructure;
use spenso::structure::concrete_index::ExpandedIndex;
use spenso::tensors::complex::RealOrComplexTensor;
use spenso::tensors::data::StorageTensor;
use spenso::tensors::parametric::to_param::ToAtom;
use spenso::tensors::parametric::{MixedTensor, ParamTensor};
use spenso_hep_lib::hep_lib_atom;
use symbolica::{
    domains::{
        dual::HyperDual,
        float::{FixedPrecision, Float as SymbolicaFloat, FloatLike as SymFloatLike},
    },
    prelude::*,
};

use statrs::function::gamma::{gamma, gamma_lr, gamma_ur};
use std::cmp::{Ord, Ordering};
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;

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 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 {
        let n = x.numerator();

        let n = match n {
            Integer::Double(f) => Float::with_val(N, f),
            Integer::Large(f) => Float::with_val(N, f),
            Integer::Single(f) => Float::with_val(N, f),
        };

        let d = x.denominator();

        let d = match d {
            Integer::Double(f) => Float::with_val(N, f),
            Integer::Large(f) => Float::with_val(N, f),
            Integer::Single(f) => Float::with_val(N, f),
        };

        let r = n / d;

        VarFloat {
            float: rug::Float::with_val(N, r),
        }
    }
}

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 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 {
        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 {
        Self {
            float: value.into_inner(),
        }
    }
}

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 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 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(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()
    }
}

impl FloatLike for ArbPrec {
    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(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>>>> {
        None
    }

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

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);
        }

        let valid = Float::parse(format!("{}", x)).unwrap();
        VarFloat {
            float: rug::Float::with_val(N, valid),
        }
    }

    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);
        }

        SymbolicaFloat::parse(
            &format!("{}", x),
            Some(DoubleFloat::default().get_precision()),
        )
        .unwrap()
        .into()
    }

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

    fn floor_via_symbolica(&self) -> Self {
        SymbolicaFloat::from(self).into_inner().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 PrecisionUpgradable for ArbPrec {
    type Higher = ArbPrec;
    type Lower = f128;

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

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

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
    +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;

    #[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};

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 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 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 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 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))
    }

    #[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())
    }

    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 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(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, and the `f128` trait impls/effective epsilon logic are written so
// flipping this alias remains a valid drop-in escape hatch.
// pub type f128 = VarFloat<113>;
pub type f128 = QuadFloat;
pub type ArbPrec = VarFloat<1000>;

/// 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)
}

/// 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 sqrt_pi = new_constant(t, &t.values[0].PI().sqrt());
    let sig = if let Some(s) = sigma {
        new_constant(t, &s)
    } else {
        new_constant(t, &F::<T>::from_f64(h_function_settings.sigma))
    };
    let power = h_function_settings.power;
    match h_function_settings.function {
        crate::settings::runtime::HFunction::Exponential => {
            (-(t.clone() * t) / (sig.clone() * sig.clone())).exp()
                * new_constant(t, &F::<T>::from_f64(2_f64))
                / (sqrt_pi * sig)
        }
        crate::settings::runtime::HFunction::PolyExponential => {
            // Result of \int_0^{\infty} dt (t/sigma)^{-p} exp(2-t^2/sigma^2-sigma^2/t^2)
            let normalisation = match power {
                None | Some(0) => {
                    sqrt_pi.clone() * &sig / new_constant(t, &F::<T>::from_f64(2_f64))
                }
                Some(1) => new_constant(t, &F::<T>::from_f64(0.841_568_215_070_771_4)) * &sig,
                Some(3) => new_constant(t, &F::<T>::from_f64(1.033_476_847_068_688_6)) * &sig,
                Some(4) => new_constant(t, &F::<T>::from_f64(1.329_340_388_179_137)) * &sig,
                Some(6) => new_constant(t, &F::<T>::from_f64(2.880_237_507_721_463_7)) * &sig,
                Some(7) => new_constant(t, &F::<T>::from_f64(4.783_566_971_347_609)) * &sig,
                Some(9) => new_constant(t, &F::<T>::from_f64(16.225_745_976_182_285)) * &sig,
                Some(10) => new_constant(t, &F::<T>::from_f64(32.735_007_058_911_25)) * &sig,
                Some(12) => new_constant(t, &F::<T>::from_f64(155.837_465_922_583_42)) * &sig,
                Some(13) => new_constant(t, &F::<T>::from_f64(364.658_500_356_566_04)) * &sig,
                Some(15) => new_constant(t, &F::<T>::from_f64(2_257.637_553_015_473)) * &sig,
                Some(16) => new_constant(t, &F::<T>::from_f64(5_939.804_418_537_864)) * &sig,
                _ => panic!(
                    "Value {} of power in poly exponential h function not supported",
                    power.unwrap()
                ),
            };
            let prefactor = match power {
                None | Some(0) => normalisation.inv(),
                Some(p) => (t.clone() / &sig).inv().pow(p as u64) / normalisation,
            };
            prefactor
                * (new_constant(t, &F::<T>::from_f64(2_f64))
                    - (t.clone() * t) / (sig.clone() * &sig)
                    - (sig.clone() * sig) / (t.clone() * t))
                    .exp()
        }
        crate::settings::runtime::HFunction::PolyLeftRightExponential => {
            // Result of \int_0^{\infty} dt (t/sigma)^{-p} exp( -((t^2/sigma^2 +1)/ (t/sigma) -2) )
            let normalisation = match power {
                None | Some(0) => new_constant(t, &F::<T>::from_f64(2.066_953_694_137_377)) * &sig,
                Some(1) => new_constant(t, &F::<T>::from_f64(1.683_136_430_141_542_8)) * &sig,
                Some(3) => new_constant(t, &F::<T>::from_f64(3.750_090_124_278_92)) * &sig,
                Some(4) => new_constant(t, &F::<T>::from_f64(9.567_133_942_695_218)) * &sig,
                Some(6) => new_constant(t, &F::<T>::from_f64(139.373_101_752_153_5)) * &sig,
                Some(7) => new_constant(t, &F::<T>::from_f64(729.317_000_713_132_1)) * &sig,
                Some(9) => new_constant(t, &F::<T>::from_f64(32_336.242_742_929_753)) * &sig,
                Some(10) => new_constant(t, &F::<T>::from_f64(263_205.217_049_469)) * &sig,
                Some(12) => new_constant(t, &F::<T>::from_f64(2.427_503_717_893_097_5e7)) * &sig,
                Some(13) => new_constant(t, &F::<T>::from_f64(2.694_265_921_644_289e8)) * &sig,
                Some(15) => new_constant(t, &F::<T>::from_f64(9.040_742_057_760_125e12)) * &sig,
                Some(16) => new_constant(t, &F::<T>::from_f64(1.452_517_480_246_491_3e14)) * &sig,
                _ => panic!(
                    "Value {} of power in poly exponential h function not supported",
                    power.unwrap()
                ),
            };

            // println!("normalisation: {}", normalisation);
            // println!("t: {}", t);
            // println!("sig: {}", sig);
            // println!("power: {:?}", power);

            let prefactor = match power {
                None | Some(0) => normalisation.inv(),
                Some(p) => (t.clone() / &sig).inv().pow(p as u64) / normalisation,
            };

            // println!("prefactor: {}", prefactor);
            prefactor
                * (new_constant(t, &F::<T>::from_f64(2_f64))
                    - ((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 sqrt_pi = t.PI().sqrt();
    let sig = if let Some(s) = sigma {
        s
    } else {
        F::<T>::from_f64(h_function_settings.sigma)
    };
    let power = h_function_settings.power;
    match h_function_settings.function {
        crate::settings::runtime::HFunction::Exponential => {
            (-(t.square()) / (sig.square())).exp() * F::<T>::from_f64(2_f64) / (sqrt_pi * &sig)
        }
        crate::settings::runtime::HFunction::PolyExponential => {
            // Result of \int_0^{\infty} dt (t/sigma)^{-p} exp(2-t^2/sigma^2-sigma^2/t^2)
            let normalisation = match power {
                None | Some(0) => sqrt_pi * &sig / F::<T>::from_f64(2_f64),
                Some(1) => F::<T>::from_f64(0.841_568_215_070_771_4) * &sig,
                Some(3) => F::<T>::from_f64(1.033_476_847_068_688_6) * &sig,
                Some(4) => F::<T>::from_f64(1.329_340_388_179_137) * &sig,
                Some(6) => F::<T>::from_f64(2.880_237_507_721_463_7) * &sig,
                Some(7) => F::<T>::from_f64(4.783_566_971_347_609) * &sig,
                Some(9) => F::<T>::from_f64(16.225_745_976_182_285) * &sig,
                Some(10) => F::<T>::from_f64(32.735_007_058_911_25) * &sig,
                Some(12) => F::<T>::from_f64(155.837_465_922_583_42) * &sig,
                Some(13) => F::<T>::from_f64(364.658_500_356_566_04) * &sig,
                Some(15) => F::<T>::from_f64(2_257.637_553_015_473) * &sig,
                Some(16) => F::<T>::from_f64(5_939.804_418_537_864) * &sig,
                _ => panic!(
                    "Value {} of power in poly exponential h function not supported",
                    power.unwrap()
                ),
            };
            let prefactor = match power {
                None | Some(0) => normalisation.inv(),
                Some(p) => (t / &sig).powi(-(p as i32)) / normalisation,
            };
            prefactor
                * (F::<T>::from_f64(2_f64)
                    - (t.square()) / (sig.square())
                    - (sig.square()) / (t.square()))
                .exp()
        }
        crate::settings::runtime::HFunction::PolyLeftRightExponential => {
            // Result of \int_0^{\infty} dt (t/sigma)^{-p} exp( -((t^2/sigma^2 +1)/ (t/sigma) -2) )
            let normalisation = match power {
                None | Some(0) => F::<T>::from_f64(2.066_953_694_137_377) * &sig,
                Some(1) => F::<T>::from_f64(1.683_136_430_141_542_8) * &sig,
                Some(3) => F::<T>::from_f64(3.750_090_124_278_92) * &sig,
                Some(4) => F::<T>::from_f64(9.567_133_942_695_218) * &sig,
                Some(6) => F::<T>::from_f64(139.373_101_752_153_5) * &sig,
                Some(7) => F::<T>::from_f64(729.317_000_713_132_1) * &sig,
                Some(9) => F::<T>::from_f64(32_336.242_742_929_753) * &sig,
                Some(10) => F::<T>::from_f64(263_205.217_049_469) * &sig,
                Some(12) => F::<T>::from_f64(2.427_503_717_893_097_5e7) * &sig,
                Some(13) => F::<T>::from_f64(2.694_265_921_644_289e8) * &sig,
                Some(15) => F::<T>::from_f64(9.040_742_057_760_125e12) * &sig,
                Some(16) => F::<T>::from_f64(1.452_517_480_246_491_3e14) * &sig,
                _ => panic!(
                    "Value {} of power in poly exponential h function not supported",
                    power.unwrap()
                ),
            };

            // println!("normalisation: {}", normalisation);
            // println!("t: {}", t);
            // println!("sig: {}", sig);
            // println!("power: {:?}", power);

            let prefactor = match power {
                None | Some(0) => normalisation.inv(),
                Some(p) => (t / &sig).powi(-(p as i32)) / normalisation,
            };

            // println!("prefactor: {}", prefactor);
            prefactor
                * (F::<T>::from_f64(2_f64) - ((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();
                        if i > 0 {
                            jac *= sin_theta.pow(i 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(),
            };
            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 = F::<T>::from_f64(2.) * zero.PI() * &x[angle_offset + 2 * i];
                jac *= F::<T>::from_f64(2.) * zero.PI();
                let cos_theta = -&one + F::<T>::from_f64(2.) * &x[angle_offset + 2 * i + 1];
                jac *= F::<T>::from_f64(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] < F::<T>::from_f64(0.5) {
                branch_x[0] = &x[0] * F::<T>::from_f64(2.0);
                ParameterizationSettings {
                    mode: ParameterizationMode::Spherical,
                    mapping: settings.mapping.clone(),
                    b: settings.b,
                    power: settings.power,
                    lmb_basis_ids: Default::default(),
                }
            } else {
                branch_x[0] = (&x[0] - F::<T>::from_f64(0.5)) * F::<T>::from_f64(2.0);
                ParameterizationSettings {
                    mode: ParameterizationMode::SphericalCommonRadial,
                    mapping: settings.mapping.clone(),
                    b: settings.b,
                    power: settings.power,
                    lmb_basis_ids: Default::default(),
                }
            };
            let (momenta, jac) = global_parameterize(&branch_x, e_cm, &branch_settings);
            (momenta, jac * F::<T>::from_f64(2.0))
        }
        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 + F::<T>::from_f64(0.5) * zero.FRAC_1_PI() * x.atan2(&y)
            } else {
                F::<T>::from_f64(0.5) * zero.FRAC_1_PI() * x.atan2(&y)
            };
            xs.push(xphi);
            inv_jac /= F::<T>::from_f64(2.) * zero.PI();

            for (i, x) in cartesian_xs[..cartesian_xs.len() - 2].iter().enumerate() {
                xs.push(F::<T>::from_f64(0.5) * (&one + x / k_r_sq.sqrt()));
                inv_jac /= F::<T>::from_f64(2.);
                if i > 0 {
                    inv_jac /= (&one - (x * x / &k_r_sq)).sqrt().powi(i as i32);
                }
                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(),
            };
            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 + F::<T>::from_f64(0.5) * zero.FRAC_1_PI() * y.atan2(x)
                } else {
                    F::<T>::from_f64(0.5) * 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 /= F::<T>::from_f64(2.) * zero.PI();
                xs.push(F::<T>::from_f64(0.5) * (&one + &mom.pz / radius_i));
                inv_jac /= F::<T>::from_f64(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(),
            };
            let common_radial_settings = ParameterizationSettings {
                mode: ParameterizationMode::SphericalCommonRadial,
                mapping: settings.mapping.clone(),
                b: settings.b,
                power: settings.power,
                lmb_basis_ids: 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 / F::<T>::from_f64(2.0);
            }
            (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 = F::<T>::from_f64(2.) * zero.PI() * &x[1];
            jac *= F::<T>::from_f64(2.) * zero.PI();

            let cos_theta = -&one + F::<T>::from_f64(2.) * &x[2]; // out of range
            jac *= F::<T>::from_f64(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::Spherical {
        panic!("Inverse mapping is only implemented for 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 + F::<T>::from_f64(0.5) * zero.FRAC_1_PI() * y.atan2(x)
    } else {
        F::<T>::from_f64(0.5) * 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 = F::<T>::from_f64(0.5) * (&one + z / &k_r);
    jac /= F::<T>::from_f64(2.) * zero.PI();
    jac /= F::<T>::from_f64(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 inverse_gamma_lr(a: f64, p: f64, n_iter: usize) -> f64 {
    // this algorithm is taken from https://dl.acm.org/doi/pdf/10.1145/22721.23109

    // get an estimate for x0 to start newton iterations.
    let q = 1.0 - p;

    if (1.0 - 1.0e-8..=1.0 + 1.0e-8).contains(&a) {
        return -q.ln();
    }

    let gamma_a = gamma(a);
    let b = q * gamma_a;
    let c = 0.577_215_664_901_532_9;

    let mut x0 = 0.5;
    if a < 1.0 {
        if b > 0.6 || (b >= 0.45 && a >= 0.3) {
            let u = if b * q > 10e-8 {
                (p * gamma(a + 1.0)).powf(a.recip())
            } else {
                (-q / a - c).exp()
            };
            x0 = u / (1.0 - u / (a + 1.0));
        } else if a < 0.3 && (0.35..=0.6).contains(&b) {
            let t = (-c - b).exp();
            let u = t * t.exp();
            x0 = t * u.exp();
        } else if (0.15..=0.35).contains(&b) || ((0.15..0.45).contains(&b) && a >= 0.3) {
            let y = -b.ln();
            let u = y - (1.0 - a) * y.ln();
            x0 = y - (1.0 - a) * y.ln() - (1.0 + (1.0 - a) / (1.0 + u)).ln();
        } else if 0.01 < b && b < 0.15 {
            let y = -b.ln();
            let u = y - (1.0 - a) * y.ln();
            x0 = y
                - (1.0 - a) * u.ln()
                - ((u * u + 2.0 * (3.0 - a) * u + (2.0 - a) * (3.0 - a))
                    / (u * u + (5.0 - a) * u + 2.0))
                    .ln();
        } else if b <= 0.01 {
            let y = -b.ln();
            let c1 = (a - 1.0) * y.ln();
            let c2 = (a - 1.0) * (1.0 + c1);
            let c3 = (a - 1.0) * (-0.5 * c1 * c1 + (a - 2.0) * c1 + (3.0 * a - 5.0) * 0.5);
            let c4 = (a - 1.0)
                * (1.0 / 3.0 * c1 * c1 * c1 - (3.0 * a - 5.0) * 0.5 * c1 * c1
                    + (a * a - 6.0 * a + 7.0) * c1
                    + (11.0 * a * a - 46.0 * a + 47.0) / 6.0);
            let c5 = (a - 1.0)
                * (-0.25 * c1 * c1 * c1 * c1
                    + (11.0 * a - 7.0) / 6.0 * c1 * c1 * c1
                    + (-3.0 * a * a - 13.0) * c1 * c1
                    + (2.0 * a * a * a - 25.0 * a * a + 72.0 * a - 61.0) * 0.5 * c1
                    + (25.0 * a * a * a - 195.0 * a * a + 477.0 * a - 379.0) / 12.0);
            x0 = y + c1 + c2 / (y) + c3 / (y * y) + c4 / (y * y * y) + c5 / (y * y * y * y);

            if b <= 1.0e-28 {
                return x0;
            }
        }
    } else {
        let (pref, tau) = if p < 0.5 { (-1.0, p) } else { (1.0, q) };
        let t = (-2.0 * tau.ln()).sqrt();

        let a_0 = 3.31125922108741;
        let a_1 = 11.6616720288968;
        let a_2 = 4.28342155967104;
        let a_3 = 0.213623493715853;

        let b_1 = 6.61053765625462;
        let b_2 = 6.40691597760039;
        let b_3 = 1.27364489782223;
        let b_4 = 3.611_708_101_884_203e-2;

        let t2 = t * t;
        let t3 = t2 * t;
        let t4 = t3 * t;

        let numerator = a_0 + a_1 * t + a_2 * t2 + a_3 * t3;
        let denominator = 1.0 + b_1 * t + b_2 * t2 + b_3 * t3 + b_4 * t4;

        let s = pref * (t - numerator / denominator);
        let s2 = s * s;
        let s3 = s * s2;
        let s4 = s * s3;
        let s5 = s * s4;

        let a_sqrt = a.sqrt();

        let w = a + s * a_sqrt + (s2 - 1.0) / 3.0 + (s3 - 7.0 * s) / (36.0 * a_sqrt)
            - (3.0 * s4 + 7.0 * s2 - 16.0) / (810.0 * a)
            + (9.0 * s5 + 256.0 * s3 - 433.0 * s) / (38880.0 * a * a_sqrt);

        if a >= 500.0 && (1.0 - w / a).abs() < 1.0e-6 {
            return w;
        } else if p > 0.5 {
            if w < 3.0 * a {
                x0 = w;
            } else {
                let d = 2f64.max(a * (a - 1.0));
                if b > 10f64.powf(-d) {
                    let u = -b.ln() + (a - 1.0) * w.ln() - (1.0 + (1.0 - a) / (1.0 + w)).ln();
                    x0 = -b.ln() + (a - 1.0) * u.ln() - (1.0 + (1.0 - a) / (1.0 + u)).ln();
                } else {
                    let y = -b.ln();
                    let c1 = (a - 1.0) * y.ln();
                    let c2 = (a - 1.0) * (1.0 + c1);
                    let c3 = (a - 1.0) * (-0.5 * c1 * c1 + (a - 2.0) * c1 + (3.0 * a - 5.0) * 0.5);
                    let c4 = (a - 1.0)
                        * (1.0 / 3.0 * c1 * c1 * c1 - (3.0 * a - 5.0) * 0.5 * c1 * c1
                            + (a * a - 6.0 * a + 7.0) * c1
                            + (11.0 * a * a - 46.0 * a + 47.0) / 6.0);
                    let c5 = (a - 1.0)
                        * (-0.25 * c1 * c1 * c1 * c1
                            + (11.0 * a - 7.0) / 6.0 * c1 * c1 * c1
                            + (-3.0 * a * a - 13.0) * c1 * c1
                            + (2.0 * a * a * a - 25.0 * a * a + 72.0 * a - 61.0) * 0.5 * c1
                            + (25.0 * a * a * a - 195.0 * a * a + 477.0 * a - 379.0) / 12.0);
                    x0 = y + c1 + c2 / (y) + c3 / (y * y) + c4 / (y * y * y) + c5 / (y * y * y * y);
                }
            }
        } else {
            // this part is heavily simplified from the paper, if any issues occur this estimate
            // will need more refinement.
            let v = (p * gamma(a + 1.0)).ln();
            x0 = ((v + w) / a).exp();
        }
    }

    // start iteration
    let mut x_n = x0;
    for _ in 0..n_iter {
        let r = x_n.powf(a - 1.0) * (-x_n).exp() / gamma_a;
        if x_n <= 0. {
            x_n = 1.0e-16;
        }
        let t_n = if p <= 0.5 {
            (gamma_lr(a, x_n) - p) / r
        } else {
            -(gamma_ur(a, x_n) - q) / r
        };
        let w_n = (a - 1.0 - x_n) / 2.0;

        let h_n = if t_n.abs() <= 0.1 && (w_n * t_n).abs() <= 0.1 {
            t_n + w_n * t_n * t_n
        } else {
            t_n
        };

        x_n -= h_n;
    }

    x_n
}

#[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 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<MixedTensor<F<f64>, ExplicitKey<Aind>>, Aind>>;
type FunLibStore =
    SymbolLib<RealOrComplexTensor<F<f64>, ShadowedStructure<Aind>>, PanicMissingConcrete>;

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

pub static FUN_LIB: LazyLock<FunLibStore> = LazyLock::new(|| {
    let mut lib = PanicMissingConcrete::new_lib();
    lib.insert(INBUILTS.conj, |a| match a {
        RealOrComplexTensor::Complex(c) => RealOrComplexTensor::Complex(c.map_data(|x| x.conj())),
        RealOrComplexTensor::Real(r) => RealOrComplexTensor::Real(r),
    });
    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 INT_FUN_LIB: LazyLock<
    SymbolLib<RealOrComplexTensor<i64, ShadowedStructure<Aind>>, PanicMissingConcrete>,
> = LazyLock::new(|| {
    let mut lib = PanicMissingConcrete::new_lib();
    lib.insert(INBUILTS.conj, |a| match a {
        RealOrComplexTensor::Complex(c) => RealOrComplexTensor::Complex(c.map_data(|x| x.conj())),
        RealOrComplexTensor::Real(r) => RealOrComplexTensor::Real(r),
    });
    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_api")]
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_api")]
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_api")]
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)
    }
}