//! Ellipse shape implementation.
//!
//! This module provides an ellipse implementation for use in Euler and Venn diagrams.
//! Ellipses are more flexible than circles and can represent elongated or rotated regions,
//! making them useful for more accurate area-proportional diagrams.
//!
//! # Representation
//!
//! An ellipse is defined by:
//! - **Center point**: The center of the ellipse
//! - **Semi-major axis** (a): Half the length of the longest diameter
//! - **Semi-minor axis** (b): Half the length of the shortest diameter  
//! - **Rotation**: Rotation angle in radians (counterclockwise from x-axis)
//!
//! # Area Calculations
//!
//! The module provides several specialized area computation methods:
//!
//! - **Sector area**: The area swept by a radius from the center through an angle θ
//! - **Segment area**: The area between the ellipse boundary and a chord connecting two points
//!
//! These primitives are essential for computing intersection areas in Euler diagrams.
//!
//! # Examples
//!
//! ```
//! use eunoia::geometry::shapes::Ellipse;
//! use eunoia::geometry::traits::Area;
//! use eunoia::geometry::primitives::Point;
//!
//! // Create an ellipse centered at origin
//! let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
//!
//! // Compute total area
//! let area = ellipse.area();
//!
//! // Compute sector area for π/4 radians
//! let sector = ellipse.sector_area(std::f64::consts::PI / 4.0);
//! ```

use std::f64::consts::PI;

use nalgebra::Matrix2;

use crate::geometry::diagram::RegionMask;
use crate::geometry::primitives::Point;
use crate::geometry::projective::Conic;
use crate::geometry::shapes::{Polygon, Rectangle};
use crate::geometry::traits::{
    Area, BoundingBox, Centroid, Closed, DiagramShape, Perimeter, Polygonize,
};

/// An ellipse defined by center, semi-major and semi-minor axes, and rotation.
///
/// Ellipses are oval-shaped closed curves that generalize circles. They are particularly
/// useful in Euler diagrams when sets have elongated or directional relationships.
///
/// # Representation
///
/// The ellipse is represented in standard form with:
/// - A center point `(h, k)`
/// - Semi-major axis length `a` (≥ semi-minor axis)
/// - Semi-minor axis length `b` (> 0)
/// - Rotation angle `φ` in radians (counterclockwise from the positive x-axis)
///
/// The canonical equation in the ellipse's local coordinate system is:
/// ```text
/// (x/a)² + (y/b)² = 1
/// ```
///
/// # Examples
///
/// ```
/// use eunoia::geometry::shapes::Ellipse;
/// use eunoia::geometry::primitives::Point;
/// use eunoia::geometry::traits::Area;
///
/// // Create an ellipse at the origin with no rotation
/// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
/// assert!((ellipse.area() - 47.12).abs() < 0.01);
///
/// // Create a rotated ellipse
/// let rotated = Ellipse::new(
///     Point::new(2.0, 3.0),
///     4.0,
///     2.0,
///     std::f64::consts::PI / 4.0  // 45 degrees
/// );
/// ```
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Ellipse {
    center: Point,
    semi_major: f64,
    semi_minor: f64,
    rotation: f64, // in radians
}

impl Ellipse {
    /// Creates a new ellipse with the specified parameters.
    ///
    /// # Arguments
    ///
    /// * `center` - The center point of the ellipse
    /// * `semi_major` - The semi-major axis length (must be > 0)
    /// * `semi_minor` - The semi-minor axis length (must be > 0)
    /// * `rotation` - Rotation angle in radians (counterclockwise from x-axis)
    ///
    /// # Panics
    ///
    /// Panics in debug builds if either axis length is <= 0.
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    ///
    /// let ellipse = Ellipse::new(Point::new(1.0, 2.0), 4.0, 3.0, 0.0);
    /// assert_eq!(ellipse.semi_major(), 4.0);
    /// assert_eq!(ellipse.semi_minor(), 3.0);
    /// ```
    pub fn new(center: Point, semi_major: f64, semi_minor: f64, rotation: f64) -> Self {
        debug_assert!(semi_major > 0.0, "Semi-major axis must be > 0");
        debug_assert!(semi_minor > 0.0, "Semi-minor axis must be > 0");
        Self {
            center,
            semi_major,
            semi_minor,
            rotation,
        }
    }

    /// Creates a new ellipse from radius and log-aspect ratio parameterization.
    ///
    /// This parameterization is optimized for numerical optimization:
    /// - `radius` controls the overall size (geometric mean of axes)
    /// - `log_aspect` controls elongation on a log scale (more uniform landscape)
    /// - `log_aspect = 0.0` gives a circle (aspect_ratio = 1.0)
    /// - `log_aspect < 0` gives elongated ellipses
    /// - `log_aspect = -0.693` gives aspect_ratio ≈ 0.5 (2:1 elongation)
    ///
    /// # Arguments
    ///
    /// * `center` - The center point of the ellipse
    /// * `radius` - Geometric mean radius: sqrt(semi_major * semi_minor)
    /// * `log_aspect` - Log of aspect ratio (semi_minor/semi_major), clamped to [-1.609, 0]
    ///   which corresponds to aspect_ratio in [0.2, 1.0]
    /// * `rotation` - Rotation angle in radians
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    ///
    /// // Circle: log_aspect = 0.0
    /// let circle = Ellipse::from_radius_ratio(Point::new(0.0, 0.0), 2.0, 0.0, 0.0);
    /// assert!((circle.semi_major() - 2.0).abs() < 1e-10);
    /// assert!((circle.semi_minor() - 2.0).abs() < 1e-10);
    ///
    /// // Elongated ellipse: log_aspect = -0.693 (aspect ≈ 0.5)
    /// let ellipse = Ellipse::from_radius_ratio(Point::new(0.0, 0.0), 2.0, -0.693147, 0.0);
    /// // radius = sqrt(a * b) = 2.0, aspect = 0.5
    /// // => a = 2.828, b = 1.414
    /// assert!((ellipse.semi_major() - 2.828).abs() < 0.01);
    /// assert!((ellipse.semi_minor() - 1.414).abs() < 0.01);
    /// ```
    pub fn from_radius_ratio(center: Point, radius: f64, log_aspect: f64, rotation: f64) -> Self {
        let radius = radius.abs();

        // Clamp log_aspect to [-ln(5), 0] which gives aspect_ratio in [0.2, 1.0]
        let log_aspect = log_aspect.clamp(-1.6094379124341003, 0.0); // -ln(5) to 0
        let aspect_ratio = log_aspect.exp();

        // Normalize rotation to [0, π) since ellipse has 180° rotational symmetry
        let rotation = rotation.rem_euclid(PI);

        // Given: r = sqrt(a * b) and aspect = b/a
        // Solving: a = r / sqrt(aspect), b = r * sqrt(aspect)
        let semi_major = radius / aspect_ratio.sqrt();
        let semi_minor = radius * aspect_ratio.sqrt();

        Self::new(center, semi_major, semi_minor, rotation)
    }

    pub fn from_conic(conic: Conic) -> Option<Self> {
        let c = conic.matrix();

        // Extract algebraic coefficients from the symmetric matrix
        // Q = [ A   B/2  D/2 ]
        //     [ B/2 C    E/2 ]
        //     [ D/2 E/2  F   ]
        let m1 = c[(0, 0)];
        let m2 = 2.0 * c[(0, 1)];
        let m3 = c[(1, 1)];
        let m4 = 2.0 * c[(0, 2)];
        let m5 = 2.0 * c[(1, 2)];
        let m6 = c[(2, 2)];

        // Check ellipse condition (negative discriminant required)
        if m2 * m2 - 4.0 * m1 * m3 >= 0.0 {
            return None;
        }

        // Quadratic form matrix
        let m = Matrix2::new(m1, m2 / 2.0, m2 / 2.0, m3);

        // Solve for center:
        // M * [xc, yc]^T = -1/2 * [D, E]^T
        let rhs = nalgebra::Vector2::new(-m4 / 2.0, -m5 / 2.0);
        let center_vec = m.lu().solve(&rhs)?;
        let xc = center_vec[0];
        let yc = center_vec[1];

        // Compute translated constant term F̄
        let f_bar = m6 + m1 * xc * xc + m2 * xc * yc + m3 * yc * yc + m4 * xc + m5 * yc;

        if f_bar >= 0.0 {
            return None; // Not an ellipse
        }

        // Eigen-decompose the quadratic part
        let eig = m.symmetric_eigen();
        let lambda1 = eig.eigenvalues[0];
        let lambda2 = eig.eigenvalues[1];
        let v = eig.eigenvectors;

        // Identify major/minor axes:
        // smaller eigenvalue ⇒ bigger radius
        let (lambda_major, lambda_minor, vec_major) = if lambda1 < lambda2 {
            (lambda1, lambda2, v.column(0))
        } else {
            (lambda2, lambda1, v.column(1))
        };

        // Radii from: a² = -F̄ / λ
        let a2 = -f_bar / lambda_major;
        let b2 = -f_bar / lambda_minor;

        if a2 <= 0.0 || b2 <= 0.0 {
            return None;
        }

        let a = a2.sqrt();
        let b = b2.sqrt();

        // Rotation angle from the major-axis eigenvector
        let vx = vec_major[0];
        let vy = vec_major[1];
        let mut phi = vy.atan2(vx);

        // Normalize angle into [0, π)
        if phi < 0.0 {
            phi += std::f64::consts::PI;
        }

        Some(Ellipse {
            center: Point::new(xc, yc),
            semi_major: a,
            semi_minor: b,
            rotation: phi,
        })
    }

