eunoia 0.4.0

//! Circle shape implementation.

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

use crate::geometry::diagram::IntersectionPoint;
use crate::geometry::primitives::point;
use crate::geometry::primitives::Point;
use crate::geometry::shapes::{Polygon, Rectangle};
use crate::geometry::traits::{
    Area, BoundingBox, Centroid, Closed, DiagramShape, Distance, Perimeter, Polygonize,
};
use argmin::core::{CostFunction, Error, Executor, State};
use argmin::solver::brent::BrentOpt;

/// A circle defined by a center point and radius.
///
/// Circles are the simplest shape for Euler and Venn diagrams and are often
/// sufficient for many use cases. They have the advantage of being rotationally
/// symmetric, which simplifies some computations.
///
/// # Examples
///
/// ```
/// use eunoia::geometry::shapes::Circle;
/// use eunoia::geometry::traits::{Area, Closed};
/// use eunoia::geometry::primitives::Point;
///
/// let c1 = Circle::new(Point::new(0.0, 0.0), 2.0);
/// let c2 = Circle::new(Point::new(3.0, 0.0), 1.0);
///
/// let area1 = c1.area();
/// let overlap = c1.intersection_area(&c2);
/// ```
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Circle {
    center: Point,
    radius: f64,
}

impl Area for Circle {
    /// Computes the area of the circle using the formula A = πr².
    fn area(&self) -> f64 {
        PI * self.radius * self.radius
    }
}

impl Centroid for Circle {
    fn centroid(&self) -> Point {
        self.center
    }
}

impl Distance for Circle {
    /// Computes the minimum distance between the boundaries of two circles.
    ///
    /// Returns 0.0 if the circles overlap or touch.
    fn distance(&self, other: &Self) -> f64 {
        let center_distance = self.center.distance(&other.center);
        let radius_sum = self.radius + other.radius;

        if center_distance > radius_sum {
            center_distance - radius_sum
        } else {
            0.0
        }
    }
}

impl Perimeter for Circle {
    /// Compute the perimeter of the circle.
    fn perimeter(&self) -> f64 {
        2.0 * PI * self.radius
    }
}

impl BoundingBox for Circle {
    fn bounding_box(&self) -> Rectangle {
        let width = 2.0 * self.radius;
        let height = 2.0 * self.radius;

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

impl Closed for Circle {
    fn contains(&self, other: &Self) -> bool {
        let center_distance = self.center.distance(&other.center);
        center_distance + other.radius <= self.radius
    }

    fn contains_point(&self, point: &Point) -> bool {
        let dist = self.center.distance(point);
        dist <= self.radius
    }

    /// Checks if two circles intersect (share any common points).
    ///
    /// Note: This implementation returns `true` if circles are separate,
    /// which appears to be inverted from the typical definition. This may
    /// need correction.
    fn intersects(&self, other: &Self) -> bool {
        let center_distance = self.center.distance(&other.center);
        // Circles intersect if distance is less than sum of radii
        // Also check for containment case (distance < |r1 - r2|)
        center_distance < self.radius + other.radius
            && center_distance > (self.radius - other.radius).abs()
    }

    /// Computes the area of intersection between two circles.
    ///
    /// Uses the standard geometric formula for circle-circle intersection:
    /// - Returns 0 if circles don't overlap
    /// - Returns area of smaller circle if one contains the other
    /// - Otherwise computes the lens-shaped intersection area
    ///
    /// # Algorithm
    ///
    /// For two circles with radii r1 and r2 separated by distance d, the
    /// intersection area is computed using the formula involving circular
    /// segments from both circles.
    fn intersection_area(&self, other: &Self) -> f64 {
        let d = self.center.distance(&other.center);

        if d >= self.radius + other.radius {
            return 0.0; // No intersection
        }

        if d <= (self.radius - other.radius).abs() {
            // One circle is completely inside the other
            let smaller_radius = self.radius.min(other.radius);
            return PI * smaller_radius * smaller_radius;
        }

        let r1 = self.radius;
        let r2 = other.radius;

        let part1 = r1 * r1 * (((d * d + r1 * r1 - r2 * r2) / (2.0 * d * r1)).acos());
        let part2 = r2 * r2 * (((d * d + r2 * r2 - r1 * r1) / (2.0 * d * r2)).acos());
        let part3 = 0.5 * ((r1 + r2 - d) * (d + r1 - r2) * (d - r1 + r2) * (d + r1 + r2)).sqrt();

        part1 + part2 - part3
    }

    /// Computes the points of intersection between two circles.
    fn intersection_points(&self, other: &Self) -> Vec<Point> {
        let d = self.center.distance(&other.center);

        if d > self.radius + other.radius || d < (self.radius - other.radius).abs() {
            return vec![]; // No intersection points
        }

        let a = (self.radius * self.radius - other.radius * other.radius + d * d) / (2.0 * d);
        let h = (self.radius * self.radius - a * a).sqrt();

        let p2_x = self.center.x() + a * (other.center.x() - self.center.x()) / d;
        let p2_y = self.center.y() + a * (other.center.y() - self.center.y()) / d;

        let rx = -(other.center.y() - self.center.y()) * (h / d);
        let ry = (other.center.x() - self.center.x()) * (h / d);

        let intersection1 = Point::new(p2_x + rx, p2_y + ry);
        let intersection2 = Point::new(p2_x - rx, p2_y - ry);

        if intersection1 == intersection2 {
            vec![intersection1]
        } else {
            vec![intersection1, intersection2]
        }
    }
}

impl DiagramShape for Circle {
    fn compute_exclusive_regions(
        shapes: &[Self],
    ) -> std::collections::HashMap<crate::geometry::diagram::RegionMask, f64> {
        crate::geometry::diagram::compute_exclusive_regions(shapes)
    }

    fn params_from_circle(x: f64, y: f64, radius: f64) -> Vec<f64> {
        vec![x, y, radius]
    }

    fn n_params() -> usize {
        3 // x, y, radius
    }

    fn from_params(params: &[f64]) -> Self {
        debug_assert_eq!(
            params.len(),
            3,
            "Circle requires 3 parameters: x, y, radius"
        );
        Circle::new(Point::new(params[0], params[1]), params[2])
    }

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

    fn compute_exclusive_regions_with_gradient(
        shapes: &[Self],
    ) -> Option<crate::geometry::traits::ExclusiveRegionsAndGradient> {
        Some(crate::geometry::diagram::compute_exclusive_regions_with_gradient_circles(shapes))
    }
}

impl Polygonize for Circle {
    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);

        for i in 0..n_vertices {
            let angle = 2.0 * PI * (i as f64) / (n_vertices as f64);
            let x = self.center.x() + self.radius * angle.cos();
            let y = self.center.y() + self.radius * angle.sin();
            vertices.push(Point::new(x, y));
        }

        Polygon::new(vertices)
    }
}

