grid1d/intervals/

operations.rs

#![deny(rustdoc::broken_intra_doc_links)]

use crate::{
    Contains, Interval, IntervalTrait,
    intervals::{
        IntervalBoundsRuntime, LowerBoundRuntime, UpperBoundRuntime,
        bounded::{IntervalFiniteLength, IntervalFinitePositiveLength},
        unbounded::{IntervalInfiniteLength, IntervalLowerUnboundedUpperUnbounded},
    },
};
use duplicate::duplicate_item;
use num_valid::RealScalar;
use serde::{Deserialize, Serialize};

//------------------------------------------------------------------------------------------------
/// Common trait for intervals with operations between different interval types.
///
/// The [`IntervalOperations`] trait provides the foundation for interval arithmetic
/// and set operations. It defines the core operations that can be performed between
/// intervals, such as intersection, while maintaining type safety and mathematical
/// correctness.
///
/// ## Design Philosophy
///
/// This trait enables:
/// - **Cross-type operations**: Operations between different interval types
/// - **Unified interface**: Common operations across all interval implementations
/// - **Type preservation**: Results maintain appropriate mathematical types
/// - **Performance**: Operations are designed for efficiency
///
/// ## Required Associated Types
///
/// - `RealType`: The scalar type used for interval bounds and operations
///
/// ## Core Operations
///
/// ### Intersection
/// The primary operation defined by this trait is interval intersection, which
/// computes the overlapping region between two intervals. The intersection
/// algorithm automatically handles:
/// - Different boundary types (open/closed combinations)
/// - Boundary precedence (most restrictive boundary wins)
/// - Result type determination (returns most appropriate interval type)
/// - Edge cases (touching boundaries, disjoint intervals)
///
/// ## Mathematical Properties
///
/// All implementations must maintain these mathematical properties:
/// - **Commutativity**: `A.intersection(B) = B.intersection(A)`
/// - **Associativity**: `(A ∩ B) ∩ C = A ∩ (B ∩ C)`
/// - **Idempotence**: `A.intersection(A) = Some(A)`
/// - **Monotonicity**: If `A ⊆ B`, then `A ∩ C ⊆ B ∩ C`
///
/// ## Generic Programming
///
/// This trait enables writing generic functions that work with any interval type:
///
/// ```rust
/// use grid1d::intervals::*;
/// use num_valid::RealScalar;
///
/// fn find_overlap<I1, I2, T>(a: &I1, b: &I2) -> Option<Interval<T>>
/// where
///     I1: IntervalOperations<RealType = T>,
///     I2: IntervalTrait<RealType = T>,
///     T: RealScalar,
/// {
///     a.intersection(b)
/// }
/// ```
///
/// ## Implementation Requirements
///
/// Types implementing this trait must:
/// - Support intersection with any other interval type
/// - Return mathematically correct results
/// - Handle boundary conditions properly
/// - Maintain performance characteristics
pub trait IntervalOperations: Contains {
    /// Computes the intersection of this interval with another interval.
    ///
    /// Returns `Some(interval)` containing the overlap if the intervals intersect,
    /// or `None` if they are disjoint. The returned interval uses the most specific
    /// type possible to represent the intersection.
    ///
    /// # Mathematical Definition
    ///
    /// For intervals A and B, their intersection A ∩ B is the set of all points
    /// that belong to both A and B. If no such points exist, the intersection is empty.
    ///
    /// # Boundary Logic
    ///
    /// The intersection algorithm follows these rules:
    /// - Lower bound: Maximum of the two lower bounds, respecting open/closed semantics
    /// - Upper bound: Minimum of the two upper bounds, respecting open/closed semantics
    /// - Result type: Most restrictive combination of boundary types
    ///
    /// # Special Cases
    ///
    /// - **Singleton result**: When intervals meet at exactly one point
    /// - **Empty result**: When intervals don't overlap (`None` is returned)
    /// - **Unbounded result**: When both intervals extend infinitely in the same direction
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    /// use try_create::New;
    ///
    /// let a = IntervalClosed::new(0.0, 2.0);        // [0., 2.]
    /// let b = IntervalClosed::new(1.0, 3.0);        // [1., 3.]
    /// let intersection = a.intersection(&b).unwrap();
    /// // Result: [1., 2.]
    /// let expected_intersection = IntervalClosed::new(1., 2.); // [1., 2.]
    /// assert!(matches!(intersection, Interval::FiniteLength(
    ///     IntervalFiniteLength::PositiveLength(
    ///         IntervalFinitePositiveLength::Closed(expected_intersection)))));
    ///
    /// let c = IntervalOpen::new(0.0, 1.0);          // (0., 1.)
    /// let d = IntervalOpen::new(2.0, 3.0);          // (2., 3.)
    /// assert!(c.intersection(&d).is_none());        // Disjoint intervals
    ///
    /// let e = IntervalClosed::new(0.0, 1.0);        // [0., 1.]
    /// let f = IntervalClosed::new(1.0, 2.0);        // [1., 2.]
    /// let singleton = e.intersection(&f).unwrap();
    /// // Result: singleton at point 1.0
    /// let expected_singleton = IntervalSingleton::new(1.); // [1.]
    /// assert!(matches!(singleton, Interval::FiniteLength(
    ///     IntervalFiniteLength::ZeroLength(expected_singleton))));
    ///
    /// let g = IntervalOpen::new(0.0, 1.0);          // (0, 1)
    /// let h = IntervalLowerClosedUpperOpen::new(0.5, 2.0); // [0.5, 2)
    /// let mixed = g.intersection(&h).unwrap();
    /// // Result: [0.5, 1) - takes most restrictive bounds
    /// let expected_intersection = IntervalLowerClosedUpperOpen::new(0.5, 1.0); // [0.5, 1.)
    /// assert!(matches!(mixed, Interval::FiniteLength(
    ///     IntervalFiniteLength::PositiveLength(
    ///         IntervalFinitePositiveLength::LowerClosedUpperOpen(expected_intersection)))));
    /// ```
    ///
    /// # Performance
    ///
    /// This operation is O(1) in computation time. The result may require allocation
    /// for the returned interval, but the intersection logic itself is constant time.
    ///
    /// # Return Type
    ///
    /// The method always returns the most general [`Interval`] enum to accommodate
    /// any possible intersection result. You can use [`TryFrom`] to convert to more
    /// specific types if needed:
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let a = IntervalClosed::new(0.0, 2.0);
    /// let b = IntervalClosed::new(1.0, 3.0);
    ///
    /// if let Some(intersection) = a.intersection(&b) {
    ///     // Try to convert to specific type
    ///     if let Ok(closed) = IntervalClosed::try_from(intersection) {
    ///         println!("Intersection is a closed interval: {:?}", closed);
    ///     }
    /// }
    /// ```
    fn intersection<IntervalType: IntervalTrait<RealType = Self::RealType>>(
        &self,
        other: &IntervalType,
    ) -> Option<Interval<Self::RealType>> {
        // get the maximum between the lower bounds of `self` and `other`
        let lb_max = self.max_lower_bound(other);

        // get the minimum between the upper bounds of `self` and `other`
        let ub_min = self.min_upper_bound(other);

        match (lb_max, ub_min) {
            (Some(lower_bound), Some(upper_bound)) => {
                // If the lower bound is less than the upper bound, we have a valid interval (i.e. a non-empty intersection).
                // If the lower and upper bounds are both closed, a valid interval (a singleton)
                // is built also when the value of the lower bound is equal to the value of the upper bound.
                IntervalFiniteLength::try_new(lower_bound, upper_bound)
                    .ok()
                    .map(|interval| interval.into())
            }
            (None, Some(upper_bound)) => {
                Some(IntervalInfiniteLength::new_lower_unbounded(upper_bound).into())
            }
            (Some(lower_bound), None) => {
                Some(IntervalInfiniteLength::new_upper_unbounded(lower_bound).into())
            }
            (None, None) => {
                // unbounded - unbounded case
                Some(IntervalLowerUnboundedUpperUnbounded::new().into())
            }
        }
    }

    /// Computes the union of this interval with another.
    ///
    /// The union of two intervals can result in either a single connected interval
    /// (if they overlap or touch) or a set of two disjoint intervals. This method
    /// returns an [`IntervalUnion`] enum to represent these two possibilities.
    ///
    /// ## Mathematical Properties
    ///
    /// - **Commutativity**: `A.union(B) == B.union(A)`
    /// - **Associativity**: `(A.union(B)).union(C) == A.union(B.union(C))`
    /// - **Identity**: `A.union(∅) = A`
    ///
    /// ## Return Value
    ///
    /// - [`IntervalUnion::SingleConnected`]: If the intervals overlap or are adjacent,
    ///   their union is a single, larger interval.
    /// - [`IntervalUnion::TwoDisjoint`]: If the intervals are separate,
    ///   the union is a collection of two disjoint intervals.
    ///
    /// ## Examples
    ///
    /// ### Overlapping Intervals (Single Connected Union)
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let i1 = IntervalClosed::new(0.0, 2.0); // [0, 2]
    /// let i2 = IntervalClosed::new(1.0, 3.0); // [1, 3]
    ///
    /// let union = i1.union(&i2);
    ///
    /// match union {
    ///     IntervalUnion::SingleConnected(result) => {
    ///         // The union is [0, 3]
    ///         let expected: Interval<f64> = IntervalClosed::new(0.0, 3.0).into();
    ///         assert_eq!(result, expected);
    ///     }
    ///     _ => panic!("Expected a single connected union"),
    /// }
    /// ```
    ///
    /// ### Disjoint Intervals
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let i1 = IntervalClosed::new(0.0, 1.0); // [0, 1]
    /// let i2 = IntervalClosed::new(2.0, 3.0); // [2, 3]
    ///
    /// let union = i1.union(&i2);
    ///
    /// match union {
    ///     IntervalUnion::TwoDisjoint{ left, right } => {
    ///         assert_eq!(left, i1.into());
    ///         assert_eq!(right, i2.into());
    ///     }
    ///     _ => panic!("Expected a disjoint union"),
    /// }
    /// ```
    ///
    /// ### Touching Intervals (Single Connected Union)
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let i1 = IntervalLowerClosedUpperOpen::new(0.0, 1.0); // [0, 1)
    /// let i2 = IntervalLowerClosedUpperUnbounded::new(1.0); // [1, +inf)
    ///
    /// let union = i1.union(&i2);
    ///
    /// match union {
    ///     IntervalUnion::SingleConnected(result) => {
    ///         // The union is [0, +inf)
    ///         let expected: Interval<f64> = IntervalLowerClosedUpperUnbounded::new(0.0).into();
    ///         assert_eq!(result, expected);
    ///     }
    ///     _ => panic!("Expected a single connected union"),
    /// }
    /// ```
    fn union<IntervalType: IntervalTrait<RealType = Self::RealType>>(
        &self,
        other: &IntervalType,
    ) -> IntervalUnion<Self::RealType>
    where
        Interval<Self::RealType>: From<Self> + From<IntervalType>,
    {
        let self_lb = self.lower_bound_runtime();
        let self_ub = self.upper_bound_runtime();

        let other_lb = other.lower_bound_runtime();
        let other_ub = other.upper_bound_runtime();

        if self_lb == other_lb && self_ub == other_ub {
            // both intervals are identical
            IntervalUnion::SingleConnected(self.clone().into())
        } else {
            // the intervals are not identical

            let self_is_on_the_left = self_lb <= other_lb;

            let a = (self_lb, self_ub);
            let b = (other_lb, other_ub);

            // order the intervals s.t. the left one is the one with lowest lower bound
            let (left, right) = if self_is_on_the_left { (a, b) } else { (b, a) };

            // check if the two intervals overlap
            let upper_bound_left = &left.1;
            let lower_bound_right = &right.0;

            let intervals_touch_or_overlap = match (upper_bound_left, lower_bound_right) {
                (None, _) | (_, None) => true,
                // when the two bounds are open, there is overlap or touch
                // only when the upper bound of the interval on the left
                // is strictly bigger than the lower bound of the interval on the right
                (Some(UpperBoundRuntime::Open(ub)), Some(LowerBoundRuntime::Open(lb))) => {
                    ub.as_ref() > lb.as_ref()
                }
                // when a bound is closed there is a touch also when the values of bounds are equal
                (Some(UpperBoundRuntime::Open(ub)), Some(LowerBoundRuntime::Closed(lb))) => {
                    ub.as_ref() >= lb.as_ref()
                }
                (Some(UpperBoundRuntime::Closed(ub)), Some(LowerBoundRuntime::Open(lb))) => {
                    ub.as_ref() >= lb.as_ref()
                }
                (Some(UpperBoundRuntime::Closed(ub)), Some(LowerBoundRuntime::Closed(lb))) => {
                    ub.as_ref() >= lb.as_ref()
                }
            };

            if intervals_touch_or_overlap {
                let lower_bound = left.0;
                let upper_bound = match (left.1, right.1) {
                    (None, _) | (_, None) => {
                        // the upper bound is unbounded
                        None
                    }
                    (Some(ub_left), Some(ub_right)) => {
                        Some(crate::bounds::max_upper_bound(ub_left, ub_right))
                    }
                };

                IntervalUnion::SingleConnected(
                    Interval::try_new(lower_bound, upper_bound)
                        .expect("Failed to create interval for union"),
                )
            } else if self_is_on_the_left {
                IntervalUnion::TwoDisjoint {
                    left: self.clone().into(),
                    right: other.clone().into(),
                }
            } else {
                IntervalUnion::TwoDisjoint {
                    left: other.clone().into(),
                    right: self.clone().into(),
                }
            }
        }
    }