    /// Returns the center point of the ellipse.
    pub fn center(&self) -> Point {
        self.center
    }

    /// Returns the semi-major axis length.
    pub fn semi_major(&self) -> f64 {
        self.semi_major
    }

    /// Returns the semi-minor axis length.
    pub fn semi_minor(&self) -> f64 {
        self.semi_minor
    }

    /// Returns the rotation angle in radians.
    pub fn rotation(&self) -> f64 {
        self.rotation
    }

    /// Computes the area of an elliptical sector from the center through angle θ.
    ///
    /// An elliptical sector is the region bounded by:
    /// - Two radii from the center to the ellipse boundary
    /// - The ellipse arc between those radii
    ///
    /// This uses the exact formula for elliptical sectors, which accounts for
    /// the non-uniform curvature of the ellipse (unlike circular sectors where
    /// area is simply ½r²θ).
    ///
    /// # Arguments
    ///
    /// * `theta` - The angle in radians from the semi-major axis
    ///
    /// # Returns
    ///
    /// The area of the sector from 0 to θ
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    /// use eunoia::geometry::traits::Area;
    ///
    /// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.0, 0.0);
    ///
    /// // Quarter sector (π/2 radians)
    /// let quarter = ellipse.sector_area(std::f64::consts::PI / 2.0);
    /// assert!((quarter - ellipse.area() / 4.0).abs() < 1e-10);
    /// ```
    pub fn sector_area(&self, theta: f64) -> f64 {
        let a = self.semi_major;
        let b = self.semi_minor;

        let num = (b - a) * (2.0 * theta).sin();
        let den = a + b + (b - a) * (2.0 * theta).cos();

        0.5 * a * b * (theta - num.atan2(den))
    }

    /// Computes the area of an elliptical sector between two angles.
    ///
    /// This is equivalent to `sector_area(theta1) - sector_area(theta0)` but
    /// handles counter-clockwise angle wrapping automatically.
    ///
    /// # Arguments
    ///
    /// * `theta0` - Starting angle in radians
    /// * `theta1` - Ending angle in radians
    ///
    /// # Returns
    ///
    /// The area of the sector from θ₀ to θ₁ (counter-clockwise)
    ///
    /// # Note
    ///
    /// If θ₁ < θ₀, the function adds 2π to θ₁ to compute the counter-clockwise
    /// sector that wraps through angle 0.
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    ///
    /// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);
    ///
    /// // Sector from 45° to 135°
    /// let sector = ellipse.sector_area_between(
    ///     std::f64::consts::PI / 4.0,
    ///     3.0 * std::f64::consts::PI / 4.0
    /// );
    /// assert!(sector > 0.0);
    /// ```
    pub fn sector_area_between(&self, theta0: f64, theta1: f64) -> f64 {
        let t0 = theta0;
        let mut t1 = theta1;

        // Ensure CCW ordering.
        if t1 < t0 {
            t1 += 2.0 * PI;
        }

        self.sector_area(t1) - self.sector_area(t0)
    }

    /// Computes the area of an ellipse segment between two boundary points.
    ///
    /// An ellipse segment is the region bounded by:
    /// - The ellipse arc between two points on the boundary
    /// - The chord (straight line) connecting those two points
    ///
    /// This method automatically handles:
    /// - Coordinate system transformation to the ellipse's local frame
    /// - Counter-clockwise angle ordering
    /// - Minor arc (≤ 180°) vs major arc (> 180°) selection
    ///
    /// # Algorithm
    ///
    /// 1. Transform points to ellipse coordinate system (centered, unrotated)
    /// 2. Compute angles θ₀ and θ₁ for each point
    /// 3. Ensure counter-clockwise ordering (add 2π if needed)
    /// 4. Calculate segment as: sector - triangle (for minor arc)
    ///    or: total_area - complementary_sector + triangle (for major arc)
    ///
    /// # Arguments
    ///
    /// * `p0` - First boundary point (in world coordinates)
    /// * `p1` - Second boundary point (in world coordinates)
    ///
    /// # Returns
    ///
    /// The area of the segment. When the angular span > π, returns the **major arc**
    /// segment (the larger of the two possible segments).
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    /// use eunoia::geometry::traits::Area;
    ///
    /// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
    ///
    /// // Segment between points on opposite sides (semicircle)
    /// let p0 = Point::new(5.0, 0.0);   // Right end of major axis
    /// let p1 = Point::new(-5.0, 0.0);  // Left end of major axis
    /// let segment = ellipse.ellipse_segment(p0, p1);
    ///
    /// // Should be approximately half the ellipse area
    /// assert!((segment - ellipse.area() / 2.0).abs() < 1e-10);
    /// ```
    pub fn ellipse_segment(&self, p0: Point, p1: Point) -> f64 {
        // 1. Move into ellipse coordinate system.
        let p0 = p0.to_ellipse_frame(self);
        let p1 = p1.to_ellipse_frame(self);

        // 2. Compute angles
        let theta0 = p0.angle_from_origin();
        let mut theta1 = p1.angle_from_origin();

        // 3. Ensure CCW ordering.
        if theta1 < theta0 {
            theta1 += 2.0 * PI;
        }

        // 4. Triangle correction (signed area of parallelogram / 2).
        let triangle = 0.5 * (p1.x() * p0.y() - p0.x() * p1.y()).abs();

        // 5. Minor or major arc?
        if (theta1 - theta0) <= PI {
            self.sector_area(theta1) - self.sector_area(theta0) - triangle
        } else {
            self.area() - (self.sector_area(theta0 + 2.0 * PI) - self.sector_area(theta1))
                + triangle
        }
    }

    /// Computes the lens-shaped intersection area between two ellipses with exactly two intersection points.
    ///
    /// Follows eulerr's polysegments: processes both arcs around the intersection perimeter.
    fn compute_lens_area(&self, other: &Self, p0: &Point, p1: &Point) -> f64 {
        // Two arcs: p1→p0 and p0→p1
        // For each arc, compute triangle + min(segment for common parents)

        // Arc 1: p1 → p0
        let seg1_10 = self.ellipse_segment(*p1, *p0);
        let seg2_10 = other.ellipse_segment(*p1, *p0);
        let triangle_10 = 0.5 * ((p0.x() + p1.x()) * (p0.y() - p1.y()));
        let arc1 = triangle_10 + seg1_10.min(seg2_10);

        // Arc 2: p0 → p1
        let seg1_01 = self.ellipse_segment(*p0, *p1);
        let seg2_01 = other.ellipse_segment(*p0, *p1);
        let triangle_01 = 0.5 * ((p1.x() + p0.x()) * (p1.y() - p0.y()));
        let arc2 = triangle_01 + seg1_01.min(seg2_01);

        arc1 + arc2
    }

    /// Computes intersection area for cases with 3 or 4 intersection points.
    ///
    /// This follows eulerr's polysegments algorithm.
    fn compute_multi_point_intersection_area(&self, other: &Self, points: &[Point]) -> f64 {
        if points.len() < 3 {
            return 0.0;
        }

        // Use the centroid of intersection points as a reference center
        let cx = points.iter().map(|p| p.x()).sum::<f64>() / points.len() as f64;
        let cy = points.iter().map(|p| p.y()).sum::<f64>() / points.len() as f64;

        // Sort points by angle around the centroid (match eulerr's atan2(x-cx, y-cy) order)
        let mut sorted_points: Vec<Point> = points.to_vec();
        sorted_points.sort_by(|a, b| {
            let angle_a = (a.x() - cx).atan2(a.y() - cy);
            let angle_b = (b.x() - cx).atan2(b.y() - cy);
            angle_a.partial_cmp(&angle_b).unwrap()
        });

        // Compute area using polysegments algorithm
        let mut total_area = 0.0;
        let n = sorted_points.len();

        for k in 0..n {
            let l = if k == 0 { n - 1 } else { k - 1 };

            let pi = sorted_points[k]; // current point (matches eulerr's i)
            let pj = sorted_points[l]; // previous point (matches eulerr's j)

            // Triangle contribution (shoelace formula) - matches eulerr line 90
            let triangle = 0.5 * ((pj.x() + pi.x()) * (pj.y() - pi.y()));

            // Segments from i to j (matches eulerr line 86: ellipse_segment(..., points[i], points[j]))
            let seg1 = self.ellipse_segment(pi, pj);
            let seg2 = other.ellipse_segment(pi, pj);

            // Add triangle + min(segment) as in eulerr line 89-91
            total_area += triangle + seg1.min(seg2);
        }

        total_area
    }
}

impl Area for Ellipse {
    /// Computes the total area of the ellipse using the formula A = πab.
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    /// use eunoia::geometry::traits::Area;
    ///
    /// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
    /// let area = ellipse.area();
    /// assert!((area - std::f64::consts::PI * 15.0).abs() < 1e-10);
    /// ```
    fn area(&self) -> f64 {
        PI * self.semi_major * self.semi_minor
    }
}