struct SeparationCost {
    r1: f64,
    r2: f64,
    target_overlap: f64,
}

impl CostFunction for SeparationCost {
    type Param = f64;
    type Output = f64;

    fn cost(&self, distance: &Self::Param) -> Result<Self::Output, Error> {
        let c1 = Circle::new(Point::new(0.0, 0.0), self.r1);
        let c2 = Circle::new(Point::new(*distance, 0.0), self.r2);

        let current_overlap = c1.intersection_area(&c2);
        let cost = (current_overlap - self.target_overlap).powi(2);

        Ok(cost)
    }
}

impl Circle {
    /// Creates a new circle with the specified center and radius.
    ///
    /// # Arguments
    ///
    /// * `center` - The center point of the circle
    /// * `radius` - The radius of the circle (must be positive)
    ///
    /// # Examples
    ///
    /// ```
    /// use eunoia::geometry::shapes::Circle;
    /// use eunoia::geometry::primitives::Point;
    ///
    /// let circle = Circle::new(Point::new(1.0, 2.0), 3.0);
    /// ```
    pub fn new(center: Point, radius: f64) -> Self {
        Circle { center, radius }
    }

    /// Returns a reference to the circle's center point.
    pub fn center(&self) -> &Point {
        &self.center
    }

    /// Returns the circle's radius.
    pub fn radius(&self) -> f64 {
        self.radius
    }

    /// Sets the center of the circle.
    pub fn set_center(&mut self, center: Point) {
        self.center = center;
    }

    /// Circular sector area given the radius and angle in radians.
    pub fn sector_area(&self, angle_rad: f64) -> f64 {
        0.5 * self.radius * self.radius * angle_rad
    }

    /// Circular segment area from radius and angle
    pub fn segment_area_from_angle(&self, angle_rad: f64) -> f64 {
        0.5 * self.radius * self.radius * (angle_rad - angle_rad.sin())
    }

    /// Circular segment area from radius and chord length
    ///
    /// # Panics
    ///
    /// Panics in debug builds if `chord_length > 2 * radius` (impossible geometry).
    /// In release builds, this validation is skipped for performance.
    pub fn segment_area_from_chord(&self, chord_length: f64) -> f64 {
        let r = self.radius;
        debug_assert!(
            chord_length <= 2.0 * r,
            "Chord length {} cannot exceed diameter {}",
            chord_length,
            2.0 * r
        );
        let theta = 2.0 * (chord_length / (2.0 * r)).asin();
        self.segment_area_from_angle(theta)
    }

    pub fn segment_area_from_points(&self, p1: &Point, p2: &Point) -> f64 {
        let chord_length = p1.distance(p2);
        self.segment_area_from_chord(chord_length)
    }
}

/// Computes the distance required between two circles to achieve a specified overlap area.
pub(crate) fn distance_for_overlap(
    r1: f64,
    r2: f64,
    overlap: f64,
    tol: Option<f64>,
    max_iter: Option<u64>,
) -> Result<f64, Error> {
    let min_distance = (r1 - r2).abs();
    let max_distance = r1 + r2;

    // If the desired overlap is zero, then the circles should
    // at most be touching.
    if overlap <= 0.0 {
        return Ok(max_distance);
    }

    let cost_fun = SeparationCost {
        r1,
        r2,
        target_overlap: overlap,
    };

    let solver = BrentOpt::new(min_distance, max_distance);

    let result = Executor::new(cost_fun, solver)
        .configure(|state| {
            state
                .max_iters(max_iter.unwrap_or(1000))
                .target_cost(tol.unwrap_or(f64::EPSILON.sqrt())) // Match eulerr: sqrt(machine epsilon)
        })
        .run()?;

    Ok(*result.state.get_best_param().unwrap())
}

#[deprecated(
    since = "0.3.1",
    note = "Returns wrong area when one circle in the mask contains the others' lens. Use `crate::geometry::diagram::compute_exclusive_regions` or the boundary-arc helpers (`region_boundary_arcs` + `area_from_boundary_arcs`) instead."
)]
pub fn multiple_overlap_areas(circles: &[Circle], points: &[IntersectionPoint]) -> f64 {
    let n_circles = circles.len();

    // Filter to only points that are in ALL circles (full intersection)
    let full_intersection_points: Vec<&IntersectionPoint> = points
        .iter()
        .filter(|ip| ip.adopters().len() == n_circles)
        .collect();

    if full_intersection_points.is_empty() {
        return 0.0;
    }

    let n_points = full_intersection_points.len();

    // Sort the points by their angles around the centroid
    let centroid = point::centroid(
        &full_intersection_points
            .iter()
            .map(|ip| *ip.point())
            .collect::<Vec<Point>>(),
    );

    let mut indices: Vec<usize> = (0..n_points).collect();
    indices.sort_by(|&i, &j| {
        full_intersection_points[i]
            .point()
            .angle_to(&centroid)
            .partial_cmp(&full_intersection_points[j].point().angle_to(&centroid))
            .unwrap_or(std::cmp::Ordering::Less)
    });

    let mut area = 0.0;

    let mut l = n_points - 1;

    for k in 0..n_points {
        let i = indices[k];
        let j = indices[l];

        let p1 = &full_intersection_points[i].point();
        let p2 = &full_intersection_points[j].point();

        // Now we need to discover which of the circles the two points are
        // coming from so that we can compute the segment area.
        // This should be the set intersection of the parents of both points.
        // In some cases, the intersection may be of length 2, in which
        // case we need to compute both segment areas and pick the
        // smaller one.
        let parents1 = &full_intersection_points[i].parents();
        let parents2 = &full_intersection_points[j].parents();

        let common_parents: Vec<usize> = vec![parents1.0, parents1.1]
            .into_iter()
            .filter(|p| *p == parents2.0 || *p == parents2.1)
            .collect();

        let mut segment_areas = Vec::with_capacity(common_parents.len());

        if common_parents.is_empty() {
            // This should not happen in a well-formed set of intersection points
            panic!("No common parent circles found for intersection points");
        }

        for &circle_index in &common_parents {
            let circle = &circles[circle_index];
            let seg_area = circle.segment_area_from_points(p1, p2);

            debug_assert!(seg_area >= 0.0, "Segment area should be non-negative");

            segment_areas.push(seg_area);
        }

        let triangle_area = 0.5 * ((p1.x() + p2.x()) * (p1.y() - p2.y()));
        // Note: triangle_area can be negative (signed area from shoelace algorithm)
        // We take abs() at the end to get the final area

        let min_segment = segment_areas
            .into_iter()
            .fold(f64::INFINITY, |a, b| a.min(b));

        area += triangle_area;
        area += min_segment;

        l = k;
    }

    area.abs()
}