    /// Computes the set difference of this interval with another interval.
    ///
    /// Returns the set of all points that are in `self` but not in `other`, denoted as `A \ B`
    /// or `A - B`. The result can be:
    /// - **`None`**: If `other` completely contains `self` (empty difference)
    /// - **`Some(Single)`**: If `other` doesn't intersect or partially overlaps with `self`
    /// - **`Some(TwoDisjoint)`**: If `other` splits `self` into two parts
    ///
    /// # Mathematical Definition
    ///
    /// For intervals A and B, the set difference A \ B is defined as:
    /// ```text
    /// A \ B = { x : x ∈ A ∧ x ∉ B }
    /// ```
    ///
    /// # Properties
    ///
    /// - **Non-commutativity**: `A \ B ≠ B \ A` in general
    /// - **Non-associativity**: `(A \ B) \ C ≠ A \ (B \ C)` in general
    /// - **Identity**: `A \ ∅ = A` (returns `Some(Single(A))`) and `A \ A = ∅` (returns `None`)
    /// - **Complement-like**: When B ⊇ A, then `A \ B = ∅` (returns `None`)
    ///
    /// # Boundary Logic
    ///
    /// The difference operation carefully handles boundary points:
    /// - If `self` includes a boundary point but `other` excludes it, the point remains
    /// - If both include or both exclude a boundary point, standard set logic applies
    /// - Boundary transitions create appropriate open/closed bounds in the result
    ///
    /// # Examples
    ///
    /// ## No Overlap (Result: Original Interval)
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let a = IntervalClosed::new(0.0, 2.0);    // [0, 2]
    /// let b = IntervalClosed::new(3.0, 5.0);    // [3, 5]
    ///
    /// if let Some(diff) = a.difference(&b) {
    ///     match diff {
    ///         IntervalDifference::SingleConnected(result) => {
    ///             // A \ B = [0, 2] since they don't overlap
    ///             let expected: Interval<f64> = IntervalClosed::new(0.0, 2.0).into();
    ///             assert_eq!(result, expected);
    ///         }
    ///         _ => panic!("Expected single interval"),
    ///     }
    /// }
    /// ```
    ///
    /// ## Partial Overlap (Result: Single Interval)
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let a = IntervalClosed::new(0.0, 3.0);    // [0, 3]
    /// let b = IntervalClosed::new(2.0, 5.0);    // [2, 5]
    ///
    /// if let Some(diff) = a.difference(&b) {
    ///     match diff {
    ///         IntervalDifference::SingleConnected(result) => {
    ///             // A \ B = [0, 2) - removes overlap [2, 3]
    ///             let expected: Interval<f64> = IntervalLowerClosedUpperOpen::new(0.0, 2.0).into();
    ///             assert_eq!(result, expected);
    ///         }
    ///         _ => panic!("Expected single interval"),
    ///     }
    /// }
    /// ```
    ///
    /// ## Complete Containment (Result: None)
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let a = IntervalClosed::new(1.0, 2.0);    // [1, 2]
    /// let b = IntervalClosed::new(0.0, 3.0);    // [0, 3]
    ///
    /// let diff = a.difference(&b);
    /// assert!(diff.is_none());  // A \ B = ∅ since B contains A
    /// ```
    ///
    /// ## Interior Removal (Result: Two Disjoint Intervals)
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let a = IntervalClosed::new(0.0, 4.0);    // [0, 4]
    /// let b = IntervalOpen::new(1.0, 3.0);      // (1, 3)
    ///
    /// if let Some(diff) = a.difference(&b) {
    ///     match diff {
    ///         IntervalDifference::TwoDisjoint { left, right } => {
    ///             // A \ B = [0, 1] ∪ [3, 4] - middle part removed
    ///             let expected_left: Interval<f64> = IntervalClosed::new(0.0, 1.0).into();
    ///             let expected_right: Interval<f64> = IntervalClosed::new(3.0, 4.0).into();
    ///             assert_eq!(left, expected_left);
    ///             assert_eq!(right, expected_right);
    ///         }
    ///         _ => panic!("Expected two disjoint intervals"),
    ///     }
    /// }
    /// ```
    ///
    /// ## Boundary Subtlety with Open/Closed Bounds
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let a = IntervalClosed::new(0.0, 2.0);        // [0, 2]
    /// let b = IntervalLowerOpenUpperClosed::new(1.0, 2.0); // (1, 2]
    ///
    /// if let Some(diff) = a.difference(&b) {
    ///     match diff {
    ///         IntervalDifference::SingleConnected(result) => {
    ///             // A \ B = [0, 1] - boundary point 1.0 is kept since b excludes it
    ///             let expected: Interval<f64> = IntervalClosed::new(0.0, 1.0).into();
    ///             assert_eq!(result, expected);
    ///         }
    ///         _ => panic!("Expected single interval"),
    ///     }
    /// }
    /// ```
    ///
    /// ## With Unbounded Intervals
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let a = IntervalLowerClosedUpperUnbounded::new(0.0); // [0, +∞)
    /// let b = IntervalClosed::new(5.0, 10.0);              // [5, 10]
    ///
    /// if let Some(diff) = a.difference(&b) {
    ///     match diff {
    ///         IntervalDifference::TwoDisjoint { left, right } => {
    ///             // A \ B = [0, 5) ∪ (10, +∞)
    ///             let expected_left: Interval<f64> =
    ///                 IntervalLowerClosedUpperOpen::new(0.0, 5.0).into();
    ///             let expected_right: Interval<f64> =
    ///                 IntervalLowerOpenUpperUnbounded::new(10.0).into();
    ///             assert_eq!(left, expected_left);
    ///             assert_eq!(right, expected_right);
    ///         }
    ///         _ => panic!("Expected two disjoint intervals"),
    ///     }
    /// }
    /// ```
    ///
    /// # Performance
    ///
    /// This operation is O(1) in computation time. The result may require allocation
    /// for the returned interval(s), but the difference logic itself is constant time.
    ///
    /// # Return Type
    ///
    /// Returns `Option<IntervalDifference>`:
    /// - `None`: No points remain after subtraction (empty difference)
    /// - `Some(IntervalDifference::SingleConnected)`: Result is a single interval
    /// - `Some(IntervalDifference::TwoDisjoint)`: Result is two separate intervals
    ///
    /// # See Also
    ///
    /// - [`intersection`](Self::intersection) - Compute overlap between intervals
    /// - [`union`](Self::union) - Combine two intervals
    /// - [`contains_interval`](Contains::contains_interval) - Test if one interval contains another
    fn difference<IntervalType: IntervalTrait<RealType = Self::RealType>>(
        &self,
        other: &IntervalType,
    ) -> Option<IntervalDifference<Self::RealType>>
    where
        Interval<Self::RealType>: From<Self> + From<IntervalType>,
    {
        // Check if there's any intersection at all
        let intersection = self.intersection(other);

        if intersection.is_none() {
            // No overlap: A \ B = A
            return Some(IntervalDifference::SingleConnected(self.clone().into()));
        }

        // Get bounds of both intervals
        let self_lb = self.lower_bound_runtime();
        let self_ub = self.upper_bound_runtime();
        let other_lb = other.lower_bound_runtime();
        let other_ub = other.upper_bound_runtime();

        // Check if other completely contains self
        let other_contains_self_lower = match (&other_lb, &self_lb) {
            (None, _) => true,  // other is unbounded below, contains any lower bound
            (_, None) => false, // self is unbounded below, cannot be contained
            (Some(other_lb_val), Some(self_lb_val)) => other_lb_val <= self_lb_val,
        };

        let other_contains_self_upper = match (&other_ub, &self_ub) {
            (None, _) => true,  // other is unbounded above, contains any upper bound
            (_, None) => false, // self is unbounded above, cannot be contained
            (Some(other_ub_val), Some(self_ub_val)) => other_ub_val >= self_ub_val,
        };

        if other_contains_self_lower && other_contains_self_upper {
            // other completely contains self: A \ B = ∅
            return None;
        }

        // Determine which parts of self remain after subtracting other
        let has_left_part = match (&self_lb, &other_lb) {
            (None, None) => false,    // Both unbounded below, no left part
            (None, Some(_)) => true,  // self unbounded, other bounded: left part exists
            (Some(_), None) => false, // other unbounded below, covers entire left
            (Some(self_lb_val), Some(other_lb_val)) => self_lb_val < other_lb_val,
        };

        let has_right_part = match (&self_ub, &other_ub) {
            (None, None) => false,    // Both unbounded above, no right part
            (None, Some(_)) => true,  // self unbounded, other bounded: right part exists
            (Some(_), None) => false, // other unbounded above, covers entire right
            (Some(self_ub_val), Some(other_ub_val)) => self_ub_val > other_ub_val,
        };

        match (has_left_part, has_right_part) {
            (false, false) => {
                unreachable!(
                    "Invariant violation: if no intersection exists, we return early at line ~1161; \
                     if other completely contains self, we return None at line ~1188. \
                     Therefore, at least one of has_left_part or has_right_part must be true."
                )
            }
            (true, false) => {
                // Only left part remains: [self.lb, other.lb) or variant
                // Invariant: self.lb < other.lb (from has_left_part=true)
                // We flip other's lower bound to create the complementary upper bound
                // If other has [b (closed), we need b) (open) to exclude what other includes
                // If other has (b (open), we need b] (closed) to include what other excludes
                let left_ub = other_lb
                    .as_ref()
                    .map(|lb| lb.clone().flip_bound_side_and_type());

                // Use Interval::try_new() to handle edge cases like singletons (e.g., [2.0, 2.0])
                Some(IntervalDifference::SingleConnected(
                    Interval::try_new(self_lb, left_ub)
                        .expect("Failed to create interval difference"),
                ))
            }
            (false, true) => {
                // Only right part remains: (other.ub, self.ub] or variant
                // Invariant: other.ub < self.ub (from has_right_part=true)
                // We flip other's upper bound to create the complementary lower bound
                // If other has b] (closed), we need (b (open) to exclude what other includes
                // If other has b) (open), we need [b (closed) to include what other excludes
                let right_lb = other_ub
                    .as_ref()
                    .map(|ub| ub.clone().flip_bound_side_and_type());

                Some(IntervalDifference::SingleConnected(
                    Interval::try_new(right_lb, self_ub)
                        .expect("Failed to create interval difference"),
                ))
            }
            (true, true) => {
                // Both left and right parts remain: [self.lb, other.lb) ∪ (other.ub, self.ub]
                // Invariants:
                // - self.lb < other.lb (from has_left_part=true)
                // - other.ub < self.ub (from has_right_part=true)

                // Flip both bounds to create complementary boundaries
                // If other has [b, we need b); if other has (b, we need b]
                let left_ub = other_lb
                    .as_ref()
                    .map(|lb| lb.clone().flip_bound_side_and_type());
                let right_lb = other_ub
                    .as_ref()
                    .map(|ub| ub.clone().flip_bound_side_and_type());

                // Both intervals are guaranteed to be valid because:
                // - has_left_part = true guarantees self.lb < other.lb, so left interval is non-degenerate
                // - has_right_part = true guarantees other.ub < self.ub, so right interval is non-degenerate
                let left = Interval::try_new(self_lb.clone(), left_ub)
                    .expect("Left interval is guaranteed valid when has_left_part = true");
                let right = Interval::try_new(right_lb, self_ub.clone())
                    .expect("Right interval is guaranteed valid when has_right_part = true");
                Some(IntervalDifference::TwoDisjoint { left, right })
            }
        }
    }
}
//------------------------------------------------------------------------------------------------

//------------------------------------------------------------------------------------------------
/// Represents the union of two intervals.
///
/// The union of two intervals `A ∪ B` contains all points that belong to either `A` or `B`
/// (or both). Depending on the relationship between the intervals, the result can be:
///
/// - **[`SingleConnected`]**: If the intervals overlap or touch, their union is a single
///   continuous interval
/// - **[`TwoDisjoint`]**: If the intervals are separated by a gap, the union consists of
///   two disjoint intervals
///
/// ## Mathematical Foundation
///
/// For intervals A and B, the union is defined as:
/// ```text
/// A ∪ B = { x : x ∈ A ∨ x ∈ B }
/// ```
///
/// ## Properties
///
/// - **Commutativity**: `A ∪ B = B ∪ A`
/// - **Associativity**: `(A ∪ B) ∪ C = A ∪ (B ∪ C)`
/// - **Identity**: `A ∪ ∅ = A`
/// - **Idempotence**: `A ∪ A = A`
///
/// ## Visual Examples
///
/// ```text
/// Case 1: Overlapping intervals → SingleConnected
///   A:   [------]
///   B:       [------]
///   A∪B: [----------]
///
/// Case 2: Touching intervals → SingleConnected
///   A:   [-----)
///   B:        [------]
///   A∪B: [----------]
///
/// Case 3: Disjoint intervals → TwoDisjoint
///   A:   [----]
///   B:            [----]
///   A∪B: [----]   [----]
/// ```
///
/// ## Usage Example
///
/// ```rust
/// use grid1d::intervals::*;
///
/// // Overlapping intervals produce a single connected union
/// let a = IntervalClosed::new(0.0, 2.0);
/// let b = IntervalClosed::new(1.0, 3.0);
/// let union = a.union(&b);
///
/// match union {
///     IntervalUnion::SingleConnected(interval) => {
///         println!("Union is a single interval: {:?}", interval);
///     }
///     IntervalUnion::TwoDisjoint { left, right } => {
///         println!("Union is two intervals: {:?} and {:?}", left, right);
///     }
/// }
/// ```
///
/// ## See Also
///
/// - [`IntervalDifference`] - Result of set difference operation
/// - [`union`](IntervalOperations::union) - Method to compute unions
/// - [`gap`](IntervalUnion::gap) - Get the gap between disjoint intervals
///
/// [`SingleConnected`]: IntervalUnion::SingleConnected
/// [`TwoDisjoint`]: IntervalUnion::TwoDisjoint
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
#[serde(bound(deserialize = "RealType: for<'a> Deserialize<'a>"))]
pub enum IntervalUnion<RealType: RealScalar> {
    /// A single connected interval resulting from overlapping or touching intervals.
    ///
    /// This variant represents the case where the union of two intervals produces
    /// one continuous interval with no gaps.
    SingleConnected(Interval<RealType>),