impl Perimeter for Ellipse {
    /// Computes an approximation of the ellipse perimeter.
    ///
    /// Uses Ramanujan's second approximation formula, which provides excellent
    /// accuracy for all ellipses:
    ///
    /// ```text
    /// P ≈ π(a + b)(1 + 3h / (10 + √(4 - 3h)))
    /// where h = ((a - b) / (a + b))²
    /// ```
    ///
    /// This approximation has a relative error of less than 0.01% for typical cases.
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    /// use eunoia::geometry::traits::Perimeter;
    ///
    /// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
    /// let perimeter = ellipse.perimeter();
    /// assert!(perimeter > 0.0);
    /// ```
    fn perimeter(&self) -> f64 {
        // Approximation using Ramanujan's second formula
        let a = self.semi_major;
        let b = self.semi_minor;
        let h = ((a - b).powi(2)) / ((a + b).powi(2));
        PI * (a + b) * (1.0 + (3.0 * h) / (10.0 + (4.0 - 3.0 * h).sqrt()))
    }
}

impl Centroid for Ellipse {
    /// Returns the centroid of the ellipse, which is its center point.
    fn centroid(&self) -> Point {
        self.center
    }
}

impl BoundingBox for Ellipse {
    /// Computes the axis-aligned bounding box that contains the ellipse.
    ///
    /// For a rotated ellipse, this computes the smallest axis-aligned rectangle
    /// that fully contains the ellipse. The calculation accounts for the rotation
    /// by projecting the ellipse axes onto the coordinate axes.
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Ellipse;
    /// use eunoia::geometry::primitives::Point;
    /// use eunoia::geometry::traits::BoundingBox;
    ///
    /// let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
    /// let bbox = ellipse.bounding_box();
    /// // For unrotated ellipse, bbox dimensions equal 2*semi-major and 2*semi-minor
    /// ```
    fn bounding_box(&self) -> Rectangle {
        let cos_theta = self.rotation.cos();
        let sin_theta = self.rotation.sin();

        let width = 2.0
            * ((self.semi_major * cos_theta).powi(2) + (self.semi_minor * sin_theta).powi(2))
                .sqrt();
        let height = 2.0
            * ((self.semi_major * sin_theta).powi(2) + (self.semi_minor * cos_theta).powi(2))
                .sqrt();

        Rectangle::new(self.center, width, height)
    }
}

impl Closed for Ellipse {
    fn contains(&self, other: &Self) -> bool {
        // Quick check: if other's center is outside self, it can't be contained
        if !self.contains_point(&other.center) {
            return false;
        }

        // If other is larger than self (by area), it can't be contained
        if other.area() > self.area() {
            return false;
        }

        // Convert both ellipses to conics
        let c1 = Conic::from_ellipse(*self);
        let c2 = Conic::from_ellipse(*other);

        // Check for boundary intersections
        let intersection_points = c1.intersect_conic(&c2);

        // If no intersections found and center is inside, other is contained
        if intersection_points.is_empty() {
            return true;
        }

        // If intersections were found, verify they're not real by checking
        // the extreme points of other (ends of major/minor axes)
        let phi = other.rotation;
        let cos_phi = phi.cos();
        let sin_phi = phi.sin();

        // Check points at the ends of the semi-major axis
        let major_offset_x = other.semi_major * cos_phi;
        let major_offset_y = other.semi_major * sin_phi;
        let p1 = Point::new(
            other.center.x() + major_offset_x,
            other.center.y() + major_offset_y,
        );
        let p2 = Point::new(
            other.center.x() - major_offset_x,
            other.center.y() - major_offset_y,
        );

        // Check points at the ends of the semi-minor axis
        let minor_offset_x = other.semi_minor * -sin_phi;
        let minor_offset_y = other.semi_minor * cos_phi;
        let p3 = Point::new(
            other.center.x() + minor_offset_x,
            other.center.y() + minor_offset_y,
        );
        let p4 = Point::new(
            other.center.x() - minor_offset_x,
            other.center.y() - minor_offset_y,
        );

        // All extreme points must be inside or on self
        self.contains_point(&p1)
            && self.contains_point(&p2)
            && self.contains_point(&p3)
            && self.contains_point(&p4)
    }

    fn contains_point(&self, point: &Point) -> bool {
        // Transform point to ellipse's local coordinate system
        let dx = point.x() - self.center.x();
        let dy = point.y() - self.center.y();

        let cos_phi = self.rotation.cos();
        let sin_phi = self.rotation.sin();

        // Rotate point to align with ellipse axes
        let x_local = dx * cos_phi + dy * sin_phi;
        // Correct inverse rotation for y': y_local = -dx*sin(phi) + dy*cos(phi)
        let y_local = -dx * sin_phi + dy * cos_phi;

        // Check if point is inside using the ellipse equation
        (x_local * x_local) / (self.semi_major * self.semi_major)
            + (y_local * y_local) / (self.semi_minor * self.semi_minor)
            <= 1.0
    }

    fn intersects(&self, other: &Self) -> bool {
        // Quick rejection: if centers are too far apart, they can't intersect
        let center_distance = self.center.distance(&other.center);
        let max_reach_self = self.semi_major;
        let max_reach_other = other.semi_major;

        if center_distance > max_reach_self + max_reach_other {
            return false;
        }

        // Check if one contains the other (no intersection in that case)
        if self.contains(other) || other.contains(self) {
            return false;
        }

        // Convert to conics and check for intersection points
        let c1 = Conic::from_ellipse(*self);
        let c2 = Conic::from_ellipse(*other);

        let intersection_points = c1.intersect_conic(&c2);

        !intersection_points.is_empty()
    }

    fn intersection_area(&self, other: &Self) -> f64 {
        // Check if one ellipse contains the other
        if self.contains(other) {
            return other.area();
        }
        if other.contains(self) {
            return self.area();
        }

        // Get intersection points
        let points = self.intersection_points(other);
        let n = points.len();

        match n {
            0 => 0.0, // No intersection
            1 => {
                // Single point of tangency - negligible area
                0.0
            }
            2 => {
                // Two intersection points - compute the lens area
                self.compute_lens_area(other, &points[0], &points[1])
            }
            _ => {
                // Three or four intersection points
                self.compute_multi_point_intersection_area(other, &points)
            }
        }
    }

    fn intersection_points(&self, other: &Self) -> Vec<Point> {
        // Convert both ellipses to conics
        let c1 = Conic::from_ellipse(*self);
        let c2 = Conic::from_ellipse(*other);

        // Get homogeneous intersection points
        let homogeneous_points = c1.intersect_conic(&c2);

        // Convert to Cartesian points
        homogeneous_points
            .into_iter()
            .map(Point::from_homogeneous)
            .collect()
    }
}

impl DiagramShape for Ellipse {
    fn compute_exclusive_regions(shapes: &[Self]) -> std::collections::HashMap<RegionMask, f64> {
        use crate::geometry::diagram::{adopters_to_mask, to_exclusive_areas, IntersectionPoint};
        use std::collections::HashMap;

        let n_sets = shapes.len();

        // 1) Collect ALL intersection points from ALL pairs (matches eulerr lines 52-63)
        let mut intersections: Vec<IntersectionPoint> = Vec::new();
        for i in 0..n_sets {
            for j in (i + 1)..n_sets {
                let pts = shapes[i].intersection_points(&shapes[j]);
                for p in pts {
                    // adopters: all shapes that contain this intersection point
                    // Use a small tolerance when checking containment to include boundary points
                    let adopters: Vec<usize> = (0..n_sets)
                        .filter(|&k| {
                            let local = p.to_ellipse_frame(&shapes[k]);
                            let val = (local.x() / shapes[k].semi_major).powi(2)
                                + (local.y() / shapes[k].semi_minor).powi(2);
                            val <= 1.0 + 1e-9
                        })
                        .collect();
                    intersections.push(IntersectionPoint::new(p, (i, j), adopters));
                }
            }
        }

        // 2) Generate all region masks (matches eulerr's id = set_index(n))
        let n_overlaps = (1 << n_sets) - 1;
        let mut overlapping_areas: HashMap<RegionMask, f64> = HashMap::new();

        for mask in 1..=n_overlaps {
            let bits: Vec<usize> = (0..n_sets).filter(|&i| (mask & (1 << i)) != 0).collect();

            if bits.len() == 1 {
                // Single set - just return its area
                overlapping_areas.insert(mask, shapes[bits[0]].area());
            } else {
                // Two or more sets - find intersection points for this region (lines 75-79)
                let region_points: Vec<usize> = intersections
                    .iter()
                    .enumerate()
                    .filter_map(|(idx, ip)| {
                        let (p1, p2) = ip.parents();
                        // Check if both parents are in the mask
                        let parents_in_mask = bits.contains(&p1) && bits.contains(&p2);
                        // Check if mask is a subset of adopters
                        let adopters_mask = adopters_to_mask(ip.adopters());
                        let mask_subset = (mask & adopters_mask) == mask;

                        if parents_in_mask && mask_subset {
                            Some(idx)
                        } else {
                            None
                        }
                    })
                    .collect();

                if region_points.is_empty() {
                    // No intersections: either disjoint or subset (line 83)
                    let area = disjoint_or_subset(shapes, &bits);
                    overlapping_areas.insert(mask, area);
                } else {
                    // Use polysegments algorithm (line 88)
                    let area = compute_polysegments(shapes, &intersections, &region_points);
                    overlapping_areas.insert(mask, area);
                }
            }
        }

        // 5) Convert overlapping to exclusive (lines 104-109)
        to_exclusive_areas(&overlapping_areas)
    }