/// Compute the area of the intersection region for a subset of circles.
///
/// This is similar to `multiple_overlap_areas` but allows specifying which circles
/// to consider for the "full intersection". This is needed for computing 3-way
/// intersections in a 4+ circle diagram where some intersection points may be
/// in more than just the 3 circles of interest.
///
/// # Arguments
/// * `circles` - All circles in the diagram
/// * `points` - Intersection points (with adopters referencing indices in `circles`)
/// * `circle_indices` - Indices of the circles that define this region
#[deprecated(
    since = "0.3.1",
    note = "Returns wrong area when one circle in the mask contains the others' lens. Use `crate::geometry::diagram::compute_exclusive_regions` or the boundary-arc helpers (`region_boundary_arcs` + `area_from_boundary_arcs`) instead."
)]
pub fn multiple_overlap_areas_with_mask(
    circles: &[Circle],
    points: &[IntersectionPoint],
    circle_indices: &[usize],
) -> f64 {
    // Filter to only points that are in ALL of the specified circles
    // A point is in the region if all circle_indices are present in its adopters
    let region_points: Vec<&IntersectionPoint> = points
        .iter()
        .filter(|ip| {
            circle_indices
                .iter()
                .all(|&idx| ip.adopters().contains(&idx))
        })
        .collect();

    if region_points.is_empty() {
        return 0.0;
    }

    let n_points = region_points.len();

    // Sort the points by their angles around the centroid
    let centroid = point::centroid(
        &region_points
            .iter()
            .map(|ip| *ip.point())
            .collect::<Vec<Point>>(),
    );

    let mut indices: Vec<usize> = (0..n_points).collect();
    indices.sort_by(|&i, &j| {
        region_points[i]
            .point()
            .angle_to(&centroid)
            .partial_cmp(&region_points[j].point().angle_to(&centroid))
            .unwrap_or(std::cmp::Ordering::Less)
    });

    let mut area = 0.0;

    let mut l = n_points - 1;

    for k in 0..n_points {
        let i = indices[k];
        let j = indices[l];

        let p1 = &region_points[i].point();
        let p2 = &region_points[j].point();

        // Find which of the region circles these points come from
        let parents1 = &region_points[i].parents();
        let parents2 = &region_points[j].parents();

        let common_parents: Vec<usize> = vec![parents1.0, parents1.1]
            .into_iter()
            .filter(|p| *p == parents2.0 || *p == parents2.1)
            .filter(|p| circle_indices.contains(p)) // Only consider circles in our region
            .collect();

        let mut segment_areas = Vec::with_capacity(common_parents.len());

        if common_parents.is_empty() {
            // Try to find any circle in the region that contains both points
            // This can happen when points come from different pairs but are connected
            // through the region
            for &circle_idx in circle_indices {
                let circle = &circles[circle_idx];
                if circle.contains_point(p1) && circle.contains_point(p2) {
                    let seg_area = circle.segment_area_from_points(p1, p2);
                    segment_areas.push(seg_area);
                }
            }

            if segment_areas.is_empty() {
                // No circle in the region connects these points - use straight line
                // This shouldn't normally happen in a well-formed region
                l = k;
                continue;
            }
        } else {
            for &circle_index in &common_parents {
                let circle = &circles[circle_index];
                let seg_area = circle.segment_area_from_points(p1, p2);

                debug_assert!(seg_area >= 0.0, "Segment area should be non-negative");

                segment_areas.push(seg_area);
            }
        }

        let triangle_area = 0.5 * ((p1.x() + p2.x()) * (p1.y() - p2.y()));

        let min_segment = segment_areas
            .into_iter()
            .fold(f64::INFINITY, |a, b| a.min(b));

        area += triangle_area;
        area += min_segment;

        l = k;
    }

    area.abs()
}

/// One circular arc on the boundary of a region, oriented to be traversed in the
/// CCW direction around the region.
///
/// `phi_start` and `phi_end` are the standard `atan2` angles (in `(-π, π]`) of
/// the arc endpoints relative to the owning circle's centre. `delta_phi` is the
/// signed short-arc angular delta from `phi_start` to `phi_end` in `(-π, π]`
/// (a full-circle arc uses `2π`); its sign matches the traversal direction
/// around the owning circle (`+1` = CCW, `−1` = CW). The pair `(phi_start,
/// delta_phi)` is sufficient to drive the boundary integral; `phi_end` is
/// stored only as a convenience for the closed-form gradient formulas.
#[derive(Debug, Clone, Copy)]
pub(crate) struct BoundaryArc {
    pub circle: usize,
    pub phi_start: f64,
    pub phi_end: f64,
    pub delta_phi: f64,
}

/// Implicit-form value `((p.x − c.x)² + (p.y − c.y)²) / r²` for point `p` and
/// circle `c`. `< 1` strictly inside, `= 1` on the boundary, `> 1` outside.
#[inline]
fn circle_implicit_value(p: &Point, c: &Circle) -> f64 {
    let dx = p.x() - c.center().x();
    let dy = p.y() - c.center().y();
    (dx * dx + dy * dy) / (c.radius() * c.radius())
}