    /// Two disjoint intervals separated by a gap.
    ///
    /// This variant represents the case where the intervals are separated and cannot
    /// be merged into a single continuous interval. Use [`gap()`](IntervalUnion::gap)
    /// to compute the interval representing the gap between them.
    TwoDisjoint {
        /// The left interval (with smaller bounds)
        left: Interval<RealType>,
        /// The right interval (with larger bounds)  
        right: Interval<RealType>,
    },
}

//------------------------------------------------------------------------------------------------
/// Represents the result of a set difference operation between two intervals.
///
/// The difference `A \ B` (read as "A minus B" or "A without B") contains all points
/// that belong to `A` but not to `B`. Depending on how the intervals overlap, the
/// result can be:
///
/// - **[`SingleConnected`]**: The result is a single continuous interval
/// - **[`TwoDisjoint`]**: Removing `B` creates a gap in `A`, resulting in two separate intervals
/// - **`None`**: If `B` completely contains `A`, the difference is empty
///
/// ## Mathematical Foundation
///
/// For intervals A and B, the set difference is defined as:
/// ```text
/// A \ B = { x : x ∈ A ∧ x ∉ B }
/// ```
///
/// ## Properties
///
/// - **Non-commutativity**: `A \ B ≠ B \ A` in general
/// - **Identity**: `A \ ∅ = A` and `A \ A = ∅`
/// - **Complement**: When `B ⊇ A`, then `A \ B = ∅` (returns `None`)
/// - **Distributivity**: `A \ (B ∪ C) = (A \ B) ∩ (A \ C)`
///
/// ## Visual Examples
///
/// ```text
/// Case 1: No overlap → SingleConnected (original interval)
///   A:     [----]
///   B:           [----]
///   A\B:   [----]
///
/// Case 2: Partial overlap → SingleConnected (truncated)
///   A:   [--------]
///   B:        [--------]
///   A\B: [----)
///
/// Case 3: Interior removal → TwoDisjoint
///   A:   [------------]
///   B:      [----]
///   A\B: [--)    (----]
///
/// Case 4: Complete containment → None (empty)
///   A:     [----]
///   B:   [----------]
///   A\B:   (empty)
/// ```
///
/// ## Usage Example
///
/// ```rust
/// use grid1d::intervals::*;
///
/// let a = IntervalClosed::new(0.0, 4.0);
/// let b = IntervalOpen::new(1.0, 3.0);
///
/// if let Some(diff) = a.difference(&b) {
///     match diff {
///         IntervalDifference::SingleConnected(interval) => {
///             println!("Difference is a single interval: {:?}", interval);
///         }
///         IntervalDifference::TwoDisjoint { left, right } => {
///             println!("Difference is two intervals: {:?} and {:?}", left, right);
///         }
///     }
/// } else {
///     println!("No points remain (empty difference)");
/// }
/// ```
///
/// ## See Also
///
/// - [`IntervalUnion`] - Result of union operation
/// - [`difference`](IntervalOperations::difference) - Method to compute differences
/// - [`to_intervals`](IntervalDifference::to_intervals) - Convert to vector of intervals
///
/// [`SingleConnected`]: IntervalDifference::SingleConnected
/// [`TwoDisjoint`]: IntervalDifference::TwoDisjoint
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
#[serde(bound(deserialize = "RealType: for<'a> Deserialize<'a>"))]
pub enum IntervalDifference<RealType: RealScalar> {
    /// A single connected interval resulting from the difference operation.
    ///
    /// This variant represents cases where removing `B` from `A` leaves a
    /// single continuous interval (including the case where `B` doesn't overlap `A` at all).
    SingleConnected(Interval<RealType>),

    /// Two disjoint intervals resulting from removing an interior portion.
    ///
    /// This variant represents the case where `B` removes a middle section of `A`,
    /// creating two separate intervals. Use [`to_intervals()`](IntervalDifference::to_intervals)
    /// to collect both intervals into a vector.
    TwoDisjoint {
        /// The left interval (with smaller bounds)
        left: Interval<RealType>,
        /// The right interval (with larger bounds)
        right: Interval<RealType>,
    },
}

// Use duplicate_item to generate common implementations for both enums
#[duplicate_item(
    EnumName             contains_doc;
    [IntervalUnion]      ["contained in this union"];
    [IntervalDifference] ["contained in this difference result"];
)]
impl<RealType: RealScalar> EnumName<RealType> {
    /// Returns `true` if the result is a single connected interval.
    pub fn is_single_connected(&self) -> bool {
        matches!(self, Self::SingleConnected(_))
    }

    /// Returns `true` if the result is two disjoint intervals.
    pub fn is_two_disjoint(&self) -> bool {
        matches!(self, Self::TwoDisjoint { .. })
    }

    /// Returns the single interval if the result is a single connected interval.
    pub fn as_single_connected(&self) -> Option<&Interval<RealType>> {
        match self {
            Self::SingleConnected(interval) => Some(interval),
            _ => None,
        }
    }

    /// Returns the two intervals if the result is two disjoint intervals.
    pub fn as_two_disjoint(&self) -> Option<(&Interval<RealType>, &Interval<RealType>)> {
        match self {
            Self::TwoDisjoint { left, right } => Some((left, right)),
            _ => None,
        }
    }

    #[doc = concat!("Checks if a point is ", contains_doc, ".")]
    pub fn contains_point(&self, x: &RealType) -> bool {
        match self {
            Self::SingleConnected(interval) => interval.contains_point(x),
            Self::TwoDisjoint { left, right } => left.contains_point(x) || right.contains_point(x),
        }
    }

    /// For disjoint intervals, returns the gap between them as an interval.
    /// Returns `None` if the union is a single connected interval.
    ///
    /// The gap represents the space between the two disjoint intervals,
    /// with boundary types "flipped" to exclude the endpoints that belong to the original intervals.
    pub fn gap(&self) -> Option<IntervalFinitePositiveLength<RealType>> {
        match self {
            Self::SingleConnected(_) => None,
            Self::TwoDisjoint { left, right } => {
                // The gap is from left.upper_bound to right.lower_bound
                let left_ub = left
                    .upper_bound_runtime()
                    .expect("The left interval in a disjoint union cannot be unbounded above");
                let right_lb = right
                    .lower_bound_runtime()
                    .expect("The right interval in a disjoint union cannot be unbounded below");

                // Flip the boundary types and sides to create the gap interval
                // The gap excludes the endpoints of the touching intervals
                // If left ends with a], gap starts with (a
                // If left ends with a), gap starts with [a
                // If right starts with [b, gap ends with b)
                // If right starts with (b, gap ends with b]
                let lb = left_ub.flip_bound_side_and_type();
                let ub = right_lb.flip_bound_side_and_type();

                // Create the appropriate interval type based on boundary types
                Some(
                    IntervalFinitePositiveLength::try_new(lb, ub)
                        .expect("Failed to create the interval representing the gap!"),
                )
            }
        }
    }

    /// Collects all intervals in the difference result into a vector.
    ///
    /// This is useful when you need to process all resulting intervals uniformly,
    /// regardless of whether the result is single connected or two disjoint.
    ///
    /// # Returns
    ///
    /// - Single-element vector for `SingleConnected`
    /// - Two-element vector for `TwoDisjoint`
    ///
    /// # Examples
    ///
    /// ```rust
    /// use grid1d::intervals::*;
    ///
    /// let a = IntervalClosed::new(0.0, 4.0);
    /// let b = IntervalOpen::new(1.0, 3.0);
    ///
    /// if let Some(diff) = a.difference(&b) {
    ///     let intervals = diff.to_intervals();
    ///     assert_eq!(intervals.len(), 2);
    ///     for interval in intervals {
    ///         println!("Interval: {:?}", interval);
    ///     }
    /// }
    /// ```
    pub fn to_intervals(self) -> Vec<Interval<RealType>> {
        match self {
            Self::SingleConnected(interval) => vec![interval],
            Self::TwoDisjoint { left, right } => vec![left, right],
        }
    }
}