    fn params_from_circle(x: f64, y: f64, radius: f64) -> Vec<f64> {
        // Convert circle to ellipse using direct semi-axis representation:
        // semi_major = semi_minor = radius (circle), rotation = 0
        vec![x, y, radius, radius, 0.0]
    }

    fn n_params() -> usize {
        5 // x, y, semi_major, semi_minor, rotation
    }

    fn from_params(params: &[f64]) -> Self {
        assert_eq!(
            params.len(),
            5,
            "Ellipse requires 5 parameters: x, y, semi_major, semi_minor, rotation"
        );

        Ellipse::new(
            Point::new(params[0], params[1]),
            params[2].abs(), // Ensure positive semi_major
            params[3].abs(), // Ensure positive semi_minor
            params[4],
        )
    }

    fn to_params(&self) -> Vec<f64> {
        vec![
            self.center.x(),
            self.center.y(),
            self.semi_major,
            self.semi_minor,
            self.rotation,
        ]
    }
}

/// Check if region is disjoint or one ellipse contains another (eulerr's disjoint_or_subset)
fn disjoint_or_subset(shapes: &[Ellipse], indices: &[usize]) -> f64 {
    // Check containment
    for &i in indices {
        for &j in indices {
            if i != j && shapes[i].contains(&shapes[j]) {
                return shapes[j].area();
            }
        }
    }
    // Disjoint
    0.0
}

/// Compute area using polysegments algorithm (eulerr's polysegments function)
fn compute_polysegments(
    shapes: &[Ellipse],
    intersections: &[crate::geometry::diagram::IntersectionPoint],
    region_points: &[usize],
) -> f64 {
    let n = region_points.len();
    if n == 0 {
        return 0.0;
    }

    // Compute centroid of the intersection points (lines 44-50)
    let cx = region_points
        .iter()
        .map(|&idx| intersections[idx].point().x())
        .sum::<f64>()
        / n as f64;
    let cy = region_points
        .iter()
        .map(|&idx| intersections[idx].point().y())
        .sum::<f64>()
        / n as f64;

    // Sort points by angle around centroid (lines 52-62)
    let mut sorted_indices: Vec<usize> = (0..n).collect();
    sorted_indices.sort_by(|&a, &b| {
        let pa = intersections[region_points[a]].point();
        let pb = intersections[region_points[b]].point();
        let angle_a = (pa.x() - cx).atan2(pa.y() - cy);
        let angle_b = (pb.x() - cx).atan2(pb.y() - cy);
        angle_a.partial_cmp(&angle_b).unwrap()
    });

    // Compute area by walking around the boundary (lines 65-100)
    let mut area = 0.0;

    for k in 0..n {
        let l = if k == 0 { n - 1 } else { k - 1 };

        let i_idx = region_points[sorted_indices[k]];
        let j_idx = region_points[sorted_indices[l]];

        let pi = intersections[i_idx].point();
        let pj = intersections[j_idx].point();

        // Find common parents (lines 72-78)
        let parents_i = intersections[i_idx].parents();
        let parents_j = intersections[j_idx].parents();

        let mut common_parents = Vec::new();
        for &p in &[parents_i.0, parents_i.1] {
            if p == parents_j.0 || p == parents_j.1 {
                common_parents.push(p);
            }
        }

        if common_parents.is_empty() {
            // Failure case - should not happen with valid input
            return 0.0;
        }

        // Triangle contribution (line 90)
        let triangle = 0.5 * ((pj.x() + pi.x()) * (pj.y() - pi.y()));

        // Compute segments from common parent ellipses (lines 85-86)
        let mut segments = Vec::new();
        for &ellipse_idx in &common_parents {
            let seg = shapes[ellipse_idx].ellipse_segment(*pi, *pj);
            segments.push(seg);
        }

        // Add triangle + min(segments) (line 89-91)
        if let Some(&min_seg) = segments.iter().min_by(|a, b| a.partial_cmp(b).unwrap()) {
            area += triangle + min_seg;
        }
    }

    area
}

impl Polygonize for Ellipse {
    fn polygonize(&self, n_vertices: usize) -> Polygon {
        use std::f64::consts::PI;

        let n_vertices = n_vertices.max(3);
        let mut vertices = Vec::with_capacity(n_vertices);

        let cos_rot = self.rotation.cos();
        let sin_rot = self.rotation.sin();

        for i in 0..n_vertices {
            let angle = 2.0 * PI * (i as f64) / (n_vertices as f64);

            // Point in local ellipse frame
            let local_x = self.semi_major * angle.cos();
            let local_y = self.semi_minor * angle.sin();

            // Rotate and translate to world frame
            let x = self.center.x() + local_x * cos_rot - local_y * sin_rot;
            let y = self.center.y() + local_x * sin_rot + local_y * cos_rot;

            vertices.push(Point::new(x, y));
        }

        Polygon::new(vertices)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::geometry::shapes::Circle;

    fn approx_eq(a: f64, b: f64) -> bool {
        (a - b).abs() < 1e-10
    }

    #[test]
    fn test_sector_area_between_quarter() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);

        // Quarter sector from 0 to π/2
        let sector_between = ellipse.sector_area_between(0.0, PI / 2.0);
        let expected = ellipse.sector_area(PI / 2.0) - ellipse.sector_area(0.0);