/// Decide whether an arc whose midpoint is on `∂C_j` is owned by `j` for
/// the purpose of region-boundary contributions. The midpoint must be inside
/// every other mask circle; when it lies on another mask circle's boundary
/// (boundaries coincide), the smaller-index circle owns the arc to avoid
/// double-counting (the identical-circles case in particular would otherwise
/// emit one full-circle arc per circle).
///
/// Tolerance is scaled per comparator circle: a geometric near-tangency of
/// `~δ` translates to an implicit-value deviation of `~2δ/r`, so we pick
/// `eps_l = 2·BOUNDARY_COINCIDENCE_GEOM_TOL / r_l` (clamped against
/// floating-point noise) so the threshold corresponds to a roughly constant
/// geometric gap regardless of circle scale.
fn arc_midpoint_owned_by_j(j: usize, mid: &Point, indices: &[usize], circles: &[Circle]) -> bool {
    for &l in indices {
        if l == j {
            continue;
        }
        let cl = &circles[l];
        let eps = (2.0 * BOUNDARY_COINCIDENCE_GEOM_TOL / cl.radius()).max(IMPLICIT_VALUE_FP_TOL);
        let v = circle_implicit_value(mid, cl);
        if v > 1.0 + eps {
            return false;
        }
        if (v - 1.0).abs() <= eps && l < j {
            return false;
        }
    }
    true
}

/// Geometric tolerance (in world-frame distance units) for treating two
/// circle boundaries as coincident at a probe point. Converted to a per-circle
/// implicit-value tolerance inside the tiebreaker.
const BOUNDARY_COINCIDENCE_GEOM_TOL: f64 = 1e-7;

/// Floor for the per-circle implicit-value tolerance, guarding against
/// floating-point noise on huge circles where the geometric tolerance would
/// otherwise map to a sub-ulp implicit-value threshold.
const IMPLICIT_VALUE_FP_TOL: f64 = 1e-12;

/// Build the CCW-oriented list of boundary arcs for an overlapping region.
///
/// For each circle in the mask, the function identifies the IPs that lie on
/// that circle's boundary and inside every other circle in the mask, sorts
/// them by angle around the circle's centre, and emits an arc for every
/// inter-IP segment whose midpoint sits inside every other mask circle. When
/// no such IPs exist on a circle, a probe point is used to decide whether the
/// circle contributes a full-circle arc (i.e. it is wholly inside every
/// other) or no arc.
///
/// All emitted arcs have `delta_phi > 0` (CCW around the owning circle), since
/// for a region that is the intersection of disks, R sits inside each
/// boundary circle and CCW-around-Cₖ on the boundary coincides with CCW
/// around R. This is robust across degenerate cases that the older
/// "polygon + min-segment" decomposition would miss — most notably 3+-way
/// regions where the IPs all come from a single pair of circles because a
/// third circle fully contains their lens.
///
/// When circle boundaries coincide (identical circles, or arcs shared between
/// two circles), an index tiebreaker — `arc_midpoint_owned_by_j` — hands the
/// arc to the smallest-index circle so the area is not multi-counted.
pub(crate) fn region_boundary_arcs(
    mask: crate::geometry::diagram::RegionMask,
    circles: &[Circle],
    intersections: &[IntersectionPoint],
    n_sets: usize,
) -> Vec<BoundaryArc> {
    use crate::geometry::diagram::{adopters_to_mask, mask_to_indices};
    let circle_count = mask.count_ones();
    if circle_count == 0 {
        return Vec::new();
    }
    if circle_count == 1 {
        let idx = mask.trailing_zeros() as usize;
        return vec![BoundaryArc {
            circle: idx,
            phi_start: 0.0,
            phi_end: 0.0,
            delta_phi: 2.0 * PI,
        }];
    }

    let indices = mask_to_indices(mask, n_sets);
    let mut arcs = Vec::new();
    let two_pi = 2.0 * PI;
    // Numerical guard for collapsing arcs at tangencies / shared IPs.
    let arc_eps = 1e-9;

    for &j in &indices {
        let cj = &circles[j];
        // IPs on ∂C_j (j as a parent) that are inside every other mask circle.
        let mut j_phis: Vec<f64> = intersections
            .iter()
            .filter(|ip| {
                let (p1, p2) = ip.parents();
                if p1 != j && p2 != j {
                    return false;
                }
                let am = adopters_to_mask(ip.adopters());
                (mask & am) == mask
            })
            .map(|ip| {
                let p = ip.point();
                (p.y() - cj.center().y()).atan2(p.x() - cj.center().x())
            })
            .collect();

        if j_phis.is_empty() {
            // Either C_j is wholly inside every other mask circle (full-circle
            // arc) or wholly outside the region (no arc). Probe at φ = 0; if
            // the probe lies on ∂C_l for some l ≠ j (boundaries coincide),
            // the index tiebreaker hands the arc to the smallest-index circle
            // so identical / coincident-boundary circles don't all emit
            // duplicate full-circle arcs.
            let probe = Point::new(cj.center().x() + cj.radius(), cj.center().y());
            if arc_midpoint_owned_by_j(j, &probe, &indices, circles) {
                arcs.push(BoundaryArc {
                    circle: j,
                    phi_start: 0.0,
                    phi_end: 0.0,
                    delta_phi: two_pi,
                });
            }
            continue;
        }

        j_phis.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        let m = j_phis.len();

        for k in 0..m {
            let phi_a = j_phis[k];
            let phi_b = j_phis[(k + 1) % m];
            // CCW arc length from phi_a to phi_b.
            let mut delta = phi_b - phi_a;
            if delta <= 0.0 {
                delta += two_pi;
            }
            if delta < arc_eps || delta > two_pi - arc_eps {
                continue;
            }
            // Midpoint of the candidate arc; include only if it sits inside
            // every other circle in the mask (i.e. on R's boundary). The
            // index tiebreaker handles boundary-coincident arcs the same way
            // as the empty-IP branch above.
            let phi_mid = phi_a + delta * 0.5;
            let mid = Point::new(
                cj.center().x() + cj.radius() * phi_mid.cos(),
                cj.center().y() + cj.radius() * phi_mid.sin(),
            );
            if !arc_midpoint_owned_by_j(j, &mid, &indices, circles) {
                continue;
            }
            // Canonicalise phi_end to (-π, π] so sin/cos in the area &
            // gradient formulas operate on raw atan2 values; the periodicity
            // of sin/cos makes phi_end == phi_a + delta (continuous) and
            // phi_end (canonical) interchangeable.
            let phi_end_canonical = wrap_angle(phi_a + delta);
            arcs.push(BoundaryArc {
                circle: j,
                phi_start: phi_a,
                phi_end: phi_end_canonical,
                delta_phi: delta,
            });
        }
    }

    arcs
}