// Additional documentation for IntervalDifference (generated above by duplicate_item)
/// ## IntervalDifference Documentation
///
/// Represents the result of a set difference operation between two intervals.
///
/// When computing `A \ B` (read as "A minus B" or "A without B"), the result can take
/// one of two forms depending on how the intervals overlap:
///
/// - **SingleConnected**: The result is a single continuous interval
/// - **TwoDisjoint**: B removes a middle section of A, creating two separate intervals
///
/// If the result is empty (B completely contains A), the `difference` method returns `None`.
///
/// ## Mathematical Foundation
///
/// For intervals A and B, the set difference is defined as:
/// ```text
/// A \ B = { x : x ∈ A ∧ x ∉ B }
/// ```
///
/// ## Visual Examples
///
/// ```text
/// Case 1: Empty (B contains A) → None
///   A:     [----]
///   B:   [----------]
///   A\B:   (empty)
///
/// Case 2: SingleConnected (partial overlap) → Some(SingleConnected)
///   A:   [--------]
///   B:        [--------]
///   A\B: [----)
///
/// Case 3: TwoDisjoint (B in middle of A) → Some(TwoDisjoint)
///   A:   [------------]
///   B:      [----]
///   A\B: [--)    (----]
/// ```
///
/// ## Properties
///
/// - **Non-commutativity**: `A \ B ≠ B \ A` in general
/// - **Identity**: `A \ ∅ = A` and `A \ A = ∅` (returns `None`)
/// - **Complement**: When B ⊇ A, then `A \ B = ∅` (returns `None`)
///
/// ## Usage Example
///
/// ```rust
/// use grid1d::intervals::*;
///
/// let a = IntervalClosed::new(0.0, 4.0);
/// let b = IntervalOpen::new(1.0, 3.0);
///
/// if let Some(result) = a.difference(&b) {
///     match result {
///         IntervalDifference::SingleConnected(interval) => {
///             println!("Result is a single interval: {:?}", interval);
///         }
///         IntervalDifference::TwoDisjoint { left, right } => {
///             println!("Result is two intervals: {:?} and {:?}", left, right);
///         }
///     }
/// } else {
///     println!("No points remain (empty difference)");
/// }
/// ```
// Specialized implementation for IntervalDifference only
impl<RealType: RealScalar> IntervalDifference<RealType> {}
//------------------------------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{
        SubIntervalInPartition,
        intervals::{
            GetLowerBoundValue, GetUpperBoundValue, Interval, IntervalClosed,
            IntervalFinitePositiveLengthTrait, IntervalLowerClosedUpperOpen,
            IntervalLowerClosedUpperUnbounded, IntervalLowerOpenUpperClosed,
            IntervalLowerOpenUpperUnbounded, IntervalLowerUnboundedUpperClosed,
            IntervalLowerUnboundedUpperOpen, IntervalLowerUnboundedUpperUnbounded, IntervalOpen,
            IntervalSingleton,
        },
    };
    use num_valid::scalars::PositiveRealScalar;
    use try_create::{New, TryNew};

    mod difference {
        use super::*;

        #[test]
        fn difference_no_overlap() {
            // A = [0, 2], B = [3, 5]
            // A \ B = [0, 2]
            let a = IntervalClosed::new(0.0, 2.0);
            let b = IntervalClosed::new(3.0, 5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(result) => {
                    let expected: Interval<f64> = IntervalClosed::new(0.0, 2.0).into();
                    assert_eq!(result, expected);
                }
                _ => panic!("Expected single interval"),
            }
        }

        #[test]
        fn difference_partial_overlap_right() {
            // A = [0, 3], B = [2, 5]
            // A \ B = [0, 2)
            let a = IntervalClosed::new(0.0, 3.0);
            let b = IntervalClosed::new(2.0, 5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(result) => {
                    let expected: Interval<f64> =
                        IntervalLowerClosedUpperOpen::new(0.0, 2.0).into();
                    assert_eq!(result, expected);
                }
                _ => panic!("Expected single interval"),
            }
        }

        #[test]
        fn difference_partial_overlap_left() {
            // A = [2, 5], B = [0, 3]
            // A \ B = (3, 5]
            let a = IntervalClosed::new(2.0, 5.0);
            let b = IntervalClosed::new(0.0, 3.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(result) => {
                    let expected: Interval<f64> =
                        IntervalLowerOpenUpperClosed::new(3.0, 5.0).into();
                    assert_eq!(result, expected);
                }
                _ => panic!("Expected single interval"),
            }
        }

        #[test]
        fn difference_complete_containment() {
            // A = [1, 2], B = [0, 3]
            // A \ B = ∅
            let a = IntervalClosed::new(1.0, 2.0);
            let b = IntervalClosed::new(0.0, 3.0);
            let diff = a.difference(&b);

            assert!(diff.is_none());
        }

        #[test]
        fn difference_interior_removal() {
            // A = [0, 4], B = (1, 3)
            // A \ B = [0, 1] ∪ [3, 4]
            let a = IntervalClosed::new(0.0, 4.0);
            let b = IntervalOpen::new(1.0, 3.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> = IntervalClosed::new(0.0, 1.0).into();
                    let expected_right: Interval<f64> = IntervalClosed::new(3.0, 4.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals"),
            }
        }

        #[test]
        fn difference_boundary_subtlety_open_closed() {
            // A = [0, 2], B = (1, 2]
            // A \ B = [0, 1] (boundary point 1.0 is kept)
            let a = IntervalClosed::new(0.0, 2.0);
            let b = IntervalLowerOpenUpperClosed::new(1.0, 2.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(result) => {
                    let expected: Interval<f64> = IntervalClosed::new(0.0, 1.0).into();
                    assert_eq!(result, expected);
                }
                _ => panic!("Expected single interval"),
            }
        }

        #[test]
        fn difference_boundary_subtlety_closed_open() {
            // A = [0, 2], B = [1, 2)
            // A \ B = [0, 1) ∪ {2} but since {2} = [2, 2], result is [0, 1) ∪ [2, 2]
            let a = IntervalClosed::new(0.0, 2.0);
            let b = IntervalLowerClosedUpperOpen::new(1.0, 2.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> =
                        IntervalLowerClosedUpperOpen::new(0.0, 1.0).into();
                    let expected_right: Interval<f64> = IntervalSingleton::new(2.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals"),
            }
        }

        #[test]
        fn difference_with_unbounded_intervals() {
            // A = [0, +∞), B = [5, 10]
            // A \ B = [0, 5) ∪ (10, +∞)
            let a = IntervalLowerClosedUpperUnbounded::new(0.0);
            let b = IntervalClosed::new(5.0, 10.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> =
                        IntervalLowerClosedUpperOpen::new(0.0, 5.0).into();
                    let expected_right: Interval<f64> =
                        IntervalLowerOpenUpperUnbounded::new(10.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals"),
            }
        }

        #[test]
        fn difference_unbounded_minus_unbounded() {
            // A = (-∞, +∞), B = [0, 10]
            // A \ B = (-∞, 0) ∪ (10, +∞)
            let a = IntervalLowerUnboundedUpperUnbounded::new();
            let b = IntervalClosed::new(0.0, 10.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> =
                        IntervalLowerUnboundedUpperOpen::new(0.0).into();
                    let expected_right: Interval<f64> =
                        IntervalLowerOpenUpperUnbounded::new(10.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals"),
            }
        }

        #[test]
        fn difference_unbounded_complete_containment() {
            // A = [5, 10], B = (-∞, +∞)
            // A \ B = ∅
            let a = IntervalClosed::new(5.0, 10.0);
            let b = IntervalLowerUnboundedUpperUnbounded::new();
            let diff = a.difference(&b);

            assert!(diff.is_none());
        }

        #[test]
        fn difference_unbounded_no_overlap() {
            // A = [0, 5], B = [10, +∞)
            // A \ B = [0, 5]
            let a = IntervalClosed::new(0.0, 5.0);
            let b = IntervalLowerClosedUpperUnbounded::new(10.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(result) => {
                    let expected: Interval<f64> = IntervalClosed::new(0.0, 5.0).into();
                    assert_eq!(result, expected);
                }
                _ => panic!("Expected single interval"),
            }
        }

        #[test]
        fn difference_with_singleton() {
            // A = [0, 2], B = {1}
            // A \ B = [0, 1) ∪ (1, 2]
            let a = IntervalClosed::new(0.0, 2.0);
            let b = IntervalSingleton::new(1.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> =
                        IntervalLowerClosedUpperOpen::new(0.0, 1.0).into();
                    let expected_right: Interval<f64> =
                        IntervalLowerOpenUpperClosed::new(1.0, 2.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals"),
            }
        }

        #[test]
        fn difference_singleton_minus_interval() {
            // A = {1}, B = [0, 2]
            // A \ B = ∅
            let a = IntervalSingleton::new(1.0);
            let b = IntervalClosed::new(0.0, 2.0);
            let diff = a.difference(&b);

            assert!(diff.is_none());
        }

        #[test]
        fn difference_singleton_minus_non_overlapping() {
            // A = {5}, B = [0, 2]
            // A \ B = {5}
            let a = IntervalSingleton::new(5.0);
            let b = IntervalClosed::new(0.0, 2.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(result) => {
                    let expected: Interval<f64> = IntervalSingleton::new(5.0).into();
                    assert_eq!(result, expected);
                }
                _ => panic!("Expected single interval"),
            }
        }

        #[test]
        fn difference_identical_intervals() {
            // A = [0, 1], B = [0, 1]
            // A \ B = ∅
            let a = IntervalClosed::new(0.0, 1.0);
            let b = IntervalClosed::new(0.0, 1.0);
            let diff = a.difference(&b);

            assert!(diff.is_none());
        }

        #[test]
        fn difference_open_vs_closed_same_values() {
            // A = [0, 2], B = (0, 2)
            // A \ B = {0} ∪ {2} = [0, 0] ∪ [2, 2]
            let a = IntervalClosed::new(0.0, 2.0);
            let b = IntervalOpen::new(0.0, 2.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> = IntervalSingleton::new(0.0).into();
                    let expected_right: Interval<f64> = IntervalSingleton::new(2.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals (singletons at boundaries)"),
            }
        }

        #[test]
        fn difference_contains_point_test() {
            // A = [0, 4], B = (1, 3)
            // A \ B = [0, 1] ∪ [3, 4]
            let a = IntervalClosed::new(0.0, 4.0);
            let b = IntervalOpen::new(1.0, 3.0);
            let diff = a.difference(&b).unwrap();

            // Points in left part
            assert!(diff.contains_point(&0.0));
            assert!(diff.contains_point(&0.5));
            assert!(diff.contains_point(&1.0));

            // Points removed by B
            assert!(!diff.contains_point(&1.5));
            assert!(!diff.contains_point(&2.0));
            assert!(!diff.contains_point(&2.5));

            // Points in right part
            assert!(diff.contains_point(&3.0));
            assert!(diff.contains_point(&3.5));
            assert!(diff.contains_point(&4.0));

            // Points outside A
            assert!(!diff.contains_point(&-1.0));
            assert!(!diff.contains_point(&5.0));
        }

        #[test]
        fn difference_to_intervals_test() {
            // Test empty (None)
            let a = IntervalClosed::new(1.0, 2.0);
            let b = IntervalClosed::new(0.0, 3.0);
            let diff = a.difference(&b);
            assert!(diff.is_none());

            // Test single
            let a = IntervalClosed::new(0.0, 2.0);
            let b = IntervalClosed::new(3.0, 5.0);
            let diff = a.difference(&b);
            assert!(diff.is_some());
            assert_eq!(diff.unwrap().to_intervals().len(), 1);

            // Test two disjoint
            let a = IntervalClosed::new(0.0, 4.0);
            let b = IntervalOpen::new(1.0, 3.0);
            let diff = a.difference(&b);
            assert!(diff.is_some());
            assert_eq!(diff.unwrap().to_intervals().len(), 2);
        }

        #[test]
        fn difference_helper_methods_test() {
            let a = IntervalClosed::new(0.0, 4.0);
            let b = IntervalOpen::new(1.0, 3.0);
            let diff = a.difference(&b).unwrap();

            assert!(!diff.is_single_connected());
            assert!(diff.is_two_disjoint());

            assert!(diff.as_single_connected().is_none());
            assert!(diff.as_two_disjoint().is_some());

            let (left, right) = diff.as_two_disjoint().unwrap();
            assert!(left.contains_point(&0.5));
            assert!(right.contains_point(&3.5));
        }

        #[test]
        fn difference_lower_unbounded_partial() {
            // A = (-∞, 5], B = [0, 3]
            // A \ B = (-∞, 0) ∪ (3, 5]
            let a = IntervalLowerUnboundedUpperClosed::new(5.0);
            let b = IntervalClosed::new(0.0, 3.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> =
                        IntervalLowerUnboundedUpperOpen::new(0.0).into();
                    let expected_right: Interval<f64> =
                        IntervalLowerOpenUpperClosed::new(3.0, 5.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals"),
            }
        }

        /// Test difference of two open intervals (no overlap case)
        #[test]
        fn difference_open_open_no_overlap() {
            // A = (0, 2), B = (3, 5)
            // A \ B = (0, 2)
            let a = IntervalOpen::new(0.0, 2.0);
            let b = IntervalOpen::new(3.0, 5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(result) => {
                    let expected: Interval<f64> = IntervalOpen::new(0.0, 2.0).into();
                    assert_eq!(result, expected);
                }
                _ => panic!("Expected single interval"),
            }
        }

        /// Test difference of two open intervals (partial overlap)
        #[test]
        fn difference_open_open_partial_overlap() {
            // A = (0, 4), B = (2, 6)
            // A \ B = (0, 2]
            let a = IntervalOpen::new(0.0, 4.0);
            let b = IntervalOpen::new(2.0, 6.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(result) => {
                    let expected: Interval<f64> =
                        IntervalLowerOpenUpperClosed::new(0.0, 2.0).into();
                    assert_eq!(result, expected);
                }
                _ => panic!("Expected single interval"),
            }
        }

        /// Test difference where result would have degenerate left interval
        /// This tests the (false, true) branch in the (true, true) case
        #[test]
        fn difference_degenerate_left_valid_right() {
            // A = [0, 4], B = [0, 2)
            // A \ B = [2, 4]
            // Left would be [0, 0) which is degenerate, only right part [2, 4] remains
            let a = IntervalClosed::new(0.0, 4.0);
            let b = IntervalLowerClosedUpperOpen::new(0.0, 2.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(result) => {
                    let expected: Interval<f64> = IntervalClosed::new(2.0, 4.0).into();
                    assert_eq!(
                        result, expected,
                        "Should only have right part, left is degenerate"
                    );
                }
                _ => panic!("Expected single interval (right part only)"),
            }
        }

        /// Test difference where result would have degenerate right interval
        /// This tests the (true, false) branch in the (true, true) case
        #[test]
        fn difference_valid_left_degenerate_right() {
            // A = [0, 4], B = (2, 4]
            // A \ B = [0, 2]
            // Right would be (4, 4] which is degenerate, only left part [0, 2] remains
            let a = IntervalClosed::new(0.0, 4.0);
            let b = IntervalLowerOpenUpperClosed::new(2.0, 4.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(result) => {
                    let expected: Interval<f64> = IntervalClosed::new(0.0, 2.0).into();
                    assert_eq!(
                        result, expected,
                        "Should only have left part, right is degenerate"
                    );
                }
                _ => panic!("Expected single interval (left part only)"),
            }
        }

        /// Test difference of two identical singletons
        #[test]
        fn difference_singleton_same_value() {
            // A = {5}, B = {5}
            // A \ B = ∅ (singleton minus itself is empty)
            let a = IntervalSingleton::new(5.0);
            let b = IntervalSingleton::new(5.0);
            let diff = a.difference(&b);

            assert!(
                diff.is_none(),
                "Expected None when subtracting identical singletons ({{5}} \\ {{5}} = ∅)"
            );
        }

        /// Test difference of two different singletons
        #[test]
        fn difference_singleton_different_values() {
            // A = {3}, B = {7}
            // A \ B = {3} (disjoint singletons)
            let a = IntervalSingleton::new(3.0);
            let b = IntervalSingleton::new(7.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> = IntervalSingleton::new(3.0).into();
                    assert_eq!(
                        interval, expected,
                        "Expected {{3}} when singletons are disjoint"
                    );
                }
                _ => panic!("Expected single singleton interval"),
            }
        }

        /// Test difference of singleton contained in interval
        #[test]
        fn difference_singleton_contained_in_interval() {
            // A = {5}, B = [0, 10]
            // A \ B = ∅ (singleton completely contained)
            let a = IntervalSingleton::new(5.0);
            let b = IntervalClosed::new(0.0, 10.0);
            let diff = a.difference(&b);

            assert!(
                diff.is_none(),
                "Expected None when singleton is contained in interval"
            );
        }

        /// Test difference of interval containing singleton at endpoint
        #[test]
        fn difference_interval_minus_singleton_at_endpoint() {
            // A = [0, 10], B = {0}
            // A \ B = (0, 10] (remove left endpoint)
            let a = IntervalClosed::new(0.0, 10.0);
            let b = IntervalSingleton::new(0.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> =
                        IntervalLowerOpenUpperClosed::new(0.0, 10.0).into();
                    assert_eq!(
                        interval, expected,
                        "Expected (0, 10] when removing left endpoint singleton"
                    );
                }
                _ => panic!("Expected single interval"),
            }
        }

        /// Test difference of interval minus singleton at midpoint
        #[test]
        fn difference_interval_minus_singleton_interior() {
            // A = [0, 10], B = {5}
            // A \ B = [0, 5) ∪ (5, 10] (split at midpoint)
            let a = IntervalClosed::new(0.0, 10.0);
            let b = IntervalSingleton::new(5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> =
                        IntervalLowerClosedUpperOpen::new(0.0, 5.0).into();
                    let expected_right: Interval<f64> =
                        IntervalLowerOpenUpperClosed::new(5.0, 10.0).into();
                    assert_eq!(left, expected_left, "Expected [0, 5) as left part");
                    assert_eq!(right, expected_right, "Expected (5, 10] as right part");
                }
                _ => panic!("Expected two disjoint intervals"),
            }
        }

        /// Test difference of open interval minus singleton at boundary
        #[test]
        fn difference_open_interval_minus_singleton_at_boundary() {
            // A = (0, 10), B = {0}
            // A \ B = (0, 10) (singleton not in open interval)
            let a = IntervalOpen::new(0.0, 10.0);
            let b = IntervalSingleton::new(0.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> = IntervalOpen::new(0.0, 10.0).into();
                    assert_eq!(
                        interval, expected,
                        "Expected (0, 10) unchanged since {{0}} is not in (0, 10)"
                    );
                }
                _ => panic!("Expected single interval"),
            }
        }

        /// Test difference of open intervals with interior removal
        #[test]
        fn difference_open_interior_removal() {
            // A = (0, 6), B = [2, 4]
            // A \ B = (0, 2) ∪ (4, 6)
            let a = IntervalOpen::new(0.0, 6.0);
            let b = IntervalClosed::new(2.0, 4.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> = IntervalOpen::new(0.0, 2.0).into();
                    let expected_right: Interval<f64> = IntervalOpen::new(4.0, 6.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint open intervals"),
            }
        }

        /// Test difference with half-open intervals (mixed boundary types)
        #[test]
        fn difference_half_open_intervals() {
            // A = [0, 5), B = (2, 4]
            // A \ B = [0, 2] ∪ (4, 5)
            let a = IntervalLowerClosedUpperOpen::new(0.0, 5.0);
            let b = IntervalLowerOpenUpperClosed::new(2.0, 4.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> = IntervalClosed::new(0.0, 2.0).into();
                    let expected_right: Interval<f64> = IntervalOpen::new(4.0, 5.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals with mixed bounds"),
            }
        }

        /// Test difference where other is unbounded on both sides
        #[test]
        fn difference_with_fully_unbounded_other() {
            // A = [0, 5], B = (-∞, +∞)
            // A \ B = ∅ (completely contained)
            let a = IntervalClosed::new(0.0, 5.0);
            let b = IntervalLowerUnboundedUpperUnbounded::new();
            let diff = a.difference(&b);

            assert!(
                diff.is_none(),
                "Finite interval minus universal interval should be empty"
            );
        }

        /// Test difference of unbounded intervals
        #[test]
        fn difference_unbounded_partial_overlap() {
            // A = [0, +∞), B = (-∞, 5]
            // A \ B = (5, +∞)
            let a = IntervalLowerClosedUpperUnbounded::new(0.0);
            let b = IntervalLowerUnboundedUpperClosed::new(5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(result) => {
                    let expected: Interval<f64> = IntervalLowerOpenUpperUnbounded::new(5.0).into();
                    assert_eq!(result, expected);
                }
                _ => panic!("Expected single unbounded interval"),
            }
        }

        #[test]
        fn difference_open_lower_unbounded_minus_finite() {
            // A = (-∞, 5), B = [0, 3]
            // A \ B = (-∞, 0) ∪ (3, 5)
            let a = IntervalLowerUnboundedUpperOpen::new(5.0);
            let b = IntervalClosed::new(0.0, 3.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> =
                        IntervalLowerUnboundedUpperOpen::new(0.0).into();
                    let expected_right: Interval<f64> = IntervalOpen::new(3.0, 5.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals"),
            }
        }

        #[test]
        fn difference_open_upper_unbounded_minus_finite() {
            // A = (0, +∞), B = [5, 10]
            // A \ B = (0, 5) ∪ (10, +∞)
            let a = IntervalLowerOpenUpperUnbounded::new(0.0);
            let b = IntervalClosed::new(5.0, 10.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> = IntervalOpen::new(0.0, 5.0).into();
                    let expected_right: Interval<f64> =
                        IntervalLowerOpenUpperUnbounded::new(10.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals"),
            }
        }

        #[test]
        fn difference_lower_unbounded_minus_open_finite() {
            // A = (-∞, 10], B = (2, 5)
            // A \ B = (-∞, 2] ∪ [5, 10]
            let a = IntervalLowerUnboundedUpperClosed::new(10.0);
            let b = IntervalOpen::new(2.0, 5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> =
                        IntervalLowerUnboundedUpperClosed::new(2.0).into();
                    let expected_right: Interval<f64> = IntervalClosed::new(5.0, 10.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals"),
            }
        }

        #[test]
        fn difference_upper_unbounded_minus_open_finite() {
            // A = [0, +∞), B = (5, 10)
            // A \ B = [0, 5] ∪ [10, +∞)
            let a = IntervalLowerClosedUpperUnbounded::new(0.0);
            let b = IntervalOpen::new(5.0, 10.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> = IntervalClosed::new(0.0, 5.0).into();
                    let expected_right: Interval<f64> =
                        IntervalLowerClosedUpperUnbounded::new(10.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals"),
            }
        }

        #[test]
        fn difference_same_direction_upper_unbounded_disjoint() {
            // A = [0, +∞), B = [10, +∞)
            // A \ B = [0, 10)
            let a = IntervalLowerClosedUpperUnbounded::new(0.0);
            let b = IntervalLowerClosedUpperUnbounded::new(10.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> =
                        IntervalLowerClosedUpperOpen::new(0.0, 10.0).into();
                    assert_eq!(interval, expected);
                }
                _ => panic!("Expected single connected interval"),
            }
        }

        #[test]
        fn difference_same_direction_upper_unbounded_overlap() {
            // A = [5, +∞), B = [0, +∞)
            // A \ B = ∅ (A ⊆ B)
            let a = IntervalLowerClosedUpperUnbounded::new(5.0);
            let b = IntervalLowerClosedUpperUnbounded::new(0.0);
            let diff = a.difference(&b);

            assert!(diff.is_none());
        }

        #[test]
        fn difference_same_direction_lower_unbounded_disjoint() {
            // A = (-∞, 10], B = (-∞, 5]
            // A \ B = (5, 10]
            let a = IntervalLowerUnboundedUpperClosed::new(10.0);
            let b = IntervalLowerUnboundedUpperClosed::new(5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> =
                        IntervalLowerOpenUpperClosed::new(5.0, 10.0).into();
                    assert_eq!(interval, expected);
                }
                _ => panic!("Expected single connected interval"),
            }
        }

        #[test]
        fn difference_same_direction_lower_unbounded_overlap() {
            // A = (-∞, 5], B = (-∞, 10]
            // A \ B = ∅ (A ⊆ B)
            let a = IntervalLowerUnboundedUpperClosed::new(5.0);
            let b = IntervalLowerUnboundedUpperClosed::new(10.0);
            let diff = a.difference(&b);

            assert!(diff.is_none());
        }

        #[test]
        fn difference_universal_minus_upper_unbounded() {
            // A = (-∞, +∞), B = [5, +∞)
            // A \ B = (-∞, 5)
            let a = IntervalLowerUnboundedUpperUnbounded::new();
            let b = IntervalLowerClosedUpperUnbounded::new(5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> = IntervalLowerUnboundedUpperOpen::new(5.0).into();
                    assert_eq!(interval, expected);
                }
                _ => panic!("Expected single connected interval"),
            }
        }

        #[test]
        fn difference_universal_minus_lower_unbounded() {
            // A = (-∞, +∞), B = (-∞, 5]
            // A \ B = (5, +∞)
            let a = IntervalLowerUnboundedUpperUnbounded::new();
            let b = IntervalLowerUnboundedUpperClosed::new(5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> = IntervalLowerOpenUpperUnbounded::new(5.0).into();
                    assert_eq!(interval, expected);
                }
                _ => panic!("Expected single connected interval"),
            }
        }

        #[test]
        fn difference_universal_minus_open_finite() {
            // A = (-∞, +∞), B = (0, 10)
            // A \ B = (-∞, 0] ∪ [10, +∞)
            let a = IntervalLowerUnboundedUpperUnbounded::new();
            let b = IntervalOpen::new(0.0, 10.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::TwoDisjoint { left, right } => {
                    let expected_left: Interval<f64> =
                        IntervalLowerUnboundedUpperClosed::new(0.0).into();
                    let expected_right: Interval<f64> =
                        IntervalLowerClosedUpperUnbounded::new(10.0).into();
                    assert_eq!(left, expected_left);
                    assert_eq!(right, expected_right);
                }
                _ => panic!("Expected two disjoint intervals"),
            }
        }

        #[test]
        fn difference_universal_minus_universal() {
            // A = (-∞, +∞), B = (-∞, +∞)
            // A \ B = ∅ (A \ A = ∅)
            let a: IntervalLowerUnboundedUpperUnbounded<f64> =
                IntervalLowerUnboundedUpperUnbounded::new();
            let b: IntervalLowerUnboundedUpperUnbounded<f64> =
                IntervalLowerUnboundedUpperUnbounded::new();
            let diff = a.difference(&b);

            assert!(diff.is_none());
        }

        #[test]
        fn difference_finite_closed_minus_lower_unbounded_disjoint() {
            // A = [10, 20], B = (-∞, 5]
            // A \ B = [10, 20] (no overlap)
            let a = IntervalClosed::new(10.0, 20.0);
            let b = IntervalLowerUnboundedUpperClosed::new(5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> = IntervalClosed::new(10.0, 20.0).into();
                    assert_eq!(interval, expected);
                }
                _ => panic!("Expected single connected interval"),
            }
        }

        #[test]
        fn difference_finite_closed_minus_lower_unbounded_overlap() {
            // A = [0, 10], B = (-∞, 5]
            // A \ B = (5, 10]
            let a = IntervalClosed::new(0.0, 10.0);
            let b = IntervalLowerUnboundedUpperClosed::new(5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> =
                        IntervalLowerOpenUpperClosed::new(5.0, 10.0).into();
                    assert_eq!(interval, expected);
                }
                _ => panic!("Expected single connected interval"),
            }
        }

        #[test]
        fn difference_finite_open_minus_upper_unbounded_disjoint() {
            // A = (0, 5), B = [10, +∞)
            // A \ B = (0, 5) (no overlap)
            let a = IntervalOpen::new(0.0, 5.0);
            let b = IntervalLowerClosedUpperUnbounded::new(10.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> = IntervalOpen::new(0.0, 5.0).into();
                    assert_eq!(interval, expected);
                }
                _ => panic!("Expected single connected interval"),
            }
        }

        #[test]
        fn difference_finite_open_minus_lower_unbounded_overlap() {
            // A = (0, 10), B = (-∞, 5]
            // A \ B = (5, 10)
            let a = IntervalOpen::new(0.0, 10.0);
            let b = IntervalLowerUnboundedUpperClosed::new(5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> = IntervalOpen::new(5.0, 10.0).into();
                    assert_eq!(interval, expected);
                }
                _ => panic!("Expected single connected interval"),
            }
        }

        #[test]
        fn difference_lower_unbounded_minus_upper_unbounded_disjoint() {
            // A = (-∞, 5], B = [10, +∞)
            // A \ B = (-∞, 5] (no overlap)
            let a = IntervalLowerUnboundedUpperClosed::new(5.0);
            let b = IntervalLowerClosedUpperUnbounded::new(10.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> =
                        IntervalLowerUnboundedUpperClosed::new(5.0).into();
                    assert_eq!(interval, expected);
                }
                _ => panic!("Expected single connected interval"),
            }
        }

        #[test]
        fn difference_lower_unbounded_minus_upper_unbounded_overlap() {
            // A = (-∞, 10], B = [5, +∞)
            // A \ B = (-∞, 5)
            let a = IntervalLowerUnboundedUpperClosed::new(10.0);
            let b = IntervalLowerClosedUpperUnbounded::new(5.0);
            let diff = a.difference(&b);

            assert!(diff.is_some());
            match diff.unwrap() {
                IntervalDifference::SingleConnected(interval) => {
                    let expected: Interval<f64> = IntervalLowerUnboundedUpperOpen::new(5.0).into();
                    assert_eq!(interval, expected);
                }
                _ => panic!("Expected single connected interval"),
            }
        }

        #[test]
        fn difference_none_half_open_contained() {
            // A = [1, 5), B = [0, 10]
            // A \ B = ∅ (A is completely contained in B)
            let a = IntervalLowerClosedUpperOpen::new(1.0, 5.0);
            let b = IntervalClosed::new(0.0, 10.0);
            let diff = a.difference(&b);

            assert!(
                diff.is_none(),
                "Expected None when interval is completely contained"
            );
        }

        #[test]
        fn difference_none_open_interval_contained() {
            // A = (2, 4), B = [2, 4]
            // A \ B = ∅ (open interval contained in closed interval)
            let a = IntervalOpen::new(2.0, 4.0);
            let b = IntervalClosed::new(2.0, 4.0);
            let diff = a.difference(&b);

            assert!(
                diff.is_none(),
                "Expected None when open interval is contained in closed interval"
            );
        }

        #[test]
        fn difference_none_identical_intervals() {
            // A = [3, 7], B = [3, 7]
            // A \ B = ∅ (identical intervals)
            let a = IntervalClosed::new(3.0, 7.0);
            let b = IntervalClosed::new(3.0, 7.0);
            let diff = a.difference(&b);

            assert!(diff.is_none(), "Expected None when A = B (A \\ A = ∅)");
        }

        #[test]
        fn difference_none_open_interval_subset() {
            // A = (1, 3), B = (0, 5)
            // A \ B = ∅ (A ⊂ B)
            let a = IntervalOpen::new(1.0, 3.0);
            let b = IntervalOpen::new(0.0, 5.0);
            let diff = a.difference(&b);

            assert!(
                diff.is_none(),
                "Expected None when open interval is subset of another"
            );
        }

        #[test]
        fn difference_none_lower_unbounded_contained() {
            // A = (-∞, 5], B = (-∞, 10]
            // A \ B = ∅ (A ⊆ B)
            let a = IntervalLowerUnboundedUpperClosed::new(5.0);
            let b = IntervalLowerUnboundedUpperClosed::new(10.0);
            let diff = a.difference(&b);

            assert!(
                diff.is_none(),
                "Expected None when lower unbounded interval is subset"
            );
        }

        #[test]
        fn difference_none_upper_unbounded_contained() {
            // A = [10, +∞), B = [5, +∞)
            // A \ B = ∅ (A ⊆ B)
            let a = IntervalLowerClosedUpperUnbounded::new(10.0);
            let b = IntervalLowerClosedUpperUnbounded::new(5.0);
            let diff = a.difference(&b);

            assert!(
                diff.is_none(),
                "Expected None when upper unbounded interval is subset"
            );
        }

        #[test]
        fn difference_none_mixed_boundaries_contained() {
            // A = (2, 8), B = [1, 9]
            // A \ B = ∅ (open interval contained in closed interval)
            let a = IntervalOpen::new(2.0, 8.0);
            let b = IntervalClosed::new(1.0, 9.0);
            let diff = a.difference(&b);

            assert!(
                diff.is_none(),
                "Expected None when interval with open boundaries is contained"
            );
        }
    }

    mod intersection {
        use super::*;

        #[test]
        fn intersection_closed_closed() {
            let a = IntervalClosed::new(0.0, 2.0);
            let b = IntervalClosed::new(1.0, 3.0);
            let intersection = a.intersection(&b).unwrap();
            let expected_intersection =
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::from(IntervalClosed::new(1.0, 2.0)),
                ));
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_closed_open() {
            let a = IntervalClosed::new(0.0, 2.0);
            let b = IntervalOpen::new(1.0, 3.0);
            let intersection = a.intersection(&b).unwrap();
            let expected_intersection =
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::from(IntervalLowerOpenUpperClosed::new(1.0, 2.0)),
                ));
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_open_open() {
            let a = IntervalOpen::new(0.0, 2.0);
            let b = IntervalOpen::new(1.0, 3.0);
            let intersection = a.intersection(&b).unwrap();
            let expected_intersection =
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::from(IntervalOpen::new(1.0, 2.0)),
                ));
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_closed_left_half_open() {
            let a = IntervalClosed::new(0.0, 2.0);
            let b = IntervalLowerOpenUpperClosed::new(1.0, 3.0);
            let intersection = a.intersection(&b).unwrap();
            let expected_intersection =
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::from(IntervalLowerOpenUpperClosed::new(1.0, 2.0)),
                ));
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_closed_right_half_open() {
            let a = IntervalClosed::new(0.0, 2.0);
            let b = IntervalLowerClosedUpperOpen::new(1.0, 3.0);
            let intersection = a.intersection(&b).unwrap();
            let expected_intersection =
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::from(IntervalClosed::new(1.0, 2.0)),
                ));
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_singleton_closed() {
            let a = IntervalSingleton::new(1.0);
            let b = IntervalClosed::new(0.0, 2.0);
            let intersection = a.intersection(&b).unwrap();
            let expected_intersection = Interval::FiniteLength(IntervalFiniteLength::from(a));
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_open_left_half_open() {
            let a = IntervalOpen::new(0.0, 2.0);
            let b = IntervalLowerOpenUpperClosed::new(1.0, 3.0);
            let intersection = a.intersection(&b).unwrap();
            let expected_intersection =
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::from(IntervalOpen::new(1.0, 2.0)),
                ));
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_left_half_open_right_half_open() {
            let a = IntervalLowerOpenUpperClosed::new(0.0, 2.0);
            let b = IntervalLowerClosedUpperOpen::new(1.0, 3.0);
            let intersection = a.intersection(&b).unwrap();
            let expected_intersection =
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::from(IntervalClosed::new(1.0, 2.0)),
                ));
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_disjoint() {
            let a = IntervalClosed::new(0.0, 1.0);
            let b = IntervalClosed::new(2.0, 3.0);
            let intersection = a.intersection(&b);
            assert!(intersection.is_none());
        }

        /// Test intersection of two disjoint open intervals
        /// This tests the open-open branch where lower_bound >= upper_bound (returns None)
        #[test]
        fn intersection_disjoint_open_open() {
            // (0.0, 1.0) ∩ (2.0, 3.0) = ∅
            let a = IntervalOpen::new(0.0, 1.0);
            let b = IntervalOpen::new(2.0, 3.0);
            let intersection = a.intersection(&b);
            assert!(
                intersection.is_none(),
                "Disjoint open intervals should have empty intersection"
            );
        }

        /// Test intersection of two open intervals that share a boundary point
        /// This also tests the open-open branch where lower_bound >= upper_bound
        #[test]
        fn intersection_open_open_touching_boundary() {
            // (0.0, 1.0) ∩ (1.0, 2.0) = ∅ (they touch at 1.0 but both exclude it)
            let a = IntervalOpen::new(0.0, 1.0);
            let b = IntervalOpen::new(1.0, 2.0);
            let intersection = a.intersection(&b);
            assert!(
                intersection.is_none(),
                "Open intervals touching at boundary should have empty intersection"
            );
        }

        #[test]
        fn intersection_singleton_open() {
            let a = IntervalSingleton::new(1.0);
            let b = IntervalOpen::new(0.0, 2.0);
            let intersection = a.intersection(&b).unwrap();
            let expected_intersection = Interval::FiniteLength(IntervalFiniteLength::from(a));
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_unbounded_closed() {
            let a = IntervalInfiniteLength::LowerUnboundedUpperUnbounded(
                IntervalLowerUnboundedUpperUnbounded::new(),
            );
            let b = IntervalClosed::new(0.0, 2.0);
            let intersection = a.intersection(&b).unwrap();
            let expected_intersection = Interval::FiniteLength(
                IntervalFiniteLength::PositiveLength(IntervalFinitePositiveLength::from(b)),
            );
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_unbounded_unbounded() {
            let a = IntervalInfiniteLength::<f64>::LowerUnboundedUpperUnbounded(
                IntervalLowerUnboundedUpperUnbounded::new(),
            );
            let b = IntervalInfiniteLength::<f64>::LowerUnboundedUpperUnbounded(
                IntervalLowerUnboundedUpperUnbounded::new(),
            );
            let intersection = a.intersection(&b).unwrap();
            let expected_intersection =
                Interval::InfiniteLength(IntervalInfiniteLength::LowerUnboundedUpperUnbounded(
                    IntervalLowerUnboundedUpperUnbounded::new(),
                ));
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_subintervals() {
            let a =
                SubIntervalInPartition::<IntervalClosed<f64>>::First(IntervalClosed::new(0.0, 2.0));
            let b = SubIntervalInPartition::<IntervalClosed<f64>>::Last(
                IntervalLowerOpenUpperClosed::new(1.0, 2.0),
            );
            let intersection = a.intersection(&b).unwrap();
            let expected_intersection =
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::from(IntervalLowerOpenUpperClosed::new(1.0, 2.0)),
                ));
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_none() {
            let i1 = IntervalClosed::new(0.0, 1.0);
            let i2 = IntervalClosed::new(2.0, 3.0);
            assert!(i1.intersection(&i2).is_none());
        }

        #[test]
        fn intersection_singleton() {
            let i1 = IntervalClosed::new(0.0, 1.0);
            let i2 = IntervalClosed::new(1.0, 2.0);
            let intersection = i1.intersection(&i2).unwrap();
            let expected_intersection =
                Interval::FiniteLength(IntervalFiniteLength::from(IntervalSingleton::new(1.0)));
            assert_eq!(intersection, expected_intersection);
        }

        #[test]
        fn intersection_01() {
            let lower_bound = 0.;
            let upper_bound_0 = 1.;
            let upper_bound_1 = 0.5;

            let closed_0 = IntervalClosed::new(lower_bound, upper_bound_0); // [0., 1.]
            let closed_1 = IntervalClosed::new(lower_bound, upper_bound_1); // [0., 0.5]

            let i0: IntervalFinitePositiveLength<_> = closed_0
                .intersection(&closed_1)
                .unwrap()
                .try_into()
                .unwrap(); // [0., 0.5]
            assert_eq!(i0.lower_bound_value(), &0.);
            assert_eq!(i0.upper_bound_value(), &0.5);

            let open_1 = IntervalOpen::new(lower_bound, upper_bound_1); // (0., 0.5)

            let i1: IntervalFinitePositiveLength<_> =
                closed_0.intersection(&open_1).unwrap().try_into().unwrap(); // [0., 0.5)
            assert_eq!(i1.lower_bound_value(), &0.);
            assert_eq!(i1.upper_bound_value(), &0.5);

            let singleton_1 = IntervalSingleton::new(lower_bound); // [0., 0.]

            let i2 = singleton_1.intersection(&closed_1).unwrap();
            match i2 {
                Interval::InfiniteLength(_) => panic!("The interval should be a singleton!"),
                Interval::FiniteLength(interval_finite_length) => match interval_finite_length {
                    IntervalFiniteLength::PositiveLength(_) => {
                        panic!("The interval should be a singleton!")
                    }
                    IntervalFiniteLength::ZeroLength(ref singleton) => {
                        assert_eq!(singleton.value(), &0.);
                    }
                },
            };

            let i3 = singleton_1.intersection(&open_1);
            assert!(i3.is_none());
        }

        #[test]
        fn intersection_02() {
            let lower_bound = 0.;
            let upper_bound = 0.5;

            let open_1 = IntervalOpen::new(lower_bound, upper_bound);
            let singleton_1 = IntervalSingleton::new(lower_bound);

            let i3 = singleton_1.intersection(&open_1);
            assert!(i3.is_none());
        }

        /// Test intersection with unbounded lower and open upper bound
        /// This tests the (None, Some(UpperBoundRuntime::Open)) branch
        #[test]
        fn intersection_unbounded_lower_open_upper() {
            // (-∞, 5.0) ∩ [0.0, 10.0] = [0.0, 5.0)
            let a = IntervalLowerUnboundedUpperOpen::new(5.0);
            let b = IntervalClosed::new(0.0, 10.0);

            let intersection = a.intersection(&b).unwrap();

            match intersection {
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::LowerClosedUpperOpen(interval),
                )) => {
                    assert_eq!(interval.lower_bound_value(), &0.0);
                    assert_eq!(interval.upper_bound_value(), &5.0);
                }
                _ => panic!("Expected [0.0, 5.0) interval"),
            }
        }

        /// Test intersection with unbounded lower and open upper bound (disjoint case)
        #[test]
        fn intersection_unbounded_lower_open_upper_disjoint() {
            // (-∞, 0.0) ∩ [1.0, 10.0] = ∅
            let a = IntervalLowerUnboundedUpperOpen::new(0.0);
            let b = IntervalClosed::new(1.0, 10.0);

            let intersection = a.intersection(&b);
            assert!(intersection.is_none());
        }

        /// Test intersection with unbounded lower and closed upper bound
        /// This tests the (None, Some(UpperBoundRuntime::Closed)) branch
        #[test]
        fn intersection_unbounded_lower_closed_upper() {
            // (-∞, 5.0] ∩ [0.0, 10.0] = [0.0, 5.0]
            let a = IntervalLowerUnboundedUpperClosed::new(5.0);
            let b = IntervalClosed::new(0.0, 10.0);

            let intersection = a.intersection(&b).unwrap();

            match intersection {
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::Closed(interval),
                )) => {
                    assert_eq!(interval.lower_bound_value(), &0.0);
                    assert_eq!(interval.upper_bound_value(), &5.0);
                }
                _ => panic!("Expected [0.0, 5.0] interval"),
            }
        }

        /// Test intersection with open lower and unbounded upper bound
        /// This tests the (Some(LowerBoundRuntime::Open), None) branch
        #[test]
        fn intersection_open_lower_unbounded_upper() {
            // (3.0, +∞) ∩ [0.0, 10.0] = (3.0, 10.0]
            let a = IntervalLowerOpenUpperUnbounded::new(3.0);
            let b = IntervalClosed::new(0.0, 10.0);

            let intersection = a.intersection(&b).unwrap();

            match intersection {
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::LowerOpenUpperClosed(interval),
                )) => {
                    assert_eq!(interval.lower_bound_value(), &3.0);
                    assert_eq!(interval.upper_bound_value(), &10.0);
                }
                _ => panic!("Expected (3.0, 10.0] interval"),
            }
        }

        /// Test intersection with open lower and unbounded upper bound (disjoint case)
        #[test]
        fn intersection_open_lower_unbounded_upper_disjoint() {
            // (10.0, +∞) ∩ [0.0, 5.0] = ∅
            let a = IntervalLowerOpenUpperUnbounded::new(10.0);
            let b = IntervalClosed::new(0.0, 5.0);

            let intersection = a.intersection(&b);
            assert!(intersection.is_none());
        }

        /// Test intersection with closed lower and unbounded upper bound
        /// This tests the (Some(LowerBoundRuntime::Closed), None) branch
        #[test]
        fn intersection_closed_lower_unbounded_upper() {
            // [3.0, +∞) ∩ [0.0, 10.0] = [3.0, 10.0]
            let a = IntervalLowerClosedUpperUnbounded::new(3.0);
            let b = IntervalClosed::new(0.0, 10.0);

            let intersection = a.intersection(&b).unwrap();

            match intersection {
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::Closed(interval),
                )) => {
                    assert_eq!(interval.lower_bound_value(), &3.0);
                    assert_eq!(interval.upper_bound_value(), &10.0);
                }
                _ => panic!("Expected [3.0, 10.0] interval"),
            }
        }

        /// Test intersection between two half-unbounded intervals
        #[test]
        fn intersection_half_unbounded_both() {
            // [3.0, +∞) ∩ (-∞, 10.0] = [3.0, 10.0]
            let a = IntervalLowerClosedUpperUnbounded::new(3.0);
            let b = IntervalLowerUnboundedUpperClosed::new(10.0);

            let intersection = a.intersection(&b).unwrap();

            match intersection {
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::Closed(interval),
                )) => {
                    assert_eq!(interval.lower_bound_value(), &3.0);
                    assert_eq!(interval.upper_bound_value(), &10.0);
                }
                _ => panic!("Expected [3.0, 10.0] interval"),
            }
        }

        /// Test intersection between two half-unbounded intervals (open bounds)
        #[test]
        fn intersection_half_unbounded_both_open() {
            // (3.0, +∞) ∩ (-∞, 10.0) = (3.0, 10.0)
            let a = IntervalLowerOpenUpperUnbounded::new(3.0);
            let b = IntervalLowerUnboundedUpperOpen::new(10.0);

            let intersection = a.intersection(&b).unwrap();

            match intersection {
                Interval::FiniteLength(IntervalFiniteLength::PositiveLength(
                    IntervalFinitePositiveLength::Open(interval),
                )) => {
                    assert_eq!(interval.lower_bound_value(), &3.0);
                    assert_eq!(interval.upper_bound_value(), &10.0);
                }
                _ => panic!("Expected (3.0, 10.0) interval"),
            }
        }

        /// Test intersection between two half-unbounded intervals (disjoint)
        #[test]
        fn intersection_half_unbounded_disjoint() {
            // [10.0, +∞) ∩ (-∞, 5.0] = ∅
            let a = IntervalLowerClosedUpperUnbounded::new(10.0);
            let b = IntervalLowerUnboundedUpperClosed::new(5.0);

            let intersection = a.intersection(&b);
            assert!(intersection.is_none());
        }

        /// Test intersection preserving unbounded lower with closed upper
        #[test]
        fn intersection_preserves_unbounded_lower_closed() {
            // (-∞, 10.0] ∩ (-∞, 5.0] = (-∞, 5.0]
            let a = IntervalLowerUnboundedUpperClosed::new(10.0);
            let b = IntervalLowerUnboundedUpperClosed::new(5.0);

            let intersection = a.intersection(&b).unwrap();

            match intersection {
                Interval::InfiniteLength(IntervalInfiniteLength::LowerUnboundedUpperClosed(
                    interval,
                )) => {
                    assert_eq!(interval.upper_bound_value(), &5.0);
                }
                _ => panic!("Expected (-∞, 5.0] interval"),
            }
        }

        /// Test intersection preserving unbounded lower with open upper
        #[test]
        fn intersection_preserves_unbounded_lower_open() {
            // (-∞, 10.0) ∩ (-∞, 5.0) = (-∞, 5.0)
            let a = IntervalLowerUnboundedUpperOpen::new(10.0);
            let b = IntervalLowerUnboundedUpperOpen::new(5.0);

            let intersection = a.intersection(&b).unwrap();

            match intersection {
                Interval::InfiniteLength(IntervalInfiniteLength::LowerUnboundedUpperOpen(
                    interval,
                )) => {
                    assert_eq!(interval.upper_bound_value(), &5.0);
                }
                _ => panic!("Expected (-∞, 5.0) interval"),
            }
        }

        /// Test intersection preserving closed lower with unbounded upper
        #[test]
        fn intersection_preserves_closed_lower_unbounded() {
            // [5.0, +∞) ∩ [3.0, +∞) = [5.0, +∞)
            let a = IntervalLowerClosedUpperUnbounded::new(5.0);
            let b = IntervalLowerClosedUpperUnbounded::new(3.0);

            let intersection = a.intersection(&b).unwrap();

            match intersection {
                Interval::InfiniteLength(IntervalInfiniteLength::LowerClosedUpperUnbounded(
                    interval,
                )) => {
                    assert_eq!(interval.lower_bound_value(), &5.0);
                }
                _ => panic!("Expected [5.0, +∞) interval"),
            }
        }

        /// Test intersection preserving open lower with unbounded upper
        #[test]
        fn intersection_preserves_open_lower_unbounded() {
            // (5.0, +∞) ∩ (3.0, +∞) = (5.0, +∞)
            let a = IntervalLowerOpenUpperUnbounded::new(5.0);
            let b = IntervalLowerOpenUpperUnbounded::new(3.0);

            let intersection = a.intersection(&b).unwrap();

            match intersection {
                Interval::InfiniteLength(IntervalInfiniteLength::LowerOpenUpperUnbounded(
                    interval,
                )) => {
                    assert_eq!(interval.lower_bound_value(), &5.0);
                }
                _ => panic!("Expected (5.0, +∞) interval"),
            }
        }
    }

    mod union {
        use super::*;
        use crate::intervals::{
            Interval, IntervalClosed, IntervalLowerClosedUpperOpen, IntervalLowerOpenUpperClosed,
            IntervalOpen, IntervalSingleton,
        };
        use try_create::New;

        #[test]
        fn union_disjoint() {
            let i1 = IntervalClosed::new(0.0, 1.0);
            let i2 = IntervalClosed::new(2.0, 3.0);
            let union = i1.union(&i2);
            if let IntervalUnion::TwoDisjoint { left, right } = union {
                assert_eq!(left, i1.into());
                assert_eq!(right, i2.into());
            } else {
                panic!("Expected disjoint union");
            }
        }

        /// Test union of two disjoint intervals where `self` is on the right
        /// This tests that the union method correctly orders intervals regardless of call order
        #[test]
        fn union_disjoint_self_on_right() {
            // i1 is [3.0, 4.0] and i2 is [0.0, 1.0]
            // When calling i1.union(&i2), self (i1) is on the right of other (i2)
            let i1 = IntervalClosed::new(3.0, 4.0); // self - right interval
            let i2 = IntervalClosed::new(0.0, 1.0); // other - left interval
            let union = i1.union(&i2);

            // Convert to Interval for comparison
            let i1_interval: Interval<_> = i1.into();
            let i2_interval: Interval<_> = i2.into();

            // The result should still have the intervals ordered correctly (left, right)
            if let IntervalUnion::TwoDisjoint {
                ref left,
                ref right,
            } = union
            {
                // left should be i2 (the interval with lower values)
                assert_eq!(left, &i2_interval, "Left should be the interval [0.0, 1.0]");
                // right should be i1 (the interval with higher values)
                assert_eq!(
                    right, &i1_interval,
                    "Right should be the interval [3.0, 4.0]"
                );
            } else {
                panic!("Expected disjoint union");
            }

            // Verify commutativity: i1.union(&i2) should equal i2.union(&i1)
            let i1 = IntervalClosed::new(3.0, 4.0);
            let i2 = IntervalClosed::new(0.0, 1.0);
            let union_reversed = i2.union(&i1);
            assert_eq!(union, union_reversed, "Union should be commutative");
        }

        /// Test union of disjoint open intervals where `self` is on the right
        #[test]
        fn union_disjoint_open_self_on_right() {
            // (5.0, 6.0) ∪ (1.0, 2.0) where self is on the right
            let i1 = IntervalOpen::new(5.0, 6.0); // self - right interval
            let i2 = IntervalOpen::new(1.0, 2.0); // other - left interval
            let union = i1.union(&i2);

            assert!(union.is_two_disjoint(), "Should be two disjoint intervals");

            // Test the gap before destructuring
            let gap = union.gap();
            assert!(gap.is_some(), "Should have a gap");
            if let Some(ref gap_interval) = gap {
                assert!(gap_interval.contains_point(&3.0), "Gap should contain 3.0");
                assert!(
                    gap_interval.contains_point(&2.0),
                    "Gap should include lower boundary"
                );
                assert!(
                    gap_interval.contains_point(&5.0),
                    "Gap should include upper boundary"
                );
            }

            // Convert to Interval for comparison
            let i1_interval: Interval<_> = i1.into();
            let i2_interval: Interval<_> = i2.into();

            if let IntervalUnion::TwoDisjoint { left, right } = union {
                // Intervals should be ordered correctly
                assert_eq!(left, i2_interval, "Left should be (1.0, 2.0)");
                assert_eq!(right, i1_interval, "Right should be (5.0, 6.0)");
            } else {
                panic!("Expected TwoDisjoint variant");
            }
        }

        #[test]
        fn union_overlapping() {
            let i1 = IntervalClosed::new(0.0, 2.0);
            let i2 = IntervalClosed::new(1.0, 3.0);
            let union = i1.union(&i2);
            let expected: Interval<f64> = IntervalClosed::new(0.0, 3.0).into();
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, expected);
            } else {
                panic!("Expected connected union");
            }
        }

        #[test]
        fn union_adjacent_closed_open() {
            let i1 = IntervalClosed::new(0.0, 1.0);
            let i2 = IntervalOpen::new(1.0, 2.0);
            let union = i1.union(&i2);
            assert!(union.is_single_connected());

            let expected: Interval<f64> = IntervalLowerClosedUpperOpen::new(0.0, 2.0).into();
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, expected);
            } else {
                panic!("Expected connected union");
            }
        }

        #[test]
        fn union_adjacent_open_closed() {
            let i1 = IntervalOpen::new(0.0, 1.0);
            let i2 = IntervalClosed::new(1.0, 2.0);
            let union = i1.union(&i2);
            assert!(union.is_single_connected());

            let expected: Interval<f64> = IntervalLowerOpenUpperClosed::new(0.0, 2.0).into();
            if let IntervalUnion::SingleConnected(ref connected) = union {
                assert_eq!(connected, &expected);
            } else {
                panic!("Expected connected union");
            }

            let connected = union.as_single_connected().unwrap();
            assert_eq!(connected, &expected);

            assert!(union.as_two_disjoint().is_none());

            assert!(union.contains_point(&1.0));
        }

        #[test]
        fn union_adjacent_open_open() {
            let i1 = IntervalOpen::new(0.0, 1.0);
            let i2 = IntervalOpen::new(1.0, 2.0);
            let union = i1.union(&i2);
            assert!(union.is_two_disjoint());

            let i1: Interval<_> = i1.into();
            let i2: Interval<_> = i2.into();

            if let IntervalUnion::TwoDisjoint {
                ref left,
                ref right,
            } = union
            {
                assert_eq!(left, &i1);
                assert_eq!(right, &i2);
            } else {
                panic!("Expected disjoint union because of the gap at 1.0");
            }

            let (left, right) = union.as_two_disjoint().unwrap();
            assert_eq!(left, &i1);
            assert_eq!(right, &i2);

            assert!(union.as_single_connected().is_none());

            assert!(!union.contains_point(&1.0));
        }

        /// Test union of two disjoint open intervals with a gap between them
        /// This is different from union_adjacent_open_open where intervals touch at boundary
        #[test]
        fn union_disjoint_open_open() {
            // (0.0, 1.0) ∪ (3.0, 4.0) = TwoDisjoint with gap [1.0, 3.0]
            let i1 = IntervalOpen::new(0.0, 1.0);
            let i2 = IntervalOpen::new(3.0, 4.0);
            let union = i1.union(&i2);

            assert!(union.is_two_disjoint(), "Expected two disjoint intervals");

            // Test the gap
            let gap = union.gap();
            assert!(
                gap.is_some(),
                "Should have a gap between disjoint intervals"
            );

            if let Some(ref gap_interval) = gap {
                // Gap should be [1.0, 3.0]
                assert!(
                    gap_interval.contains_point(&2.0),
                    "Gap should contain point 2.0"
                );
                assert!(
                    gap_interval.contains_point(&1.0),
                    "Gap should include lower boundary 1.0"
                );
                assert!(
                    gap_interval.contains_point(&3.0),
                    "Gap should include upper boundary 3.0"
                );
                assert!(
                    !gap_interval.contains_point(&0.5),
                    "Gap should not contain points in left interval"
                );
                assert!(
                    !gap_interval.contains_point(&3.5),
                    "Gap should not contain points in right interval"
                );
            }

            // Test the disjoint intervals
            let i1_interval: Interval<_> = i1.into();
            let i2_interval: Interval<_> = i2.into();

            if let IntervalUnion::TwoDisjoint { left, right } = union {
                assert_eq!(left, i1_interval);
                assert_eq!(right, i2_interval);
            } else {
                panic!("Expected TwoDisjoint variant");
            }

            // Test point containment on new union
            let i1 = IntervalOpen::new(0.0, 1.0);
            let i2 = IntervalOpen::new(3.0, 4.0);
            let union = i1.union(&i2);
            assert!(
                union.contains_point(&0.5),
                "Should contain point in first interval"
            );
            assert!(
                union.contains_point(&3.5),
                "Should contain point in second interval"
            );
            assert!(
                !union.contains_point(&2.0),
                "Should not contain point in gap"
            );
            assert!(
                !union.contains_point(&1.0),
                "Should not contain excluded boundary"
            );
            assert!(
                !union.contains_point(&3.0),
                "Should not contain excluded boundary"
            );
        }

        /// Test union of disjoint open intervals with half-open intervals
        #[test]
        fn union_disjoint_open_half_open() {
            // (0.0, 1.0) ∪ [2.0, 3.0) with clear gap
            let i1 = IntervalOpen::new(0.0, 1.0);
            let i2 = IntervalLowerClosedUpperOpen::new(2.0, 3.0);
            let union = i1.union(&i2);

            assert!(union.is_two_disjoint());

            if let IntervalUnion::TwoDisjoint { left, right } = union {
                assert_eq!(left, i1.into());
                assert_eq!(right, i2.into());
            } else {
                panic!("Expected disjoint union");
            }
        }

        #[test]
        fn union_one_contains_other() {
            let i1 = IntervalClosed::new(0.0, 3.0);
            let i2 = IntervalClosed::new(1.0, 2.0);
            let union = i1.union(&i2);
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, i1.into());
            } else {
                panic!("Expected connected union");
            }
        }

        #[test]
        fn union_with_singleton_inside() {
            let i1 = IntervalClosed::new(0.0, 2.0);
            let i2 = IntervalSingleton::new(1.0);
            let union = i1.union(&i2);
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, i1.into());
            } else {
                panic!("Expected connected union");
            }
        }

        #[test]
        fn union_with_singleton_at_boundary_closed() {
            let i1 = IntervalClosed::new(0.0, 1.0);
            let i2 = IntervalSingleton::new(1.0);
            let union = i1.union(&i2);
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, i1.into());
            } else {
                panic!("Expected connected union");
            }
        }

        #[test]
        fn union_with_singleton_at_boundary_open() {
            let i1 = IntervalLowerClosedUpperOpen::new(0.0, 1.0);
            let i2 = IntervalSingleton::new(1.0);
            let union = i1.union(&i2);
            let expected: Interval<f64> = IntervalClosed::new(0.0, 1.0).into();
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, expected);
            } else {
                panic!("Expected connected union");
            }
        }

        #[test]
        fn union_with_singleton_outside() {
            let i1 = IntervalClosed::new(0.0, 1.0);
            let i2 = IntervalSingleton::new(2.0);
            let union = i1.union(&i2);
            if let IntervalUnion::TwoDisjoint { left, right } = union {
                assert_eq!(left, i1.into());
                assert_eq!(right, i2.into());
            } else {
                panic!("Expected disjoint union");
            }
        }

        #[test]
        fn union_unbounded_overlapping() {
            let i1 = IntervalLowerClosedUpperUnbounded::new(0.0);
            let i2 = IntervalLowerClosedUpperUnbounded::new(1.0);
            let union = i1.union(&i2);
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, i1.into());
            } else {
                panic!("Expected connected union");
            }
        }

        #[test]
        fn union_lower_unbounded_with_upper_unbounded() {
            let i1 = IntervalLowerUnboundedUpperClosed::new(0.0);
            let i2 = IntervalLowerClosedUpperUnbounded::new(1.0);
            let union = i1.union(&i2);
            if let IntervalUnion::TwoDisjoint { left, right } = union {
                assert_eq!(left, i1.into());
                assert_eq!(right, i2.into());
            } else {
                panic!("Expected disjoint union for (-inf, 0] U [1, inf)");
            }
        }

        #[test]
        fn union_lower_unbounded_with_upper_unbounded_overlapping() {
            let i1 = IntervalLowerUnboundedUpperClosed::new(2.0);
            let i2 = IntervalLowerClosedUpperUnbounded::new(1.0);
            let union = i1.union(&i2);
            let expected: Interval<f64> = IntervalLowerUnboundedUpperUnbounded::new().into();
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, expected);
            } else {
                panic!("Expected connected union resulting in (-inf, inf)");
            }
        }

        #[test]
        fn union_finite_with_unbounded() {
            let i1 = IntervalClosed::new(0.0, 5.0);
            let i2 = IntervalLowerClosedUpperUnbounded::new(3.0);
            let union = i1.union(&i2);
            let expected: Interval<f64> = IntervalLowerClosedUpperUnbounded::new(0.0).into();
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, expected);
            } else {
                panic!("Expected connected union");
            }
        }

        #[test]
        fn union_identical_intervals() {
            let i1 = IntervalClosed::new(0.0, 1.0);
            let i2 = IntervalClosed::new(0.0, 1.0);
            let union = i1.union(&i2);
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, i1.into());
            } else {
                panic!("Expected connected union");
            }
        }

        #[test]
        fn union_different_bound_types_contained() {
            let i1 = IntervalClosed::new(0.0, 2.0);
            let i2 = IntervalOpen::new(0.0, 2.0);
            let union = i1.union(&i2);
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, i1.into());
            } else {
                panic!("Expected connected union");
            }
        }

        #[test]
        fn union_different_bound_types_overlapping() {
            let i1 = IntervalLowerClosedUpperOpen::new(0.0, 2.0); // [0, 2)
            let i2 = IntervalLowerOpenUpperClosed::new(1.0, 3.0); // (1, 3]
            let union = i1.union(&i2);
            let expected: Interval<f64> = IntervalClosed::new(0.0, 3.0).into();
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, expected);
            } else {
                panic!("Expected connected union");
            }
        }

        #[test]
        fn union_with_fully_unbounded() {
            let i1 = IntervalClosed::new(0.0, 1.0);
            let i2 = IntervalLowerUnboundedUpperUnbounded::new();
            let union = i1.union(&i2);
            if let IntervalUnion::SingleConnected(connected) = union {
                assert_eq!(connected, i2.into());
            } else {
                panic!("Expected connected union");
            }
        }

        mod gap {
            use super::*;
            use num_valid::{RealNative64StrictFinite, functions::Abs};

            type Real = RealNative64StrictFinite;

            #[test]
            fn gap_with_disjoint_intervals_positive_gap() {
                // Test case: [1, 2] and [4, 5] with gap (2, 4)
                let interval1 =
                    IntervalClosed::new(Real::try_new(1.0).unwrap(), Real::try_new(2.0).unwrap());
                let interval2 =
                    IntervalClosed::new(Real::try_new(4.0).unwrap(), Real::try_new(5.0).unwrap());

                let union = interval1.union(&interval2);

                let gap = union.gap().unwrap();

                // Gap should be the open interval (2, 4)
                match gap {
                    IntervalFinitePositiveLength::Open(open_gap) => {
                        assert_eq!(open_gap.lower_bound_value(), &2.0);
                        assert_eq!(open_gap.upper_bound_value(), &4.0);
                        assert_eq!(open_gap.length().as_ref(), &2.0);
                    }
                    _ => panic!("Expected open interval for gap"),
                }
            }

            #[test]
            fn gap_with_disjoint_intervals_small_gap() {
                // Test case: [0, 1] and [1.5, 2] with gap (1, 1.5)
                let interval1 =
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap());
                let interval2 =
                    IntervalClosed::new(Real::try_new(1.5).unwrap(), Real::try_new(2.0).unwrap());

                let union = interval1.union(&interval2);

                let gap = union.gap().unwrap();

                match gap {
                    IntervalFinitePositiveLength::Open(open_gap) => {
                        assert_eq!(open_gap.lower_bound_value(), &1.0);
                        assert_eq!(open_gap.upper_bound_value(), &1.5);
                        assert_eq!(open_gap.length().as_ref(), &0.5);
                    }
                    _ => panic!("Expected open interval for gap"),
                }
            }

            #[test]
            fn gap_with_disjoint_half_open_intervals() {
                // Test case: [0, 1) and (2, 3] with gap [1, 2]
                let interval1 = IntervalLowerClosedUpperOpen::new(
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                );
                let interval2 = IntervalLowerOpenUpperClosed::new(
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                );

                let union = interval1.union(&interval2);

                let gap = union.gap().unwrap();

                match gap {
                    IntervalFinitePositiveLength::Closed(closed_gap) => {
                        assert_eq!(closed_gap.lower_bound_value(), &1.0);
                        assert_eq!(closed_gap.upper_bound_value(), &2.0);
                        assert_eq!(closed_gap.length().as_ref(), &1.0);
                    }
                    _ => panic!("Expected closed interval for gap"),
                }
            }

            #[test]
            fn gap_with_touching_intervals_no_gap() {
                // Test case: [1, 2] and [2, 3] - touching at point 2, no gap
                let interval1 =
                    IntervalClosed::new(Real::try_new(1.0).unwrap(), Real::try_new(2.0).unwrap());
                let interval2 =
                    IntervalClosed::new(Real::try_new(2.0).unwrap(), Real::try_new(3.0).unwrap());

                let union = interval1.union(&interval2);

                assert!(union.gap().is_none());
            }

            #[test]
            fn gap_with_overlapping_intervals_no_gap() {
                // Test case: [1, 3] and [2, 4] - overlapping, no gap
                let interval1 =
                    IntervalClosed::new(Real::try_new(1.0).unwrap(), Real::try_new(3.0).unwrap());
                let interval2 =
                    IntervalClosed::new(Real::try_new(2.0).unwrap(), Real::try_new(4.0).unwrap());

                let union = interval1.union(&interval2);

                assert!(union.gap().is_none());
            }

            #[test]
            fn gap_with_open_intervals_disjoint() {
                // Test case: (0, 1) and (2, 3) with gap [1, 2]
                let interval1 =
                    IntervalOpen::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap());
                let interval2 =
                    IntervalOpen::new(Real::try_new(2.0).unwrap(), Real::try_new(3.0).unwrap());

                let union = interval1.union(&interval2);

                let gap = union.gap().unwrap();

                match gap {
                    IntervalFinitePositiveLength::Closed(closed_gap) => {
                        assert_eq!(closed_gap.lower_bound_value(), &1.0);
                        assert_eq!(closed_gap.upper_bound_value(), &2.0);
                        assert_eq!(closed_gap.length().as_ref(), &1.0);
                    }
                    _ => panic!("Expected closed interval for gap"),
                }
            }

            #[test]
            fn gap_with_single_connected_interval() {
                // Test case: Single interval [1, 5] - no gap possible
                let interval =
                    IntervalClosed::new(Real::try_new(1.0).unwrap(), Real::try_new(5.0).unwrap());

                let union = interval.union(&interval);

                assert!(union.gap().is_none());
            }

            #[test]
            fn gap_with_very_small_gap() {
                // Test case: [0, 1] and [1 + ε, 2] where ε is very small
                let epsilon = PositiveRealScalar::try_new(Real::try_new(1e-10).unwrap()).unwrap();
                let interval1 =
                    IntervalClosed::new(Real::try_new(0.0).unwrap(), Real::try_new(1.0).unwrap());
                let interval2 = IntervalClosed::new(
                    Real::try_new(1.0).unwrap() + epsilon.as_ref(),
                    Real::try_new(2.0).unwrap(),
                );

                let union = interval1.union(&interval2);

                let gap = union.gap().unwrap();

                match gap {
                    IntervalFinitePositiveLength::Open(open_gap) => {
                        assert_eq!(open_gap.lower_bound_value(), &1.0);
                        assert_eq!(
                            open_gap.upper_bound_value(),
                            &(Real::try_new(1.0).unwrap() + epsilon.as_ref())
                        );
                        more_asserts::assert_lt!(
                            (open_gap.length().into_inner() - epsilon.as_ref()).abs(),
                            1.0e-16
                        );
                    }
                    _ => panic!("Expected open interval for very small gap"),
                }
            }

            #[test]
            fn gap_with_negative_intervals() {
                // Test case: [-5, -3] and [-1, 1] with gap (-3, -1)
                let interval1 =
                    IntervalClosed::new(Real::try_new(-5.0).unwrap(), Real::try_new(-3.0).unwrap());
                let interval2 =
                    IntervalClosed::new(Real::try_new(-1.0).unwrap(), Real::try_new(1.0).unwrap());

                let union = IntervalUnion::TwoDisjoint {
                    left: interval1.into(),
                    right: interval2.into(),
                };

                let gap = union.gap().unwrap();

                match gap {
                    IntervalFinitePositiveLength::Open(open_gap) => {
                        assert_eq!(open_gap.lower_bound_value(), &-3.0);
                        assert_eq!(open_gap.upper_bound_value(), &-1.0);
                        assert_eq!(open_gap.length().as_ref(), &2.0);
                    }
                    _ => panic!("Expected open interval for gap"),
                }
            }

            #[test]
            fn gap_with_mixed_boundary_types() {
                // Test case: (0, 1] and [3, 4) with gap (1, 3)
                let interval1 = IntervalLowerOpenUpperClosed::new(
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                );
                let interval2 = IntervalLowerClosedUpperOpen::new(
                    Real::try_new(3.0).unwrap(),
                    Real::try_new(4.0).unwrap(),
                );

                let union = IntervalUnion::TwoDisjoint {
                    left: interval1.into(),
                    right: interval2.into(),
                };

                let gap = union.gap().unwrap();

                match gap {
                    IntervalFinitePositiveLength::Open(open_gap) => {
                        assert_eq!(open_gap.lower_bound_value(), &1.0);
                        assert_eq!(open_gap.upper_bound_value(), &3.0);
                        assert_eq!(open_gap.length().as_ref(), &2.0);
                    }
                    _ => panic!("Expected open interval for gap"),
                }
            }

            #[test]
            fn gap_preserves_boundary_semantics() {
                // Test that the gap correctly reflects which boundaries should be included
                // based on the original intervals' boundary types

                // Case 1: [0, 1) and (2, 3] should have gap [1, 2]
                let interval1 = IntervalLowerClosedUpperOpen::new(
                    Real::try_new(0.0).unwrap(),
                    Real::try_new(1.0).unwrap(),
                );
                let interval2 = IntervalLowerOpenUpperClosed::new(
                    Real::try_new(2.0).unwrap(),
                    Real::try_new(3.0).unwrap(),
                );

                let union = IntervalUnion::TwoDisjoint {
                    left: interval1.into(),
                    right: interval2.into(),
                };

                let gap = union.gap().unwrap();

                // Gap should include both 1 and 2 since they're excluded from original intervals
                match gap {
                    IntervalFinitePositiveLength::Closed(closed_gap) => {
                        assert_eq!(closed_gap.lower_bound_value(), &1.0);
                        assert_eq!(closed_gap.upper_bound_value(), &2.0);
                        assert!(closed_gap.contains_point(&Real::try_new(1.0).unwrap()));
                        assert!(closed_gap.contains_point(&Real::try_new(2.0).unwrap()));
                    }
                    _ => panic!("Expected closed interval for gap"),
                }
            }

            #[test]
            fn gap_with_singleton_intervals() {
                // Test gap between two singleton intervals
                let singleton1 = IntervalSingleton::new(Real::try_new(1.0).unwrap());
                let singleton2 = IntervalSingleton::new(Real::try_new(3.0).unwrap());

                let union = IntervalUnion::TwoDisjoint {
                    left: singleton1.into(),
                    right: singleton2.into(),
                };

                let gap = union.gap().unwrap();

                match gap {
                    IntervalFinitePositiveLength::Open(open_gap) => {
                        assert_eq!(open_gap.lower_bound_value(), &1.0);
                        assert_eq!(open_gap.upper_bound_value(), &3.0);
                        assert_eq!(open_gap.length().as_ref(), &2.0);
                        // Gap should not contain the singleton endpoints
                        assert!(!open_gap.contains_point(&Real::try_new(1.0).unwrap()));
                        assert!(!open_gap.contains_point(&Real::try_new(3.0).unwrap()));
                    }
                    _ => panic!("Expected open interval for gap between singletons"),
                }
            }
        }
    }
}