        assert!(approx_eq(sector_between, expected));
        assert!(approx_eq(sector_between, ellipse.area() / 4.0));
    }

    #[test]
    fn test_sector_area_between_with_ccw_adjustment() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.0, 0.0);

        // When theta1 < theta0, it should add 2π to theta1 for CCW ordering
        // Going from 3π/4 to π/4 CCW (wrapping around through 0)
        let theta0 = 3.0 * PI / 4.0;
        let theta1 = PI / 4.0;

        let sector_between = ellipse.sector_area_between(theta0, theta1);

        // This should be equivalent to going from 3π/4 to (π/4 + 2π)
        let expected = ellipse.sector_area(theta1 + 2.0 * PI) - ellipse.sector_area(theta0);

        assert!(
            approx_eq(sector_between, expected),
            "Expected {}, got {}",
            expected,
            sector_between
        );

        // The sector should be positive and cover most of the ellipse
        assert!(sector_between > 0.0);
        assert!(sector_between > ellipse.area() / 2.0);
    }

    #[test]
    fn test_sector_area_between_same_angle() {
        let ellipse = Ellipse::new(Point::new(1.0, 2.0), 3.0, 2.5, 0.3);

        // Same angle should give zero sector area
        let angle = PI / 3.0;
        let sector = ellipse.sector_area_between(angle, angle);

        assert!(approx_eq(sector, 0.0));
    }

    #[test]
    fn test_sector_area_between_full_rotation() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);

        // Full rotation: from 0 to 2π
        let sector = ellipse.sector_area_between(0.0, 2.0 * PI);

        assert!(approx_eq(sector, ellipse.area()));
    }

    #[test]
    fn test_sector_full_ellipse() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);
        let full_angle = 2.0 * PI;
        let sector = ellipse.sector_area(full_angle);
        assert!(approx_eq(sector, ellipse.area()));
    }

    #[test]
    fn test_sector_half_ellipse() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 4.0, 3.0, 0.0);
        let half_angle = PI;
        let sector = ellipse.sector_area(half_angle);
        assert!(approx_eq(sector, ellipse.area() / 2.0));
    }

    #[test]
    fn test_sector_quarter_ellipse() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let quarter_angle = PI / 2.0;
        let sector = ellipse.sector_area(quarter_angle);
        assert!(approx_eq(sector, ellipse.area() / 4.0));
    }

    #[test]
    fn test_sector_zero_angle() {
        let ellipse = Ellipse::new(Point::new(1.0, 1.0), 3.0, 2.0, 0.5);
        let sector = ellipse.sector_area(0.0);
        assert!(approx_eq(sector, 0.0));
    }

    #[test]
    fn test_sector_circular_ellipse_matches_circle() {
        // Create a circle-like ellipse (semi_major == semi_minor)
        let radius = 3.0;
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), radius, radius, 0.0);
        let circle = Circle::new(Point::new(0.0, 0.0), radius);

        let angle = PI / 3.0;
        let ellipse_sector = ellipse.sector_area(angle);
        let circle_sector = circle.sector_area(angle);

        assert!(approx_eq(ellipse_sector, circle_sector));
    }

    #[test]
    fn test_sector_circular_ellipse_multiple_angles() {
        let radius = 2.0;
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), radius, radius, 0.0);
        let circle = Circle::new(Point::new(0.0, 0.0), radius);

        let test_angles = vec![
            0.0,
            PI / 6.0,
            PI / 4.0,
            PI / 2.0,
            PI,
            3.0 * PI / 2.0,
            2.0 * PI,
        ];

        for angle in test_angles {
            let ellipse_sector = ellipse.sector_area(angle);
            let circle_sector = circle.sector_area(angle);
            assert!(
                approx_eq(ellipse_sector, circle_sector),
                "Mismatch at angle {}: ellipse={}, circle={}",
                angle,
                ellipse_sector,
                circle_sector
            );
        }
    }

    #[test]
    fn test_ellipse_segment_semicircle() {
        let radius = 2.0;
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), radius, radius, 0.0);

        // Points at opposite ends of diameter (along x-axis)
        let p0 = Point::new(-radius, 0.0);
        let p1 = Point::new(radius, 0.0);

        let segment = ellipse.ellipse_segment(p0, p1);
        // Semicircle segment equals half the circle area
        assert!(approx_eq(segment, ellipse.area() / 2.0));
    }

    #[test]
    fn test_ellipse_segment_circular_matches_circle_segment() {
        let radius = 3.0;
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), radius, radius, 0.0);
        let circle = Circle::new(Point::new(0.0, 0.0), radius);

        // Test angle
        let angle = PI / 4.0;

        // Points on circle boundary
        let p0 = Point::new(radius, 0.0);
        let p1 = Point::new(radius * angle.cos(), radius * angle.sin());

        let ellipse_segment = ellipse.ellipse_segment(p0, p1);
        let circle_segment = circle.segment_area_from_angle(angle);

        assert!(
            approx_eq(ellipse_segment, circle_segment),
            "Circular ellipse segment should match circle segment: ellipse={}, circle={}",
            ellipse_segment,
            circle_segment
        );
    }

    #[test]
    fn test_ellipse_segment_small_angle() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);

        let angle = PI / 12.0; // 15 degrees
        let p0 = Point::new(ellipse.semi_major(), 0.0);
        let p1 = Point::new(
            ellipse.semi_major() * angle.cos(),
            ellipse.semi_minor() * angle.sin(),
        );

        let segment = ellipse.ellipse_segment(p0, p1);

        // Segment should be positive and less than the sector
        assert!(segment > 0.0);
        let sector = ellipse.sector_area(angle);
        assert!(segment < sector);
    }

    #[test]
    fn test_ellipse_segment_zero_angle() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.5, 0.0);

        // Same point twice should give zero area
        let p = Point::new(ellipse.semi_major(), 0.0);
        let segment = ellipse.ellipse_segment(p, p);

        assert!(approx_eq(segment, 0.0));
    }

    #[test]
    fn test_ellipse_segment_with_known_values() {
        // Reference values computed for ellipse with semi_major=5.0, semi_minor=3.0
        // Total area: 47.12388980384689

        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);

        // Quarter ellipse segment (from θ=0 to θ=π/2)
        let p0 = Point::new(ellipse.semi_major(), 0.0);
        let p1 = Point::new(0.0, ellipse.semi_minor());
        let segment_quarter = ellipse.ellipse_segment(p0, p1);
        let expected_quarter = 4.280972450961723;

        assert!(
            approx_eq(segment_quarter, expected_quarter),
            "Quarter segment: expected {}, got {}",
            expected_quarter,
            segment_quarter
        );

        // Half ellipse segment (from θ=0 to θ=π)
        let p2 = Point::new(-ellipse.semi_major(), 0.0);
        let segment_half = ellipse.ellipse_segment(p0, p2);
        let expected_half = 23.561944901923447;

        assert!(
            approx_eq(segment_half, expected_half),
            "Half segment: expected {}, got {}",
            expected_half,
            segment_half
        );

        // Half segment should also equal half the ellipse area
        assert!(approx_eq(segment_half, ellipse.area() / 2.0));
    }

    #[test]
    fn test_ellipse_segment_additional_angles() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);

        // 30 degree segment (π/6)
        let angle = PI / 6.0;
        let p0 = Point::new(ellipse.semi_major(), 0.0);
        let p1 = Point::new(
            ellipse.semi_major() * angle.cos(),
            ellipse.semi_minor() * angle.sin(),
        );
        let segment_30deg = ellipse.ellipse_segment(p0, p1);
        let expected_30deg = 0.17699081698724095;

        assert!(
            approx_eq(segment_30deg, expected_30deg),
            "30° segment: expected {}, got {}",
            expected_30deg,
            segment_30deg
        );

        // 60 degree segment (π/3)
        let angle = PI / 3.0;
        let p2 = Point::new(
            ellipse.semi_major() * angle.cos(),
            ellipse.semi_minor() * angle.sin(),
        );
        let segment_60deg = ellipse.ellipse_segment(p0, p2);
        let expected_60deg = 1.3587911055911919;

        assert!(
            approx_eq(segment_60deg, expected_60deg),
            "60° segment: expected {}, got {}",
            expected_60deg,
            segment_60deg
        );

        // Additional sanity checks
        assert!(segment_30deg < segment_60deg);
        assert!(segment_60deg < ellipse.area() / 4.0);
    }

    #[test]
    fn test_ellipse_segment_270_degrees() {
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);

        // 270 degree segment (3π/2)
        // When theta1 - theta0 > π, the algorithm computes the MAJOR arc segment.
        // This is the correct behavior: going CCW from 0° to 270° covers 270° of arc,
        // which is the larger of the two possible segments between these points.
        let angle = 3.0 * PI / 2.0;
        let p0 = Point::new(ellipse.semi_major(), 0.0);
        let p1 = Point::new(
            ellipse.semi_major() * angle.cos(),
            ellipse.semi_minor() * angle.sin(),
        );
        let segment_270deg = ellipse.ellipse_segment(p0, p1);
        let expected_270deg = 42.84291735288517;

        assert!(
            approx_eq(segment_270deg, expected_270deg),
            "270° segment: expected {}, got {}",
            expected_270deg,
            segment_270deg
        );

        // Verify this is indeed the major arc (> half the ellipse area)
        assert!(segment_270deg > ellipse.area() / 2.0);

        // The complementary segment (minor arc) would be the small 90° piece
        let complementary = ellipse.area() - segment_270deg;
        assert!(complementary > 0.0);
        assert!(complementary < ellipse.area() / 4.0);
    }

    #[test]
    fn test_contains_point_center_and_boundary() {
        let ellipse = Ellipse::new(Point::new(2.0, 3.0), 5.0, 3.0, 0.0);

        // Center should be inside
        assert!(ellipse.contains_point(&Point::new(2.0, 3.0)));

        // Points on the semi-major axis (horizontal)
        assert!(ellipse.contains_point(&Point::new(7.0, 3.0))); // right edge
        assert!(ellipse.contains_point(&Point::new(-3.0, 3.0))); // left edge

        // Points on the semi-minor axis (vertical)
        assert!(ellipse.contains_point(&Point::new(2.0, 6.0))); // top edge
        assert!(ellipse.contains_point(&Point::new(2.0, 0.0))); // bottom edge

        // Point clearly outside
        assert!(!ellipse.contains_point(&Point::new(10.0, 10.0)));
    }

    #[test]
    fn test_contains_point_rotated_ellipse() {
        // Create an ellipse rotated 45 degrees
        let ellipse = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.0, PI / 4.0);

        // Center should be inside
        assert!(ellipse.contains_point(&Point::new(0.0, 0.0)));

        // Point on the rotated semi-major axis
        let dist = 4.0 / (2.0_f64).sqrt(); // 4 * cos(45°) = 4 * sin(45°)
        assert!(ellipse.contains_point(&Point::new(dist, dist)));
        assert!(ellipse.contains_point(&Point::new(-dist, -dist)));

        // Point clearly outside
        assert!(!ellipse.contains_point(&Point::new(5.0, 0.0)));
        assert!(!ellipse.contains_point(&Point::new(0.0, 5.0)));

        // Point inside but not on axis
        assert!(ellipse.contains_point(&Point::new(0.5, 0.5)));
    }

    #[test]
    fn test_contains_ellipse_concentric() {
        // Larger ellipse should contain smaller concentric one
        let outer = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let inner = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);

        assert!(outer.contains(&inner));
        assert!(!inner.contains(&outer));
    }

    #[test]
    fn test_contains_ellipse_offset_centers() {
        // Large ellipse at origin
        let outer = Ellipse::new(Point::new(0.0, 0.0), 10.0, 8.0, 0.0);

        // Small ellipse offset but still inside
        let inner = Ellipse::new(Point::new(2.0, 1.0), 2.0, 1.5, 0.0);

        assert!(outer.contains(&inner));

        // Small ellipse offset too far (outside)
        let outside = Ellipse::new(Point::new(9.0, 0.0), 2.0, 1.5, 0.0);

        assert!(!outer.contains(&outside));
    }

    #[test]
    fn test_contains_ellipse_rotated() {
        // Test with rotated ellipses
        let outer = Ellipse::new(Point::new(0.0, 0.0), 6.0, 4.0, PI / 6.0);
        let inner = Ellipse::new(Point::new(0.0, 0.0), 2.0, 1.5, PI / 4.0);

        assert!(outer.contains(&inner));
        assert!(!inner.contains(&outer));
    }

    #[test]
    fn test_contains_ellipse_intersecting() {
        // Overlapping but not containing
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 0.0), 4.0, 2.5, 0.0);

        assert!(!e1.contains(&e2));
        assert!(!e2.contains(&e1));
    }

    #[test]
    fn test_intersects_overlapping() {
        // Two ellipses that intersect
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 0.0), 4.0, 2.5, 0.0);

        assert!(e1.intersects(&e2));
        assert!(e2.intersects(&e1));
    }

    #[test]
    fn test_intersects_separate() {
        // Two ellipses that don't touch
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 2.0, 1.5, 0.0);
        let e2 = Ellipse::new(Point::new(10.0, 0.0), 3.0, 2.0, 0.0);

        assert!(!e1.intersects(&e2));
        assert!(!e2.intersects(&e1));
    }

    #[test]
    fn test_intersects_contained() {
        // One ellipse contains the other - no intersection
        let outer = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let inner = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);

        assert!(!outer.intersects(&inner));
        assert!(!inner.intersects(&outer));
    }

    #[test]
    fn test_intersects_touching() {
        // Two ellipses that barely touch (externally tangent)
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);
        let e2 = Ellipse::new(Point::new(6.0, 0.0), 3.0, 2.0, 0.0);

        // They should intersect (tangent point counts as intersection)
        assert!(e1.intersects(&e2));
    }

    #[test]
    fn test_intersection_points_overlapping() {
        // Two ellipses that intersect at multiple points
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 0.0), 4.0, 2.5, 0.0);

        let points = e1.intersection_points(&e2);

        // Should have 2 or 4 intersection points (depending on overlap)
        assert!(!points.is_empty());
        assert!(points.len() <= 4);

        // Note: Due to numerical precision, points on the boundary may not satisfy
        // contains_point with the <= 1.0 check, so we just verify they exist
    }

    #[test]
    fn test_intersection_points_separate() {
        // Two ellipses that don't intersect
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 2.0, 1.5, 0.0);
        let e2 = Ellipse::new(Point::new(10.0, 0.0), 3.0, 2.0, 0.0);

        let points = e1.intersection_points(&e2);

        // Note: The conic intersection algorithm may return spurious points
        // for widely separated ellipses due to numerical issues
        // In practice, these should be filtered or we should check intersects() first
        assert!(points.len() <= 4);
    }

    #[test]
    fn test_intersection_points_contained() {
        // One ellipse contains the other - no boundary intersection
        let outer = Ellipse::new(Point::new(0.0, 0.0), 5.0, 3.0, 0.0);
        let inner = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);

        let points = outer.intersection_points(&inner);

        // Should have no intersection points (or handle numerical artifacts)
        // The contains check would return true, so we expect 0 or duplicate points
        assert!(points.len() <= 4);
    }

    #[test]
    fn test_intersection_points_tangent() {
        // Two circles (special case of ellipses) that are externally tangent
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 2.0, 2.0, 0.0); // circle
        let e2 = Ellipse::new(Point::new(4.0, 0.0), 2.0, 2.0, 0.0); // circle

        let points = e1.intersection_points(&e2);

        // Should have 1 or 2 intersection points (tangent at one point, possibly duplicated)
        assert!(!points.is_empty());

        // The tangent point should be at (2.0, 0.0)
        if !points.is_empty() {
            let first_point = &points[0];
            assert!((first_point.x() - 2.0).abs() < 1e-6);
            assert!(first_point.y().abs() < 1e-6);
        }
    }

    #[test]
    fn test_intersection_area_no_overlap() {
        // Two ellipses that don't overlap
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 2.0, 1.5, 0.0);
        let e2 = Ellipse::new(Point::new(10.0, 0.0), 2.0, 1.5, 0.0);

        let area = e1.intersection_area(&e2);
        assert!(approx_eq(area, 0.0), "Expected 0, got {}", area);
    }

    #[test]
    fn test_intersection_area_one_contains_other() {
        // Larger ellipse contains smaller one
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 5.0, 4.0, 0.0);
        let e2 = Ellipse::new(Point::new(0.0, 0.0), 2.0, 1.5, 0.0);

        let area = e1.intersection_area(&e2);
        let expected = e2.area();

        assert!(
            (area - expected).abs() < 1e-8,
            "Expected {}, got {}",
            expected,
            area
        );

        // Symmetric test
        let area_reverse = e2.intersection_area(&e1);
        assert!(
            (area_reverse - expected).abs() < 1e-8,
            "Expected {}, got {}",
            expected,
            area_reverse
        );
    }

    #[test]
    fn test_intersection_area_identical_ellipses() {
        let e1 = Ellipse::new(Point::new(1.0, 2.0), 3.0, 2.0, PI / 4.0);
        let e2 = Ellipse::new(Point::new(1.0, 2.0), 3.0, 2.0, PI / 4.0);

        let area = e1.intersection_area(&e2);
        let expected = e1.area();

        assert!(
            (area - expected).abs() < 1e-8,
            "Expected {}, got {}",
            expected,
            area
        );
    }

    #[test]
    fn test_intersection_area_circles_partial_overlap() {
        // Two identical circles with partial overlap
        let radius = 3.0;
        let e1 = Ellipse::new(Point::new(0.0, 0.0), radius, radius, 0.0);
        let e2 = Ellipse::new(Point::new(4.0, 0.0), radius, radius, 0.0);

        let area = e1.intersection_area(&e2);

        // For circles, we can use the analytical formula to verify
        let c1 = Circle::new(Point::new(0.0, 0.0), radius);
        let c2 = Circle::new(Point::new(4.0, 0.0), radius);
        let expected = c1.intersection_area(&c2);

        assert!(
            (area - expected).abs() < 1e-6,
            "Expected {}, got {}",
            expected,
            area
        );
    }

    #[test]
    fn test_intersection_area_tangent_ellipses() {
        // Two circles that are externally tangent (touch at one point)
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 2.0, 2.0, 0.0);
        let e2 = Ellipse::new(Point::new(4.0, 0.0), 2.0, 2.0, 0.0);

        let area = e1.intersection_area(&e2);

        // Tangent circles should have zero intersection area
        assert!(area < 1e-6, "Expected ~0, got {}", area);
    }

    #[test]
    fn test_intersection_area_symmetric() {
        // Test that intersection_area is symmetric
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 4.0, 3.0, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 2.0), 3.5, 2.5, PI / 6.0);

        let area1 = e1.intersection_area(&e2);
        let area2 = e2.intersection_area(&e1);

        assert!(
            (area1 - area2).abs() < 1e-8,
            "Areas should be symmetric: {} vs {}",
            area1,
            area2
        );
    }

    #[test]
    fn test_intersection_area_bounds() {
        // Intersection area should be bounded by the smaller ellipse area
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 5.0, 4.0, 0.0);
        let e2 = Ellipse::new(Point::new(2.0, 1.0), 3.0, 2.5, PI / 8.0);

        let area = e1.intersection_area(&e2);
        let min_area = e1.area().min(e2.area());

        assert!(
            area >= 0.0,
            "Intersection area should be non-negative: {}",
            area
        );
        assert!(
            area <= min_area + 1e-6,
            "Intersection area {} should not exceed smaller ellipse area {}",
            area,
            min_area
        );
    }

    #[test]
    fn test_intersection_area_rotated_ellipses() {
        // Two rotated ellipses with partial overlap
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.5, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 0.0), 4.0, 2.5, PI / 2.0);

        let area = e1.intersection_area(&e2);

        // Should have some intersection
        assert!(area > 0.0, "Expected positive area, got {}", area);
        assert!(
            area < e1.area() && area < e2.area(),
            "Area {} should be less than both ellipse areas",
            area
        );
    }

    #[test]
    fn test_intersection_area_nearly_identical() {
        // Two nearly identical ellipses (slightly offset)
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 3.0, 2.0, 0.0);
        let e2 = Ellipse::new(Point::new(0.01, 0.01), 3.0, 2.0, 0.0);

        let area = e1.intersection_area(&e2);
        let expected = e1.area();

        // Should be very close to full area
        assert!(
            (area - expected).abs() < 0.1,
            "Expected ~{}, got {}",
            expected,
            area
        );
    }

    // Helper function for Monte Carlo validation
    fn monte_carlo_intersection_area(
        e1: &Ellipse,
        e2: &Ellipse,
        n_samples: usize,
        seed: u64,
    ) -> f64 {
        use rand::rngs::StdRng;
        use rand::{Rng, SeedableRng};

        let mut rng = StdRng::seed_from_u64(seed);
        let mut in_both = 0;

        // Bounding box for sampling - intersection of both ellipse bounding boxes
        let bbox1 = e1.bounding_box();
        let bbox2 = e2.bounding_box();
        let (min1, max1) = bbox1.to_points();
        let (min2, max2) = bbox2.to_points();

        let x_min = min1.x().max(min2.x());
        let x_max = max1.x().min(max2.x());
        let y_min = min1.y().max(min2.y());
        let y_max = max1.y().min(max2.y());

        if x_min >= x_max || y_min >= y_max {
            return 0.0; // No bounding box overlap
        }

        let bbox_area = (x_max - x_min) * (y_max - y_min);

        for _ in 0..n_samples {
            let x = x_min + (x_max - x_min) * rng.random::<f64>();
            let y = y_min + (y_max - y_min) * rng.random::<f64>();
            let p = Point::new(x, y);

            if e1.contains_point(&p) && e2.contains_point(&p) {
                in_both += 1;
            }
        }

        (in_both as f64 / n_samples as f64) * bbox_area
    }

    #[test]
    fn test_intersection_area_monte_carlo_circles() {
        // Validate circle intersection against Monte Carlo sampling
        let radius = 3.0;
        let e1 = Ellipse::new(Point::new(0.0, 0.0), radius, radius, 0.0);
        let e2 = Ellipse::new(Point::new(4.0, 0.0), radius, radius, 0.0);

        let exact = e1.intersection_area(&e2);
        let mc_estimate = monte_carlo_intersection_area(&e1, &e2, 100_000, 42);

        let error = (exact - mc_estimate).abs() / exact;
        assert!(
            error < 0.05,
            "Monte Carlo and exact should agree within 5%: exact={}, mc={}, error={:.1}%",
            exact,
            mc_estimate,
            error * 100.0
        );
    }

    #[test]
    fn test_intersection_area_monte_carlo_rotated_ellipses() {
        // Validate rotated ellipse intersection against Monte Carlo
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.5, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 0.0), 4.0, 2.5, PI / 2.0);

        let exact = e1.intersection_area(&e2);
        let mc_estimate = monte_carlo_intersection_area(&e1, &e2, 100_000, 42);

        let error = (exact - mc_estimate).abs() / exact.max(mc_estimate);
        assert!(
            error < 0.05,
            "Monte Carlo and exact should agree within 5%: exact={}, mc={}, error={:.1}%",
            exact,
            mc_estimate,
            error * 100.0
        );
    }

    #[test]
    fn test_intersection_area_monte_carlo_general_ellipses() {
        // Validate general ellipse intersection
        let e1 = Ellipse::new(Point::new(1.0, 0.5), 3.5, 2.0, PI / 6.0);
        let e2 = Ellipse::new(Point::new(3.0, 1.5), 3.0, 2.5, PI / 4.0);

        let exact = e1.intersection_area(&e2);
        let mc_estimate = monte_carlo_intersection_area(&e1, &e2, 150_000, 42);

        // For smaller intersections, allow slightly larger relative error
        if exact > 0.5 && mc_estimate > 0.5 {
            let error = (exact - mc_estimate).abs() / exact.max(mc_estimate);
            assert!(
                error < 0.06,
                "Monte Carlo and exact should agree within 6%: exact={}, mc={}, error={:.1}%",
                exact,
                mc_estimate,
                error * 100.0
            );
        } else {
            // Both should be small
            assert!(
                exact < 1.0 && mc_estimate < 1.0,
                "Both should be small: exact={}, mc={}",
                exact,
                mc_estimate
            );
        }
    }

    #[test]
    fn test_intersection_area_monte_carlo_contained() {
        // Validate containment case
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 5.0, 4.0, 0.0);
        let e2 = Ellipse::new(Point::new(0.5, 0.5), 2.0, 1.5, PI / 8.0);

        let exact = e1.intersection_area(&e2);
        let expected = e2.area(); // Should be the smaller ellipse

        // Verify containment
        assert!(e1.contains(&e2), "e1 should contain e2");
        assert!(
            (exact - expected).abs() < 1e-6,
            "Intersection should equal smaller ellipse area: exact={}, expected={}",
            exact,
            expected
        );

        // Monte Carlo validation
        let mc_estimate = monte_carlo_intersection_area(&e1, &e2, 100_000, 42);

        let error = (exact - mc_estimate).abs() / exact;
        assert!(
            error < 0.03,
            "Monte Carlo should agree with exact (contained case): exact={}, mc={}, error={:.1}%",
            exact,
            mc_estimate,
            error * 100.0
        );
    }

    // ------------------------
    // Radius+ratio parameterization tests
    // ------------------------

    #[test]
    fn test_from_radius_ratio_circle() {
        // aspect_ratio = 1.0 should give a circle
        let e = Ellipse::from_radius_ratio(Point::new(1.0, 2.0), 3.0, 1.0, 0.5);
        assert!((e.semi_major() - 3.0).abs() < 1e-10);
        assert!((e.semi_minor() - 3.0).abs() < 1e-10);
        assert_eq!(e.center(), Point::new(1.0, 2.0));
        assert_eq!(e.rotation(), 0.5);
    }

    #[test]
    fn test_from_radius_ratio_elongated() {
        // log_aspect = ln(0.5) ≈ -0.693 means aspect_ratio = 0.5
        let radius = 2.0;
        let log_aspect = -std::f64::consts::LN_2; // ln(0.5)
        let e = Ellipse::from_radius_ratio(Point::new(0.0, 0.0), radius, log_aspect, 0.0);

        // radius = sqrt(a * b), aspect = exp(log_aspect) = 0.5
        // => a = radius / sqrt(aspect), b = radius * sqrt(aspect)
        let aspect = log_aspect.exp();
        let expected_a = radius / aspect.sqrt();
        let expected_b = radius * aspect.sqrt();

        assert!((e.semi_major() - expected_a).abs() < 1e-10);
        assert!((e.semi_minor() - expected_b).abs() < 1e-10);

        // Verify geometric mean
        let geom_mean = (e.semi_major() * e.semi_minor()).sqrt();
        assert!((geom_mean - radius).abs() < 1e-10);
    }

    #[test]
    fn test_from_radius_ratio_preserves_area() {
        // For fixed radius, changing log_aspect changes shape but preserves area
        let radius = 3.0;
        let expected_area = PI * radius * radius;

        for log_aspect in [-1.6, -0.693, -0.223, 0.0] {
            let e = Ellipse::from_radius_ratio(Point::new(0.0, 0.0), radius, log_aspect, 0.0);
            assert!((e.area() - expected_area).abs() < 1e-9);
        }
    }

    #[test]
    fn test_from_radius_ratio_clamping() {
        // log_aspect should be clamped to [-ln(5), 0] which is [-1.609, 0]
        let e_low = Ellipse::from_radius_ratio(Point::new(0.0, 0.0), 2.0, -5.0, 0.0);
        let e_high = Ellipse::from_radius_ratio(Point::new(0.0, 0.0), 2.0, 2.0, 0.0);

        // Both should be valid ellipses (clamped to valid range)
        assert!(e_low.semi_major() > 0.0);
        assert!(e_low.semi_minor() > 0.0);
        assert!(e_high.semi_major() > 0.0);
        assert!(e_high.semi_minor() > 0.0);

        // e_low should be at minimum aspect (most elongated)
        // log_aspect = -1.609 => aspect = 0.2
        assert!((e_low.semi_minor() / e_low.semi_major() - 0.2).abs() < 0.01);

        // e_high should be clamped to 0 (circle)
        assert!((e_high.semi_minor() / e_high.semi_major() - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_diagram_shape_params_from_circle() {
        use crate::geometry::traits::DiagramShape;

        // params_from_circle should give direct semi-axis parameterization
        let params = Ellipse::params_from_circle(1.0, 2.0, 3.0);
        assert_eq!(params.len(), 5);
        assert_eq!(params[0], 1.0); // x
        assert_eq!(params[1], 2.0); // y
        assert_eq!(params[2], 3.0); // semi_major (same as radius for circle)
        assert_eq!(params[3], 3.0); // semi_minor (same as radius for circle)
        assert_eq!(params[4], 0.0); // rotation

        // from_params should reconstruct a circle
        let e = Ellipse::from_params(&params);
        assert!((e.semi_major() - 3.0).abs() < 1e-10);
        assert!((e.semi_minor() - 3.0).abs() < 1e-10);
    }

    // ------------------------
    // Exclusive region tests
    // ------------------------

    #[test]
    fn test_exclusive_regions_two_ellipses_partial_overlap() {
        use crate::geometry::traits::DiagramShape;

        let e1 = Ellipse::new(Point::new(0.0, 0.0), 4.0, 2.5, 0.0);
        let e2 = Ellipse::new(Point::new(3.0, 1.0), 3.5, 2.0, PI / 6.0);

        let areas = Ellipse::compute_exclusive_regions(&[e1, e2]);
        let a_only = areas.get(&(1usize << 0)).copied().unwrap_or(0.0);
        let b_only = areas.get(&(1usize << 1)).copied().unwrap_or(0.0);
        let both = areas
            .get(&((1usize << 0) | (1usize << 1)))
            .copied()
            .unwrap_or(0.0);

        // Basic sanity: non-negative and sums to union
        assert!(a_only >= 0.0 && b_only >= 0.0 && both >= 0.0);
        let union = a_only + b_only + both;
        // Union must be <= sum of individual ellipse areas
        assert!(union <= e1.area() + e2.area() + 1e-6);

        // Intersection area should be close to direct computation
        let exact = e1.intersection_area(&e2);
        let error = (both - exact).abs() / exact.max(1e-6);
        assert!(
            error < 0.1,
            "Error too large: both={}, exact={}, error={:.2}%",
            both,
            exact,
            error * 100.0
        );
    }

    #[test]
    fn test_exclusive_regions_two_ellipses_contained() {
        use crate::geometry::traits::DiagramShape;

        let outer = Ellipse::new(Point::new(0.0, 0.0), 6.0, 4.0, 0.2);
        let inner = Ellipse::new(Point::new(0.5, 0.3), 2.0, 1.5, -0.3);

        let areas = Ellipse::compute_exclusive_regions(&[outer, inner]);
        let inner_only = areas.get(&(1usize << 1)).copied().unwrap_or(0.0);
        let both = areas
            .get(&((1usize << 0) | (1usize << 1)))
            .copied()
            .unwrap_or(0.0);

        // When contained, exclusive intersection should equal inner area
        assert!((both - inner.area()).abs() < 1e-3);
        // Inner exclusive should be ~0
        assert!(inner_only < 1e-6);
    }

    #[test]
    fn test_exclusive_regions_three_ellipses_basic() {
        use crate::geometry::traits::DiagramShape;

        let e1 = Ellipse::new(Point::new(-2.0, 0.0), 3.0, 2.0, 0.2);
        let e2 = Ellipse::new(Point::new(2.0, 0.5), 3.2, 2.1, -0.1);
        let e3 = Ellipse::new(Point::new(0.0, 1.5), 2.8, 1.8, 0.5);

        let areas = Ellipse::compute_exclusive_regions(&[e1, e2, e3]);
        let mask_all = (1usize << 0) | (1usize << 1) | (1usize << 2);
        let all = areas.get(&mask_all).copied().unwrap_or(0.0);

        // Non-negative and bounded by smallest ellipse area
        let min_area = e1.area().min(e2.area()).min(e3.area());
        assert!(all >= 0.0 && all <= min_area + 1e-6);

        // Union should not exceed sum of individuals
        let union: f64 = areas.values().sum();
        assert!(union <= e1.area() + e2.area() + e3.area() + 1e-3);
    }

    // Helper function for Monte Carlo validation of pairwise intersection
    fn monte_carlo_pairwise_intersection(
        e1: &Ellipse,
        e2: &Ellipse,
        n_samples: usize,
        seed: u64,
    ) -> f64 {
        use rand::rngs::StdRng;
        use rand::{Rng, SeedableRng};

        let mut rng = StdRng::seed_from_u64(seed);
        let mut in_both = 0;

        // Bounding box for sampling - intersection of both ellipse bounding boxes
        let bbox1 = e1.bounding_box();
        let bbox2 = e2.bounding_box();
        let (min1, max1) = bbox1.to_points();
        let (min2, max2) = bbox2.to_points();

        let x_min = min1.x().max(min2.x());
        let x_max = max1.x().min(max2.x());
        let y_min = min1.y().max(min2.y());
        let y_max = max1.y().min(max2.y());

        if x_min >= x_max || y_min >= y_max {
            return 0.0; // No bounding box overlap
        }

        let bbox_area = (x_max - x_min) * (y_max - y_min);

        for _ in 0..n_samples {
            let x = x_min + (x_max - x_min) * rng.random::<f64>();
            let y = y_min + (y_max - y_min) * rng.random::<f64>();
            let p = Point::new(x, y);

            if e1.contains_point(&p) && e2.contains_point(&p) {
                in_both += 1;
            }
        }

        (in_both as f64 / n_samples as f64) * bbox_area
    }

    #[test]
    #[ignore = "slow stochastic validation"]
    fn test_random_ellipse_intersections_monte_carlo() {
        use rand::rngs::StdRng;
        use rand::{Rng, SeedableRng};

        const N_TESTS: usize = 100;
        const MC_SAMPLES: usize = 50_000;
        const MAX_RELATIVE_ERROR: f64 = 0.10; // 10% tolerance

        let mut failures = Vec::new();

        for seed in 0..N_TESTS {
            let mut rng = StdRng::seed_from_u64(seed as u64);

            // Random ellipse 1
            let x1 = rng.random_range(-5.0..5.0);
            let y1 = rng.random_range(-5.0..5.0);
            let a1 = rng.random_range(0.5..3.0);
            let b1 = rng.random_range(0.5..3.0);
            let rot1 = rng.random_range(0.0..PI);

            // Random ellipse 2 (overlapping region)
            let x2 = x1 + rng.random_range(-2.0..2.0);
            let y2 = y1 + rng.random_range(-2.0..2.0);
            let a2 = rng.random_range(0.5..3.0);
            let b2 = rng.random_range(0.5..3.0);
            let rot2 = rng.random_range(0.0..PI);

            let e1 = Ellipse::new(Point::new(x1, y1), a1, b1, rot1);
            let e2 = Ellipse::new(Point::new(x2, y2), a2, b2, rot2);

            // Compute exact intersection area
            let exact = e1.intersection_area(&e2);

            // Compute Monte Carlo estimate
            let mc_estimate = monte_carlo_pairwise_intersection(&e1, &e2, MC_SAMPLES, seed as u64);

            // Skip cases with very small intersection (harder to validate)
            if exact < 0.01 && mc_estimate < 0.01 {
                continue;
            }

            // Compute relative error
            let max_val = exact.max(mc_estimate);
            let error = (exact - mc_estimate).abs() / max_val;

            if error > MAX_RELATIVE_ERROR {
                let n_points = e1.intersection_points(&e2).len();
                failures.push((seed, exact, mc_estimate, error, n_points));
            }
        }

        assert!(
            failures.is_empty(),
            "Monte Carlo validation failed for {} of {} cases: {:?}",
            failures.len(),
            N_TESTS,
            failures.iter().take(10).collect::<Vec<_>>()
        );
    }

    #[test]
    fn test_three_ellipse_complete_overlap() {
        use crate::geometry::traits::DiagramShape;

        // Test case: A=B=C=1, A&B&C=1 (all three completely overlap)
        // Three identical circles
        let e1 = Ellipse::new(Point::new(0.0, 0.0), 0.564, 0.564, 0.0); // area ≈ 1.0
        let e2 = Ellipse::new(Point::new(0.0, 0.0), 0.564, 0.564, 0.0);
        let e3 = Ellipse::new(Point::new(0.0, 0.0), 0.564, 0.564, 0.0);

        let areas = Ellipse::compute_exclusive_regions(&[e1, e2, e3]);

        let mask_all = (1usize << 0) | (1usize << 1) | (1usize << 2);
        let all_three = areas.get(&mask_all).copied().unwrap_or(0.0);

        // When all three are identical and overlapping, the 3-way intersection should be ~1.0
        assert!(
            (all_three - 1.0).abs() < 0.01,
            "Complete overlap should give area ~1.0, got {}",
            all_three
        );
    }

    #[test]
    fn test_three_ellipse_partial_overlap_monte_carlo() {
        use rand::rngs::StdRng;
        use rand::{Rng, SeedableRng};

        // Test a specific 3-ellipse configuration with Monte Carlo
        let seed = 42u64;
        let mut rng = StdRng::seed_from_u64(seed);

        let x1 = 0.0;
        let y1 = 0.0;
        let a1 = 1.5;
        let b1 = 1.0;
        let rot1 = 0.0;

        let x2 = 1.0;
        let y2 = 0.0;
        let a2 = 1.5;
        let b2 = 1.0;
        let rot2 = PI / 6.0;

        let x3 = 0.5;
        let y3 = 1.0;
        let a3 = 1.5;
        let b3 = 1.0;
        let rot3 = -PI / 6.0;

        let e1 = Ellipse::new(Point::new(x1, y1), a1, b1, rot1);
        let e2 = Ellipse::new(Point::new(x2, y2), a2, b2, rot2);
        let e3 = Ellipse::new(Point::new(x3, y3), a3, b3, rot3);

        // Monte Carlo estimate for 3-way intersection
        let n_samples = 100_000;
        let mut in_all_three = 0;

        // Find bounding box
        let bbox1 = e1.bounding_box();
        let bbox2 = e2.bounding_box();
        let bbox3 = e3.bounding_box();

        let (min1, max1) = bbox1.to_points();
        let (min2, max2) = bbox2.to_points();
        let (min3, max3) = bbox3.to_points();

        let x_min = min1.x().max(min2.x()).max(min3.x());
        let x_max = max1.x().min(max2.x()).min(max3.x());
        let y_min = min1.y().max(min2.y()).max(min3.y());
        let y_max = max1.y().min(max2.y()).min(max3.y());

        if x_min >= x_max || y_min >= y_max {
            return;
        }

        let bbox_area = (x_max - x_min) * (y_max - y_min);

        for _ in 0..n_samples {
            let x = x_min + (x_max - x_min) * rng.random::<f64>();
            let y = y_min + (y_max - y_min) * rng.random::<f64>();
            let p = Point::new(x, y);

            if e1.contains_point(&p) && e2.contains_point(&p) && e3.contains_point(&p) {
                in_all_three += 1;
            }
        }

        let mc_area = (in_all_three as f64 / n_samples as f64) * bbox_area;

        // Compute using compute_exclusive_regions
        use crate::geometry::traits::DiagramShape;
        let areas = Ellipse::compute_exclusive_regions(&[e1, e2, e3]);
        let mask_all = (1usize << 0) | (1usize << 1) | (1usize << 2);
        let exact_area = areas.get(&mask_all).copied().unwrap_or(0.0);

        if mc_area < 0.01 && exact_area < 0.01 {
            return;
        }

        let error = (exact_area - mc_area).abs() / mc_area.max(exact_area);

        assert!(
            error < 0.15,
            "Three-ellipse Monte Carlo error too large: {:.1}%",
            error * 100.0
        );
    }
}