/// Compute the area enclosed by a CCW boundary arc list via the line-integral
/// `A = (1/2) ∮ (x dy − y dx)`. For an arc on circle k parametrised by φ:
/// ```text
/// x dy − y dx = (xₖ rₖ cos φ + yₖ rₖ sin φ + rₖ²) dφ
/// ```
/// integrating from `phi_start` to `phi_start + delta_phi` (continuous).
pub(crate) fn area_from_boundary_arcs(arcs: &[BoundaryArc], circles: &[Circle]) -> f64 {
    let mut total = 0.0;
    for arc in arcs {
        let cj = &circles[arc.circle];
        let xk = cj.center().x();
        let yk = cj.center().y();
        let r = cj.radius();
        let phi_a = arc.phi_start;
        let phi_b = phi_a + arc.delta_phi;
        total += xk * r * (phi_b.sin() - phi_a.sin());
        total += yk * r * (phi_a.cos() - phi_b.cos());
        total += r * r * arc.delta_phi;
    }
    0.5 * total
}

/// Wrap an angle into `(-π, π]`.
#[inline]
pub(crate) fn wrap_angle(x: f64) -> f64 {
    let two_pi = 2.0 * PI;
    let y = x.rem_euclid(two_pi);
    if y > PI {
        y - two_pi
    } else {
        y
    }
}

/// Accumulate the gradient of an overlapping region's area into `grad`, where
/// `grad` is a length-`3 · n_sets` vector laid out as `[x₀, y₀, r₀, x₁, …]`.
///
/// Each boundary arc on circle `k` contributes:
/// ```text
/// ∂A/∂xₖ += sign · rₖ · (sin φ_end − sin φ_start)
/// ∂A/∂yₖ −= sign · rₖ · (cos φ_end − cos φ_start)
/// ∂A/∂rₖ += rₖ · |delta_phi|
/// ```
/// where `sign = sign(delta_phi)`, derived from the boundary-velocity identity
/// `dA/dθ = ∮_∂R (v_θ · n) ds` with `θ ∈ {xₖ, yₖ, rₖ}`.
pub(crate) fn accumulate_region_overlap_gradient(
    arcs: &[BoundaryArc],
    circles: &[Circle],
    grad: &mut [f64],
) {
    for arc in arcs {
        let r = circles[arc.circle].radius();
        let sign = arc.delta_phi.signum();
        let off = arc.circle * 3;
        grad[off] += sign * r * (arc.phi_end.sin() - arc.phi_start.sin());
        grad[off + 1] -= sign * r * (arc.phi_end.cos() - arc.phi_start.cos());
        grad[off + 2] += r * arc.delta_phi.abs();
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    const EPSILON: f64 = 1e-10;

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

    #[test]
    fn test_circle_new() {
        let center = Point::new(1.0, 2.0);
        let circle = Circle::new(center, 5.0);
        assert_eq!(circle.radius(), 5.0);
        assert_eq!(circle.center().x(), 1.0);
        assert_eq!(circle.center().y(), 2.0);
    }

    #[test]
    fn test_circle_area() {
        let circle = Circle::new(Point::new(0.0, 0.0), 1.0);
        assert!(approx_eq(circle.area(), PI));

        let circle2 = Circle::new(Point::new(0.0, 0.0), 2.0);
        assert!(approx_eq(circle2.area(), 4.0 * PI));

        let circle3 = Circle::new(Point::new(5.0, 5.0), 3.0);
        assert!(approx_eq(circle3.area(), 9.0 * PI));
    }

    #[test]
    fn test_circle_distance_no_overlap() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 1.0);
        let circle2 = Circle::new(Point::new(5.0, 0.0), 1.0);
        assert_eq!(circle1.distance(&circle2), 3.0);
    }

    #[test]
    fn test_circle_distance_touching() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 1.0);
        let circle2 = Circle::new(Point::new(2.0, 0.0), 1.0);
        assert_eq!(circle1.distance(&circle2), 0.0);
    }

    #[test]
    fn test_circle_distance_overlapping() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 2.0);
        let circle2 = Circle::new(Point::new(1.0, 0.0), 2.0);
        assert_eq!(circle1.distance(&circle2), 0.0);
    }

    #[test]
    fn test_circle_contains_smaller() {
        let large = Circle::new(Point::new(0.0, 0.0), 5.0);
        let small = Circle::new(Point::new(1.0, 1.0), 2.0);
        assert!(large.contains(&small));
    }

    #[test]
    fn test_circle_contains_self() {
        let circle = Circle::new(Point::new(0.0, 0.0), 3.0);
        assert!(circle.contains(&circle));
    }

    #[test]
    fn test_circle_not_contains() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 2.0);
        let circle2 = Circle::new(Point::new(5.0, 0.0), 2.0);
        assert!(!circle1.contains(&circle2));
    }

    #[test]
    fn test_circle_not_contains_partial_overlap() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 3.0);
        let circle2 = Circle::new(Point::new(2.0, 0.0), 2.0);
        assert!(!circle1.contains(&circle2));
    }

    #[test]
    fn test_circle_intersects_separate() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 1.0);
        let circle2 = Circle::new(Point::new(5.0, 0.0), 1.0);
        // Distance = 5, sum of radii = 2, circles are separated
        assert!(!circle1.intersects(&circle2));
    }

    #[test]
    fn test_circle_intersects_touching() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 1.0);
        let circle2 = Circle::new(Point::new(2.0, 0.0), 1.0);
        // Distance = 2, sum of radii = 2, circles touch at exactly one point
        // This is a boundary case - could be either true or false depending on definition
        // We treat tangent as NOT intersecting (no area overlap)
        assert!(!circle1.intersects(&circle2));
    }

    #[test]
    fn test_circle_intersects_overlapping() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 2.0);
        let circle2 = Circle::new(Point::new(1.0, 0.0), 2.0);
        // Distance = 1, sum of radii = 4, circles overlap
        assert!(circle1.intersects(&circle2));
    }

    #[test]
    fn test_intersection_area_no_overlap() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 1.0);
        let circle2 = Circle::new(Point::new(10.0, 0.0), 1.0);
        assert_eq!(circle1.intersection_area(&circle2), 0.0);
    }

    #[test]
    fn test_intersection_area_touching() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 1.0);
        let circle2 = Circle::new(Point::new(2.0, 0.0), 1.0);
        let area = circle1.intersection_area(&circle2);
        assert!(approx_eq(area, 0.0));
    }

    #[test]
    fn test_intersection_area_complete_overlap_same_size() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 2.0);
        let circle2 = Circle::new(Point::new(0.0, 0.0), 2.0);
        let expected = PI * 4.0;
        assert!(approx_eq(circle1.intersection_area(&circle2), expected));
    }

    #[test]
    fn test_intersection_area_one_inside_other() {
        let large = Circle::new(Point::new(0.0, 0.0), 5.0);
        let small = Circle::new(Point::new(1.0, 0.0), 2.0);
        let expected = PI * 4.0; // Area of smaller circle
        assert!(approx_eq(large.intersection_area(&small), expected));
        assert!(approx_eq(small.intersection_area(&large), expected));
    }

    #[test]
    fn test_intersection_area_partial_overlap() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 1.0);
        let circle2 = Circle::new(Point::new(1.0, 0.0), 1.0);
        let area = circle1.intersection_area(&circle2);

        // For two unit circles with centers 1 apart, there's a known formula
        // The intersection area should be positive and less than π
        assert!(area > 0.0);
        assert!(area < PI);
    }

    #[test]
    fn test_intersection_area_symmetric() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 2.0);
        let circle2 = Circle::new(Point::new(1.5, 0.0), 1.5);
        let area1 = circle1.intersection_area(&circle2);
        let area2 = circle2.intersection_area(&circle1);
        assert!(approx_eq(area1, area2));
    }

    #[test]
    fn test_intersection_area_different_sizes() {
        let circle1 = Circle::new(Point::new(0.0, 0.0), 3.0);
        let circle2 = Circle::new(Point::new(2.0, 0.0), 1.0);
        let area = circle1.intersection_area(&circle2);

        // Should be positive and at most the smaller circle's area
        assert!(area > 0.0);
        assert!(area <= PI * 1.0 * 1.0);
    }

    #[test]
    fn test_distance_for_overlap_zero_overlap() {
        let r1 = 2.0;
        let r2 = 1.5;
        let overlap = 0.0;

        let distance = distance_for_overlap(r1, r2, overlap, None, None).unwrap();

        // Should return the sum of radii (circles just touching)
        assert!(approx_eq(distance, r1 + r2));
    }

    #[test]
    fn test_distance_for_overlap_negative_overlap() {
        let r1 = 2.0;
        let r2 = 1.5;
        let overlap = -1.0;

        let distance = distance_for_overlap(r1, r2, overlap, None, None).unwrap();

        // Should return the sum of radii (circles separated)
        assert!(approx_eq(distance, r1 + r2));
    }

    #[test]
    fn test_distance_for_overlap_full_overlap() {
        let r1 = 3.0;
        let r2 = 2.0;
        let overlap = PI * r2 * r2; // Full area of smaller circle

        let distance = distance_for_overlap(r1, r2, overlap, None, None).unwrap();

        // Distance should be close to |r1 - r2| (one inside the other)
        assert!(distance <= (r1 - r2).abs() + 0.1); // Allow small tolerance
    }

    #[test]
    fn test_distance_for_overlap_partial_overlap_equal_radii() {
        let r1 = 2.0;
        let r2 = 2.0;
        let target_overlap = 2.0; // Some specific overlap area

        let distance = distance_for_overlap(r1, r2, target_overlap, None, None).unwrap();

        // Verify the result by computing the actual overlap at this distance
        let c1 = Circle::new(Point::new(0.0, 0.0), r1);
        let c2 = Circle::new(Point::new(distance, 0.0), r2);
        let actual_overlap = c1.intersection_area(&c2);

        // Should match target within tolerance (relaxed for optimization convergence)
        assert!((actual_overlap - target_overlap).abs() < 1e-2);
    }

    #[test]
    fn test_distance_for_overlap_partial_overlap_different_radii() {
        let r1 = 3.0;
        let r2 = 1.5;
        let target_overlap = 1.0;

        let distance = distance_for_overlap(r1, r2, target_overlap, None, None).unwrap();

        // Verify the result
        let c1 = Circle::new(Point::new(0.0, 0.0), r1);
        let c2 = Circle::new(Point::new(distance, 0.0), r2);
        let actual_overlap = c1.intersection_area(&c2);

        assert!((actual_overlap - target_overlap).abs() < 1e-2);
    }

    #[test]
    fn test_distance_for_overlap_custom_tolerance() {
        let r1 = 2.0;
        let r2 = 1.0;
        let target_overlap = 0.5;
        let custom_tol = 1e-8;

        let distance =
            distance_for_overlap(r1, r2, target_overlap, Some(custom_tol), None).unwrap();

        let c1 = Circle::new(Point::new(0.0, 0.0), r1);
        let c2 = Circle::new(Point::new(distance, 0.0), r2);
        let actual_overlap = c1.intersection_area(&c2);

        // Should be very close to target with custom tolerance
        assert!((actual_overlap - target_overlap).abs() < 1e-3);
    }

    #[test]
    fn test_distance_for_overlap_custom_max_iter() {
        let r1 = 2.0;
        let r2 = 1.5;
        let target_overlap = 1.0;
        let max_iter = 50;

        let distance = distance_for_overlap(r1, r2, target_overlap, None, Some(max_iter)).unwrap();

        // Should still converge within fewer iterations
        let c1 = Circle::new(Point::new(0.0, 0.0), r1);
        let c2 = Circle::new(Point::new(distance, 0.0), r2);
        let actual_overlap = c1.intersection_area(&c2);

        assert!((actual_overlap - target_overlap).abs() < 1e-2);
    }

    #[test]
    fn test_distance_for_overlap_small_circles() {
        let r1 = 0.5;
        let r2 = 0.3;
        let target_overlap = 0.1;

        let distance = distance_for_overlap(r1, r2, target_overlap, None, None).unwrap();

        let c1 = Circle::new(Point::new(0.0, 0.0), r1);
        let c2 = Circle::new(Point::new(distance, 0.0), r2);
        let actual_overlap = c1.intersection_area(&c2);

        assert!((actual_overlap - target_overlap).abs() < 1e-4);
    }

    #[test]
    fn test_distance_for_overlap_large_circles() {
        let r1 = 100.0;
        let r2 = 75.0;
        let target_overlap = 500.0;

        let distance = distance_for_overlap(r1, r2, target_overlap, None, None).unwrap();

        let c1 = Circle::new(Point::new(0.0, 0.0), r1);
        let c2 = Circle::new(Point::new(distance, 0.0), r2);
        let actual_overlap = c1.intersection_area(&c2);

        assert!((actual_overlap - target_overlap).abs() < 1.0); // Larger tolerance for large circles
    }

    #[test]
    fn test_distance_for_overlap_bounds() {
        let r1 = 2.0;
        let r2 = 1.5;
        let target_overlap = 1.0;

        let distance = distance_for_overlap(r1, r2, target_overlap, None, None).unwrap();

        // Distance should be between min (one inside other) and max (circles touching)
        let min_distance = (r1 - r2).abs();
        let max_distance = r1 + r2;

        assert!(distance >= min_distance);
        assert!(distance <= max_distance);
    }

    #[test]
    fn test_intersection_points_no_intersection() {
        // Circles too far apart
        let c1 = Circle::new(Point::new(0.0, 0.0), 1.0);
        let c2 = Circle::new(Point::new(5.0, 0.0), 1.0);

        let points = c1.intersection_points(&c2);
        assert_eq!(points.len(), 0);
    }

    #[test]
    fn test_intersection_points_one_inside_other() {
        // One circle completely inside the other
        let c1 = Circle::new(Point::new(0.0, 0.0), 3.0);
        let c2 = Circle::new(Point::new(0.0, 0.0), 1.0);

        let points = c1.intersection_points(&c2);
        assert_eq!(points.len(), 0);
    }

    #[test]
    fn test_intersection_points_touching_externally() {
        // Circles touch at exactly one point (externally)
        let c1 = Circle::new(Point::new(0.0, 0.0), 2.0);
        let c2 = Circle::new(Point::new(4.0, 0.0), 2.0);

        let points = c1.intersection_points(&c2);
        assert_eq!(points.len(), 1);

        // The touching point should be at (2.0, 0.0)
        assert!(approx_eq(points[0].x(), 2.0));
        assert!(approx_eq(points[0].y(), 0.0));
    }

    #[test]
    fn test_intersection_points_two_points() {
        // Circles intersect at two points
        let c1 = Circle::new(Point::new(0.0, 0.0), 2.0);
        let c2 = Circle::new(Point::new(2.0, 0.0), 2.0);

        let points = c1.intersection_points(&c2);
        assert_eq!(points.len(), 2);

        // Both points should be on the circles
        for point in &points {
            let dist_to_c1 = c1.center.distance(point);
            let dist_to_c2 = c2.center.distance(point);
            assert!(approx_eq(dist_to_c1, c1.radius));
            assert!(approx_eq(dist_to_c2, c2.radius));
        }

        // Points should be at (1.0, sqrt(3)) and (1.0, -sqrt(3))
        let expected_x = 1.0;
        let expected_y = 3.0_f64.sqrt();

        // Check one point is at (1.0, sqrt(3))
        let found_positive = points
            .iter()
            .any(|p| approx_eq(p.x(), expected_x) && approx_eq(p.y(), expected_y));
        // Check other point is at (1.0, -sqrt(3))
        let found_negative = points
            .iter()
            .any(|p| approx_eq(p.x(), expected_x) && approx_eq(p.y(), -expected_y));

        assert!(found_positive);
        assert!(found_negative);
    }

    #[test]
    fn test_intersection_points_vertical_alignment() {
        // Test with circles aligned vertically
        let c1 = Circle::new(Point::new(0.0, 0.0), 1.5);
        let c2 = Circle::new(Point::new(0.0, 2.0), 1.5);

        let points = c1.intersection_points(&c2);
        assert_eq!(points.len(), 2);

        // Both points should be equidistant from both centers
        for point in &points {
            let dist_to_c1 = c1.center.distance(point);
            let dist_to_c2 = c2.center.distance(point);
            assert!(approx_eq(dist_to_c1, c1.radius));
            assert!(approx_eq(dist_to_c2, c2.radius));
        }
    }

    #[test]
    fn test_intersection_points_equal_radii_partial_overlap() {
        // Two circles with equal radii, partially overlapping
        let c1 = Circle::new(Point::new(0.0, 0.0), 1.0);
        let c2 = Circle::new(Point::new(1.0, 0.0), 1.0);

        let points = c1.intersection_points(&c2);
        assert_eq!(points.len(), 2);

        // The intersection points should be symmetric about the line connecting centers
        // They should both have x-coordinate = 0.5
        for point in &points {
            assert!(approx_eq(point.x(), 0.5));
        }
    }

    #[test]
    fn test_intersection_points_different_radii() {
        // Two circles with different radii
        let c1 = Circle::new(Point::new(0.0, 0.0), 3.0);
        let c2 = Circle::new(Point::new(2.0, 0.0), 1.5);

        let points = c1.intersection_points(&c2);
        assert_eq!(points.len(), 2);

        // Verify both points lie on both circles
        for point in &points {
            let dist_to_c1 = c1.center.distance(point);
            let dist_to_c2 = c2.center.distance(point);
            assert!(approx_eq(dist_to_c1, c1.radius));
            assert!(approx_eq(dist_to_c2, c2.radius));
        }
    }

    #[test]
    fn test_perimeter() {
        let c = Circle::new(Point::new(0.0, 0.0), 1.0);
        let perimeter = c.perimeter();
        assert!(approx_eq(perimeter, 2.0 * PI));
    }

    #[test]
    fn test_perimeter_various_radii() {
        let test_cases = vec![(0.5, PI), (2.0, 4.0 * PI), (10.0, 20.0 * PI)];

        for (radius, expected) in test_cases {
            let c = Circle::new(Point::new(0.0, 0.0), radius);
            assert!(approx_eq(c.perimeter(), expected));
        }
    }

    #[test]
    fn test_sector_area_full_circle() {
        let c = Circle::new(Point::new(0.0, 0.0), 2.0);
        let full_circle_angle = 2.0 * PI;
        let sector = c.sector_area(full_circle_angle);
        assert!(approx_eq(sector, c.area()));
    }

    #[test]
    fn test_sector_area_half_circle() {
        let c = Circle::new(Point::new(0.0, 0.0), 3.0);
        let half_circle_angle = PI;
        let sector = c.sector_area(half_circle_angle);
        assert!(approx_eq(sector, c.area() / 2.0));
    }

    #[test]
    fn test_sector_area_quarter_circle() {
        let c = Circle::new(Point::new(0.0, 0.0), 4.0);
        let quarter_circle_angle = PI / 2.0;
        let sector = c.sector_area(quarter_circle_angle);
        assert!(approx_eq(sector, c.area() / 4.0));
    }

    #[test]
    fn test_sector_area_zero() {
        let c = Circle::new(Point::new(0.0, 0.0), 5.0);
        let sector = c.sector_area(0.0);
        assert!(approx_eq(sector, 0.0));
    }

    #[test]
    fn test_segment_area_from_angle_semicircle() {
        let c = Circle::new(Point::new(0.0, 0.0), 2.0);
        let angle = PI;
        let segment = c.segment_area_from_angle(angle);
        // For a semicircle, segment area equals sector area (half circle)
        assert!(approx_eq(segment, c.area() / 2.0));
    }

    #[test]
    fn test_segment_area_from_angle_zero() {
        let c = Circle::new(Point::new(0.0, 0.0), 3.0);
        let segment = c.segment_area_from_angle(0.0);
        assert!(approx_eq(segment, 0.0));
    }

    #[test]
    fn test_segment_area_from_angle_small() {
        let c = Circle::new(Point::new(0.0, 0.0), 1.0);
        let angle = PI / 6.0; // 30 degrees
        let segment = c.segment_area_from_angle(angle);

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

    #[test]
    fn test_segment_area_from_chord_diameter() {
        let radius = 2.0;
        let c = Circle::new(Point::new(0.0, 0.0), radius);
        let chord_length = 2.0 * radius; // Diameter
        let segment = c.segment_area_from_chord(chord_length);

        // A chord equal to diameter creates a semicircle
        assert!(approx_eq(segment, c.area() / 2.0));
    }

    #[test]
    fn test_segment_area_from_chord_small() {
        let c = Circle::new(Point::new(0.0, 0.0), 5.0);
        let chord_length = 2.0;
        let segment = c.segment_area_from_chord(chord_length);

        // Small chord should create small segment
        assert!(segment > 0.0);
        assert!(segment < c.area() / 4.0);
    }

    #[test]
    fn test_segment_area_from_chord_vs_angle_consistency() {
        let c = Circle::new(Point::new(0.0, 0.0), 3.0);
        let angle = PI / 3.0; // 60 degrees

        // Calculate chord length from angle
        let chord_length = 2.0 * c.radius() * (angle / 2.0).sin();

        let segment_from_angle = c.segment_area_from_angle(angle);
        let segment_from_chord = c.segment_area_from_chord(chord_length);

        // Both methods should give same result
        assert!(approx_eq(segment_from_angle, segment_from_chord));
    }

    #[test]
    fn test_segment_area_relationships() {
        let c = Circle::new(Point::new(0.0, 0.0), 1.0);
        let angle = PI / 4.0; // 45 degrees

        let sector = c.sector_area(angle);
        let segment = c.segment_area_from_angle(angle);

        // Segment area should be less than sector area (triangle is subtracted)
        assert!(segment < sector);

        // For small angles, segment should be much smaller than sector
        let small_angle = 0.1;
        let small_sector = c.sector_area(small_angle);
        let small_segment = c.segment_area_from_angle(small_angle);
        assert!(small_segment < small_sector / 2.0);
    }

    #[test]
    #[should_panic(expected = "Chord length")]
    fn test_segment_area_from_chord_invalid() {
        let c = Circle::new(Point::new(0.0, 0.0), 2.0);
        let chord_length = 5.0; // Impossible: longer than diameter
        c.segment_area_from_chord(chord_length);
    }

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

        // Three identical unit circles. Without the index tiebreaker each
        // circle would emit a full-circle arc and A∩B∩C would return ~3π.
        let c = Circle::new(Point::new(0.0, 0.0), 1.0);
        let areas = Circle::compute_exclusive_regions(&[c, c, c]);

        let mask_all = 0b111;
        let all_three = areas.get(&mask_all).copied().unwrap_or(0.0);

        assert!(
            (all_three - PI).abs() < 1e-6,
            "Complete overlap should give area ~π, got {}",
            all_three
        );
    }

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

        // Two coincident unit circles plus a smaller circle of radius 0.5
        // sharing the same centre. Expected A∩B∩C = π·0.25.
        // Without the tiebreaker, both unit circles emit a full-circle arc
        // (each one sits on the other's boundary), so the inclusion-exclusion
        // pipeline overshoots.
        let big = Circle::new(Point::new(0.0, 0.0), 1.0);
        let small = Circle::new(Point::new(0.0, 0.0), 0.5);
        let areas = Circle::compute_exclusive_regions(&[big, big, small]);

        let mask_all = 0b111;
        let all_three = areas.get(&mask_all).copied().unwrap_or(0.0);
        let expected = PI * 0.25;

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