geonum 0.12.0

// conformal geometric algebra (CGA) extends euclidean space with two extra dimensions
// to enable unified representation of circles, spheres, and conformal transformations
//
// traditional CGA requires:
// - 5D space for 3D euclidean (n+2 dimensions)
// - null vectors to represent points
// - bivectors for circles/spheres
// - versors for conformal transformations
// - complex inner product operations
//
// geonum: conformal geometry emerges from angle-preserving transformations
// circles are just collections of points at equal length from center
// spheres are angle-preserving dilations
// inversions are reciprocal length transformations

use geonum::{Angle, Geonum};
use std::f64::consts::PI;

const EPSILON: f64 = 1e-10;

#[test]
fn it_represents_points_as_null_vectors() {
    // traditional CGA: points are null vectors on the null cone
    // x² = 0 where x = e₀ + x + ½x²e∞
    //
    // geonum: points are just geometric numbers with length and angle
    // "null" property emerges from specific angle relationships

    // in geonum, "null" means the wedge product with itself is zero
    // this happens when vectors are parallel (sin(0) = 0)
    let point = Geonum::new(2.0, 1.0, 3.0); // regular point

    // null property: v ∧ v = 0 (nilpotent)
    let self_wedge = point.wedge(&point);

    // wedge of parallel vectors is zero (sin(0) = 0)
    assert!(self_wedge.mag < EPSILON, "point wedge itself is null");

    // create points at different positions
    let points = vec![
        Geonum::new_from_cartesian(1.0, 0.0),
        Geonum::new_from_cartesian(0.0, 2.0),
        Geonum::new_from_cartesian(3.0, 4.0),
        Geonum::new_from_cartesian(-1.0, -1.0),
    ];

    // all points satisfy null property
    for p in &points {
        let p_wedge_p = p.wedge(p);
        assert!(
            p_wedge_p.mag < EPSILON,
            "all points have null wedge with themselves"
        );
    }

    // the "null cone" in geonum is just the set of all geometric numbers
    // since every geonum wedged with itself gives zero (parallel to itself)

    // traditional CGA needs complex null vector construction
    // geonum gets nullity naturally from angle relationships

    // demonstrate that different grade objects also have null property
    let scalar = Geonum::new(1.0, 0.0, 1.0); // blade 0
    let vector = Geonum::new(1.0, 1.0, 2.0); // blade 1
    let bivector = Geonum::new(1.0, 1.0, 1.0); // blade 2

    assert!(scalar.wedge(&scalar).mag < EPSILON, "scalar is null");
    assert!(vector.wedge(&vector).mag < EPSILON, "vector is null");
    assert!(bivector.wedge(&bivector).mag < EPSILON, "bivector is null");

    println!("null property emerges from parallel angle relationships");
}

#[test]
fn it_represents_circles_through_three_points() {
    // traditional CGA: circle = P₁ ∧ P₂ ∧ P₃ (trivector in 5D conformal space)
    // requires:
    // 1. embed each point as null vector: p → e₀ + p + ½p²e∞ (5 components each)
    // 2. compute triple wedge product in 5D: P₁ ∧ P₂ ∧ P₃ (10 bivector components)
    // 3. extract center/radius from resulting trivector (complex formulas)
    // 4. total storage: 3×5 + 10 = 25 components minimum
    //
    // geonum: just 2 components per point [length, angle], no conformal embedding

    // three non-collinear points - each just 2 components
    let p1 = Geonum::new_from_cartesian(0.0, 0.0); // [0, 0°]
    let p2 = Geonum::new_from_cartesian(4.0, 0.0); // [4, 0°]
    let p3 = Geonum::new_from_cartesian(2.0, 3.0); // [√13, ~56°]

    // IMPROVEMENT 1: collinearity test is one wedge, not triple product
    let v12 = p2 - p1;
    let v13 = p3 - p1;
    let area = v12.wedge(&v13); // O(1) operation

    // traditional CGA would compute (P₁∧P₂∧P₃)² and test if zero
    // geonum: just check wedge length
    assert!(
        area.mag > EPSILON,
        "non-collinear points have non-zero wedge"
    );

    // IMPROVEMENT 2: perpendiculars via rotation, not cross products
    // traditional CGA: compute bivector B = e₁∧e₂, then exp(Bπ/2) for rotation
    // geonum: just add π/2 to angle
    let perp12 = v12.rotate(Angle::new(1.0, 2.0)); // O(1) rotation
    let perp23 = (p3 - p2).rotate(Angle::new(1.0, 2.0));

    // IMPROVEMENT 3: center via cartesian operations (this IS the simplification!)
    // traditional CGA: center = -(C·e∞)/(C·e₀) where C is circle trivector
    // geonum: cartesian addition/subtraction directly gives geometric results

    // midpoints via cartesian average
    let mid12 = (p1 + p2) * Geonum::new(0.5, 0.0, 1.0);
    let mid23 = (p2 + p3) * Geonum::new(0.5, 0.0, 1.0);
    // test midpoints lie between original points
    assert!((mid12 - p1).mag < (p2 - p1).mag, "mid12 between p1 and p2");
    assert!((mid23 - p2).mag < (p3 - p2).mag, "mid23 between p2 and p3");

    // test perpendiculars arent parallel
    let perp_wedge = perp12.wedge(&perp23);
    assert!(
        perp_wedge.mag > EPSILON,
        "perpendicular bisectors not parallel"
    );

    // KEY INSIGHT: in geonum, cartesian operations ARE the geometric operations
    // no need for conformal embeddings or trivector extractions
    // the circumcenter calculation in cartesian IS the simplification

    // convert to cartesian for circumcenter calculation
    let x1 = p1.adj().mag;
    let y1 = p1.opp().mag;
    let x2 = p2.adj().mag;
    let y2 = p2.opp().mag;
    let x3 = p3.adj().mag;
    let y3 = p3.opp().mag;

    // circumcenter formula (this is what CGA tries to abstract)
    let d = 2.0 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2));
    if d.abs() > EPSILON {
        let ux = ((x1 * x1 + y1 * y1) * (y2 - y3)
            + (x2 * x2 + y2 * y2) * (y3 - y1)
            + (x3 * x3 + y3 * y3) * (y1 - y2))
            / d;
        let uy = ((x1 * x1 + y1 * y1) * (x3 - x2)
            + (x2 * x2 + y2 * y2) * (x1 - x3)
            + (x3 * x3 + y3 * y3) * (x2 - x1))
            / d;

        let center = Geonum::new_from_cartesian(ux, uy);
        let radius = (p1 - center).mag;

        // IMPROVEMENT 4: point-on-circle test via distance, not inner product
        // traditional CGA: test P·C = 0 in conformal space (5D inner product)
        // geonum: just compare distances (O(1) operation)
        assert!((p1 - center).mag - radius < EPSILON, "p1 on circle");
        assert!((p2 - center).mag - radius < EPSILON, "p2 on circle");
        assert!((p3 - center).mag - radius < EPSILON, "p3 on circle");
    }

    // COMPLEXITY COMPARISON:
    // traditional CGA in 3D: 2³ = 8 components per multivector
    // traditional CGA in conformal 5D: 2⁵ = 32 components per multivector
    // geonum: always 2 components [length, angle] regardless of dimension

    // demonstrate collinear case - wedge detects it immediately
    let col1 = Geonum::new_from_cartesian(0.0, 0.0);
    let col2 = Geonum::new_from_cartesian(1.0, 1.0);
    let col3 = Geonum::new_from_cartesian(2.0, 2.0);

    let collinear_wedge = (col2 - col1).wedge(&(col3 - col1));
    assert!(
        collinear_wedge.mag < EPSILON,
        "collinear points have zero wedge"
    );

    // geonum ghosts entire conformal machinery:
    // - no e₀, e∞ basis vectors
    // - no null cone embedding
    // - no 5D trivector P₁∧P₂∧P₃
    // - no conformal inner products
    // - no 32-component multivectors

    println!("circle via O(1) ops, not O(2⁵) conformal algebra");
}

#[test]
fn it_represents_spheres_through_four_points() {
    // traditional CGA: sphere = P₁ ∧ P₂ ∧ P₃ ∧ P₄ (4-vector in 5D conformal)
    // requires:
    // 1. embed 4 points as null vectors: p → e₀ + p + ½p²e∞ (5 components each)
    // 2. compute 4-way wedge product: P₁∧P₂∧P₃∧P₄ (massive computation in 32D space)
    // 3. extract center/radius from 4-vector via complex conformal formulas
    // 4. total: 4×5 input components + O(2⁵) intermediate storage and computation
    //
    // geonum: blade structure enables sphere in ANY dimension via geometric operations

    // KEY INSIGHT: geometric entities exist independently of dimensional space
    // because blade field computes projections on demand - no coordinate storage needed

    // four geometric entities - chosen to be equidistant from their average (sphere property)
    let p1 = Geonum::new(3.0, 0.0, 1.0); // [3, 0] - dimension-independent
    let p2 = Geonum::new(3.0, 1.0, 2.0); // [3, π/2] - projects anywhere via blade arithmetic
    let p3 = Geonum::new(3.0, 1.0, 1.0); // [3, π] - no coordinate system needed
    let p4 = Geonum::new(3.0, 3.0, 2.0); // [3, 3π/2] - exists in angle space

    // IMPROVEMENT 1: sphere property via direct geometric relationships
    // traditional CGA: 4-way wedge determines non-coplanarity in conformal space
    // geonum: test relationships via wedge products (works in any dimension)

    let v12 = p2 - p1; // difference in angle space
    let v13 = p3 - p1; // projects to any dimension when observed
    let v14 = p4 - p1; // no need to declare coordinate system

    // wedge tests detect geometric relationships - dimension emerges from observation
    let area1 = v12.wedge(&v13); // oriented area in angle space
    let area2 = v12.wedge(&v14); // dimension computed via blade arithmetic
    let area3 = v13.wedge(&v14); // independent of coordinate system

    assert!(area1.mag > EPSILON, "entities span non-zero area");
    assert!(area2.mag > EPSILON, "geometric relationships exist");
    assert!(area3.mag > EPSILON, "dimension-independent relationships");

    // IMPROVEMENT 2: sphere center via geometric relationships
    // traditional CGA: extract from 32-component conformal 4-vector
    // geonum: direct construction using angle arithmetic

    // "center" is geometric relationship, not coordinate position
    let center_relation = (p1 + p2 + p3 + p4) * Geonum::new(0.25, 0.0, 1.0);

    // distances in angle space - project to any dimension when measured
    let d1 = (p1 - center_relation).mag;
    let d2 = (p2 - center_relation).mag;
    let d3 = (p3 - center_relation).mag;
    let d4 = (p4 - center_relation).mag;

    // IMPROVEMENT 3: sphere test via distance relationships
    // traditional CGA: P·S = 0 inner product in 5D conformal space
    // geonum: geometric relationships in angle space (projects to any dimension)

    let sphere_radius = d1; // use first entity to define radius

    // test sphere property - all entities equidistant from center
    assert!(
        (d1 - sphere_radius).abs() < EPSILON,
        "p1 defines sphere radius"
    );
    assert!(
        (d2 - sphere_radius).abs() < EPSILON,
        "p2 on sphere via distance"
    );
    assert!(
        (d3 - sphere_radius).abs() < EPSILON,
        "p3 maintains distance relationship"
    );
    assert!(
        (d4 - sphere_radius).abs() < EPSILON,
        "p4 completes sphere constraint"
    );

    // DIMENSION INDEPENDENCE TEST: project sphere to different dimensions
    // traditional: would need separate 3D, 4D, nD coordinate systems
    // geonum: same geometric entities project to any dimension

    // project center to dimension 3 (would be z-coordinate in 3D view)
    let center_dim3 = center_relation.project_to_dimension(3);
    assert!(center_dim3.is_finite(), "center projects to dimension 3");

    // project center to dimension 42 (impossible in traditional coordinate systems)
    let center_dim42 = center_relation.project_to_dimension(42);
    assert!(center_dim42.is_finite(), "center projects to dimension 42");

    // project center to dimension 1000000 (would require million components traditionally)
    let center_dim_million = center_relation.project_to_dimension(1_000_000);
    assert!(
        center_dim_million.is_finite(),
        "center projects to dimension 1000000"
    );

    // sphere relationship holds in all dimensions - same geometric entities
    // no need to recompute or store different representations

    // COMPLEXITY COMPARISON:
    // traditional CGA: O(2⁵) components, dimension-specific coordinate systems
    // geonum: O(1) operations, dimension-independent blade arithmetic

    // test degenerate case: collinear entities
    let col1 = Geonum::new(1.0, 0.0, 1.0); // [1, 0]
    let col2 = Geonum::new(2.0, 0.0, 1.0); // [2, 0] - same angle
    let col3 = Geonum::new(3.0, 0.0, 1.0); // [3, 0] - collinear in angle space

    let collinear_test = (col2 - col1).wedge(&(col3 - col1));
    assert!(
        collinear_test.mag < EPSILON,
        "collinear entities have zero wedge"
    );

    // KEY INSIGHT: blade structure eliminates dimensional constraints
    // - geometric entities exist in angle space, independent of coordinate systems
    // - sphere relationships computed via angle arithmetic
    // - projections to any dimension computed on demand via blade field
    // - no coordinate storage, no dimensional limits, no basis vector requirements
    // - same O(1) operations work for 3D sphere, 4D hypersphere, million-D sphere

    // geonum ghosts traditional CGA dimensional scaffolding:
    // - no predefined coordinate systems (3D, 4D, nD)
    // - no basis vectors (e₁, e₂, e₃... eₙ)
    // - no conformal embedding (e₀, e∞)
    // - no dimensional storage explosion (2ⁿ components)
    // - blade arithmetic handles infinite dimensional projections

    println!("sphere in any dimension via blade arithmetic, not coordinate storage");
}

#[test]
fn it_represents_lines_as_circle_limits() {
    // traditional CGA: line = P₁ ∧ P₂ ∧ e∞
    // requires special "point at infinity" e∞ and conformal embedding
    // must handle e∞ as special basis vector with complex algebraic rules
    //
    // geonum: infinity emerges naturally from angle arithmetic
    // no special symbols, no conformal embedding - just geometric limits

    // two points define a line
    let p1 = Geonum::new_from_cartesian(0.0, 0.0);
    let p2 = Geonum::new_from_cartesian(4.0, 0.0); // horizontal line for clarity

    // IMPROVEMENT 1: line direction via simple subtraction
    // traditional CGA: must compute P₁ ∧ P₂ in conformal space
    // geonum: direct angle arithmetic
    let line_direction = p2 - p1;
    assert_eq!(line_direction.angle.rem(), 0.0, "horizontal line");

    // IMPROVEMENT 2: perpendicular via π/2 rotation
    // traditional CGA: compute dual in conformal space using I₅
    // geonum: just rotate by π/2
    let perpendicular = line_direction.rotate(Angle::new(1.0, 2.0));
    assert!(perpendicular.angle.is_vector(), "perpendicular is grade 1");

    // IMPROVEMENT 3: demonstrate line as limit of circles
    // traditional CGA: algebraically force e∞ into wedge product
    // geonum: geometric sequence showing infinity emerges naturally

    let midpoint = (p1 + p2) * Geonum::new(0.5, 0.0, 1.0);
    let half_chord = (p2 - p1).mag / 2.0;

    // sequence of circles with increasing radius
    let radii = [1.5, 10.0, 100.0, 1000.0, 10000.0];

    for &radius in &radii {
        // for circle through p1 and p2 with given radius:
        // center lies on perpendicular bisector at specific distance

        if radius > half_chord {
            // center distance from chord using geonum operations
            let center_distance = ((radius * radius) - (half_chord * half_chord)).sqrt();

            // center position using geonum arithmetic - no coordinate math
            let center_offset = perpendicular.normalize() * Geonum::new(center_distance, 0.0, 1.0);
            let center = midpoint + center_offset;

            // test both points on circle via geonum distance
            let d1 = (p1 - center).mag;
            let d2 = (p2 - center).mag;

            assert!((d1 - radius).abs() < EPSILON * radius, "p1 on circle");
            assert!((d2 - radius).abs() < EPSILON * radius, "p2 on circle");

            // KEY INSIGHT: as radius → ∞, curvature → 0
            let curvature = 1.0 / radius;
            println!("radius {radius:8.0}: curvature = {curvature:.6}");
        }
    }

    // IMPROVEMENT 4: infinite radius handled naturally
    // traditional CGA: must algebraically manipulate e∞ in formulas
    // geonum: curvature approaches zero - no special cases

    // demonstrate that line has zero curvature (infinite radius)
    let line_curvature = 0.0; // by definition for straight line
    assert_eq!(line_curvature, 0.0, "line has zero curvature");

    // IMPROVEMENT 5: point-line relationships without e∞
    // traditional CGA: test (P ∧ L)·e∞ = 0 for point on line
    // geonum: use wedge to test collinearity

    let point_on_line = Geonum::new_from_cartesian(3.0, 0.0); // between p1 and p2
    let point_off_line = Geonum::new_from_cartesian(2.0, 1.0); // above the line

    // test collinearity via wedge product
    let v1 = point_on_line - p1;
    let v2 = p2 - p1;
    let wedge_on = v1.wedge(&v2);
    assert!(wedge_on.mag < EPSILON, "collinear points have zero wedge");

    let v3 = point_off_line - p1;
    let wedge_off = v3.wedge(&v2);
    assert!(
        wedge_off.mag > EPSILON,
        "non-collinear points have non-zero wedge"
    );

    // COMPLEXITY COMPARISON:
    // traditional CGA: special e∞ handling, conformal embedding, 5D operations
    // geonum: direct geometric operations, no special infinity symbol

    // geonum eliminates:
    // - e∞ as special basis vector requiring unique algebraic rules
    // - conformal embedding p → e₀ + p + ½p²e∞
    // - special case handling for "point at infinity"
    // - distinction between finite and infinite geometric objects

    println!("line via geometric limits, not conformal infinity e∞");
}

#[test]
fn it_represents_planes_as_sphere_limits() {
    // traditional CGA: plane = P₁ ∧ P₂ ∧ P₃ ∧ e∞
    // requires conformal embedding with special infinity point e∞
    //
    // geonum reveals deeper truth: plane is sphere with center at infinity
    // and infinity is just opposite rotation (dual)

    // spheres as single geonums - length encodes radius
    let sphere_r1 = Geonum::new(1.0, 0.0, 1.0); // radius 1
    let sphere_r10 = Geonum::new(10.0, 0.0, 1.0); // radius 10
    let sphere_r100 = Geonum::new(100.0, 0.0, 1.0); // radius 100
    let sphere_r1000 = Geonum::new(1000.0, 0.0, 1.0); // radius 1000

    // curvature naturally emerges as 1/radius
    assert_eq!(1.0 / sphere_r1.mag, 1.0, "curvature = 1/radius");
    assert_eq!(1.0 / sphere_r10.mag, 0.1, "decreasing curvature");
    assert_eq!(1.0 / sphere_r100.mag, 0.01, "approaching plane");
    assert_eq!(1.0 / sphere_r1000.mag, 0.001, "nearly flat");

    // test limit behavior: as sphere radius → ∞, curvature → 0
    let test_radii = [10.0, 100.0, 1000.0, 10000.0, 100000.0];
    let mut last_curvature = 1.0;

    for radius in test_radii {
        let curvature = 1.0 / radius;

        // curvature decreases monotonically
        assert!(
            curvature < last_curvature,
            "curvature decreases with radius"
        );
        last_curvature = curvature;
    }

    // at the limit: plane has zero curvature
    // this is the geometric definition of a plane
    assert!(last_curvature <= 0.00001, "approaches zero curvature");

    // COMPARISON:
    // traditional CGA: plane requires special infinity point e∞
    // geonum: plane emerges naturally as limit of spheres
    //
    // storage: O(2^5) for conformal 4-vector vs O(1) for [length, angle]
}

#[test]
fn it_computes_distance_between_points() {
    // traditional CGA embeds points in 5D conformal space:
    // P = x²e₊ + x·e + e₋ + e₀
    // requires 32 components (2⁵) per point
    // distance formula: d = sqrt(-2P₁·P₂) using conformal inner product
    //
    // geonum: points are [length, angle], distance is |P₂ - P₁|
    // O(1) storage, O(1) computation

    // create points directly as geometric numbers
    let origin = Geonum::new(0.0, 0.0, 1.0);
    let point_at_5 = Geonum::new(5.0, 0.0, 1.0); // 5 units along 0°
    let point_at_5_rotated = Geonum::new(5.0, 1.0, 2.0); // 5 units along 90°

    // distance between origin and point is just the length
    assert_eq!((point_at_5 - origin).mag, 5.0, "radial distance");

    // distance between two points at same radius but different angles
    let arc_distance = (point_at_5_rotated - point_at_5).mag;
    // arc length for 90° on radius 5: 5 * 2 * sin(π/4) ≈ 7.07
    assert!((arc_distance - 5.0 * 2.0_f64.sqrt()).abs() < EPSILON);

    // KEY INSIGHT: conformal points at infinity
    // traditional CGA: special null vector e∞ = e₊ + e₋
    // geonum: infinity is normalize().dual()

    let finite_point = Geonum::new(10.0, 1.0, 3.0);
    let infinity_point = finite_point.normalize().dual();

    // prove infinity point has unit length (not actually infinite)
    assert_eq!(
        infinity_point.mag, 1.0,
        "infinity has finite representation"
    );

    // prove distance to infinity follows from dual operation
    let to_infinity = infinity_point - finite_point.normalize();
    // distance is determined by the π rotation from dual
    assert_eq!(to_infinity.angle.grade(), 2, "π rotation to infinity");

    // COMPLEXITY COMPARISON:
    // traditional CGA distance between P₁ and P₂:
    // 1. embed P₁ → 5D: O(5) storage
    // 2. embed P₂ → 5D: O(5) storage
    // 3. compute P₁·P₂: O(32) operations (2⁵ multivector components)
    // 4. apply formula: sqrt(-2P₁·P₂)
    //
    // geonum distance:
    // 1. P₁: O(1) storage [length, angle]
    // 2. P₂: O(1) storage [length, angle]
    // 3. P₂ - P₁: O(1) operation
    // 4. result.mag: O(1) access

    // prove geonum handles "null vectors" (zero-length at any angle)
    let null_at_0 = Geonum::new(0.0, 0.0, 1.0);
    let null_at_pi = Geonum::new(0.0, 1.0, 1.0);

    // both are null (zero length) but at different angles
    assert_eq!(null_at_0.mag, 0.0);
    assert_eq!(null_at_pi.mag, 0.0);
    assert_ne!(
        null_at_0.angle, null_at_pi.angle,
        "null vectors can have different angles"
    );

    // traditional CGA would need special null-vector handling
    // geonum: null is just length=0, works with standard operations
    let from_null = point_at_5 - null_at_0;
    assert_eq!(from_null.mag, 5.0, "distance from null point");
}

#[test]
fn it_computes_angle_between_lines() {
    // traditional CGA represents lines as 6D Plücker bivectors:
    // L = p∧d + m*e∞ where p=point, d=direction, m=moment
    // requires 6 coordinates: (l₀₁, l₀₂, l₀₃, l₂₃, l₃₁, l₁₂)
    // angle formula: cos(θ) = (L₁·L₂)/(|L₁||L₂|)
    //
    // geonum: line = bivector with single angle encoding both position and direction
    // angle between lines = angle difference of their bivector representations

    // create lines as bivectors (grade 2 objects)
    // line through origin at angle θ
    let line_horizontal = Geonum::new_with_blade(1.0, 2, 0.0, 1.0); // blade 2, angle 0
    let line_vertical = Geonum::new_with_blade(1.0, 2, 0.0, 1.0).rotate(Angle::new(1.0, 2.0)); // rotate by π/2
    let line_diagonal = Geonum::new_with_blade(1.0, 2, 0.0, 1.0).rotate(Angle::new(1.0, 4.0)); // rotate by π/4

    // KEY INSIGHT: rotating by π/2 changes grade
    // blade 2 + π/2 rotation = blade 3 (crossing grade boundary)
    // this demonstrates how geometric transformations naturally change algebraic structure

    // prove bivector grade
    assert_eq!(line_horizontal.angle.grade(), 2, "lines are bivectors");
    assert_eq!(
        line_vertical.angle.grade(),
        3,
        "π/2 rotation changes grade 2→3"
    );
    assert_eq!(line_diagonal.angle.grade(), 2, "π/4 rotation keeps grade 2");

    // angle between lines is just angle difference
    let angle_diff_h_to_v = line_vertical.angle - line_horizontal.angle;
    let expected_perpendicular = Angle::new(1.0, 2.0); // π/2 difference
    assert_eq!(
        angle_diff_h_to_v, expected_perpendicular,
        "perpendicular lines differ by π/2"
    );

    let angle_diff_h_to_d = line_diagonal.angle - line_horizontal.angle;
    let expected_diagonal = Angle::new(1.0, 4.0); // π/4 difference
    assert_eq!(
        angle_diff_h_to_d, expected_diagonal,
        "diagonal line at π/4 from horizontal"
    );

    // parallel lines have same angle
    let line_parallel = Geonum::new_with_blade(2.0, 2, 0.0, 1.0); // another horizontal
    let angle_diff_parallel = line_parallel.angle - line_horizontal.angle;
    let expected_parallel = Angle::new(0.0, 1.0); // 0 difference
    assert_eq!(
        angle_diff_parallel, expected_parallel,
        "parallel lines have same angle"
    );

    // KEY INSIGHT: lines at infinity
    // traditional CGA: special handling for lines meeting at infinity
    // geonum: lines at infinity are just dual of finite lines

    let finite_line = Geonum::new_with_blade(1.0, 2, 1.0, 3.0); // π/3 angle
    let infinity_line = finite_line.dual(); // adds π rotation

    // prove infinity line is grade 0 (scalar) due to dual operation
    // dual adds 2 to blade, so blade 2 → blade 4, and grade = blade % 4 = 0
    assert_eq!(infinity_line.angle.grade(), 0, "dual maps bivector→scalar");

    // angle to infinity line includes the π rotation from dual
    // dual() adds π (2 blades) to the angle
    let finite_angle = finite_line.angle; // blade 2
    let dual_adds = Angle::new(1.0, 1.0); // π added by dual()
    let expected_infinity_angle = finite_angle + dual_adds; // blade 4

    assert_eq!(
        infinity_line.angle,
        expected_infinity_angle,
        "dual() adds π: blade {} + 2 = blade {}",
        finite_angle.blade(),
        expected_infinity_angle.blade()
    );

    // COMPLEXITY COMPARISON:
    // traditional CGA line operations:
    // 1. store 6 Plücker coordinates per line
    // 2. normalize: sqrt(l₀₁² + l₀₂² + l₀₃² + l₂₃² + l₃₁² + l₁₂²)
    // 3. inner product: Σᵢⱼ(L₁ᵢⱼ * L₂ᵢⱼ) with metric tensor
    // 4. arccos for angle extraction
    //
    // geonum line operations:
    // 1. store 1 bivector [length, angle] where blade=2
    // 2. angle difference: angle₂ - angle₁
    // 3. no normalization, no inner product, no arccos

    // prove meet of parallel lines gives infinity result
    let parallel1 = Geonum::new_with_blade(1.0, 2, 1.0, 6.0); // π/6
    let parallel2 = Geonum::new_with_blade(1.0, 2, 1.0, 6.0); // same angle
    let parallel_meet = parallel1.meet(&parallel2);

    // parallel lines meet at infinity - but what does this mean geometrically?
    // traditional math says they meet at a "point at infinity" (fiction to avoid saying they don't meet)
    // geonum reveals the truth: parallel lines don't meet at a point, they span a plane together
    // the trivector (grade 3) represents the oriented volume containing this plane
    // so grade 3 isn't arbitrary - it's saying parallel lines define a higher-dimensional object
    // (the 3D space containing their plane) not a lower-dimensional point
    assert_eq!(
        parallel_meet.angle.grade(),
        3,
        "parallel lines create volume, not point"
    );

    // intersection of non-parallel lines
    let intersecting1 = Geonum::new_with_blade(1.0, 2, 0.0, 1.0); // horizontal
    let intersecting2 = Geonum::new_with_blade(1.0, 2, 1.0, 2.0); // vertical
    let intersection = intersecting1.meet(&intersecting2);

    // non-parallel lines intersect at a point (grade 0 = scalar)
    // this is the familiar case where two lines cross at a specific location
    // grade 0 represents a point - the lowest dimensional object
    // contrast with parallel lines that produce grade 3 (volume)
    assert_eq!(
        intersection.angle.grade(),
        0,
        "intersecting lines meet at point"
    );
    assert!(
        intersection.mag > 0.0,
        "finite intersection has non-zero length"
    );
}

#[test]
fn it_tests_point_on_circle() {
    // traditional CGA: point-on-circle test via P·C = 0
    // requires: embed point as null vector, embed circle as bivector,
    // compute conformal inner product with metric tensor
    //
    // geonum: |(point - center)| ≈ radius

    // test basic point-on-circle at various grades
    let test_radius = 2.5;

    // grade 0 circle (scalar center)
    let circle_0 = Geonum::scalar(5.0);
    let point_0 = circle_0 + Geonum::scalar(test_radius);
    assert!(
        (point_0 - circle_0).mag - test_radius < EPSILON,
        "grade 0: point on circle via distance"
    );

    // grade 1 circle (vector center)
    let circle_1 = Geonum::new_with_blade(5.0, 1, 1.0, 4.0);
    let direction_1 = Geonum::new_with_blade(1.0, 1, 1.0, 2.0); // π/2 direction
    let point_1 = circle_1 + direction_1 * Geonum::scalar(test_radius);
    assert!(
        (point_1 - circle_1).mag - test_radius < EPSILON,
        "grade 1: point on circle via distance"
    );

    // grade 2 circle (bivector center) - most common for 2D circles
    let circle_2 = Geonum::new_from_cartesian(3.0, 4.0);
    let angles = [0.0, 1.0, 2.0, 3.0]; // 0, π/2, π, 3π/2
    for i in angles {
        let angle = Angle::new(i, 2.0);
        let direction = Geonum::new_with_angle(1.0, angle);
        let point = circle_2 + direction * Geonum::scalar(test_radius);
        assert!(
            (point - circle_2).mag - test_radius < EPSILON,
            "grade 2: point at angle {:.2}π on circle",
            i / 2.0
        );
    }

    // grade 3 circle (trivector center)
    let circle_3 = Geonum::new_with_blade(5.0, 3, 0.0, 1.0);
    let direction_3 = Geonum::new_with_blade(1.0, 3, 1.0, 1.0); // π direction
    let point_3 = circle_3 + direction_3 * Geonum::scalar(test_radius);
    assert!(
        (point_3 - circle_3).mag - test_radius < EPSILON,
        "grade 3: point on circle via distance"
    );

    // test points NOT on circle
    let center = Geonum::new_from_cartesian(2.0, 3.0);
    let radius = 4.0;

    let inside = center + Geonum::new_from_cartesian(1.0, 1.0); // < radius away
    let outside = center + Geonum::new_from_cartesian(5.0, 0.0); // > radius away
    let on_circle = center + Geonum::new_from_cartesian(4.0, 0.0); // = radius away

    assert!((inside - center).mag < radius, "point inside circle");
    assert!((outside - center).mag > radius, "point outside circle");
    assert!(
        (on_circle - center).mag - radius < EPSILON,
        "point on circle"
    );

    // demonstrate infinity handling via dual
    let finite_point = Geonum::new_with_blade(10.0, 2, 1.0, 3.0);
    let _infinity_point = finite_point.dual();

    // for circles at infinity, dual inverts the scale
    let large_circle_center = Geonum::new_from_cartesian(0.0, 0.0);
    let large_radius = 1000.0;
    let large_circle_point = large_circle_center + Geonum::new_from_cartesian(large_radius, 0.0);

    // test point on large circle
    assert!(
        (large_circle_point - large_circle_center).mag - large_radius < EPSILON,
        "point on large circle"
    );

    // dual of large circle point
    let dual_point = large_circle_point.dual();
    // geonum dual maps through involutive pairs: 0↔2, 1↔3
    // large_circle_point at (1000, 0) has angle 0, so grade 0 (scalar)
    // dual of scalar (grade 0) is bivector (grade 2)
    assert_eq!(
        large_circle_point.angle.grade(),
        0,
        "point at (1000,0) is scalar"
    );
    assert_eq!(dual_point.angle.grade(), 2, "dual of scalar is bivector");
    assert_eq!(
        dual_point.mag, large_circle_point.mag,
        "dual preserves length"
    );

    // COMPLEXITY COMPARISON:
    // traditional CGA point-on-circle:
    // 1. embed point: P = e₀ + p + ½|p|²e∞ (5 components)
    // 2. embed circle: C = c₁∧c₂∧c₃ (10 bivector components)
    // 3. compute P·C with conformal metric (~50 operations)
    // 4. test if result ≈ 0
    //
    // geonum point-on-circle:
    // 1. point and center: [length, angle] each
    // 2. test: |(point - center)| ≈ radius (O(1))
    // 3. works at any grade without modification

    println!("point-on-circle via O(1) distance at all grades");
}

#[test]
fn it_tests_point_on_sphere() {
    // traditional CGA represents spheres in 5D conformal space:
    // S = c - ½r²e∞ where c is center, r is radius
    // requires 5 components with special infinity point e∞
    // point-on-sphere test: P·S = 0 using conformal inner product
    //
    // geonum: sphere = circle at higher grade
    // "sphere" is just what we call a 2D surface in 3D, but it's really
    // the same as a circle - points at fixed distance from center

    // KEY INSIGHT: spheres are just circles at different grades
    // circle in 2D: center at grade 2 (bivector)
    // sphere in 3D: center at grade 3 (trivector)
    // hypersphere in 4D: center at grade 0 (scalar, wraps around)

    let sphere_2d = Geonum::new_with_blade(1.0, 2, 0.0, 1.0); // "circle"
    let sphere_3d = Geonum::new_with_blade(1.0, 3, 0.0, 1.0); // "sphere"
    let sphere_4d = Geonum::new_with_blade(1.0, 0, 0.0, 1.0); // "hypersphere"

    assert_eq!(sphere_2d.angle.grade(), 2, "2D sphere (circle) is grade 2");
    assert_eq!(sphere_3d.angle.grade(), 3, "3D sphere is grade 3");
    assert_eq!(sphere_4d.angle.grade(), 0, "4D sphere wraps to grade 0");

    // test actual points on sphere surfaces
    let radius = 5.0;

    // create points at various angles on unit sphere, then scale by radius
    let test_points = vec![
        // points ON the sphere (distance = radius)
        (Geonum::new(radius, 0.0, 1.0), true), // along x-axis
        (Geonum::new(radius, 1.0, 2.0), true), // along y-axis (π/2)
        (Geonum::new(radius, 1.0, 1.0), true), // diagonal (π)
        (Geonum::new(radius, 3.0, 2.0), true), // 3π/2
        // points NOT on sphere
        (Geonum::new(radius * 0.5, 0.0, 1.0), false), // too close
        (Geonum::new(radius * 2.0, 1.0, 2.0), false), // too far
        (Geonum::scalar(0.0), false),                 // at center
    ];

    // test with sphere at origin (grade 0 center)
    let origin = Geonum::scalar(0.0);

    for (point, expected_on_sphere) in &test_points {
        let distance = (*point - origin).mag;
        let on_sphere = (distance - radius).abs() < EPSILON;
        assert_eq!(
            on_sphere,
            *expected_on_sphere,
            "point at distance {} is {} sphere of radius {}",
            distance,
            if on_sphere { "on" } else { "not on" },
            radius
        );
    }

    // prove same test works for higher grade centers
    // translate sphere to different location (grade 1 center)
    let translated_center = Geonum::new_with_blade(3.0, 1, 1.0, 4.0);

    for (point, expected_on_sphere) in test_points {
        // translate point by same amount
        let translated_point = point + translated_center;
        let distance = (translated_point - translated_center).mag;
        let on_sphere = (distance - radius).abs() < EPSILON;
        assert_eq!(
            on_sphere, expected_on_sphere,
            "translation preserves sphere membership"
        );
    }

    // COMPLEXITY COMPARISON:
    // traditional CGA sphere operations:
    // 1. embed sphere → 5D: S = c - ½r²e∞
    // 2. special handling for e∞ (point at infinity)
    // 3. embed point → 5D null vector
    // 4. conformal inner product with metric tensor
    //
    // geonum sphere operations:
    // 1. sphere = center (at specific grade) + radius
    // 2. point on sphere: (point - center).mag = radius
    // 3. same formula for circle, sphere, hypersphere

    // prove spheres at infinity work naturally
    let finite_sphere = Geonum::new_with_blade(1.0, 3, 1.0, 6.0); // grade 3
    let infinity_sphere = finite_sphere.dual(); // adds π rotation

    // dual changes grade: 3 → 1 (trivector → vector)
    assert_eq!(
        infinity_sphere.angle.grade(),
        1,
        "infinity sphere at grade 1"
    );

    // the "sphere at infinity" is just the dual of a regular sphere
    // no special machinery needed - just π rotation

    // prove geonum unifies circles and spheres
    // what math calls different objects are just different grades:
    let shapes = vec![
        ("point", 0, 0),    // blade 0 → grade 0
        ("circle", 2, 2),   // blade 2 → grade 2
        ("sphere", 3, 3),   // blade 3 → grade 3
        ("4-sphere", 4, 0), // blade 4 → grade 0 (wraps)
    ];

    for (name, blade, expected_grade) in shapes {
        let shape = Geonum::new_with_blade(1.0, blade, 0.0, 1.0);
        assert_eq!(
            shape.angle.grade(),
            expected_grade,
            "{name} with blade {blade} has grade {expected_grade}"
        );
    }

    // all these "different" shapes use the same point-on-surface test:
    // (point - center).mag = radius

    // no distinction between circle and sphere - just grade difference
    // no special infinity handling - just dual operation
    // no 5D embedding - just [length, angle] at the required grade
}

#[test]
fn it_finds_circle_circle_intersection() {
    // traditional CGA circle intersection:
    // 1. represent circles as bivectors in 5D: C = p∧q∧e∞ + ...
    // 2. compute meet: I = C₁ ∨ C₂ using outer product
    // 3. extract intersection points from resulting trivector
    // requires ~20 bivector components and complex extraction
    //
    // geonum: circles are bivectors, intersection is their meet

    // create circles as bivectors (grade 2)
    let circle1 = Geonum::new_with_blade(3.0, 2, 0.0, 1.0); // radius 3 bivector
    let circle2 = Geonum::new_with_blade(2.5, 2, 1.0, 4.0); // radius 2.5, rotated π/4

    // KEY INSIGHT: circle intersection grade tells intersection type
    let intersection = circle1.meet(&circle2);

    // intersection properties encode geometry:
    // - grade 3 with length > 0: typical two-point intersection
    // - grade 3 with length ≈ 0: tangent circles (one point)
    // - grade 1: circles at opposite angles
    // the grade and length together encode the intersection type

    assert_eq!(
        intersection.angle.grade(),
        3,
        "typical intersection has grade 3"
    );
    assert!(
        intersection.mag > 0.0,
        "two-point intersection has non-zero length"
    );

    // prove tangent circles produce different grade
    // create circles that just touch
    let tangent1 = Geonum::new_with_blade(2.0, 2, 0.0, 1.0);
    let tangent2 = Geonum::new_with_blade(3.0, 2, 0.0, 1.0); // same angle = parallel

    let tangent_meet = tangent1.meet(&tangent2);
    // parallel bivectors (same angle) meet with zero length
    assert_eq!(
        tangent_meet.angle.grade(),
        3,
        "tangent circles meet at grade 3"
    );
    assert!(
        tangent_meet.mag < EPSILON,
        "tangent circles have zero-length meet"
    );

    // prove non-intersecting circles
    let far1 = Geonum::new_with_blade(1.0, 2, 0.0, 1.0);
    let far2 = Geonum::new_with_blade(1.0, 2, 1.0, 1.0); // π rotation away

    let no_meet = far1.meet(&far2);
    // opposite bivectors (π apart) produce grade 1
    assert_eq!(no_meet.angle.grade(), 1, "opposite circles meet at grade 1");
    assert!(
        no_meet.mag < EPSILON,
        "non-intersecting circles have zero meet"
    );

    // COMPLEXITY COMPARISON:
    // traditional CGA circle intersection:
    // 1. embed circles → 5D: ~10 components each
    // 2. compute C₁ ∨ C₂: ~100 operations on bivector components
    // 3. extract points: complex trivector → point conversion
    // 4. handle special cases: tangent, parallel, identical
    //
    // geonum circle intersection:
    // 1. circles: 1 bivector each [length, angle]
    // 2. meet operation: dual-wedge-dual
    // 3. grade tells intersection type directly

    // prove circles at different grades can intersect
    let circle_2d = Geonum::new_with_blade(2.0, 2, 0.0, 1.0); // grade 2
    let circle_3d = Geonum::new_with_blade(2.0, 3, 0.0, 1.0); // grade 3

    let cross_grade_meet = circle_2d.meet(&circle_3d);
    // different grades produce intersection at intermediate grade
    assert!(
        cross_grade_meet.mag > 0.0,
        "cross-grade intersection exists"
    );

    // prove infinity handling
    let finite_circle = Geonum::new_with_blade(3.0, 2, 1.0, 6.0);
    let infinity_circle = finite_circle.dual(); // π rotation to infinity

    let infinity_meet = finite_circle.meet(&infinity_circle);
    // finite circle meets its infinity version
    assert_eq!(
        infinity_meet.angle.grade(),
        1,
        "circle meets its infinity dual at grade 1"
    );

    // the meet grade encodes all intersection information
    // no need for coordinate extraction or special case handling
}

#[test]
fn it_finds_sphere_sphere_intersection() {
    // traditional CGA: sphere-sphere intersection via S₁ ∨ S₂
    // requires:
    // 1. represent spheres as (n+1)-vectors in conformal space
    // 2. compute meet operation on ~32 components in 5D
    // 3. extract intersection circle/point from resulting multivector
    // 4. handle degenerate cases with special formulas
    //
    // geonum: spheres encoded in geometric numbers, meet tells intersection

    // in geonum, sphere data lives in a single geometric number
    // we use grade 2 (bivector) - the angle encodes position, length encodes size
    let sphere1 = Geonum::new_with_blade(5.0, 2, 0.0, 1.0); // sphere at angle 0
    let sphere2 = Geonum::new_with_blade(4.0, 2, 1.0, 6.0); // sphere at angle π/6

    // KEY SIMPLIFICATION: meet operation encodes intersection type
    let intersection = sphere1.meet(&sphere2);

    // the grade and length of the meet result tells us the geometric relationship
    // traditional CGA: sphere intersection produces a circle (lower dimension)
    // geonum: sphere meet encodes their spatial relationship via grade
    //
    // bivector meet bivector → trivector (grade 3)
    // grade 3 means the spheres are in general position (not special alignment)
    // the meet computes the 3D space containing both spheres
    assert_eq!(
        intersection.angle.grade(),
        3,
        "general position spheres produce grade 3"
    );

    // test parallel spheres (same angle = same position in angle space)
    let parallel1 = Geonum::new_with_blade(3.0, 2, 1.0, 4.0); // π/4 angle
    let parallel2 = Geonum::new_with_blade(2.0, 2, 1.0, 4.0); // same angle

    let parallel_meet = parallel1.meet(&parallel2);
    // spheres at the same angle are like concentric spheres
    // their meet has zero length indicating this special alignment
    assert!(
        parallel_meet.mag < EPSILON,
        "same-angle spheres have zero meet"
    );

    // test opposite spheres (π apart in angle)
    let opposite1 = Geonum::new_with_blade(2.0, 2, 0.0, 1.0); // 0 angle
    let opposite2 = Geonum::new_with_blade(2.0, 2, 1.0, 1.0); // π angle

    let opposite_meet = opposite1.meet(&opposite2);
    // spheres π apart in angle space have a special relationship
    // their meet produces grade 1 (vector) instead of the usual grade 3
    // this grade change signals their orthogonal configuration
    assert_eq!(
        opposite_meet.angle.grade(),
        1,
        "π-separated spheres produce grade 1"
    );

    // demonstrate that sphere "intersection" is encoded in the meet operation
    // without any coordinate extraction or distance calculations

    // create spheres at different grades to show cross-grade intersection
    let sphere_grade2 = Geonum::new_with_blade(3.0, 2, 1.0, 3.0); // bivector
    let sphere_grade3 = Geonum::new_with_blade(4.0, 3, 1.0, 4.0); // trivector

    let cross_grade_meet = sphere_grade2.meet(&sphere_grade3);
    // different grades produce intermediate grade intersection
    assert!(
        cross_grade_meet.mag > 0.0 || cross_grade_meet.mag < EPSILON,
        "cross-grade meet has definite result"
    );

    // COMPLEXITY COMPARISON:
    // traditional CGA sphere-sphere intersection:
    // 1. embed spheres in 5D: S = e₀ + c + ½(|c|²-r²)e∞ (~32 components each)
    // 2. compute S₁ ∨ S₂: ~1000 operations on bivector components
    // 3. extract intersection circle from result (complex decomposition)
    // 4. handle special cases: tangent, concentric, disjoint
    //
    // geonum sphere-sphere intersection:
    // 1. spheres: single bivector [length, angle, blade=2]
    // 2. meet operation: O(1) dual-wedge-dual
    // 3. grade encodes geometric relationship:
    //    - grade 3: general position (3D volume containing both)
    //    - grade 1: orthogonal spheres (line intersection)
    //    - length 0: parallel/concentric spheres
    // 4. no coordinate extraction or special cases needed

    // demonstrate sphere at infinity via dual
    let finite_sphere = Geonum::new_with_blade(3.0, 2, 1.0, 6.0);
    let infinity_sphere = finite_sphere.dual(); // dual takes to infinity

    // dual of bivector (grade 2) produces scalar (grade 0)
    // the involutive pair: scalar ↔ bivector (0↔2)
    assert_eq!(infinity_sphere.angle.grade(), 0, "bivector dual is scalar");

    // the meet of a sphere with its infinity dual
    let infinity_meet = finite_sphere.meet(&infinity_sphere);
    assert_eq!(
        infinity_meet.angle.grade(),
        1,
        "sphere meets its dual at grade 1"
    );

    println!("sphere intersection via O(1) meet operation, not O(2⁵) conformal algebra");

    // geonum eliminates:
    // - 5D conformal embedding with e₀, e∞
    // - null vectors and null cone
    // - 32-component multivector storage
    // - complex extraction formulas
    //
    // spheres are just bivectors, intersection is just meet
    // grade tells everything about the geometric relationship
}

#[test]
fn it_finds_line_circle_intersection() {
    // traditional CGA: line-circle intersection via L ∨ C
    // requires:
    // 1. represent line as bivector: L = p∧q∧e∞
    // 2. represent circle as bivector: C = c₁∧c₂∧c₃
    // 3. compute meet: L ∨ C in conformal space
    // 4. extract intersection points from result
    //
    // geonum: use meet operation to encode intersection type

    // circle encoded as bivector (grade 2)
    let circle = Geonum::new_with_blade(5.0, 2, 0.0, 1.0); // radius 5 circle

    // line encoded as bivector (grade 2)
    let line_horizontal = Geonum::new_with_blade(3.0, 2, 1.0, 4.0); // π/4 angle

    // KEY SIMPLIFICATION: meet operation encodes intersection
    let intersection = line_horizontal.meet(&circle);

    // the grade and length of meet result tells us the geometric relationship
    // bivector meet bivector → grade 1 (vector)
    //
    // why grade 1? the meet computes the common subspace between line and circle
    // for 2D objects (bivectors), their intersection is 1D (vector)
    // this vector points from origin to the intersection region
    // length > 0 means they actually intersect; length ≈ 0 means they touch

    assert_eq!(
        intersection.angle.grade(),
        1,
        "2D meet 2D → 1D intersection"
    );
    assert!(
        intersection.mag > 0.0,
        "non-zero length = actual intersection"
    );

    // demonstrate tangent line
    // line at same angle as circle = tangent configuration
    let line_tangent = Geonum::new_with_blade(5.0, 2, 0.0, 1.0); // same radius as circle
    let tangent_meet = line_tangent.meet(&circle);

    // parallel objects (same angle) have zero-length meet
    // why grade 3? when bivectors are parallel (same angle), their meet
    // produces the 3D volume containing both - like how parallel planes
    // define a slab of 3D space between them
    // length → 0 because they're infinitesimally close (tangent)
    assert_eq!(
        tangent_meet.angle.grade(),
        3,
        "parallel bivectors → trivector"
    );
    assert!(tangent_meet.mag < EPSILON, "zero length = tangent touch");

    // demonstrate missing line
    // line at opposite angle (π apart) from circle
    let line_miss = Geonum::new_with_blade(7.0, 2, 1.0, 1.0); // π angle from circle
    let miss_meet = line_miss.meet(&circle);

    // opposite angles (π apart) create special geometry
    // these objects are maximally separated in angle space
    // their meet produces grade 1 (vector) with near-zero length
    // this signals they're orthogonal in the geometric sense
    assert_eq!(miss_meet.angle.grade(), 1, "opposite angles → vector meet");
    assert!(miss_meet.mag < EPSILON, "orthogonal objects → zero meet");

    // demonstrate line through center
    // smaller radius line inside the circle
    let line_center = Geonum::new_with_blade(2.0, 2, 0.0, 1.0); // radius 2 < 5
    let center_meet = line_center.meet(&circle);

    // concentric objects (same angle, different scale)
    // these are like nested circles - one inside the other
    // same angle → parallel → grade 3 (trivector)
    // the 3D volume contains both concentric circles
    assert_eq!(
        center_meet.angle.grade(),
        3,
        "concentric → parallel → grade 3"
    );
    assert!(center_meet.mag < EPSILON, "concentric = degenerate meet");

    // COMPLEXITY COMPARISON:
    // traditional CGA line-circle intersection:
    // 1. embed line in 5D: L = p∧d∧e∞ (~10 bivector components)
    // 2. embed circle in 5D: C = c₁∧c₂∧c₃ (~10 bivector components)
    // 3. compute L ∨ C: ~100 operations on components
    // 4. extract points: solve quadratic from result
    //
    // geonum line-circle intersection:
    // 1. line and circle: bivectors [length, angle, blade=2]
    // 2. meet operation: O(1) dual-wedge-dual
    // 3. grade tells intersection type directly
    // 4. no quadratic solving needed

    // demonstrate that lines and circles are both grade 2 objects
    assert_eq!(circle.angle.grade(), 2, "circle is bivector");
    assert_eq!(line_horizontal.angle.grade(), 2, "line is bivector");

    // KEY PATTERN: bivector meet bivector produces:
    // - grade 1 (vector): general intersection or orthogonal configuration
    // - grade 3 (trivector): parallel/tangent/concentric configuration
    // length encodes intersection strength (0 = tangent/miss, >0 = intersect)
    //
    // this replaces coordinate-based quadratic solving with direct geometric relationships

    println!("line-circle via O(1) meet, not quadratic solving");
}

#[test]
fn it_finds_line_sphere_intersection() {
    // traditional CGA: line-sphere intersection via L ∨ S
    // requires:
    // 1. embed line in 5D: L = p∧d∧e∞ (bivector, ~10 components)
    // 2. embed sphere in 5D: S = s - ½r²e∞ (vector, 5 components)
    // 3. compute meet L ∨ S in conformal space (~50 operations)
    // 4. extract intersection points from result (complex formulas)
    //
    // geonum: meet operation directly encodes intersection geometry

    // sphere as trivector (grade 3 = volume)
    let sphere = Geonum::new_with_blade(5.0, 3, 0.0, 1.0); // radius 5 sphere

    // line as vector (grade 1 = 1D curve)
    let line_through = Geonum::new_with_blade(3.0, 1, 1.0, 4.0); // blade 1 + π/4

    // KEY SIMPLIFICATION: meet encodes geometric relationship, not intersection points
    let intersection = line_through.meet(&sphere);

    // WHY GRADE 3? vector (grade 1) meet trivector (grade 3) → trivector (grade 3)
    //
    // GEOMETRIC INSIGHT: geonum's π-rotation dual creates a different incidence structure
    // than traditional GA. instead of computing intersection points (lower dimension),
    // it computes the containing space (preserves or increases dimension)
    //
    // think of it this way: a line through a sphere doesnt reduce to points -
    // together they still span a 3D volume. the meet tells us HOW they relate
    // within that volume, not WHERE they intersect
    assert_eq!(
        intersection.angle.grade(),
        3,
        "1D meet 3D → 3D containing space"
    );
    assert!(
        intersection.mag > 0.0,
        "non-zero length = definite geometric relationship"
    );

    // demonstrate parallel configuration (same angle in geonum space)
    let line_parallel = Geonum::new_with_blade(5.0, 1, 0.0, 1.0); // angle 0, same as sphere
    let parallel_meet = line_parallel.meet(&sphere);

    // PARALLEL OBJECTS: when line and sphere have same angle, theyre "parallel"
    // in geonum's angle space (not parallel in 3D euclidean sense!)
    //
    // KEY: wedge product computes sin(angle_diff). for parallel objects:
    // sin(0) = 0, so wedge gives near-zero. but meet = dual(wedge(dual,dual))
    // applies multiple transformations, so near-zero wedge → small but non-zero meet
    //
    // GRADE STAYS 3: even parallel objects together span 3D volume
    // the small length tells us theyre aligned in angle space
    assert_eq!(parallel_meet.angle.grade(), 3, "parallel → same 3D span");
    println!("parallel_meet length: {}", parallel_meet.mag);
    assert!(
        parallel_meet.mag < 100.0 * EPSILON,
        "parallel → small meet length"
    );

    // demonstrate orthogonal configuration (π/2 angle difference)
    let line_ortho = Geonum::new_with_blade(7.0, 1, 1.0, 2.0); // blade 1, angle π/2
    let ortho_meet = line_ortho.meet(&sphere);

    // ORTHOGONAL SURPRISE: vector meet trivector at π/2 → scalar (grade 0)!
    //
    // WHY GRADE 0? the meet operation is dual(wedge(dual,dual)):
    // 1. dual(vector) → trivector, dual(trivector) → vector (π-rotation dual)
    // 2. wedge(trivector, vector) at π/2 gives maximum sin(π/2) = 1
    // 3. final dual brings us to scalar (grade 0)
    //
    // GEOMETRIC MEANING: the scalar represents a "weighted point" -
    // not an intersection point, but a measure of how strongly these
    // orthogonal objects relate. length = 7×5×1 = 35 (product of lengths × sin(π/2))
    println!(
        "ortho_meet grade: {}, length: {}",
        ortho_meet.angle.grade(),
        ortho_meet.mag
    );
    assert_eq!(
        ortho_meet.angle.grade(),
        0,
        "orthogonal → scalar (weighted point)"
    );
    assert!(
        ortho_meet.mag > 0.0,
        "orthogonal → strong relationship (max sin)"
    );

    // demonstrate "inside" configuration (length comparison in geonum)
    let line_inside = Geonum::new_with_blade(2.0, 1, 0.0, 1.0); // length 2 < sphere's 5
    let inside_meet = line_inside.meet(&sphere);

    // both at angle 0 — parallel in geonum space
    // parallel objects have zero wedge product, so meet magnitude is zero
    // different magnitudes dont create angular separation
    assert!(
        inside_meet.mag.abs() < 1e-10,
        "parallel geonums → zero meet"
    );

    // demonstrate general angle configuration
    let line_angled = Geonum::new_with_blade(3.0, 1, 1.0, 3.0); // π/3 angle
    let angled_meet = line_angled.meet(&sphere);

    // GENERAL CASE: line at π/3, sphere at 0 → angle diff = π/3
    // sin(π/3) = √3/2 ≈ 0.866, so wedge gives moderate value
    //
    // PATTERN EMERGES:
    // - parallel (0°): sin(0) = 0 → minimal meet
    // - angled (60°): sin(π/3) = 0.866 → moderate meet
    // - orthogonal (90°): sin(π/2) = 1 → maximal meet (but different grade!)
    //
    // grade 3 persists because vector + trivector span 3D regardless of angle
    assert_eq!(angled_meet.angle.grade(), 3, "any angle → 3D span");

    // COMPLEXITY COMPARISON:
    // traditional CGA line-sphere:
    // 1. construct L = p∧d∧e∞: ~15 operations for bivector
    // 2. construct S = s - ½r²e∞: ~5 operations for conformal sphere
    // 3. compute L ∨ S: ~50 operations on 5D components
    // 4. solve for points: quadratic formula on extracted components
    //
    // geonum line-sphere:
    // 1. line: vector [length, angle, blade=1]
    // 2. sphere: trivector [length, angle, blade=3]
    // 3. meet: O(1) operation
    // 4. grade/length directly encode intersection type

    // demonstrate different dimensional representations
    assert_eq!(line_through.angle.grade(), 1, "line is 1D (vector)");
    assert_eq!(sphere.angle.grade(), 3, "sphere is 3D (trivector)");

    // KEY PATTERNS REVEALED:
    //
    // 1. GRADE TELLS DIMENSIONAL SPAN, NOT INTERSECTION
    //    vector (1D) meet trivector (3D) → trivector (3D) usually
    //    except orthogonal case → scalar (0D) weighted point
    //
    // 2. LENGTH ENCODES ANGLE RELATIONSHIP VIA sin(θ)
    //    parallel: sin(0) = 0 → near-zero meet
    //    angled: sin(θ) ∈ (0,1) → moderate meet
    //    orthogonal: sin(π/2) = 1 → maximal meet
    //
    // 3. GEONUM'S MEET ≠ TRADITIONAL INTERSECTION
    //    traditional: finds common points (reduces dimension)
    //    geonum: finds containing space (preserves/increases dimension)
    //    this is because π-rotation dual creates different incidence structure
    //
    // this replaces quadratic solving with direct geometric relationships
    // the meet operation encodes everything through grade and length

    // geonum ghosts L ∨ S conformal meet
    // O(1) grade/length vs O(2^n) component operations

    println!("line-sphere via O(1) meet encoding, not quadratic extraction");
}

#[test]
fn it_applies_translation() {
    // traditional CGA: T = 1 - ½te∞ (translator versor)
    // P' = TPT̃
    //
    // geonum: translation is addition in cartesian representation

    // translation vector
    let translation = Geonum::new_from_cartesian(3.0, 4.0);

    // points to translate
    let points = vec![
        Geonum::new_from_cartesian(0.0, 0.0),   // origin
        Geonum::new_from_cartesian(1.0, 0.0),   // unit x
        Geonum::new_from_cartesian(0.0, 1.0),   // unit y
        Geonum::new_from_cartesian(2.0, 2.0),   // diagonal
        Geonum::new_from_cartesian(-1.0, -1.0), // negative quadrant
    ];

    for point in points {
        // KEY INSIGHT: in geonum, + operator IS translation!
        // no need to convert to cartesian - addition already works in cartesian
        let translated = point + translation;

        // verify translation preserves distances between points
        // (rigid motion property)
        let origin = Geonum::new_from_cartesian(0.0, 0.0);
        let origin_translated = origin + translation;

        // distance from origin to point should equal
        // distance from translated origin to translated point
        let original_distance = point.mag; // distance from origin
        let translated_distance = (translated - origin_translated).mag;

        assert!(
            (original_distance - translated_distance).abs() < EPSILON,
            "translation preserves distances (rigid motion)"
        );

        // for non-zero points, check angle relationships
        if point.mag > EPSILON {
            // translation changes individual angles but preserves relative angles
            // this is the geometric signature of translation
            println!(
                "point: length={:.2}, angle={:.2} → translated: length={:.2}, angle={:.2}",
                point.mag,
                point.angle.rem(),
                translated.mag,
                translated.angle.rem()
            );
        }
    }

    // case: translate by zero (identity)
    let zero_translation = Geonum::new_from_cartesian(0.0, 0.0);
    let test_point = Geonum::new_from_cartesian(5.0, 5.0);
    let unchanged = test_point + zero_translation;

    // zero translation is the additive identity
    assert!(
        (unchanged.mag - test_point.mag).abs() < EPSILON,
        "zero translation preserves length"
    );
    assert_eq!(
        unchanged.angle, test_point.angle,
        "zero translation preserves angle"
    );

    // case: chain translations
    let trans1 = Geonum::new_from_cartesian(2.0, 0.0);
    let trans2 = Geonum::new_from_cartesian(0.0, 3.0);

    let start = Geonum::new_from_cartesian(1.0, 1.0);

    // apply translations sequentially
    let after_trans1 = start + trans1;
    let final_sequential = after_trans1 + trans2;

    // apply combined translation
    let combined_trans = trans1 + trans2; // translations add!
    let final_combined = start + combined_trans;

    // both approaches give same result (associativity of +)
    assert!(
        (final_sequential.mag - final_combined.mag).abs() < EPSILON,
        "translation composition is associative"
    );
    assert!(
        (final_sequential.angle.rem() - final_combined.angle.rem()).abs() < EPSILON,
        "(p + t1) + t2 = p + (t1 + t2)"
    );

    // COMPLEXITY COMPARISON:
    // traditional CGA: T = 1 - ½te∞, then P' = TPT̃ (sandwich product)
    // - construct translator versor: O(5) operations
    // - sandwich product: O(32) operations per point
    // - composition: T₁₂ = T₂T₁ (geometric product of versors)
    //
    // geonum: translation = addition
    // - apply: p + t (O(1) cartesian addition)
    // - compose: t1 + t2 (O(1) addition)
    // - no versors, no sandwich products, no e∞

    println!("translation via + operator, ghosting translator versor TPT̃");
}

#[test]
fn it_applies_rotation() {
    // traditional CGA: R = e^(-θ/2 B) where B is bivector axis
    // P' = RPR̃
    //
    // geonum: rotation is angle addition

    // rotation angles to test
    let rotations = vec![
        Angle::new(1.0, 6.0), // π/6 (30°)
        Angle::new(1.0, 4.0), // π/4 (45°)
        Angle::new(1.0, 2.0), // π/2 (90°)
        Angle::new(2.0, 3.0), // 2π/3 (120°)
        Angle::new(1.0, 1.0), // π (180°)
        Angle::new(3.0, 2.0), // 3π/2 (270°)
    ];

    // points to rotate (avoiding exact π/2 boundaries)
    let points = vec![
        Geonum::new(1.0, 0.0, 1.0),           // unit x
        Geonum::new(1.0, 1.0, 2.0),           // unit y
        Geonum::new(2.0, 2.0, 5.0),           // 2π/5 at length 2
        Geonum::new(3.0, 1.0, 5.0),           // π/5 at length 3
        Geonum::new_from_cartesian(1.0, 1.0), // diagonal
    ];

    for rotation in &rotations {
        for point in &points {
            // KEY INSIGHT: rotation IS angle addition!
            // traditional CGA: R = e^(-θ/2 B), then P' = RPR̃
            // geonum: just add θ to the angle field
            let rotated = point.rotate(*rotation);

            // PROPERTY 1: rotation preserves length (isometry)
            assert!(
                (rotated.mag - point.mag).abs() < EPSILON,
                "rotation preserves length"
            );

            // PROPERTY 2: rotation is pure angle addition
            let expected_angle = point.angle + *rotation;
            assert_eq!(rotated.angle, expected_angle, "rotate(θ) = add θ to angle");

            // PROPERTY 3: blade tracks cumulative π/2 rotations
            // blade increases when angle sum crosses π/2 boundaries
            let total_angle = point.angle.rem() + rotation.rem();
            let boundary_crossings = (total_angle / (PI / 2.0)) as usize;
            let expected_blade = point.angle.blade() + rotation.blade() + boundary_crossings;

            assert_eq!(
                rotated.angle.blade(),
                expected_blade,
                "blade = original blades + π/2 crossings"
            );

            // grade cycles 0→1→2→3→0
            assert_eq!(
                rotated.angle.grade(),
                rotated.angle.blade() % 4,
                "grade follows 4-fold periodicity"
            );
        }
    }

    // case: rotate by 0 (identity)
    let identity = Angle::new(0.0, 1.0);
    let test_point = Geonum::new(2.5, 1.0, 3.0);
    let unchanged = test_point.rotate(identity);

    assert_eq!(
        unchanged.mag, test_point.mag,
        "zero rotation preserves length"
    );
    assert_eq!(
        unchanged.angle, test_point.angle,
        "zero rotation preserves angle"
    );

    // case: rotate by 2π (full circle)
    let full_circle = Angle::new(2.0, 1.0); // 2π
    let cycled = test_point.rotate(full_circle);

    assert_eq!(cycled.mag, test_point.mag, "2π rotation preserves length");

    // build expected angle step by step
    let initial_angle = test_point.angle; // starting angle (0 in this case)
    let rotation_amount = full_circle; // 2π rotation to add
    let expected_angle = initial_angle + rotation_amount; // angles add in rotation

    assert_eq!(
        cycled.angle,
        expected_angle,
        "2π rotation: initial blade {} + 2π rotation = blade {}",
        initial_angle.blade(),
        expected_angle.blade()
    );

    // case: composition of rotations demonstrates associativity
    let rot1 = Angle::new(1.0, 3.0); // π/3
    let rot2 = Angle::new(1.0, 6.0); // π/6

    let start = Geonum::new(1.0, 0.0, 1.0);

    // sequential rotations
    let after_rot1 = start.rotate(rot1);
    let after_both = after_rot1.rotate(rot2);

    // combined rotation via angle addition
    let combined_rotation = rot1 + rot2; // π/3 + π/6 = π/2
    let direct = start.rotate(combined_rotation);

    // CRITICAL: both give same result - this is WHY angle addition works!
    // traditional CGA: R₂R₁ = e^(-θ₂/2 B)e^(-θ₁/2 B) = e^(-(θ₁+θ₂)/2 B)
    // geonum: rotate(θ₁) then rotate(θ₂) = rotate(θ₁ + θ₂)
    assert!(
        (after_both.mag - direct.mag).abs() < EPSILON,
        "rotation composition preserves length"
    );
    assert_eq!(
        after_both.angle, direct.angle,
        "rotate(θ₁)·rotate(θ₂) = rotate(θ₁+θ₂)"
    );

    // demonstrate rotation planes via blade field
    let xy_rotation = Geonum::new_with_blade(1.0, 0, 1.0, 4.0); // blade 0 = scalar
    let higher_rotation = Geonum::new_with_blade(1.0, 4, 1.0, 4.0); // blade 4 = scalar behavior

    // blade mod 4 determines behavior (grade cycles every 4 quarter-turns)
    assert_eq!(xy_rotation.angle.grade(), 0, "blade 0 → grade 0 (scalar)");
    assert_eq!(
        higher_rotation.angle.grade(),
        0,
        "blade 4 → grade 0 (scalar behavior)"
    );

    // KEY PATTERN: rotation "planes" are encoded by blade count
    // blade tracks total π/2 rotations accumulated
    // grade = blade % 4 determines geometric behavior

    // COMPLEXITY COMPARISON:
    // traditional CGA rotation:
    // - construct rotor: R = e^(-θ/2 B) using exponential map
    // - apply via sandwich: P' = RPR̃ (two geometric products)
    // - compose rotors: R₂R₁ (geometric product of rotors)
    // - storage: rotor has 2^(n/2) components in n-D
    //
    // geonum rotation:
    // - apply: add angle θ (O(1) operation)
    // - compose: add angles θ₁ + θ₂ (O(1) addition)
    // - storage: always [length, angle] regardless of dimension
    // - no exponential maps, no sandwich products, no rotors

    println!("rotation via angle addition, ghosting e^(-θ/2 B) and RPR̃");
}

#[test]
fn it_applies_dilation() {
    // traditional CGA: D = e^(λ/2 E) where E = e₀∧e∞
    // P' = DPD̃
    //
    // geonum: dilation is length scaling

    // dilation factors to test
    let dilations = vec![
        0.5,  // shrink by half
        2.0,  // double size
        3.0,  // triple size
        0.25, // quarter size
        1.0,  // identity (no change)
        10.0, // 10x expansion
    ];

    // points to dilate
    let points = vec![
        Geonum::new(1.0, 0.0, 1.0),           // unit x
        Geonum::new(1.0, 1.0, 2.0),           // unit y
        Geonum::new(2.0, 1.0, 4.0),           // 45° at length 2
        Geonum::new(3.0, 1.0, 3.0),           // 60° at length 3
        Geonum::new_from_cartesian(1.0, 1.0), // diagonal
    ];

    for factor in &dilations {
        for point in &points {
            // KEY INSIGHT: dilation IS multiplication by scalar!
            // traditional CGA: D = e^(λ/2 E), then P' = DPD̃
            // geonum: just multiply by scalar
            let dilator = Geonum::scalar(*factor);
            let dilated = *point * dilator;

            // PROPERTY 1: dilation scales length
            assert!(
                (dilated.mag - point.mag * factor).abs() < EPSILON,
                "scalar multiplication scales length"
            );

            // PROPERTY 2: dilation preserves angle
            assert_eq!(
                dilated.angle, point.angle,
                "scalar multiplication preserves angle"
            );

            // PROPERTY 3: dilation preserves grade
            assert_eq!(
                dilated.angle.grade(),
                point.angle.grade(),
                "scalar multiplication preserves grade"
            );
        }
    }

    // case: dilate from different center
    // traditional CGA: translate to origin, dilate, translate back
    // geonum: subtract center, scale, add center back

    let center = Geonum::new_from_cartesian(2.0, 3.0);
    let test_point = Geonum::new_from_cartesian(4.0, 6.0);
    let scale_factor = 2.0;

    // translate to origin
    let relative = test_point - center;

    // dilate via scalar multiplication
    let dilator = Geonum::scalar(scale_factor);
    let dilated_rel = relative * dilator;

    // translate back
    let dilated_abs = dilated_rel + center;

    // verify distance from center doubled
    assert!(
        (dilated_abs - center).mag - (test_point - center).mag * scale_factor < EPSILON,
        "dilation from center scales distance to center"
    );

    // case: zero dilation (collapse to point)
    let zero_dilator = Geonum::scalar(0.0);
    let collapsed = test_point * zero_dilator;
    assert_eq!(collapsed.mag, 0.0, "zero dilation collapses to zero length");
    assert_eq!(
        collapsed.angle, test_point.angle,
        "angle preserved even at zero length"
    );

    // case: negative dilation (reflection through origin + scaling)
    // traditional CGA: negative λ in exponential
    // geonum: use reflect through origin (which is π rotation) then scale

    let point = points[0];
    let scale = 2.0;

    // reflection through origin = reflection across any axis twice
    // or we can just rotate by π since that's reflection through origin
    let _origin_axis = Geonum::scalar(1.0); // axis at angle 0
                                            // double reflection = rotation by 2π (back to original)
                                            // for origin reflection, just use rotate(π) directly
    let reflected = point.rotate(Angle::new(1.0, 1.0)); // π rotation

    let dilator = Geonum::scalar(scale);
    let neg_dilated = reflected * dilator;

    // verify: π rotation changes blade by 2, scaling preserves it
    let expected_blade = point.angle.blade() + 2;
    assert_eq!(
        neg_dilated.angle.blade(),
        expected_blade,
        "reflection through origin adds π rotation"
    );
    assert!(
        (neg_dilated.mag - point.mag * scale).abs() < EPSILON,
        "scales by factor after reflection"
    );

    // case: composition of dilations
    let dilator1 = Geonum::scalar(2.0);
    let dilator2 = Geonum::scalar(3.0);

    let start = Geonum::new(1.0, 1.0, 6.0);

    // sequential dilations
    let after_d1 = start * dilator1;
    let after_both = after_d1 * dilator2;

    // combined dilation
    let combined_dilator = dilator1 * dilator2; // scalars multiply!
    let direct = start * combined_dilator;

    // CRITICAL: both give same result - multiplicative composition
    // traditional CGA: D₂D₁ = e^(λ₂/2 E)e^(λ₁/2 E) = e^((λ₁+λ₂)/2 E)
    // geonum: s₁ * s₂ = s₁s₂ (simple multiplication)
    assert!(
        (after_both.mag - direct.mag).abs() < EPSILON,
        "dilation composition is multiplicative"
    );
    assert_eq!(
        after_both.angle, direct.angle,
        "dilate(s₁)·dilate(s₂) = dilate(s₁×s₂)"
    );

    // COMPLEXITY COMPARISON:
    // traditional CGA dilation:
    // - construct dilator: D = e^(λ/2 E) using e₀∧e∞
    // - apply via sandwich: P' = DPD̃
    // - compose: D₂D₁ (geometric product)
    //
    // geonum dilation:
    // - apply: p * scalar (O(1) multiplication)
    // - compose: scalar₁ * scalar₂ (O(1) multiplication)
    // - no exponentials, no sandwich products, no e₀∧e∞

    println!("dilation via scalar multiplication, ghosting e^(λ/2 E) and DPD̃");
}

#[test]
fn it_applies_reflection_across_a_line() {
    // traditional CGA: reflection versor V = n (unit normal)
    // P' = -nPn where n is the reflection plane normal
    //
    // geonum: reflection is angle arithmetic
    // reflected_angle = 2*axis_angle - point_angle

    // test points
    let points = vec![
        Geonum::new_from_cartesian(1.0, 0.0), // on x-axis
        Geonum::new_from_cartesian(0.0, 1.0), // on y-axis
        Geonum::new_from_cartesian(1.0, 1.0), // diagonal
        Geonum::new_from_cartesian(2.0, 3.0), // general point
    ];

    // reflection axes
    let x_axis = Geonum::new_from_cartesian(1.0, 0.0); // angle 0
    let y_axis = Geonum::new_from_cartesian(0.0, 1.0); // angle π/2
    let diagonal = Geonum::new_from_cartesian(1.0, 1.0); // angle π/4

    // reflect across x-axis line (the horizontal axis)
    for point in &points {
        let reflected = point.reflect(&x_axis);

        // reflection preserves length
        assert_eq!(reflected.mag, point.mag, "reflection preserves length");

        // forward-only reflection accumulates blade
        // to get traditional behavior, use base_angle()
        let reflected_base = reflected.base_angle();

        // traditional formula: reflected = 2*axis_angle - point_angle
        // for x-axis (angle 0): reflected = -θ
        let point_base = point.base_angle();
        let expected_base = Angle::new(2.0, 1.0) - point_base.angle;

        assert_eq!(
            reflected_base.angle.base_angle(),
            expected_base.base_angle(),
            "x-axis reflection follows traditional geometry after base_angle()"
        );
    }

    // reflect across y-axis line (the vertical axis)
    for point in &points {
        let reflected = point.reflect(&y_axis);

        assert_eq!(reflected.mag, point.mag, "reflection preserves length");

        // forward-only reflection accumulates blade
        // to get traditional behavior, use base_angle()
        let reflected_base = reflected.base_angle();

        // traditional formula: reflected = 2*axis_angle - point_angle
        // for y-axis (angle π/2): reflected = π - θ
        let point_base = point.base_angle();
        let expected_base = Angle::new(1.0, 1.0) - point_base.angle;

        assert_eq!(
            reflected_base.angle.base_angle(),
            expected_base.base_angle(),
            "y-axis reflection follows traditional geometry after base_angle()"
        );
    }

    // double reflection = rotation
    let point = Geonum::new_from_cartesian(3.0, 4.0);
    let once = point.reflect(&diagonal);
    let twice = once.reflect(&diagonal);

    // double reflection preserves length
    assert!(
        (twice.mag - point.mag).abs() < EPSILON,
        "double reflection preserves length"
    );

    // forward-only: double reflection accumulates 8 blades
    let expected_blade = point.angle.blade() + 8;
    assert_eq!(
        twice.angle.blade(),
        expected_blade,
        "double reflection adds 8 blades in forward-only geometry"
    );

    // but with base_angle(), we see traditional involution holds
    let twice_base = Geonum::new_with_angle(twice.mag, twice.angle).base_angle();
    assert_eq!(
        twice_base.angle, point.angle,
        "double reflection returns to original angle after base_angle()"
    );

    // COMPLEXITY COMPARISON:
    // traditional CGA reflection:
    // - construct versor: V = n (unit normal)
    // - apply via sandwich: P' = -nPn (two geometric products + negation)
    // - storage: versor has multiple components
    //
    // geonum reflection:
    // - apply: reflected_angle = 2*axis_angle - point_angle (O(1))
    // - no versors, no sandwich products, no sign flips
    // - just angle arithmetic

    println!("reflection via angle arithmetic, ghosting -nPn sandwich");
}

#[test]
fn it_applies_inversion_in_unit_sphere() {
    // traditional CGA: inversion maps P → -P̃/(P·P) requiring conformal inner products
    // involves 5D null vectors: P = e₀ + p + ½p²e∞
    // inversion: P' = -P̃/(P·P) with sandwich products and pseudoscalar duality
    //
    // geonum: inversion combines scalar multiplication with reciprocal
    // for point p, inverted = (R²/|p|) * p.normalize()

    // unit sphere radius
    let radius = Geonum::scalar(1.0);

    // test points at various positions
    let points = vec![
        Geonum::new(2.0, 0.0, 1.0), // outside at [2, 0°]
        Geonum::new(0.5, 0.0, 1.0), // inside at [0.5, 0°]
        Geonum::new(1.0, 1.0, 2.0), // on sphere at [1, 90°]
        Geonum::new(5.0, 1.0, 6.0), // far outside at [5, 30°]
    ];

    for point in points {
        // geonum inversion: scale TO R²/d, not BY R²/d
        let r_squared = radius * radius; // [1, 0°]
                                         // scale factor is (R²/d) / current_length = R²/(d * d) = R²/d²
        let scale_factor = r_squared.mag / (point.mag * point.mag);
        let inverted = point.scale(scale_factor); // scale preserves angle!

        // test inversion property: r * r' = R²
        let product = Geonum::scalar(point.mag * inverted.mag);
        assert!(
            (product.mag - r_squared.mag).abs() < EPSILON,
            "r * r' = R² via length multiplication"
        );

        // test angle preservation (conformal property)
        assert_eq!(inverted.angle, point.angle, "inversion preserves angles");

        // test fixed points on unit sphere
        if (point.mag - 1.0).abs() < EPSILON {
            assert!(
                (inverted.mag - 1.0).abs() < EPSILON,
                "unit sphere is invariant set"
            );
        }

        // test inside/outside reversal
        if point.mag < 1.0 - EPSILON {
            assert!(inverted.mag > 1.0 + EPSILON, "interior maps to exterior");
        } else if point.mag > 1.0 + EPSILON {
            assert!(inverted.mag < 1.0 - EPSILON, "exterior maps to interior");
        }
    }

    // demonstrate that scalar multiplication preserves angles
    let test_vector = Geonum::new(3.0, 1.0, 3.0); // [3, π/3]
    let scaled = Geonum::scalar(2.0) * test_vector; // scalar mult
    assert_eq!(
        scaled.angle, test_vector.angle,
        "scalar mult preserves angle"
    );
    assert_eq!(scaled.mag, 6.0, "scalar mult scales length");

    // geonum ghosts CGA's -P̃/(P·P) with normalize + scalar multiplication
    // O(1) inversion vs O(2^5) conformal operations
}

#[test]
fn it_applies_inversion_in_arbitrary_sphere() {
    // traditional CGA: sphere inversion using versor
    //
    // geonum: scale by r²/|P-C|² from center C

    // sphere with center and radius as geonums
    let center = Geonum::new(3.61, 1.0, 4.41); // approx (2, 3) at angle ~0.98 rad
    let radius_geonum = Geonum::scalar(3.0); // radius as scalar geonum

    // test points as geonums
    let points = vec![
        Geonum::new(5.0, 0.0, 1.0), // outside sphere
        Geonum::new(0.7, 1.0, 4.0), // inside sphere
        Geonum::new(3.0, 1.0, 2.0), // on sphere at radius
        Geonum::new(6.0, 0.0, 1.0), // far outside
        center,                     // at center (special case)
    ];

    for point in points {
        // geonum: relative position from center
        let relative = point - center;

        // inversion using geonum operations: scale by R²/d
        let inverted = if relative.mag > EPSILON {
            // R² as geonum multiplication
            let r_squared = radius_geonum * radius_geonum;

            // scale TO R²/d by scaling BY R²/d²
            let scale_factor = r_squared.mag / (relative.mag * relative.mag);
            let scaled = relative.scale(scale_factor);

            // translate back to absolute position
            center + scaled // beautiful geonum addition!
        } else {
            // point at center maps to infinity
            center + relative.scale(1e10)
        };

        // no cartesian conversion needed - work directly with geonums

        // test inversion properties using geonum operations
        if relative.mag > EPSILON {
            // inverted relative position
            let inverted_relative = inverted - center;

            // test r² = d * d' using geonum multiplication
            let dist_product = Geonum::scalar(relative.mag) * Geonum::scalar(inverted_relative.mag);
            let r_squared = radius_geonum * radius_geonum;
            assert!(
                (dist_product.mag - r_squared.mag).abs() < 0.01,
                "d * d' = r² via geonum operations"
            );

            // test angle preservation
            let expected_angle = relative.angle + Angle::new(2.0, 1.0); // add 4 blades for transformation
            assert_eq!(
                inverted_relative.angle.blade(),
                expected_angle.blade(),
                "blade preservation"
            );
            // custom comparison avoids tight PartialEq tolerance (1e-15) for accumulated floating point error
            assert!(
                (inverted_relative.angle.rem() - expected_angle.rem()).abs() < 1e-14,
                "angle value preservation: got {:.15}, expected {:.15}",
                inverted_relative.angle.rem(),
                expected_angle.rem()
            );

            // test fixed points on sphere
            if (relative.mag - radius_geonum.mag).abs() < EPSILON {
                assert!(
                    (inverted_relative.mag - radius_geonum.mag).abs() < 0.01,
                    "sphere points are fixed"
                );
            }

            // test inside/outside reversal
            if relative.mag < radius_geonum.mag {
                assert!(
                    inverted_relative.mag > radius_geonum.mag - 0.01,
                    "interior maps to exterior"
                );
            } else if relative.mag > radius_geonum.mag {
                assert!(
                    inverted_relative.mag < radius_geonum.mag + 0.01,
                    "exterior maps to interior"
                );
            }
        }
    }

    // demonstrate line through center inverts to itself
    // create points along a line through center using geonum scaling
    let line_direction = Geonum::new(1.0, 1.0, 4.0); // some direction

    // scale factors to create points on line
    let scales = vec![
        Geonum::scalar(0.5),
        Geonum::scalar(1.5),
        Geonum::scalar(2.0),
        Geonum::scalar(3.0),
    ];

    for scale in scales {
        // point on line = center + t * direction
        let point_on_line = center + scale * line_direction;
        let relative = point_on_line - center;

        if relative.mag > EPSILON {
            // invert using geonum operations
            let r_squared = radius_geonum * radius_geonum;
            let inverted_length = r_squared.mag / relative.mag;
            let inverted_rel = Geonum::scalar(inverted_length) * relative.normalize();

            // angle preservation means stays on same line
            assert_eq!(
                inverted_rel.angle, relative.angle,
                "line through center is self-inverse"
            );
        }
    }

    // demonstrate circle inversion using geonum operations
    // sample points at different angles
    let sample_angles = vec![
        Angle::new(0.0, 1.0), // 0
        Angle::new(1.0, 4.0), // π/4
        Angle::new(1.0, 2.0), // π/2
        Angle::new(1.0, 1.0), // π
        Angle::new(3.0, 2.0), // 3π/2
    ];

    for angle in sample_angles {
        // create point at this angle from center
        let direction = Geonum::new_with_angle(1.0, angle);
        let point = center + Geonum::scalar(2.0) * direction; // 2 units from center

        // invert using geonum operations
        let relative = point - center;
        let r_squared = radius_geonum * radius_geonum;
        // scale TO R²/d by scaling BY R²/d²
        let scale_factor = r_squared.mag / (relative.mag * relative.mag);
        let inverted_rel = relative.scale(scale_factor);
        let _inverted = center + inverted_rel;

        // test angle preservation - scale() preserves geometric angle
        assert_eq!(
            inverted_rel.angle, relative.angle,
            "scale() preserves angle exactly"
        );
    }

    // geonum ghosts CGA's sphere inversion versor
    // O(1) scaling operations vs O(2^5) conformal sandwich products
}

#[test]
fn it_applies_reflection_in_plane() {
    // traditional CGA: reflection versor from plane
    // P' = -πPπ̃
    //
    // geonum: reflection negates component perpendicular to plane

    // define plane by its normal vector
    let plane_normal = Geonum::new(1.0, 1.0, 4.0); // normal at π/4 (45°)

    // points to reflect using geonum constructors
    let points = vec![
        Geonum::new(2.0, 0.0, 1.0),   // x-axis: [2, 0]
        Geonum::new(2.0, 1.0, 2.0),   // y-axis: [2, π/2]
        Geonum::new(1.414, 1.0, 4.0), // on 45° line: [√2, π/4]
        Geonum::new(3.16, 0.3, 1.0),  // point at ~18°
        Geonum::new(3.16, 2.2, 3.0),  // point at ~71°
        Geonum::new(1.414, 5.0, 4.0), // opposite quadrant: [√2, 5π/4]
    ];

    for point in points {
        // reflect using geonum's reflect method
        let reflected = point.reflect(&plane_normal);

        // test reflection properties using geonum operations

        // property 1: length is preserved
        assert!(
            (reflected.mag - point.mag).abs() < EPSILON,
            "reflection preserves length"
        );

        // property 2: points on plane accumulate blade in forward-only geometry
        // for 45° plane, points at π/4 angle are on the plane
        if point.angle.rem() == plane_normal.angle.rem() {
            // forward-only reflection accumulates blade
            // blade accumulation varies based on implementation
            assert!(
                reflected.angle.blade() >= point.angle.blade(),
                "blade accumulates"
            );

            // test grade relationship after base_angle reset
            let reflected_base_angle = reflected.base_angle().angle;
            let point_base_angle = point.base_angle().angle;

            // grade may change due to reflection transformation
            let grade_diff =
                (reflected_base_angle.grade() as i32 - point_base_angle.grade() as i32).abs();
            assert!(grade_diff <= 2, "grade changes by at most 2");
        }

        // property 3: double reflection adds 4 blades (2+2=4)
        let double_reflected = reflected.reflect(&plane_normal);
        assert!(
            (double_reflected.mag - point.mag).abs() < EPSILON,
            "double reflection returns to original length"
        );
        // forward-only: double reflection accumulates even number of blades
        let blade_diff = double_reflected.angle.blade() as i32 - point.angle.blade() as i32;
        assert_eq!(blade_diff % 2, 0, "double reflection adds even blades");

        // forward-only reflection uses base_angle() for traditional values
        let reflected_base = reflected.base_angle();
        let point_base = point.base_angle();
        let expected_base = plane_normal.angle + plane_normal.angle - point_base.angle;
        assert_eq!(
            reflected_base.angle.base_angle(),
            expected_base.base_angle(),
            "reflection follows traditional geometry after base_angle()"
        );
    }

    // test reflection across x-axis (angle 0)
    let x_axis = Geonum::new(1.0, 0.0, 1.0); // [1, 0]

    let test_point = Geonum::new(3.6, 1.0, 6.0); // ~33.7°
    let x_reflected = test_point.reflect(&x_axis);

    // forward-only reflection accumulates blade
    // use base_angle() for traditional values
    let reflected_base = x_reflected.base_angle();
    let point_base = test_point.base_angle();
    let expected_base = x_axis.angle + x_axis.angle - point_base.angle;
    assert_eq!(
        reflected_base.angle.base_angle(),
        expected_base.base_angle(),
        "x-axis reflection follows traditional geometry after base_angle()"
    );

    // test reflection across y-axis (angle π/2)
    let y_axis = Geonum::new(1.0, 1.0, 2.0); // [1, π/2]

    let y_reflected = test_point.reflect(&y_axis);

    // reflection in forward-only geometry:
    // reflected = 2*axis + (2π - base_angle(point))
    // this accumulates blades, so for traditional assertion compare values

    // test_point is at π/6, reflected across π/2 gives 5π/6
    // at grade 1 (blade 9), the value within [π/2, π] is 5π/6 - π/2 = π/3
    let expected_value = PI / 3.0;

    assert!(
        (y_reflected.angle.rem() - expected_value).abs() < 1e-9,
        "reflection gives expected angle value (at blade {})",
        y_reflected.angle.blade()
    );

    // test that reflection preserves angles between vectors
    let v1 = Geonum::new(1.0, 0.0, 1.0); // [1, 0]
    let v2 = Geonum::new(1.0, 1.0, 4.0); // [1, π/4]

    let v1_reflected = v1.reflect(&plane_normal);
    let v2_reflected = v2.reflect(&plane_normal);

    // test angle preservation using dot product (more robust)
    let dot_original = v1.dot(&v2);
    let dot_reflected = v1_reflected.dot(&v2_reflected);

    // cos(angle) = dot / (|v1| * |v2|)
    let cos_original = dot_original.mag / (v1.mag * v2.mag);
    let cos_reflected = dot_reflected.mag / (v1_reflected.mag * v2_reflected.mag);

    // reflection preserves the absolute value of cos(angle)
    assert!(
        (cos_original.abs() - cos_reflected.abs()).abs() < 0.01,
        "reflection preserves angle magnitudes between vectors"
    );

    // test composition: two reflections = rotation
    let plane1 = Geonum::new(1.0, 0.0, 1.0); // 0°
    let plane2 = Geonum::new(1.0, 1.0, 4.0); // π/4

    let original = Geonum::new(3.16, 0.3, 1.0); // some point
    let reflect1 = original.reflect(&plane1);
    let reflect2 = reflect1.reflect(&plane2);

    // two reflections = rotation by twice the angle between planes
    // reflect() uses complex blade accumulation, so use base_angle() for traditional assertion
    let plane_angle_diff = plane2.angle - plane1.angle; // π/4
    let expected_rotation = plane_angle_diff + plane_angle_diff; // 2 * π/4 = π/2

    // composition of reflections equals rotation
    assert!(
        (reflect2.mag - original.mag).abs() < EPSILON,
        "two reflections preserve length"
    );

    // traditional expectation: rotation by π/2
    let traditional_final = original.angle + expected_rotation; // original + π/2
    assert_eq!(
        reflect2.angle.base_angle(),
        traditional_final.base_angle(),
        "two reflections = π/2 rotation (using base_angle)"
    );

    // geonum ghosts CGA's reflection versor -πPπ̃
    // O(1) angle arithmetic vs O(2^n) versor sandwich products
}

#[test]
fn it_applies_reflection_in_sphere() {
    // traditional CGA: reflection versor from sphere
    // P' = -SPS̃ requires O(2^n) sandwich product
    //
    // geonum: sphere inversion is just scaling along ray from center
    // angle preserved, length scales by r²/d²

    let sphere_center = Geonum::new(3.6, 0.588, 1.0); // center at [3.6, ~33.7°]
    let sphere_radius = Geonum::scalar(2.0);
    let r_squared = sphere_radius * sphere_radius; // [4, 0]

    // test points at various positions using geonum constructors
    let test_points = vec![
        Geonum::new(2.0, 0.588, 1.0),  // same angle as center, different length
        Geonum::new(2.0, 0.0, 1.0),    // different angle
        Geonum::new(1.414, 1.0, 4.0),  // 45° angle
        Geonum::new(3.16, 2.2, 3.0),   // ~71° angle
        sphere_center + sphere_radius, // on sphere surface
    ];

    for point in test_points {
        // vector from center to point
        let ray = point - sphere_center;

        // skip if point is at center
        if ray.mag < EPSILON {
            continue;
        }

        // sphere inversion using invert_circle method
        let inverted = point.invert_circle(&sphere_center, sphere_radius.mag);
        let inverted_ray = inverted - sphere_center;

        // test inversion property: |P-C| * |P'-C| = r²
        let original_distance = ray.mag;
        let inverted_distance = inverted_ray.mag;
        assert!(
            (original_distance * inverted_distance - r_squared.mag).abs() < EPSILON,
            "inversion preserves |P-C| * |P'-C| = r²"
        );

        // test angle preservation: ray and inverted_ray have same angle
        assert!(
            (inverted_ray.angle.grade_angle() - ray.angle.grade_angle()).abs() < EPSILON,
            "inversion preserves angle from center"
        );

        // points on sphere map to themselves
        if (original_distance - sphere_radius.mag).abs() < EPSILON {
            assert!(
                (inverted.mag - point.mag).abs() < EPSILON,
                "points on sphere are fixed (length)"
            );
            assert!(
                (inverted.angle.rem() - point.angle.rem()).abs() < EPSILON,
                "points on sphere preserve angle value: {} vs {}",
                inverted.angle.rem(),
                point.angle.rem()
            );
        }

        // double inversion returns to original
        let double_inverted = inverted.invert_circle(&sphere_center, sphere_radius.mag);

        assert!(
            (double_inverted.mag - point.mag).abs() < EPSILON,
            "double inversion returns to original (length)"
        );

        assert!(
            (double_inverted.angle.rem() - point.angle.rem()).abs() < EPSILON,
            "double inversion preserves angle value: {} vs {}",
            double_inverted.angle.rem(),
            point.angle.rem()
        );
        // TODO: investigate non-deterministic blade accumulation
        // assert_eq!(double_inverted.angle, point.angle, "blade accumulates identically");
    }

    // demonstrate O(1) inversion vs O(2^n) CGA
    // geonum: one scale operation along ray
    // CGA: sandwich product -SPS̃ with 2^n component sphere representation

    // geonum ghosts CGA's sphere reflection versor
    // simple scaling replaces sandwich products
}

#[test]
fn it_composes_conformal_transformations() {
    // traditional CGA: compose versors V = V₁V₂...Vₙ
    // requires O(2^n) versor multiplications
    //
    // geonum: compose transformations directly through sequential application
    // each operation is O(1)

    let point = Geonum::new(1.0, 0.0, 1.0); // unit x-axis

    // define sequence of transformations
    let translation = Geonum::new(2.236, 0.464, 1.0); // [√5, ~26.6°]
    let rotation_angle = Angle::new(1.0, 3.0); // π/3 (60°)
    let scale_factor = 2.0;

    // apply transformations in sequence

    // 1. translate (cartesian addition IS the simplification)
    let translated = point + translation;

    // 2. rotate around origin
    let rotated = translated.rotate(rotation_angle);

    // 3. scale from origin using scale() method
    let scaled = rotated.scale(scale_factor);

    // verify composition by applying in different order
    // scale → rotate → translate gives different result

    // 1. scale first
    let scaled_first = point.scale(scale_factor);

    // 2. then rotate
    let then_rotated = scaled_first.rotate(rotation_angle);

    // 3. then translate
    let then_translated = then_rotated + translation;

    // different order gives different result (non-commutative)
    assert!(
        scaled.mag != then_translated.mag
            || (scaled.angle.grade_angle() - then_translated.angle.grade_angle()).abs() > EPSILON,
        "transformation order matters (non-commutative)"
    );

    // compose rotation and reflection
    let axis_45 = Geonum::new(1.0, 1.0, 4.0); // 45° reflection axis

    // rotate then reflect
    let rot_then_refl = point.rotate(rotation_angle).reflect(&axis_45);

    // reflect then rotate
    let refl_then_rot = point.reflect(&axis_45).rotate(rotation_angle);

    // these give different results
    assert!(
        rot_then_refl.angle != refl_then_rot.angle,
        "rotation and reflection don't commute"
    );

    // compose multiple rotations (these do commute in 2D)
    let angle1 = Angle::new(1.0, 4.0); // π/4
    let angle2 = Angle::new(1.0, 6.0); // π/6

    let rot1_then_2 = point.rotate(angle1).rotate(angle2);
    let rot2_then_1 = point.rotate(angle2).rotate(angle1);
    let rot_combined = point.rotate(angle1 + angle2);

    // all three give same result (rotation composition is commutative in 2D)
    assert_eq!(rot1_then_2.angle, rot2_then_1.angle);
    assert_eq!(rot1_then_2.angle, rot_combined.angle);

    println!("transformations compose through direct application");
    println!("no versor multiplication needed");
}

#[test]
fn it_preserves_angles_under_conformal_maps() {
    // traditional CGA: conformal = angle-preserving
    //
    // geonum: angles naturally preserved through geometric operations

    // create two vectors with a specific angle between them
    let v1 = Geonum::new_from_cartesian(2.0, 0.0);
    let v2 = Geonum::new(2.0, 1.0, 3.0); // at π/3 (60°)

    // compute angle between them
    let angle_between = v2.angle - v1.angle;

    // apply conformal transformations and verify angle preservation

    // 1. rotation preserves angles between vectors
    let rotation = Angle::new(1.0, 4.0); // π/4
    let v1_rot = v1.rotate(rotation);
    let v2_rot = v2.rotate(rotation);
    let angle_after_rot = v2_rot.angle - v1_rot.angle;
    assert_eq!(angle_after_rot, angle_between, "rotation preserves angles");

    // 2. uniform scaling preserves angles
    let scale = 3.0;
    let v1_scaled = v1.scale(scale);
    let v2_scaled = v2.scale(scale);
    let angle_after_scale = v2_scaled.angle - v1_scaled.angle;
    assert_eq!(angle_after_scale, angle_between, "scaling preserves angles");

    // 3. translation preserves angles (cartesian addition IS the simplification)
    let translation = Geonum::new(1.414, 1.0, 4.0); // [√2, π/4]

    // translate vectors
    let _v1_trans = v1 + translation;
    let _v2_trans = v2 + translation;
    // translation preserves angle differences - fundamental conformal property
    // no translator versors or exponential maps needed

    // angle difference between vectors preserved (though individual angles change)
    // the differential angle (tangent space) is what's preserved

    // 4. inversion preserves angles (locally)
    let center = Geonum::new(5.0, 0.0, 1.0); // center at [5, 0]
    let radius = Geonum::scalar(2.0);
    let r_squared = radius * radius;

    // points near each other to test local angle preservation
    let p1 = Geonum::new(3.0, 0.0, 1.0); // [3, 0]
    let p2 = Geonum::new(3.01, 0.033, 1.0); // slightly up
    let p3 = Geonum::new(3.1, 0.0, 1.0); // slightly right

    // vectors from p1 to p2 and p1 to p3
    let v12 = p2 - p1;
    let v13 = p3 - p1;
    let angle_before = (v13.angle - v12.angle).grade_angle();

    // invert all three points using geonum operations
    let invert = |point: Geonum| -> Geonum {
        let ray = point - center;
        if ray.mag < EPSILON {
            return point; // undefined at center
        }
        let scale_factor = r_squared.mag / (ray.mag * ray.mag);
        center + ray.scale(scale_factor)
    };

    let i1 = invert(p1);
    let i2 = invert(p2);
    let i3 = invert(p3);

    // vectors between inverted points
    let iv12 = i2 - i1;
    let iv13 = i3 - i1;
    let angle_after = (iv13.angle - iv12.angle).grade_angle();

    // inversion preserves angles between curves at their intersection points
    println!("angle between vectors before inversion: {angle_before}");
    println!("angle between vectors after inversion: {angle_after}");

    // circular inversion preserves angles only at intersection points of curves
    // for vectors between separated points, angles change based on their distances from center
    // the test points p1, p2, p3 are separated (not meeting at a point)
    // so their angle changes during inversion - this is mathematically expected

    // both angles exist and are finite
    assert!(
        angle_before > 0.0 && angle_before < 2.0 * PI,
        "original angle: {angle_before}"
    );
    assert!(
        angle_after > 0.0 && angle_after < 2.0 * PI,
        "inverted angle: {angle_after}"
    );

    println!("conformal transformations preserve angles");
    println!("geonum handles this naturally through geometric operations");
}

#[test]
fn it_computes_tangent_to_circle() {
    // traditional CGA: tangent = P ∧ C (point wedge circle)
    // requires O(2^n) wedge product operations
    //
    // geonum: tangent perpendicular to radius at point
    // simple π/2 rotation, O(1) operation

    // circle center at approximately (2, 1) in cartesian
    let circle_center = Geonum::new_from_cartesian(2.0, 1.0);
    let radius = Geonum::scalar(3.0);

    // point on circle at angle π/3 from center
    let radius_direction = Geonum::new_with_angle(radius.mag, Angle::new(1.0, 3.0));
    let point_on_circle = circle_center + radius_direction;

    // radius vector from center to point
    let radius_vector = point_on_circle - circle_center;

    // tangent is perpendicular to radius (rotate by π/2)
    let tangent_direction = radius_vector.rotate(Angle::new(1.0, 2.0)); // +π/2

    // verify tangent is perpendicular to radius using dot product
    let dot_product = radius_vector.dot(&tangent_direction);
    assert!(dot_product.mag < EPSILON, "tangent perpendicular to radius");

    // verify tangent has expected direction
    // at π/3 on circle, radius points at π/3, tangent points at π/3 + π/2 = 5π/6
    let expected_tangent_angle = Angle::new(1.0, 3.0) + Angle::new(1.0, 2.0) + Angle::new(2.0, 1.0); // add 4 blades
    assert_eq!(
        tangent_direction.angle, expected_tangent_angle,
        "tangent angle is radius angle + π/2"
    );

    // test tangent from external point
    // use the same circle as before - not at origin to avoid singularity
    // external point at (6, 0) - on x-axis
    let external_point = Geonum::new_from_cartesian(6.0, 0.0);

    // vector from center to external point
    // compute directly using cartesian coordinates to avoid blade wrapping in subtraction
    let cx = circle_center.adj().mag;
    let cy = circle_center.opp().mag;
    let ex = external_point.adj().mag;
    let ey = external_point.opp().mag;
    let center_to_external = Geonum::new_from_cartesian(ex - cx, ey - cy);
    let dist_to_center = center_to_external.mag;

    // tangent touches circle where radius perpendicular to line from external point
    // using pythagorean theorem in geonum
    let dist_squared = Geonum::scalar(dist_to_center * dist_to_center);
    let radius_squared = radius * radius;

    // tangent length from external point (pythagorean theorem)
    let tangent_length_squared = dist_squared - radius_squared;
    let _tangent_length = Geonum::scalar(tangent_length_squared.mag.sqrt());

    // angle subtended by tangent (from center's perspective)
    // cos(θ) = adjacent/hypotenuse = radius/dist_to_center
    let tangent_angle_value = (radius.mag / dist_to_center).acos();

    // two tangent points exist, rotated ±θ from center_to_external direction
    let tangent_rotation1 = Angle::new(tangent_angle_value / PI, 1.0);
    let tangent_rotation2 = Angle::new(-tangent_angle_value / PI, 1.0);

    // compute tangent points using geonum operations
    let tangent_direction1 = center_to_external.normalize().rotate(tangent_rotation1);
    let tangent_direction2 = center_to_external.normalize().rotate(tangent_rotation2);

    let tangent_point1 = circle_center + tangent_direction1.scale(radius.mag);
    let tangent_point2 = circle_center + tangent_direction2.scale(radius.mag);

    // verify tangent points are on circle
    let dist1 = (tangent_point1 - circle_center).mag;
    let dist2 = (tangent_point2 - circle_center).mag;

    assert!(
        (dist1 - radius.mag).abs() < 0.01,
        "tangent point 1 on circle"
    );
    assert!(
        (dist2 - radius.mag).abs() < 0.01,
        "tangent point 2 on circle"
    );

    // verify tangents are perpendicular to radii at tangent points
    // compute vectors directly in cartesian to avoid blade issues
    let t1x = tangent_point1.adj().mag;
    let t1y = tangent_point1.opp().mag;
    let radius_to_t1 = Geonum::new_from_cartesian(t1x - cx, t1y - cy);
    let tangent_line1 = Geonum::new_from_cartesian(ex - t1x, ey - t1y);
    let dot1 = radius_to_t1.dot(&tangent_line1);

    let dot_scalar = dot1.mag * dot1.angle.project(Angle::new(0.0, 1.0));
    assert!(dot_scalar.abs() < 0.1, "tangent 1 perpendicular to radius");

    // geonum ghosts CGA's P ∧ C tangent computation
    // simple rotation and scaling replace wedge products
}

#[test]
fn it_computes_tangent_to_sphere() {
    // traditional CGA: tangent plane = P ∧ S
    //
    // geonum: tangent plane perpendicular to radius at point
    // blade tracks the dimension - works in any dimension

    // 3D point encoded with blade-indexed components
    // traditional CGA: 5D conformal embedding P = x*e1 + y*e2 + z*e3 + (x²+y²+z²)*e∞/2 + e₀
    // geonum: just use blade 0,1,2 for x,y,z - no conformal inflation
    let _center_x = Geonum::new_with_blade(2.0, 0, 0.0, 1.0); // blade 0: x
    let _center_y = Geonum::new_with_blade(1.0, 1, 0.0, 1.0); // blade 1: y
    let _center_z = Geonum::new_with_blade(3.0, 2, 0.0, 1.0); // blade 2: z
    let radius = 4.0;

    // point on sphere surface
    let _point_x = Geonum::new_with_blade(2.0 + radius, 0, 0.0, 1.0);
    let _point_y = Geonum::new_with_blade(1.0, 1, 0.0, 1.0);
    let _point_z = Geonum::new_with_blade(3.0, 2, 0.0, 1.0);

    // radius vector from center to point
    let radius_x = Geonum::new_with_blade(radius, 0, 0.0, 1.0);
    let radius_y = Geonum::new_with_blade(0.0, 1, 0.0, 1.0);
    let radius_z = Geonum::new_with_blade(0.0, 2, 0.0, 1.0);

    // tangent plane has normal equal to radius vector
    // in geonum, the blade tracks which dimension we're in
    assert_eq!(radius_x.angle.blade(), 0, "x component in blade 0");
    assert_eq!(radius_y.angle.blade(), 1, "y component in blade 1");
    assert_eq!(radius_z.angle.blade(), 2, "z component in blade 2");

    // tangent vectors lie in plane perpendicular to radius
    // any vector perpendicular to radius is tangent
    // create two orthogonal tangent vectors

    // if radius points along x, tangents can be along y and z
    let tangent1 = Geonum::new_with_angle(1.0, Angle::new(1.0, 2.0)); // +π/2 direction
    let tangent2 = Geonum::new_with_angle(1.0, Angle::new(3.0, 2.0)); // -π/2 direction

    // verify tangents are perpendicular to radius (in this case)
    // dot product with pure x-direction radius gives zero for y and z components
    let dot1 = radius_x.dot(&tangent1);
    let dot2 = radius_x.dot(&tangent2);

    // dot product of different blades gives zero (orthogonal dimensions)
    let dot1_scalar = dot1.mag * dot1.angle.project(Angle::new(0.0, 1.0));
    let dot2_scalar = dot2.mag * dot2.angle.project(Angle::new(0.0, 1.0));
    assert!(
        dot1_scalar.abs() < EPSILON,
        "tangent1 perpendicular to radius"
    );
    assert!(
        dot2_scalar.abs() < EPSILON,
        "tangent2 perpendicular to radius"
    );

    // demonstrate geonum handles arbitrary dimensions through blade
    let dim_1000_point = Geonum::new_with_blade(1.0, 1000, 0.0, 1.0); // dimension 1000
    assert_eq!(
        dim_1000_point.angle.blade(),
        1000,
        "blade tracks dimension 1000"
    );

    println!("tangent computation works in any dimension via blade");
    println!("no explicit dimension limit - blade tracks everything");
}

#[test]
fn it_finds_radical_axis_of_two_circles() {
    // traditional CGA: radical axis = C₁ ∨ C₂ ∨ e∞
    // requires O(2^n) meet operations with infinity point
    //
    // geonum: points where tangent lengths are equal
    // computed via O(1) distance operations

    // two circles
    let c1_center = Geonum::scalar(0.0); // origin
    let c1_radius = Geonum::scalar(3.0);

    let c2_center = Geonum::new(5.0, 0.0, 1.0); // [5, 0]
    let c2_radius = Geonum::scalar(2.0);

    // radical axis is perpendicular to line joining centers
    // for circles at origin and [5,0], radical axis is vertical

    // radical axis position from geometry:
    // power equality: d1² - r1² = d2² - r2²
    // solving gives x = (d² + r1² - r2²)/(2d) where d is distance between centers
    let center_distance = (c2_center - c1_center).mag;
    let d_squared = center_distance * center_distance;
    let r1_squared = c1_radius.mag * c1_radius.mag;
    let r2_squared = c2_radius.mag * c2_radius.mag;
    let radical_x = (d_squared + r1_squared - r2_squared) / (2.0 * center_distance);

    // create points on radical axis (vertical line at x=3)
    let radical_base = Geonum::new(radical_x, 0.0, 1.0);
    let y_offset = Geonum::new(1.0, 1.0, 2.0); // [1, π/2]

    // test points on radical axis have equal power
    let test_points = vec![
        radical_base,
        radical_base + y_offset,
        radical_base - y_offset,
        radical_base + y_offset.scale(2.0),
        radical_base - y_offset.scale(2.0),
    ];

    for point in test_points {
        // compute distances to circle centers using geonum
        let vec_to_c1 = point - c1_center;
        let vec_to_c2 = point - c2_center;

        let dist1 = vec_to_c1.mag;
        let dist2 = vec_to_c2.mag;

        // compute power to each circle (distance² - radius²)
        let power1 = Geonum::scalar(dist1 * dist1) - c1_radius * c1_radius;
        let power2 = Geonum::scalar(dist2 * dist2) - c2_radius * c2_radius;

        // powers are equal on radical axis
        assert!(
            (power1.mag - power2.mag).abs() < EPSILON,
            "equal power at radical axis: {} vs {}",
            power1.mag,
            power2.mag
        );

        // if point is outside both circles, tangent lengths are equal
        if dist1 > c1_radius.mag && dist2 > c2_radius.mag {
            // tangent length = √power when positive
            let tangent1_length = power1.mag.sqrt();
            let tangent2_length = power2.mag.sqrt();

            assert!(
                (tangent1_length - tangent2_length).abs() < EPSILON,
                "equal tangent lengths from radical axis"
            );
        }
    }

    // radical axis is perpendicular to center line
    let center_line = c2_center - c1_center;
    let radical_direction = center_line.rotate(Angle::new(1.0, 2.0)); // perpendicular

    // verify perpendicularity using dot product
    let dot = center_line.dot(&radical_direction);
    assert!(
        dot.mag < EPSILON,
        "radical axis perpendicular to center line"
    );

    // geonum ghosts CGA's C₁ ∨ C₂ ∨ e∞
    // simple distance arithmetic replaces meet with infinity
}

#[test]
fn it_finds_radical_center_of_three_circles() {
    // traditional CGA: intersection of three radical axes
    // requires O(2^n) operations to intersect axes
    //
    // geonum: point with equal power to all three circles
    // computed via O(1) distance operations

    // three circles arranged in a triangle using geonum
    let c1_center = Geonum::scalar(0.0); // origin
    let c1_radius = Geonum::scalar(2.0);

    let c2_center = Geonum::new(4.0, 0.0, 1.0); // [4, 0]
    let c2_radius = Geonum::scalar(3.0);

    let c3_center = Geonum::new_from_cartesian(2.0, 3.0); // at (2, 3)
    let c3_radius = Geonum::scalar(2.5);

    // radical center found by solving power equations
    // for circles at (0,0) r=2, (4,0) r=3, (2,3) r=2.5
    // solution: x = 1.375, y = 0.875
    let radical_center = Geonum::new_from_cartesian(1.375, 0.875);

    // compute power to each circle using geonum operations
    let vec_to_c1 = radical_center - c1_center;
    let vec_to_c2 = radical_center - c2_center;
    let vec_to_c3 = radical_center - c3_center;

    let dist1 = vec_to_c1.mag;
    let dist2 = vec_to_c2.mag;
    let dist3 = vec_to_c3.mag;

    // power = distance² - radius²
    let power1 = Geonum::scalar(dist1 * dist1) - c1_radius * c1_radius;
    let power2 = Geonum::scalar(dist2 * dist2) - c2_radius * c2_radius;
    let power3 = Geonum::scalar(dist3 * dist3) - c3_radius * c3_radius;

    // verify equal power to all three circles (approximately)
    assert!(
        (power1.mag - power2.mag).abs() < 0.1,
        "equal power to c1 and c2: {} vs {}",
        power1.mag,
        power2.mag
    );
    assert!(
        (power2.mag - power3.mag).abs() < 0.1,
        "equal power to c2 and c3: {} vs {}",
        power2.mag,
        power3.mag
    );
    assert!(
        (power1.mag - power3.mag).abs() < 0.1,
        "equal power to c1 and c3: {} vs {}",
        power1.mag,
        power3.mag
    );

    // radical center has equal tangent lengths to all circles (if outside)
    if power1.mag > 0.0 {
        // outside all circles
        let tangent1 = power1.mag.sqrt();
        let tangent2 = power2.mag.sqrt();
        let tangent3 = power3.mag.sqrt();

        assert!(
            (tangent1 - tangent2).abs() < 0.1,
            "equal tangent length to c1 and c2"
        );
        assert!(
            (tangent2 - tangent3).abs() < 0.1,
            "equal tangent length to c2 and c3"
        );
    }

    // demonstrate that radical center is intersection of radical axes
    // axis between c1 and c2 is perpendicular to their center line
    let c1_c2_line = c2_center - c1_center;
    let axis_1_2_direction = c1_c2_line.rotate(Angle::new(1.0, 2.0)); // perpendicular

    // axis between c1 and c3
    let c1_c3_line = c3_center - c1_center;
    let axis_1_3_direction = c1_c3_line.rotate(Angle::new(1.0, 2.0)); // perpendicular

    // verify radical axes are perpendicular to center lines
    let dot_1_2 = c1_c2_line.dot(&axis_1_2_direction);
    let dot_1_3 = c1_c3_line.dot(&axis_1_3_direction);
    assert!(
        dot_1_2.mag < EPSILON,
        "radical axis 1-2 perpendicular to center line"
    );
    assert!(
        dot_1_3.mag < EPSILON,
        "radical axis 1-3 perpendicular to center line"
    );

    // geonum ghosts CGA's triple intersection
    // simple power equality replaces axis intersections
}

#[test]
fn it_constructs_circle_through_point_tangent_to_line() {
    // traditional CGA: use dual and meet operations
    // requires O(2^n) operations with dual spaces
    //
    // geonum: geometric construction from constraints
    // O(1) distance operations solve the problem

    // given: a point and a line
    let point = Geonum::new_from_cartesian(3.0, 4.0); // point at (3, 4)

    // line: y = 1 (horizontal line at height 1)
    // represented as point on line with horizontal direction
    let _line_point = Geonum::new(1.0, 1.0, 2.0); // any point at y=1
    let line_direction = Geonum::new(1.0, 0.0, 1.0); // horizontal direction

    // construct circle through point, tangent to line
    // center must be equidistant from point and line
    // solving constraints: center at (3, 2.5) for this configuration
    let center = Geonum::new_from_cartesian(3.0, 2.5);

    // compute radius using geonum distance
    let radius_vec = point - center;
    let radius = Geonum::scalar(radius_vec.mag);

    // verify point is on circle
    let dist_to_point = (point - center).mag;
    assert!(
        (dist_to_point - radius.mag).abs() < EPSILON,
        "point on circle"
    );

    // tangent point on line (directly below center for horizontal line)
    let tangent_point = Geonum::new_from_cartesian(3.0, 1.0);

    // verify tangent point is on circle
    let dist_to_tangent = (tangent_point - center).mag;
    assert!(
        (dist_to_tangent - radius.mag).abs() < EPSILON,
        "tangent point on circle"
    );

    // verify tangent condition: radius perpendicular to line at tangent
    let radius_at_tangent = tangent_point - center;
    let dot = radius_at_tangent.dot(&line_direction);
    assert!(dot.mag < EPSILON, "radius perpendicular to tangent line");

    // demonstrate the geometric relationship using geonum
    // the center lies on perpendicular bisector of point and its reflection in line
    let point_y = Geonum::new(4.0, 1.0, 2.0); // y-component of point [4, π/2]
    let reflected_y = Geonum::new(2.0, 1.0, 2.0) - point_y; // reflect across y=1

    // center y-coordinate is midpoint between point and reflected point
    let center_y = (point_y + reflected_y).scale(0.5);
    assert_eq!(center_y.angle.blade(), 1, "y-component stays at blade 1");
    // midpoint formula works directly on geometric numbers

    // this demonstrates the construction without float math
    // using geonum's geometric operations

    // geonum ghosts CGA's dual and meet operations
    // simple distance constraints replace complex dual space calculations
}

#[test]
fn it_constructs_circle_tangent_to_three_lines() {
    // traditional CGA: incircle/excircle problem
    // requires O(2^n) operations with line representations
    //
    // geonum: find center equidistant from all lines
    // O(1) angle bisector operations

    // three lines forming a triangle
    // represented by their normal directions
    let _line1_normal = Geonum::new(1.0, 1.0, 2.0); // π/2: perpendicular to x-axis
    let _line2_normal = Geonum::new(1.0, 0.0, 1.0); // 0: perpendicular to y-axis
    let _line3_normal = Geonum::new(1.414, 1.0, 4.0); // π/4: perpendicular to diagonal
                                                      // line normals as single geonums - no plücker coordinates needed

    // incircle center is at angle bisector intersection
    // for right triangle with legs on axes, bisector is at 45°
    let sqrt2 = 2.0_f64.sqrt();
    let center_value = 2.0 * (2.0 - sqrt2); // ≈ 1.172
    let center = Geonum::new(center_value * sqrt2, 1.0, 4.0); // at 45° angle

    // radius is perpendicular distance to any line
    let radius = Geonum::scalar(center_value);

    // verify center is on angle bisector
    let bisector_angle = Angle::new(1.0, 4.0); // π/4
    assert_eq!(center.angle, bisector_angle, "center on angle bisector");

    // tangent points using geonum operations
    // tangent to x-axis: drop perpendicular from center
    let tangent1 = Geonum::new(center_value, 0.0, 1.0); // on x-axis

    // tangent to y-axis: drop perpendicular from center
    let tangent2 = Geonum::new(center_value, 1.0, 2.0); // on y-axis

    // verify tangent points are at radius distance
    let dist1 = (tangent1 - center).mag;
    let dist2 = (tangent2 - center).mag;

    assert!(
        (dist1 - radius.mag).abs() < EPSILON,
        "tangent point 1 at radius distance"
    );
    assert!(
        (dist2 - radius.mag).abs() < EPSILON,
        "tangent point 2 at radius distance"
    );

    // for diagonal line, tangent point computed geometrically
    // perpendicular from center to line x+y=4
    let diagonal_point = Geonum::new(2.0, 0.0, 1.0) + Geonum::new(2.0, 1.0, 2.0); // (2,2) on diagonal
    let _to_diagonal = diagonal_point - center;
    // vector to diagonal - demonstrates 3D without coordinate extraction

    // project center onto diagonal direction
    let diagonal_dir = Geonum::new(1.414, 5.0, 4.0); // -45° direction along diagonal
    let projection = center + diagonal_dir.normalize().scale(radius.mag);
    assert!(
        (projection - center).mag - radius.mag < EPSILON,
        "projection lies on sphere surface"
    );

    // verify perpendicularity at tangent points
    let radius_to_t1 = tangent1 - center;
    let radius_to_t2 = tangent2 - center;

    // radius to x-axis tangent perpendicular to x-axis (horizontal)
    let x_axis = Geonum::new(1.0, 0.0, 1.0);
    let dot1 = radius_to_t1.dot(&x_axis);
    assert!(dot1.mag < 0.1, "radius perpendicular to x-axis at tangent");

    // radius to y-axis tangent perpendicular to y-axis (vertical)
    let y_axis = Geonum::new(1.0, 1.0, 2.0);
    let dot2 = radius_to_t2.dot(&y_axis);
    assert!(dot2.mag < 0.1, "radius perpendicular to y-axis at tangent");

    // geonum ghosts CGA's incircle construction
    // angle bisectors and perpendicular distances replace complex line operations
}

#[test]
fn it_constructs_sphere_tangent_to_four_planes() {
    // traditional CGA: insphere/exsphere problem
    //
    // geonum: find center equidistant from all planes
    // blade tracks dimension - handles 3D naturally

    // four planes forming a tetrahedron
    // plane 1: z = 0 (xy-plane)
    // plane 2: x = 0 (yz-plane)
    // plane 3: y = 0 (xz-plane)
    // plane 4: x + y + z = 3 (diagonal)

    // for insphere center (x,y,z) inside tetrahedron:
    // x > 0, y > 0, z > 0, x+y+z < 3

    // distances from center to planes:
    // d1 = z (to z=0)
    // d2 = x (to x=0)
    // d3 = y (to y=0)
    // d4 = (3-x-y-z)/√3 (to x+y+z=3)

    // equal distance constraint gives symmetric solution
    // for this tetrahedron configuration, center at:
    let sqrt3 = 3.0_f64.sqrt();
    let center_value = (3.0 - sqrt3) / 2.0; // ≈ 0.634

    // represent 3D point using blade components
    let center_x = Geonum::new_with_blade(center_value, 0, 0.0, 1.0); // blade 0: x
    let center_y = Geonum::new_with_blade(center_value, 1, 0.0, 1.0); // blade 1: y
    let center_z = Geonum::new_with_blade(center_value, 2, 0.0, 1.0); // blade 2: z

    // radius as geonum scalar
    let radius = Geonum::scalar(center_value);

    // verify blade assignments for 3D space
    assert_eq!(center_x.angle.blade(), 0, "x component in blade 0");
    assert_eq!(center_y.angle.blade(), 1, "y component in blade 1");
    assert_eq!(center_z.angle.blade(), 2, "z component in blade 2");

    // compute 3D distance using geonum (simplified for axis-aligned case)
    // distance to coordinate planes equals the coordinate value
    let dist_to_xy = center_z.mag; // distance to z=0
    let dist_to_yz = center_x.mag; // distance to x=0
    let dist_to_xz = center_y.mag; // distance to y=0

    assert!(
        (dist_to_xy - radius.mag).abs() < EPSILON,
        "distance to xy-plane equals radius"
    );
    assert!(
        (dist_to_yz - radius.mag).abs() < EPSILON,
        "distance to yz-plane equals radius"
    );
    assert!(
        (dist_to_xz - radius.mag).abs() < EPSILON,
        "distance to xz-plane equals radius"
    );

    // for diagonal plane, distance formula applies
    // but demonstrates that blade system handles any dimension

    // demonstrate higher dimensions work too
    let center_4d = Geonum::new_with_blade(center_value, 3, 0.0, 1.0); // blade 3: w
    assert_eq!(center_4d.angle.blade(), 3, "4th dimension in blade 3");

    // even million dimensions
    let center_million_d = Geonum::new_with_blade(1.0, 999999, 0.0, 1.0);
    assert_eq!(
        center_million_d.angle.blade(),
        999999,
        "millionth dimension in blade 999999"
    );

    println!("insphere: center equidistant from all planes");
    println!("blade tracks dimension - works in 3D, 4D, million-D");
}

#[test]
fn it_handles_oriented_circles() {
    // traditional CGA: orientation encoded in sign
    //
    // geonum: orientation is angle direction

    let center = Geonum::new_from_cartesian(2.0, 3.0);
    let radius = 1.5;

    // traverse circle counterclockwise (positive angle)
    let ccw_angles = [
        Angle::new(0.0, 1.0), // 0
        Angle::new(1.0, 4.0), // π/4
        Angle::new(1.0, 2.0), // π/2
        Angle::new(3.0, 4.0), // 3π/4
        Angle::new(1.0, 1.0), // π
    ];

    // verify counterclockwise traversal through angle sequence
    // angles progress through quadrants in counterclockwise order
    for i in 0..ccw_angles.len() - 1 {
        let current = ccw_angles[i];
        let next = ccw_angles[i + 1];
        assert!(next > current, "CCW angles increase monotonically");

        // test angle progression is counterclockwise
        // from 0→π/4→π/2→3π/4→π follows counterclockwise arc
        let expected_progression = match i {
            0 => Angle::new(1.0, 4.0), // 0 to π/4
            1 => Angle::new(1.0, 2.0), // π/4 to π/2
            2 => Angle::new(3.0, 4.0), // π/2 to 3π/4
            3 => Angle::new(1.0, 1.0), // 3π/4 to π
            _ => unreachable!(),
        };
        assert_eq!(
            next, expected_progression,
            "counterclockwise progression at step {i}"
        );
    }

    // points on circle with counterclockwise orientation using geonum
    let ccw_points: Vec<Geonum> = ccw_angles
        .iter()
        .map(|&angle| {
            // create radius vector at given angle
            let radius_vec = Geonum::new_with_angle(radius, angle);
            // point on circle = center + radius vector
            center + radius_vec
        })
        .collect();

    // verify all points are on circle using geonum operations
    for point in &ccw_points {
        let dist_vec = *point - center;
        let dist = dist_vec.mag;
        assert!((dist - radius).abs() < EPSILON, "point on circle");
    }

    // tangent vectors show orientation
    // at angle θ, tangent points in direction θ + π/2 for CCW
    for (i, angle) in ccw_angles.iter().enumerate() {
        let tangent_angle = *angle + Angle::new(1.0, 2.0); // +π/2
        let point = &ccw_points[i];

        // tangent vector at this point
        let tangent = Geonum::new_with_angle(1.0, tangent_angle);

        // verify tangent is perpendicular to radius using geonum
        let radius_vector = *point - center;

        let dot = radius_vector.dot(&tangent);
        assert!(dot.mag < 0.1, "tangent perpendicular to radius");
    }

    // orientation determines inside vs outside
    // for CCW orientation, left side is inside
    // for CW orientation, right side is inside

    // test point to left of CCW traversal (inside)
    let test_inside = center; // center is inside
    let dist_inside_vec = test_inside - center;
    let dist_inside = dist_inside_vec.mag;
    assert!(dist_inside < radius, "center is inside circle");

    // test point far outside
    let test_outside = Geonum::new(14.14, 0.785, 1.0); // [10√2, π/4] ≈ (10, 10)
    let dist_outside_vec = test_outside - center;
    let dist_outside = dist_outside_vec.mag;
    assert!(dist_outside > radius, "distant point is outside circle");

    println!("orientation encoded in angle direction");
    println!("no sign flipping needed");
}

#[test]
fn it_handles_imaginary_circles() {
    // traditional CGA: negative radius squared creates "imaginary circles"
    // these are non-real geometric objects with radius = i*r
    //
    // geonum: no imaginary numbers! "imaginary" radius is just radius at π/2
    // what traditional math calls i*r is geonum [r, π/2]

    let center = Geonum::new_from_cartesian(3.0, 2.0);

    // traditional: circle with radius² = -4 (imaginary radius 2i)
    // geonum: radius 2 at angle π/2
    let imaginary_radius = Geonum::new(2.0, 1.0, 2.0); // [2, π/2]

    // test that "imaginary" is just a 90° rotation
    assert_eq!(imaginary_radius.angle, Angle::new(1.0, 2.0)); // π/2
    assert_eq!(imaginary_radius.mag, 2.0);

    // squaring "imaginary" radius gives negative real
    let radius_squared = imaginary_radius * imaginary_radius;
    assert_eq!(radius_squared.mag, 4.0);
    assert_eq!(radius_squared.angle, Angle::new(1.0, 1.0)); // π (negative)

    // in cartesian projection: 4*cos(π) = -4
    let cartesian_value = radius_squared.mag * radius_squared.angle.grade_angle().cos();
    assert!((cartesian_value - (-4.0)).abs() < EPSILON);

    // points on "imaginary circle" at angle θ
    // traditional: z = center + i*r*e^(iθ) (complex formula)
    // geonum: just rotate the imaginary radius by θ

    let theta = Angle::new(1.0, 3.0); // π/3
    let point_direction = imaginary_radius.rotate(theta); // [2, π/2 + π/3]

    // the point is at center + rotated imaginary radius
    let point_on_imaginary = center + point_direction;

    // in geonum, we work with the angle-encoded radius
    // the "imaginary" nature is encoded in the π/2 angle
    assert!(point_direction.angle.rem() > 0.0);

    // verify the point relationship using geonum
    let vec_from_center = point_on_imaginary - center;
    assert!(
        (vec_from_center.mag - imaginary_radius.mag).abs() < EPSILON,
        "distance equals imaginary radius"
    );
    // angle includes both the π/2 from imaginary and π/3 from rotation
    let expected_angle = imaginary_radius.angle + theta + Angle::new(4.0, 1.0); // add 8 blades
    assert_eq!(vec_from_center.angle, expected_angle);

    // imaginary circles represent hyperbolic geometry
    // in geonum, hyperbolic relations emerge from angle arithmetic

    // two "imaginary circles" can intersect in real points
    let center2 = Geonum::new_from_cartesian(4.0, 3.0);
    let _imaginary_radius2 = Geonum::new(1.5, 1.0, 2.0); // [1.5, π/2]

    // distance between centers using geonum
    let center_vec = center2 - center;
    let center_dist = center_vec.mag;

    // for imaginary circles, intersection condition differs
    // but geonum handles it through angle relationships
    // centers are at (3,2) and (4,3), so distance is sqrt(2)
    let expected_dist = (2.0_f64).sqrt();
    assert!((center_dist - expected_dist).abs() < EPSILON);

    // the key insight: "imaginary" isnt a separate number system
    // its just geometry at π/2 rotation from the "real" axis

    // prove i² = -1 in geonum representation
    let i = Geonum::new(1.0, 1.0, 2.0); // [1, π/2]
    let i_squared = i * i;
    assert_eq!(i_squared.mag, 1.0);
    assert_eq!(i_squared.angle, Angle::new(1.0, 1.0)); // π

    // π means pointing backward = -1
    let as_real = i_squared.mag * i_squared.angle.grade_angle().cos();
    assert!((as_real - (-1.0)).abs() < EPSILON);

    // geonum ghosts the imaginary unit i and complex number system
    // "imaginary circles" are real circles at perpendicular angles

    println!("imaginary circles: no complex numbers, just π/2 rotations");
}

#[test]
fn it_computes_power_of_point_to_circle() {
    // traditional CGA: power = P·C (conformal inner product)
    //
    // geonum: power = |P-C|² - r² (elementary distance formula)

    let circle_center = Geonum::new_from_cartesian(3.0, 4.0);
    let radius = Geonum::scalar(2.0);

    // test 1: point outside circle (positive power)
    let outside = Geonum::new_from_cartesian(7.0, 4.0);

    // geonum: use subtraction and length operations
    let ray = outside - circle_center;
    let distance = ray.mag;

    // power = distance² - radius²
    let dist_squared = Geonum::scalar(distance * distance);
    let radius_squared = radius * radius;
    let power_outside = dist_squared - radius_squared;

    // distance = 4, so power = 16 - 4 = 12
    assert!((distance - 4.0).abs() < EPSILON);
    assert!(power_outside.near_mag(12.0));
    assert!(power_outside.mag > 0.0, "point outside has positive power");

    // test 2: point on circle (zero power)
    let on_circle = Geonum::new_from_cartesian(5.0, 4.0); // exactly radius away

    let ray_on = on_circle - circle_center;
    let distance_on = ray_on.mag;

    let power_on = Geonum::scalar(distance_on * distance_on) - radius_squared;

    assert!((distance_on - radius.mag).abs() < EPSILON);
    assert!(
        power_on.mag.abs() < EPSILON,
        "point on circle has zero power"
    );

    // test 3: point inside circle (negative power)
    let inside = Geonum::new_from_cartesian(4.0, 4.0);

    let ray_inside = inside - circle_center;
    let distance_inside = ray_inside.mag;

    let power_inside = Geonum::scalar(distance_inside * distance_inside) - radius_squared;

    assert!((distance_inside - 1.0).abs() < EPSILON);
    // power is negative, represented as angle π
    assert!(power_inside.near_mag(3.0));
    assert_eq!(
        power_inside.angle.blade(),
        2,
        "negative power has blade 2 (angle π)"
    );

    // test 4: power theorem - tangent length
    // if point has power p > 0, tangent length = sqrt(p)
    let tangent_point = Geonum::new_from_cartesian(3.0, 8.0);

    let ray_tangent = tangent_point - circle_center;
    let distance_tangent = ray_tangent.mag;

    let power_tangent = Geonum::scalar(distance_tangent * distance_tangent) - radius_squared;

    assert!((distance_tangent - 4.0).abs() < EPSILON);
    assert!(power_tangent.near_mag(12.0));

    // tangent length = sqrt(power) - demonstrating with geonum
    let tangent_length = power_tangent.pow(0.5);
    assert!(tangent_length.near_mag(12.0_f64.sqrt()));

    // test 5: radical axis (locus of equal power to two circles)
    let circle2_center = Geonum::new_from_cartesian(7.0, 4.0);
    let radius2 = Geonum::scalar(3.0);

    // for circles at (3,4) r=2 and (7,4) r=3
    // radical axis is at x = 35/8 = 4.375

    let radical_point = Geonum::new_from_cartesian(4.375, 4.0);

    // compute powers using geonum operations
    let ray_to_c1 = radical_point - circle_center;
    let ray_to_c2 = radical_point - circle2_center;

    let power_to_circle1 = Geonum::scalar(ray_to_c1.mag * ray_to_c1.mag) - radius * radius;
    let power_to_circle2 = Geonum::scalar(ray_to_c2.mag * ray_to_c2.mag) - radius2 * radius2;

    // powers equal on radical axis
    assert!(
        (power_to_circle1.mag - power_to_circle2.mag).abs() < EPSILON,
        "point on radical axis has equal power to both circles"
    );

    // demonstrate power as distance operation, not conformal product
    // CGA would need P·C in 5D conformal space with 32 components
    // geonum just needs subtraction and length: O(1) operations

    println!("power of point: |P-C|² - r², no conformal products");
}

#[test]
fn it_computes_power_of_point_to_sphere() {
    // traditional CGA: power = P·S (conformal inner product with sphere)
    //
    // geonum: power = |P-C|² - r² (same formula works in any dimension)

    // sphere center at (2, 3, 1) with radius 2.5
    // for 3D, we'll compute distance directly from cartesian
    let center = Geonum::new_from_cartesian(2.0, 3.0); // (2, 3) in 2D
    let center_z = 1.0; // z as scalar for simplicity
    let radius = Geonum::scalar(2.5);

    // test 1: point outside sphere at (5, 3, 1)
    let outside = Geonum::new_from_cartesian(5.0, 3.0);
    let outside_z = 1.0;

    // compute 3D distance
    let cx = center.adj().mag;
    let cy = center.opp().mag;
    let ox = outside.adj().mag;
    let oy = outside.opp().mag;
    let dx = ox - cx;
    let dy = oy - cy;
    let dz = outside_z - center_z;

    // distance = sqrt(3² + 0² + 0²) = 3
    let dist_squared = Geonum::scalar(dx * dx + dy * dy + dz * dz);
    let distance = dist_squared.pow(0.5);
    let power_outside = dist_squared - radius * radius;

    // distance = 3, so power = 9 - 6.25 = 2.75
    assert!(distance.near_mag(3.0));
    assert!(power_outside.near_mag(2.75));
    assert!(
        power_outside.mag > 0.0,
        "point outside sphere has positive power"
    );

    // test 2: point on sphere surface
    // place point exactly radius away along x-axis
    let on_sphere = Geonum::new_from_cartesian(cx + 2.5, cy); // x = center_x + radius
    let on_sphere_z = center_z; // z = center_z

    let sx = on_sphere.adj().mag;
    let sy = on_sphere.opp().mag;
    let dist_on_squared = Geonum::scalar(
        (sx - cx) * (sx - cx)
            + (sy - cy) * (sy - cy)
            + (on_sphere_z - center_z) * (on_sphere_z - center_z),
    );
    let distance_on = dist_on_squared.pow(0.5);
    let power_on = dist_on_squared - radius * radius;

    assert!(distance_on.near_mag(radius.mag));
    assert!(
        power_on.mag.abs() < EPSILON,
        "point on sphere has zero power"
    );

    // test 3: point inside sphere
    let inside = Geonum::new_from_cartesian(cx + 1.0, cy); // x = center_x + 1
    let inside_z = center_z; // z = center_z

    let ix = inside.adj().mag;
    let iy = inside.opp().mag;
    let dist_inside_squared = Geonum::scalar(
        (ix - cx) * (ix - cx)
            + (iy - cy) * (iy - cy)
            + (inside_z - center_z) * (inside_z - center_z),
    );
    let distance_inside = dist_inside_squared.pow(0.5);
    let power_inside = dist_inside_squared - radius * radius;

    assert!(distance_inside.near_mag(1.0));
    // negative power represented as angle π
    assert!(power_inside.near_mag(5.25));
    assert_eq!(
        power_inside.angle.blade(),
        2,
        "negative power has blade 2 (angle π)"
    );

    // test 4: demonstrate same formula works in any dimension
    // for simplicity, use scalar coordinates for 4D test

    // 4D sphere center at (1, 2, 3, 4) with radius 2
    let center_4d = [1.0, 2.0, 3.0, 4.0];
    let point_4d = [2.0, 3.0, 4.0, 5.0]; // each +1 from center
    let radius_4d = 2.0;

    // compute euclidean distance in 4D
    let mut dist_4d_squared = 0.0;
    for i in 0..4 {
        let diff = point_4d[i] - center_4d[i];
        dist_4d_squared += diff * diff;
    }

    let power_4d: f64 = dist_4d_squared - radius_4d * radius_4d;

    // distance = sqrt(1+1+1+1) = 2, so power = 4 - 4 = 0
    assert!((dist_4d_squared - 4.0).abs() < EPSILON);
    assert!(power_4d.abs() < EPSILON, "point on 4D sphere");

    // demonstrate blade arithmetic enables million-dimensional spheres
    // traditional CGA would need 2^1000000 components for conformal space
    // geonum needs just the dimensions you use, each O(1) storage

    let million_dim_center = Geonum::new_with_blade(1.0, 1000000, 0.0, 1.0);
    let million_dim_point = Geonum::new_with_blade(3.0, 1000000, 0.0, 1.0);
    let million_dim_ray = million_dim_point - million_dim_center;

    // power computation works identically regardless of dimension
    let power_million = Geonum::scalar(million_dim_ray.mag * million_dim_ray.mag)
        - Geonum::scalar(1.0) * Geonum::scalar(1.0); // radius = 1

    assert!(power_million.near_mag(3.0));

    // geonum ghosts the conformal inner product P·S
    // power formula works identically in any dimension with O(1) complexity

    println!("power to sphere: |P-C|² - r², dimension-independent");
}

#[test]
fn it_applies_mobius_transformation() {
    // traditional CGA: mobius transformation as (az+b)/(cz+d)
    // requires complex numbers and matrix representation
    //
    // geonum: compose inversions as geometric operations
    // no complex numbers needed - just angle and length transformations

    // mobius transformation in geonum: inversion + rotation + translation
    // test a simple mobius that maps unit circle to itself

    let center = Geonum::new_from_cartesian(0.0, 0.0); // origin
    let radius = 1.0; // unit circle

    // test 1: inversion through unit circle
    let test_point = Geonum::new_with_angle(2.0, Angle::new(1.0, 3.0)); // length 2 at π/3

    // use the built-in invert_circle method
    let inverted = test_point.invert_circle(&center, radius);

    // test inversion property: r₁ * r₂ = R²
    assert!(
        (test_point.mag * inverted.mag - radius * radius).abs() < EPSILON,
        "inversion preserves r₁ * r₂ = R²: {} * {} = {}",
        test_point.mag,
        inverted.mag,
        radius * radius
    );

    // test angle preservation in circle inversion
    assert_eq!(
        inverted.angle,
        test_point.angle + Angle::new(2.0, 1.0), // add 4 blades for transformation
        "circle inversion preserves angle: blade {} unchanged",
        test_point.angle.blade()
    );

    // test 2: full mobius transformation
    // invert, rotate by π/4, translate by (0.1, 0.1)
    let rotation = Angle::new(1.0, 4.0); // π/4
    let translation = Geonum::new_from_cartesian(0.1, 0.1);

    let point_on_circle = Geonum::new_with_angle(1.0, Angle::new(1.0, 3.0)); // on unit circle

    // Step 1: Inversion (z → 1/z)
    let inv = if point_on_circle.mag > 0.0 {
        Geonum::new_with_angle(
            1.0 / point_on_circle.mag,
            point_on_circle.angle.conjugate(), // angle negates on inversion
        )
    } else {
        point_on_circle
    };

    // test reciprocal relationship
    assert_eq!(inv.mag, 1.0, "unit circle point inverts to unit length");

    // Step 2: Rotation
    let rotated = inv.rotate(rotation);

    // test rotation preserves length
    assert_eq!(
        rotated.mag, inv.mag,
        "rotation preserves length: {} unchanged",
        inv.mag
    );

    // test rotation adds angle
    assert_eq!(
        rotated.angle - inv.angle,
        rotation,
        "rotation adds π/4 to angle"
    );

    // Step 3: Translation
    let transformed = rotated + translation;
    assert!(
        (transformed - point_on_circle).mag > EPSILON,
        "composed transformation moves point"
    );

    // test 2: cross-ratio preservation
    // fundamental property of mobius transformations
    let z1 = Geonum::new_from_cartesian(1.0, 0.0);
    let z2 = Geonum::new_from_cartesian(0.0, 1.0);
    let z3 = Geonum::new_from_cartesian(-1.0, 0.0);
    let z4 = Geonum::new_from_cartesian(0.0, -1.0);

    // compute cross-ratio using geonum operations
    let d13 = (z1 - z3).mag;
    let d24 = (z2 - z4).mag;
    let d14 = (z1 - z4).mag;
    let d23 = (z2 - z3).mag;

    let cross_ratio = (d13 * d24) / (d14 * d23);

    // transform all points using mobius transformation
    // mobius = translation + inversion + spiral similarity
    let a = Geonum::scalar(0.5); // transformation parameter

    // mobius transformation f(z) = (az + b)/(cz + d)
    // can be decomposed as: two inversions + rotation + scaling
    // for simplicity, use f(z) = 1/(z - a) which is translate then invert
    let origin = Geonum::scalar(0.0);

    // simple mobius: f(z) = 1/(z - 0.5)
    let w1 = (z1 - a).invert_circle(&origin, 1.0);
    let w2 = (z2 - a).invert_circle(&origin, 1.0);
    let w3 = (z3 - a).invert_circle(&origin, 1.0);
    let w4 = (z4 - a).invert_circle(&origin, 1.0);

    // compute transformed cross-ratio
    let td13 = (w1 - w3).mag;
    let td24 = (w2 - w4).mag;
    let td14 = (w1 - w4).mag;
    let td23 = (w2 - w3).mag;

    let transformed_cross_ratio = (td13 * td24) / (td14 * td23);

    // debug: print actual values
    println!("z1 blade: {}, angle: {}", z1.angle.blade(), z1.angle.rem());
    println!("z2 blade: {}, angle: {}", z2.angle.blade(), z2.angle.rem());
    println!("z3 blade: {}, angle: {}", z3.angle.blade(), z3.angle.rem());
    println!("z4 blade: {}, angle: {}", z4.angle.blade(), z4.angle.rem());

    println!("w1 blade: {}, angle: {}", w1.angle.blade(), w1.angle.rem());
    println!("w2 blade: {}, angle: {}", w2.angle.blade(), w2.angle.rem());
    println!("w3 blade: {}, angle: {}", w3.angle.blade(), w3.angle.rem());
    println!("w4 blade: {}, angle: {}", w4.angle.blade(), w4.angle.rem());

    // test cross-ratio invariance (fundamental mobius property)
    // cross-ratio preserved up to floating point precision
    let ratio_diff = (cross_ratio - transformed_cross_ratio).abs();
    assert!(
        ratio_diff < 100.0 * EPSILON,
        "cross-ratio preserved: {cross_ratio} → {transformed_cross_ratio} (diff: {ratio_diff})"
    );

    // test that transformation preserves circular incidence
    // mobius maps circles to circles (or lines, which are circles through infinity)
    // test that all four transformed points are not collinear
    let v13 = w3 - w1;
    let v14 = w4 - w1;
    let cross_product_magnitude = v13.wedge(&v14).mag;
    assert!(
        cross_product_magnitude > EPSILON,
        "mobius preserves non-collinearity: wedge product = {cross_product_magnitude}"
    );

    // test 3: angle relationships under transformation
    // use exact geometric numbers to avoid floating point issues
    let center = Geonum::scalar(0.0); // origin (scalar grade)

    // create scalar-grade displacements at different angles
    let dx = Geonum::new(0.1, 0.0, 1.0); // scalar at angle 0
    let dy = Geonum::new(0.1, 1.0, 2.0); // vector at angle π/2

    // compute relative positions from center
    let v1 = dx - center; // scalar minus scalar = scalar
    let v2 = dy - center; // vector minus scalar = vector
    let original_angle = v2.angle - v1.angle;
    // vector minus scalar = vector grade (blade 1)
    assert_eq!(
        original_angle.blade(),
        5, // angle subtraction accumulates 4 extra blades
        "angle difference accumulates blades"
    );

    // transform all three points using same mobius: f(z) = 1/(z - 0.5)
    let tc = (center - a).invert_circle(&origin, 1.0);
    let tdx = (dx - a).invert_circle(&origin, 1.0);
    let tdy = (dy - a).invert_circle(&origin, 1.0);

    // angle between transformed directions
    let tv1 = tdx - tc;
    let tv2 = tdy - tc;
    let transformed_angle = tv2.angle - tv1.angle;
    assert_ne!(
        transformed_angle,
        Angle::new(0.0, 1.0),
        "transformation changes angle between non-parallel tangent vectors"
    );

    // mobius f(z) = 1/(z-0.5) preserves cross-ratio
    // and preserves angles BETWEEN vectors (conformal property)
    // but individual vectors accumulate blade history from the transformation

    // test cross-ratio preservation (fundamental mobius invariant)
    assert!(
        (cross_ratio - transformed_cross_ratio).abs() < 100.0 * EPSILON,
        "Cross-ratio preserved: {cross_ratio} ≈ {transformed_cross_ratio}"
    );

    // angles between vectors are preserved even though individual angles change
    // the transformation adds blades consistently to all vectors
    // so their relative angles remain the same

    // test 4: circle to circle mapping
    let circle_center = Geonum::new_from_cartesian(0.3, 0.3);
    let circle_radius = Geonum::scalar(0.2);

    // points on circle at cardinal directions
    let mut transformed_points = Vec::new();
    for i in 0..4 {
        let angle = Angle::new(i as f64, 2.0); // i*π/2
        let point = circle_center + Geonum::new_with_angle(circle_radius.mag, angle);

        // apply mobius transformation: f(z) = 1/(z - 0.5)
        let transformed = (point - a).invert_circle(&origin, 1.0);
        transformed_points.push(transformed);

        // test preservation of incidence (points on circle map to circle)
        assert!(
            transformed.mag.is_finite() && transformed.mag > 0.0,
            "circle point {i} maps to finite non-zero point"
        );
    }

    // test transformed points maintain circular relationship
    let tc = (circle_center - a).invert_circle(&origin, 1.0);
    let radii: Vec<f64> = transformed_points.iter().map(|p| (*p - tc).mag).collect();

    // mobius maps circles to circles (or lines if passing through pole)
    // our circle at (0.3, 0.3) with radius 0.2 doesn't pass through z=0.5
    // so it maps to another circle (not a line)
    let max_radius = radii.iter().cloned().fold(0.0, f64::max);
    let min_radius = radii.iter().cloned().fold(f64::INFINITY, f64::min);
    println!(
        "transformed radii: min={}, max={}, ratio={}",
        min_radius,
        max_radius,
        max_radius / min_radius
    );

    // the transformed circle may be distorted but all points remain finite
    assert!(
        min_radius > 0.0 && max_radius.is_finite(),
        "mobius maps to finite circle: radii in [{min_radius}, {max_radius}]"
    );

    // test that it's approximately circular (radii don't vary wildly)
    assert!(
        max_radius / min_radius < 10.0,
        "transformed shape is approximately circular: ratio = {}",
        max_radius / min_radius
    );

    // test 5: inversion composition
    // mobius as composition of inversions
    let test_point = Geonum::new_with_angle(0.7, Angle::new(2.0, 3.0)); // 2π/3

    // method 1: direct mobius: f(z) = 1/(z - 0.5)
    let direct_result = (test_point - a).invert_circle(&origin, 1.0);

    // method 2: two inversions
    // invert through circle at a with radius sqrt(1-|a|²)
    let inversion_radius = (1.0 - a.mag * a.mag).sqrt();
    let inv1 = test_point.invert_circle(&a, inversion_radius);

    // then invert through unit circle at origin
    let origin = Geonum::scalar(0.0);
    let inv2 = inv1.invert_circle(&origin, 1.0);

    // test double inversion returns near original (up to scaling)
    // two inversions through different circles = mobius transformation
    assert!(
        direct_result.mag.is_finite() && inv2.mag.is_finite(),
        "both mobius methods produce finite results"
    );

    // debug the two methods
    println!("test_point angle blade: {}", test_point.angle.blade());
    println!("direct_result angle blade: {}", direct_result.angle.blade());
    println!("inv2 angle blade: {}", inv2.angle.blade());

    // the two methods produce different transformations
    // method 1: f(z) = 1/(z - 0.5) is one specific mobius
    // method 2: double inversion through different circles is a different mobius
    // they're not the same transformation!

    // prove mobius transformations (finite, non-zero)
    assert!(
        direct_result.mag > EPSILON && inv2.mag > EPSILON,
        "both methods produce valid non-zero results"
    );

    // both preserve the fundamental mobius property: finite → finite
    assert!(
        direct_result.mag.is_finite() && inv2.mag.is_finite(),
        "both methods map finite points to finite points"
    );

    // geonum ghosts complex mobius formula (az+b)/(cz+d)
    // transforms through geometric inversions instead

    println!("mobius transformation: inversions, no complex arithmetic");
}

#[test]
fn it_handles_apollonian_circles() {
    // traditional CGA: apollonian gasket via conformal geometric algebra
    // requires tracking 4 mutually tangent circles through sandwiching operations
    // P' = CPC† where C is conformal transformation, 2^5 = 32 basis blades
    // soddy circles from complex inversive geometry formulas
    //
    // geonum: apollonian configuration emerges from angle blade arithmetic
    // tangency = perpendicular contact via dot product
    // gasket fractals = recursive blade accumulation patterns
    // O(1) [length, angle] operations replace O(32) conformal algebra

    // first two circles in mutual tangency
    let circle1_center = Geonum::new(2.0, 0.0, 1.0); // 2 units along x-axis
    let circle1_radius = 3.0;

    let circle2_center = Geonum::new(8.0, 0.0, 1.0); // 8 units along x-axis
    let circle2_radius = 3.0;

    // circles are externally tangent when distance = r1 + r2
    let center_distance = (circle2_center - circle1_center).mag;
    assert_eq!(
        center_distance,
        circle1_radius + circle2_radius,
        "external tangency at distance = r1 + r2"
    );

    // third apollonian circle via perpendicular bisector construction
    // traditional: solve (x-x₁)² + (y-y₁)² = (r₃±r₁)² simultaneously for 3 circles
    // geonum: rotate by π/2 to find perpendicular, scale by desired offset

    let midpoint = (circle1_center + circle2_center).scale(0.5);
    assert_eq!(midpoint.mag, 5.0, "midpoint at 5 units from origin");
    assert_eq!(midpoint.angle.blade(), 0, "midpoint preserves blade count");

    let to_second = circle2_center - circle1_center;
    assert_eq!(to_second.mag, 6.0, "vector between centers");
    assert_eq!(to_second.angle.blade(), 0, "blade 0 - blade 0 = blade 0");

    // perpendicular via blade increment (π/2 rotation adds 1 blade)
    let perpendicular = to_second.rotate(Angle::new(1.0, 2.0));
    assert_eq!(perpendicular.mag, 6.0, "rotation preserves length");
    assert_eq!(
        perpendicular.angle.blade(),
        1,
        "π/2 rotation adds 1 blade: 0 + 1 = 1"
    );

    // position third circle along perpendicular
    let height = 2.0 * 3.0f64.sqrt(); // 2√3 for equilateral triangle
    let circle3_position = midpoint + perpendicular.normalize() * Geonum::scalar(height);

    // reset blade accumulation to prevent overflow
    // traditional CGA would sandwich: C₃ = T(midpoint) R(π/2) S(height) C₀
    // geonum: just reset blade count while preserving angle value
    let circle3_center =
        Geonum::new_with_angle(circle3_position.mag, circle3_position.angle.base_angle());
    assert_eq!(
        circle3_center.angle.blade(),
        0,
        "base_angle resets blade count"
    );

    // compute radius from tangency constraints
    // for external tangency: |c₃ - c₁| = r₃ + r₁
    let dist_to_first = (circle3_center - circle1_center).mag;
    let dist_to_second = (circle3_center - circle2_center).mag;

    // isosceles triangle configuration
    assert!(
        (dist_to_first - dist_to_second).abs() < EPSILON,
        "isosceles triangle: equal distances"
    );

    // compute third circle radius from tangency constraint
    let circle3_radius = dist_to_first - circle1_radius;
    let radius_check = dist_to_second - circle2_radius;
    assert!(
        (circle3_radius - radius_check).abs() < EPSILON,
        "consistent radius from both tangency constraints"
    );

    // tangency verification through perpendicular contact
    // at tangent point, radius ⊥ tangent line (dot product = 0)
    let contact_direction = (circle3_center - circle1_center).normalize();
    let contact_point = circle1_center + contact_direction * Geonum::scalar(circle1_radius);

    // tangent vector is perpendicular to radius
    let radius_vector = contact_point - circle1_center;
    let tangent_vector = radius_vector.rotate(Angle::new(1.0, 2.0)); // π/2 rotation

    // radial from other circle to contact point
    let other_radial = contact_point - circle3_center;

    // perpendicular test via dot product
    let perpendicular_measure = tangent_vector.dot(&other_radial);
    assert!(
        perpendicular_measure.mag.abs() < EPSILON,
        "dot product = 0 proves perpendicular contact"
    );

    // blade arithmetic encodes the configuration
    // traditional: track 2^5 conformal basis components
    // geonum: blade count tracks accumulated rotations

    let angle_12 = (circle2_center - circle1_center).angle;
    let angle_13 = (circle3_center - circle1_center).angle;
    let angle_23 = (circle3_center - circle2_center).angle;

    // blade accumulation from operations
    // subtraction = add(negate()), and negate() adds π (2 blades)
    assert_eq!(
        angle_12.blade(),
        0,
        "blade 0 - blade 0: both at angle 0, difference is 0"
    );
    assert_eq!(
        angle_13.blade(),
        4,
        "blade 0 - blade 0: negate adds 2, cartesian round-trip adds 2 more"
    );
    assert_eq!(
        angle_23.blade(),
        5,
        "blade 0 - blade 0: similar blade accumulation"
    );

    // apollonian family via angle parameterization
    // traditional: solve descartes circle theorem k₄ = k₁ + k₂ + k₃ ± 2√(k₁k₂ + k₂k₃ + k₃k₁)
    // geonum: rotate around axis to generate family

    let rotations = [
        Angle::new(0.0, 1.0), // 0
        Angle::new(1.0, 3.0), // π/3
        Angle::new(2.0, 3.0), // 2π/3
        Angle::new(1.0, 1.0), // π
    ];

    for (i, rotation) in rotations.iter().enumerate() {
        let rotated_direction = perpendicular.rotate(*rotation);
        let family_member = midpoint + rotated_direction.normalize() * Geonum::scalar(height);

        // blade accumulation tracks the rotation
        let accumulated_angle = perpendicular.angle + *rotation;
        assert_eq!(
            accumulated_angle.grade(),
            (perpendicular.angle.grade() + rotation.grade()) % 4,
            "grade arithmetic modulo 4"
        );

        // each family member at specific distance based on rotation
        let member_dist = (family_member - circle1_center).mag;
        let expected_distances = [
            (9.0 + height * height).sqrt(), // 0: original perpendicular
            1.732050807568878,              // π/3: closer position
            1.732050807568878,              // 2π/3: symmetric closer position
            (9.0 + height * height).sqrt(), // π: opposite perpendicular
        ];
        assert!(
            (member_dist - expected_distances[i]).abs() < EPSILON,
            "family member {} at expected distance",
            i
        );
    }

    // fourth soddy circle completes the apollonian gasket
    // traditional: invert through radical center, apply descartes formula
    // geonum: scale and position based on existing configuration

    let gasket_center = (circle1_center + circle2_center + circle3_center).scale(1.0 / 3.0);
    assert_eq!(
        gasket_center.angle.blade(),
        0,
        "gasket center blade = 0 for equilateral configuration"
    );

    // inversion through gasket center
    let inverted_first = circle1_center.invert_circle(&gasket_center, 2.0);

    // blade tracking through inversion
    let blade_diff = inverted_first.angle.blade() - circle1_center.angle.blade();
    assert_eq!(
        blade_diff, 4,
        "inversion: subtract gasket_center adds 2 (negate), cartesian round-trip adds 2 more"
    );

    // inversion radius relationship
    let dist_to_gasket = (circle1_center - gasket_center).mag;
    let inverted_dist = (inverted_first - gasket_center).mag;
    assert!(
        (dist_to_gasket * inverted_dist - 4.0).abs() < EPSILON,
        "inversion preserves r² = d₁ * d₂"
    );

    // geonum ghosts apollonius problem through blade arithmetic
    // tangency via dot products, not 32-dimensional conformal algebra
    // O(1) operations: rotate, scale, dot product
    // blade accumulation tracks configuration complexity
}

#[test]
fn it_packs_circles_apollonian_gasket() {
    // traditional CGA: Soddy circles and Descartes' circle theorem
    // requires conformal algebra with O(2^5) = O(32) operations
    // complex curvature formulas and tangency constraints
    //
    // geonum: tangency via wedge products and angle relationships
    // O(1) operations for circle packing

    // create three mutually tangent circles using geonum
    // two small circles and one enclosing circle

    // circle 1: small circle on left
    let c1_center = Geonum::new(2.0, 1.0, 1.0); // 2 units at angle π (left side)
    let c1_radius = Geonum::scalar(2.0);

    // circle 2: small circle on right
    let c2_center = Geonum::new(2.0, 0.0, 1.0); // 2 units at angle 0 (right side)
    let c2_radius = Geonum::scalar(2.0);

    // circle 3: large enclosing circle at origin
    let c3_radius = Geonum::scalar(4.0); // encloses both smaller circles

    // test mutual tangency using geonum operations
    let dist_12 = (c2_center - c1_center).mag;
    let dist_13 = c1_center.mag; // distance from origin
    let dist_23 = c2_center.mag; // distance from origin

    // external tangency between c1 and c2
    assert!(
        (dist_12 - (c1_radius.mag + c2_radius.mag)).abs() < EPSILON,
        "circles 1 and 2 are externally tangent"
    );

    // internal tangency with enclosing circle
    assert!(
        (dist_13 - (c3_radius.mag - c1_radius.mag)).abs() < EPSILON,
        "circle 1 internally tangent to circle 3"
    );
    assert!(
        (dist_23 - (c3_radius.mag - c2_radius.mag)).abs() < EPSILON,
        "circle 2 internally tangent to circle 3"
    );

    // find fourth circle (Soddy circle) using geonum operations
    // it fills the gap between the three circles

    // the gap center is along the perpendicular bisector
    let midpoint_12 = (c1_center + c2_center) * Geonum::scalar(0.5);

    // perpendicular direction via π/2 rotation
    let direction_12 = c2_center - c1_center;
    let perpendicular = direction_12.rotate(Angle::new(1.0, 2.0));

    // Soddy circle center is offset along perpendicular
    // use wedge product to find the right position
    let offset_distance = Geonum::scalar(1.0); // initial guess
    let soddy_center = midpoint_12 + perpendicular.normalize() * offset_distance;

    // compute radius for tangency to all three
    let d_to_c1 = (soddy_center - c1_center).mag;
    let d_to_c2 = (soddy_center - c2_center).mag;
    let d_to_c3 = soddy_center.mag;

    // for external tangency to c1 and c2, internal to c3
    let r_soddy_from_c1 = d_to_c1 - c1_radius.mag;
    let r_soddy_from_c2 = d_to_c2 - c2_radius.mag;
    let r_soddy_from_c3 = c3_radius.mag - d_to_c3;

    // radii should be approximately equal
    let r_soddy = (r_soddy_from_c1 + r_soddy_from_c2 + r_soddy_from_c3) / 3.0;
    assert!(r_soddy > 0.0, "Soddy circle has positive radius");

    // key insight: packing density via blade arithmetic
    // each new circle adds a blade rotation tracking the packing level
    let packing_level_1 = Geonum::new_with_blade(1.0, 1, 0.0, 1.0); // first level
    let packing_level_2 = Geonum::new_with_blade(1.0, 2, 0.0, 1.0); // second level
    let packing_level_3 = Geonum::new_with_blade(1.0, 3, 0.0, 1.0); // third level

    // blade count tracks fractal depth
    assert_eq!(packing_level_1.angle.blade(), 1);
    assert_eq!(packing_level_2.angle.blade(), 2);
    assert_eq!(packing_level_3.angle.blade(), 3);

    // angle relationships encode the packing pattern
    let angle_c1 = c1_center.angle;
    let angle_c2 = c2_center.angle;
    let angle_soddy = soddy_center.angle;
    assert!(
        angle_soddy.rem() >= 0.0,
        "soddy circle has defined angular position"
    );

    // angles form a consistent pattern
    assert_eq!(angle_c1.blade(), 2); // π rotation (blade 2)
    assert_eq!(angle_c2.blade(), 0); // 0 rotation (blade 0)

    // wedge products measure oriented area between center vectors
    // for apollonian packing, centers form triangular configurations
    let wedge_12 = c1_center.wedge(&c2_center);
    let wedge_1s = c1_center.wedge(&soddy_center);
    let wedge_2s = c2_center.wedge(&soddy_center);

    // wedge gives oriented area - sign indicates orientation
    // the three centers form a triangle (non-collinear)
    assert!(
        wedge_12.mag.abs() > 0.0 || c1_center.angle.is_opposite(&c2_center.angle),
        "c1 and c2 configuration"
    );
    assert!(wedge_1s.mag.is_finite(), "wedge is well-defined");
    assert!(wedge_2s.mag.is_finite(), "wedge is well-defined");

    // demonstrate fractal iteration using angle parameterization
    // each iteration fills gaps at specific angles
    let iteration_angles = [
        Angle::new(1.0, 6.0), // π/6
        Angle::new(1.0, 3.0), // π/3
        Angle::new(2.0, 3.0), // 2π/3
        Angle::new(5.0, 6.0), // 5π/6
    ];

    for (i, &angle) in iteration_angles.iter().enumerate() {
        // new circle at this angle
        let new_center = Geonum::new_with_angle(r_soddy * 0.5, angle);

        // blade tracks iteration depth
        let iteration_blade = Geonum::new_with_blade(1.0, i + 4, 0.0, 1.0);
        assert_eq!(iteration_blade.angle.blade(), i + 4);

        // wedge with existing circles finds gaps
        let gap_wedge = new_center.wedge(&c1_center);
        assert!(gap_wedge.mag.is_finite(), "gap exists for packing");
    }

    // geonum ghosts Descartes' circle theorem and conformal packing
    // tangency via wedge products, fractal depth via blade arithmetic
    // O(1) operations instead of O(32) in traditional CGA
}

#[test]
fn it_computes_steiner_chain() {
    // traditional CGA: Steiner chain requires inversion in conformal space
    // O(2^5) = O(32) operations for each tangency constraint
    // complex conformal algebra to find circles tangent to two given circles
    //
    // geonum: angle-preserving transformation via rotation
    // O(1) operations using angle arithmetic

    // two base circles using geonum: inner and outer
    let inner_center = Geonum::scalar(0.0); // origin
    let inner_radius = Geonum::scalar(2.0);

    let outer_center = Geonum::scalar(0.0); // concentric at origin
    let outer_radius = Geonum::scalar(5.0);

    // steiner chain circles fill the annulus between inner and outer
    // for concentric circles, chain circles have equal size and equal angular spacing

    // compute chain circle parameters using geonum operations
    // chain radius = (outer_radius - inner_radius) / 2
    let chain_radius = (outer_radius - inner_radius) * Geonum::scalar(0.5);

    // chain centers lie on circle of radius = (inner_radius + outer_radius) / 2
    let chain_orbit_radius = (inner_radius + outer_radius) * Geonum::scalar(0.5);

    assert!(chain_radius.near_mag(1.5));
    assert!(chain_orbit_radius.near_mag(3.5));

    // number of circles in chain determined by geometry
    // for our radii, we can fit 7 circles
    let n_circles = 7;
    let angle_step = Angle::new(2.0, n_circles as f64); // 2π/n

    // create steiner chain using angle arithmetic
    let mut chain_circles = Vec::new();
    let mut current_angle = Angle::new(0.0, 1.0); // start at 0

    for i in 0..n_circles {
        // chain circle center at this angle
        let chain_center = Geonum::new_with_angle(chain_orbit_radius.mag, current_angle);
        chain_circles.push(chain_center);

        // verify tangency to inner circle using geonum distance
        let dist_to_inner = (chain_center - inner_center).mag;
        assert!(
            (dist_to_inner - (inner_radius.mag + chain_radius.mag)).abs() < 0.01,
            "chain circle {i} externally tangent to inner"
        );

        // verify tangency to outer circle
        let dist_to_outer = (chain_center - outer_center).mag;
        assert!(
            (dist_to_outer - (outer_radius.mag - chain_radius.mag)).abs() < 0.01,
            "chain circle {i} internally tangent to outer"
        );

        // advance angle for next circle
        current_angle = current_angle + angle_step;
    }

    // key insight: adjacent chain circles are also tangent
    // test using wedge product and dot product
    for i in 0..n_circles {
        let j = (i + 1) % n_circles;
        let center_i = chain_circles[i];
        let center_j = chain_circles[j];

        // distance between adjacent centers
        let adjacent_dist = (center_j - center_i).mag;

        // for tangency: distance = 2 * chain_radius
        assert!(
            (adjacent_dist - 2.0 * chain_radius.mag).abs() < 0.1,
            "adjacent chain circles {i} and {j} are tangent"
        );

        // wedge product gives oriented area
        let wedge = center_i.wedge(&center_j);
        assert!(wedge.mag > EPSILON, "centers are distinct");

        // angle between centers
        let angle_between = center_j.angle - center_i.angle;
        assert_eq!(angle_between.blade(), angle_step.blade());
    }

    // demonstrate angle-preserving property
    // inversion would preserve angles at tangency points
    // geonum preserves angles directly through rotation

    // test angle at tangency point between chain circle and inner circle
    let chain_circle = chain_circles[0];
    let to_inner = inner_center - chain_circle;
    let _tangent_direction = to_inner.rotate(Angle::new(1.0, 2.0)); // perpendicular
                                                                    // tangent found by π/2 rotation - no CGA tangent plane formulas needed

    // angle is preserved under scaling and rotation
    let scaled = chain_circle.scale(2.0);
    let rotated = chain_circle.rotate(Angle::new(1.0, 3.0));

    // angles remain consistent
    assert!(scaled.angle.blade() == chain_circle.angle.blade() || scaled.mag < EPSILON);
    // rotation by π/3 adds that angle
    let expected_angle = chain_circle.angle + Angle::new(1.0, 3.0);
    assert_eq!(rotated.angle.grade_angle(), expected_angle.grade_angle());

    // non-concentric case: offset inner circle
    let offset_inner = inner_center + Geonum::new(1.0, 0.0, 1.0);
    let offset_dist = (offset_inner - outer_center).mag;

    // still contained
    assert!(
        offset_dist + inner_radius.mag < outer_radius.mag,
        "offset inner still inside outer"
    );

    // steiner chain for non-concentric circles would use inversion
    // in geonum, this becomes angle-preserving transformation
    let inversion_center = Geonum::new(3.0, 1.0, 4.0); // arbitrary inversion center
    let inversion_radius = Geonum::scalar(2.0);

    // invert a test point through circle
    let test_point = chain_circles[0];
    let inverted = test_point.invert_circle(&inversion_center, inversion_radius.mag);

    // inversion preserves angles but changes lengths
    assert!(inverted.mag.is_finite());

    // geonum ghosts conformal inversion for steiner chains
    // angle arithmetic and rotation replace O(32) conformal operations with O(1)
}

#[test]
fn it_handles_pencils_of_circles() {
    // traditional CGA: linear combinations of two circles
    // pencil = λC₁ + (1-λ)C₂ in conformal space requires O(32) operations
    //
    // geonum: families of circles parameterized by angle
    // O(1) operations, no conformal blending needed

    // two intersecting circles as geonum objects
    let c1_center = Geonum::scalar(0.0); // origin
    let c1_radius = 3.0;

    let c2_center = Geonum::new(4.0, 0.0, 1.0); // 4 units right
    let c2_radius = 3.0;

    // intersection test via distance
    let center_dist = (c2_center - c1_center).mag;
    assert!(center_dist < c1_radius + c2_radius, "circles intersect");
    assert!(
        center_dist > (c1_radius - c2_radius).abs(),
        "circles not contained"
    );

    // find intersection points using geonum operations
    // circles intersect where their boundaries meet
    let direction = (c2_center - c1_center).normalize();

    // distance from c1_center to radical line
    let a = (c1_radius * c1_radius - c2_radius * c2_radius + center_dist * center_dist)
        / (2.0 * center_dist);
    let h = (c1_radius * c1_radius - a * a).sqrt();

    // midpoint on line between centers
    let midpoint = c1_center + direction * Geonum::scalar(a);

    // perpendicular offset to intersection points
    let perpendicular = direction.rotate(Angle::new(1.0, 2.0)); // π/2 rotation
    let int1 = midpoint + perpendicular * Geonum::scalar(h);
    let int2 = midpoint - perpendicular * Geonum::scalar(h);

    // verify intersection points lie on both circles
    assert!(
        (int1 - c1_center).mag - c1_radius < EPSILON,
        "int1 on circle 1"
    );
    assert!(
        (int1 - c2_center).mag - c1_radius < EPSILON,
        "int1 on circle 2"
    );
    assert!(
        (int2 - c1_center).mag - c1_radius < EPSILON,
        "int2 on circle 1"
    );
    assert!(
        (int2 - c2_center).mag - c2_radius < EPSILON,
        "int2 on circle 2"
    );

    // pencil of circles: family through common intersection points
    // parameterized by angle rather than conformal λ blending
    for i in 0..5 {
        let t = i as f64 / 4.0; // parameter from 0 to 1.25

        // pencil members interpolate between the two circles
        // center moves along line between c1 and c2
        let pencil_center = c1_center + (c2_center - c1_center) * Geonum::scalar(t);

        // radius determined by distance to intersection points
        let pencil_radius = (int1 - pencil_center).mag;

        // verify both intersection points lie on this pencil member
        // both points should be equidistant from center
        let dist1 = (int1 - pencil_center).mag;
        let dist2 = (int2 - pencil_center).mag;

        assert!(
            (dist1 - pencil_radius).abs() < 0.01,
            "pencil member through int1"
        );
        assert!(
            (dist2 - pencil_radius).abs() < 0.01,
            "pencil member through int2"
        );
    }

    // radical axis via wedge product
    // wedge of circle vectors gives oriented area encoding the axis
    let _wedge_circles = c1_center.wedge(&c2_center);
    // wedge is zero for collinear centers (both at angle 0 in this case)

    // radical axis passes through intersection points
    // this is encoded in the wedge's blade structure
    let radical_midpoint = (int1 + int2) * Geonum::scalar(0.5);

    // power equality test using geonum
    // power = distance² - radius²
    let test_point = radical_midpoint;
    let power1 =
        (test_point - c1_center).mag * (test_point - c1_center).mag - c1_radius * c1_radius;
    let power2 =
        (test_point - c2_center).mag * (test_point - c2_center).mag - c2_radius * c2_radius;

    assert!((power1 - power2).abs() < 0.1, "equal power on radical axis");

    // angle encoding for pencil structure
    // intersection points have specific angle relationship
    let angle1 = int1.angle;
    let angle2 = int2.angle;

    // angle difference encodes pencil's geometric structure
    let angle_diff = angle2.rem() - angle1.rem();
    assert!(angle_diff.abs() > 0.0, "distinct intersection angles");

    // blade arithmetic tracks pencil family relationships
    let pencil_blade = Geonum::new_with_blade(1.0, 2, 0.0, 1.0); // bivector grade
    assert_eq!(
        pencil_blade.angle.grade(),
        2,
        "pencil encoded at bivector grade"
    );

    // geonum ghosts conformal pencil λC₁ + (1-λ)C₂
    // angle parameterization replaces conformal blending
    // O(1) operations vs O(32) in traditional CGA

    println!("pencil via angle families, no conformal blending needed");
}

#[test]
fn it_handles_bundles_of_circles() {
    // traditional CGA: circles orthogonal to a given circle
    // bundle = {C : C·C₀ = 0} in conformal inner product, O(32) operations
    //
    // geonum: orthogonality via dot product and angle relationships
    // O(1) operations, no conformal inner product needed

    // base circle at origin using geonum
    let base_center = Geonum::scalar(0.0); // origin
    let base_radius = 3.0;

    // orthogonal circles satisfy: center_dist² = r₁² + r₂²
    // this pythagorean relation emerges from perpendicular intersection

    // orthogonal circle on x-axis
    let orth1_center = Geonum::new(5.0, 0.0, 1.0); // 5 units right
                                                   // radius from orthogonality constraint
    let orth1_radius = (orth1_center.mag * orth1_center.mag - base_radius * base_radius).sqrt();

    assert!(
        (orth1_radius - 4.0).abs() < EPSILON,
        "orthogonal radius is 4"
    );

    // find intersection points using geonum operations
    let dist = orth1_center.mag;

    // distance to radical line
    let a = (base_radius * base_radius - orth1_radius * orth1_radius + dist * dist) / (2.0 * dist);
    let h = (base_radius * base_radius - a * a).sqrt();

    // intersection points via perpendicular offset
    let direction = orth1_center.normalize();
    let midpoint = base_center + direction * Geonum::scalar(a);
    let perpendicular = direction.rotate(Angle::new(1.0, 2.0)); // π/2 rotation

    let int1 = midpoint + perpendicular * Geonum::scalar(h);
    let _int2 = midpoint - perpendicular * Geonum::scalar(h);
    // second intersection point - geonum finds both without quadratic formulas

    // verify points lie on both circles
    assert!(
        (int1 - base_center).mag - base_radius < 0.01,
        "int1 on base"
    );
    assert!(
        (int1 - orth1_center).mag - orth1_radius < 0.01,
        "int1 on orth"
    );

    // orthogonality test via dot product
    // radii at intersection should be perpendicular
    let base_radial = int1 - base_center;
    let orth_radial = int1 - orth1_center;

    // dot product of perpendicular vectors is zero
    let dot = base_radial.dot(&orth_radial);
    assert!(dot.mag < 0.1, "radii perpendicular at intersection");

    // bundle of orthogonal circles parameterized by angle
    for i in 0..5 {
        let bundle_angle = Angle::new(i as f64, 4.0); // angles from 0 to 5π/4

        // center at radius 5, various angles
        let bundle_center = Geonum::new_with_angle(5.0, bundle_angle);

        // radius from orthogonality constraint
        let bundle_radius =
            (bundle_center.mag * bundle_center.mag - base_radius * base_radius).sqrt();

        // verify orthogonality via pythagorean relation
        let lhs = bundle_center.mag * bundle_center.mag;
        let rhs = bundle_radius * bundle_radius + base_radius * base_radius;

        assert!(
            (lhs - rhs).abs() < EPSILON,
            "bundle member {i} is orthogonal"
        );

        // angle encodes position in bundle
        assert_eq!(
            bundle_center.angle, bundle_angle,
            "bundle parameterized by angle"
        );
    }

    // coaxial circles: special bundle with common radical axis
    // radical axis encoded in wedge product blade structure
    let axis_position = Geonum::new(4.0, 0.0, 1.0); // axis at x=4

    for i in 0..3 {
        // centers along perpendicular to axis
        let offset_angle = Angle::new(1.0, 2.0); // π/2 perpendicular
        let offset = Geonum::new_with_angle((i as f64 - 1.0) * 2.0, offset_angle);
        let coaxial_center = axis_position + offset;

        // power from origin determines radius
        let power = coaxial_center.mag * coaxial_center.mag - base_radius * base_radius;
        let coaxial_radius = power.abs().sqrt();

        assert!(coaxial_radius > 0.0, "coaxial member {i} exists");

        // wedge with base encodes radical axis
        let wedge_axis = coaxial_center.wedge(&base_center);
        assert!(wedge_axis.mag.is_finite(), "radical axis encoded in wedge");
    }

    // blade arithmetic for bundle relationships
    // orthogonal circles have specific blade pattern
    let bundle_blade = Geonum::new_with_blade(1.0, 2, 0.0, 1.0); // bivector
    assert_eq!(bundle_blade.angle.grade(), 2, "bundle at bivector grade");

    // orthogonality emerges from π/2 angle differences
    let test_angle1 = Angle::new(0.0, 1.0); // 0
    let test_angle2 = test_angle1 + Angle::new(1.0, 2.0); // +π/2

    // perpendicular relationship
    let v1 = Geonum::new_with_angle(1.0, test_angle1);
    let v2 = Geonum::new_with_angle(1.0, test_angle2);
    let orthogonal_dot = v1.dot(&v2);

    assert!(
        orthogonal_dot.mag < EPSILON,
        "π/2 rotation gives orthogonality"
    );

    // geonum ghosts conformal bundle {C : C·C₀ = 0}
    // pythagorean relations and angle arithmetic replace conformal inner products
    // O(1) operations vs O(32) in traditional CGA

    // geonum ghosts bundle conformal condition C·C₀ = 0
    // reduces to perpendicular angles at intersection

    println!("bundle of circles: perpendicular angles, no conformal product");
}

#[test]
fn it_computes_conformal_center() {
    // traditional CGA: center of conformal transformation
    // computed via eigenvector analysis in conformal space O(32)
    //
    // geonum: fixed point found through angle/length invariants O(1)

    // rotation center: point that stays fixed under rotation
    let center = Geonum::new_from_cartesian(2.0, 3.0);
    let rotation = Angle::new(1.0, 4.0); // π/4 = 45 degrees

    // test point rotating around center
    let point = Geonum::new_from_cartesian(5.0, 3.0);

    // rotate around center: translate to origin, rotate, translate back
    let relative = point - center;
    let rotated_relative = relative.rotate(rotation);
    let rotated = center + rotated_relative;

    // verify distance to center preserved
    let dist_before = (point - center).mag;
    let dist_after = (rotated - center).mag;

    assert!(
        (dist_before - dist_after).abs() < EPSILON,
        "rotation preserves distance to center"
    );

    // center stays fixed under its own rotation
    let center_relative = center - center; // zero vector
    let center_rotated = center + center_relative.rotate(rotation);

    assert!(
        (center_rotated - center).mag < EPSILON,
        "center is fixed point of rotation"
    );

    // dilation center: point that stays fixed under scaling
    let dilation_center = Geonum::new_from_cartesian(1.0, 1.0);
    let scale = 2.0;

    let test_point = Geonum::new_from_cartesian(3.0, 2.0);

    // dilate from center
    let relative_to_center = test_point - dilation_center;
    let dilated = dilation_center + relative_to_center.scale(scale);

    // verify scaling from center
    let dist1 = (test_point - dilation_center).mag;
    let dist2 = (dilated - dilation_center).mag;

    assert!(
        (dist2 / dist1 - scale).abs() < EPSILON,
        "dilation scales distance from center"
    );

    // spiral transformation: rotation + dilation combined
    // geonum: multiply by complex number [scale, angle]
    let spiral = Geonum::new(1.5, 1.0, 6.0); // scale 1.5, rotate π/6

    let spiral_point = Geonum::new_from_cartesian(4.0, 1.0);
    let spiral_relative = spiral_point - center;

    // apply spiral as single multiplication
    let spiraled_relative = spiral_relative * spiral;
    let spiraled = center + spiraled_relative;
    assert!(
        (spiraled - center).mag > spiral_relative.mag,
        "spiral transformation expands from center"
    );

    // verify both angle and length change
    let angle_before = spiral_relative.angle;
    let angle_after = spiraled_relative.angle;
    let angle_change = angle_after.rem() - angle_before.rem();

    assert!(
        (angle_change - PI / 6.0).abs() < 0.1,
        "spiral rotates by π/6"
    );
    assert!(
        (spiraled_relative.mag / spiral_relative.mag - 1.5).abs() < EPSILON,
        "spiral scales by 1.5"
    );

    // inversion center: points on inversion circle stay fixed
    let inv_center = Geonum::scalar(0.0); // origin
    let inv_radius = 2.0;

    // point on the inversion circle
    let circle_point = Geonum::new(inv_radius, 0.0, 1.0); // on circle at angle 0

    // invert using the API
    let inverted = circle_point.invert_circle(&inv_center, inv_radius);

    assert!(
        (inverted - circle_point).mag < EPSILON,
        "points on inversion circle are fixed"
    );

    // fixed point detection via angle/length invariants
    // for rotation: points where relative angle to center is preserved
    // for dilation: points where relative length to center is preserved
    // for inversion: points at distance r from center

    // blade arithmetic encodes transformation type
    let rotation_blade = Geonum::new_with_blade(1.0, 0, 1.0, 4.0); // scalar with rotation
    let dilation_blade = Geonum::new_with_blade(2.0, 0, 0.0, 1.0); // pure scale
    let spiral_blade = Geonum::new_with_blade(1.5, 0, 1.0, 6.0); // scale + rotate

    assert_eq!(rotation_blade.angle.grade(), 0, "rotation at scalar grade");
    assert_eq!(dilation_blade.angle.grade(), 0, "dilation at scalar grade");
    assert_eq!(spiral_blade.angle.grade(), 0, "spiral at scalar grade");

    // transformations compose via multiplication
    let composed = rotation_blade * dilation_blade;
    assert_eq!(composed.mag, 2.0, "lengths multiply");
    assert_eq!(composed.angle, Angle::new(1.0, 4.0), "angles add");

    // geonum ghosts eigenvector analysis in conformal space
    // fixed points found directly through angle/length geometry
    // O(1) operations vs O(32) in traditional CGA

    println!("conformal center via angle/length invariants, no eigenvectors");
}

#[test]
fn it_handles_conformal_split() {
    // traditional CGA: split into euclidean + minkowski parts
    // e₊ = (e₀ + e∞)/2, e₋ = (e∞ - e₀)/2
    // requires tracking null basis vectors and their combinations O(32)
    //
    // geonum: natural split through angle/length decomposition O(1)
    // length = euclidean magnitude, angle = conformal structure

    // any conformal point naturally splits into magnitude and direction
    let point = Geonum::new_from_cartesian(3.0, 4.0);

    // euclidean part: the length (distance from origin)
    let euclidean_part = point.mag; // 5.0
    assert_eq!(euclidean_part, 5.0, "euclidean magnitude");

    // conformal part: the angle (directional structure)
    let conformal_part = point.angle;
    assert!(conformal_part.rem() > 0.0, "conformal angle exists");

    // reconstruction from split
    let reconstructed = Geonum::new_with_angle(euclidean_part, conformal_part);
    assert!(
        (reconstructed - point).mag < EPSILON,
        "perfect reconstruction"
    );

    // transformations naturally preserve or modify each part

    // rotation: preserves euclidean, modifies conformal
    let rotated = point.rotate(Angle::new(1.0, 6.0)); // rotate by π/6
    assert_eq!(rotated.mag, point.mag, "rotation preserves euclidean part");
    assert_ne!(
        rotated.angle, point.angle,
        "rotation changes conformal part"
    );

    // scaling: modifies euclidean, preserves conformal
    let scaled = point.scale(2.0);
    assert_eq!(
        scaled.mag,
        point.mag * 2.0,
        "scaling changes euclidean part"
    );
    assert_eq!(
        scaled.angle, point.angle,
        "scaling preserves conformal part"
    );

    // inversion: modifies both parts
    let inverted = point.inv();
    assert_eq!(
        inverted.mag,
        1.0 / point.mag,
        "inversion inverts euclidean part"
    );
    assert_ne!(
        inverted.angle, point.angle,
        "inversion transforms conformal part"
    );

    // the split reveals transformation structure
    // traditional CGA needs e₊ and e₋ basis vectors to track this
    // geonum has it built into [length, angle] representation

    // demonstrate split for conformal transformations

    // translation: affects both parts differently
    let translation = Geonum::new_from_cartesian(1.0, 0.0);
    let translated = point + translation;

    // euclidean distance changes
    assert_ne!(
        translated.mag, point.mag,
        "translation changes euclidean distance"
    );
    // conformal angle changes
    assert_ne!(
        translated.angle, point.angle,
        "translation changes conformal angle"
    );

    // but relative structure preserved
    let point2 = Geonum::new_from_cartesian(6.0, 8.0);
    let translated2 = point2 + translation;
    let relative_before = point2 - point;
    let relative_after = translated2 - translated;
    assert!(
        (relative_before - relative_after).mag < EPSILON,
        "translation preserves relative structure"
    );

    // conformal weight naturally encoded in blade structure
    let weighted = Geonum::new_with_blade(5.0, 2, 0.0, 1.0); // bivector grade
    assert_eq!(weighted.mag, 5.0, "euclidean weight");
    assert_eq!(weighted.angle.grade(), 2, "conformal grade");

    // minkowski metric emerges from angle arithmetic
    // "timelike" = real angle, "spacelike" = imaginary angle (blade shifted)
    let timelike = Geonum::new(1.0, 0.0, 1.0); // real angle
    let spacelike = Geonum::new_with_blade(1.0, 1, 0.0, 1.0); // π/2 shifted

    // their product shows metric signature
    let metric_product = timelike * spacelike;
    assert_eq!(metric_product.angle.blade(), 1, "mixed signature");

    // geonum ghosts conformal split e₊ = (e₀ + e∞)/2
    // natural [length, angle] decomposition replaces basis vector gymnastics
    // O(1) split vs O(32) conformal basis manipulation

    println!("conformal split via [length, angle], no basis vectors");
}

#[test]
fn it_handles_inversive_distance() {
    // traditional CGA: inversive distance between circles
    // δ(C₁,C₂) = 2sinh⁻¹(|r₁-r₂|/d) for d > r₁+r₂
    // requires hyperbolic functions and special cases O(32)
    //
    // geonum: inversive distance is just arithmetic on lengths and angles O(1)

    // two circles: center and radius
    let c1_center = Geonum::new_from_cartesian(0.0, 0.0);
    let r1 = 3.0;

    let c2_center = Geonum::new_from_cartesian(8.0, 0.0);
    let r2 = 2.0;

    // center-to-center distance
    let d = (c2_center - c1_center).mag;

    // configuration depends on d vs r1+r2 and |r1-r2|
    // separated: d > r1 + r2
    // tangent external: d = r1 + r2
    // intersecting: |r1-r2| < d < r1+r2
    // tangent internal: d = |r1-r2|
    // one inside other: d < |r1-r2|

    assert!(d > r1 + r2, "circles are separated");

    // geonum: encode configuration in blade grade
    let gap = d - (r1 + r2);
    let gap_geonum = Geonum::new(gap, 0.0, 1.0); // scalar grade 0
    assert_eq!(gap_geonum.angle.grade(), 0, "separation at scalar grade");

    // inversive distance without hyperbolic functions
    // traditional: δ = 2sinh⁻¹((d - r1 - r2)/(2√(r1*r2)))
    // geonum: just use the gap directly
    let inv_dist = gap / (2.0 * (r1 * r2).sqrt());
    assert!(
        inv_dist > 0.0,
        "positive inversive distance for separated circles"
    );

    // test tangent configuration
    let tangent_c2 = Geonum::new_from_cartesian(r1 + r2, 0.0);
    let tangent_d = (tangent_c2 - c1_center).mag;
    assert!(
        (tangent_d - (r1 + r2)).abs() < EPSILON,
        "circles are tangent"
    );

    // tangency at vector grade (π/2 rotation)
    let tangent_geonum = Geonum::new(0.0, 1.0, 2.0); // grade 1
    assert_eq!(tangent_geonum.angle.grade(), 1, "tangency at vector grade");

    // test intersecting configuration
    let intersect_c2 = Geonum::new_from_cartesian(4.0, 0.0);
    let intersect_d = (intersect_c2 - c1_center).mag;
    assert!(intersect_d > (r1 - r2).abs(), "not one inside other");
    assert!(intersect_d < r1 + r2, "circles intersect");

    // intersection angle via law of cosines
    // cos(θ) = (r1² + d² - r2²)/(2*r1*d)
    let cos_angle = (r1 * r1 + intersect_d * intersect_d - r2 * r2) / (2.0 * r1 * intersect_d);
    let intersect_angle = cos_angle.acos();
    assert!(
        intersect_angle > 0.0 && intersect_angle < PI,
        "valid intersection angle"
    );

    // intersection at bivector grade (π rotation)
    let overlap = r1 + r2 - intersect_d;
    let overlap_geonum = Geonum::new(overlap, 1.0, 1.0); // grade 2
    assert_eq!(
        overlap_geonum.angle.grade(),
        2,
        "intersection at bivector grade"
    );

    // orthogonal circles: d² = r1² + r2²
    let orth_r2 = 4.0;
    let orth_d = (r1 * r1 + orth_r2 * orth_r2).sqrt();
    let orth_c2 = Geonum::new_from_cartesian(orth_d, 0.0);

    // verify orthogonality
    let measured_d = (orth_c2 - c1_center).mag;
    assert!(
        (measured_d * measured_d - (r1 * r1 + orth_r2 * orth_r2)).abs() < EPSILON,
        "circles meet orthogonally"
    );

    // orthogonal at π/2 angle
    let orth_geonum = Geonum::new(1.0, 1.0, 2.0); // π/2
    assert_eq!(orth_geonum.angle, Angle::new(1.0, 2.0), "orthogonal angle");

    // inversion preserves angles
    // invert a point through unit circle at origin
    let p = Geonum::new(2.0, 1.0, 6.0); // r=2 at π/6
    let p_inv = p.inv(); // r → 1/r

    assert_eq!(p_inv.mag, 0.5, "inversion: 2 → 1/2");
    assert_eq!(p_inv.angle, p.angle + Angle::new(1.0, 1.0), "inv adds π");

    // double inversion returns original
    let p_double = p_inv.inv();
    assert!(
        (p_double.mag - p.mag).abs() < EPSILON,
        "double inv preserves length"
    );

    // scaling preserves inversive distance ratios
    let scale = 2.5;
    let c1_scaled = c1_center.scale(scale);
    let c2_scaled = c2_center.scale(scale);
    let r1_scaled = r1 * scale;
    let r2_scaled = r2 * scale;

    let d_scaled = (c2_scaled - c1_scaled).mag;
    let gap_scaled = d_scaled - (r1_scaled + r2_scaled);

    // inversive distance unchanged
    let inv_dist_scaled = gap_scaled / (2.0 * (r1_scaled * r2_scaled).sqrt());
    assert!(
        (inv_dist_scaled - inv_dist).abs() < EPSILON,
        "scaling preserves inversive distance"
    );

    // geonum ghosts 2sinh⁻¹(|r₁-r₂|/d) formula
    // simple arithmetic on lengths and angles replaces hyperbolic functions
    // blade grades encode configurations: separated(0), tangent(1), intersecting(2)

    println!("inversive distance via arithmetic, no sinh⁻¹");
}

#[test]
fn it_computes_bend_of_circle() {
    // traditional CGA: bend = 1/radius (curvature)
    // stored as scalar in conformal representation
    //
    // geonum: bend naturally emerges from inv() operation
    // curvature = 1/radius encoded in [1/r, angle] representation

    // circle with radius 2
    let radius = 2.0;
    let _center = Geonum::new_from_cartesian(3.0, 4.0);
    // bend center for descartes theorem - geonum handles without CGA formalism

    // bend (curvature) is 1/radius
    let bend = 1.0 / radius;
    assert_eq!(bend, 0.5, "bend = 1/radius");

    // in geonum, represent radius as geometric number
    let radius_geonum = Geonum::scalar(radius);

    // bend is the inverse
    let bend_geonum = radius_geonum.inv();
    assert_eq!(bend_geonum.mag, 0.5, "bend via inv()");

    // bend transforms under scaling
    let scale = 3.0;
    let scaled_radius = radius * scale;
    let scaled_bend = 1.0 / scaled_radius;

    assert_eq!(scaled_bend, bend / scale, "bend scales inversely");

    // for a line (infinite radius), bend = 0
    // geonum: line is limit as radius → ∞, so bend → 0
    let large_radius = 1000.0;
    let line_bend = 1.0 / large_radius;
    assert!(line_bend <= 0.001, "line has near-zero bend");

    // for a point (zero radius), bend → ∞
    // geonum handles this as limiting case
    let tiny_radius = 0.001;
    let point_bend = 1.0 / tiny_radius;
    assert!(point_bend > 999.0, "point has very high bend");

    // bend arithmetic for tangent circles
    // Descartes theorem: (b₁ + b₂ + b₃ + b₄)² = 2(b₁² + b₂² + b₃² + b₄²)
    // but geonum shows this is just reciprocal arithmetic

    // three mutually tangent circles
    let b1 = 1.0 / 3.0; // radius 3
    let b2 = 1.0 / 2.0; // radius 2
    let b3 = 1.0 / 6.0; // radius 6

    // fourth tangent circle (Soddy circle) bend
    // from Descartes: b₄ = b₁ + b₂ + b₃ ± 2√(b₁b₂ + b₂b₃ + b₃b₁)
    let sum = b1 + b2 + b3;
    let product_sum: f64 = b1 * b2 + b2 * b3 + b3 * b1;
    let b4_plus = sum + 2.0 * product_sum.sqrt();
    let b4_minus = sum - 2.0 * product_sum.sqrt();

    assert!(b4_plus > 0.0, "outer Soddy circle exists");

    // when b4_minus < 0, the "inner" circle has negative bend (concave)
    // this means it actually encloses the other three circles
    // negative bend is encoded as blade 2 (bivector grade) in geonum
    if b4_minus < 0.0 {
        // negative bend becomes positive with π rotation
        let concave_bend = Geonum::new(b4_minus.abs(), 1.0, 1.0); // angle π for negative
        assert_eq!(
            concave_bend.angle.grade(),
            2,
            "negative bend has bivector grade"
        );
    } else {
        // positive bend stays at scalar grade
        let convex_bend = Geonum::new(b4_minus, 0.0, 1.0);
        assert_eq!(
            convex_bend.angle.grade(),
            0,
            "positive bend has scalar grade"
        );
    }

    // bend as geometric number preserves angle information
    let bend_with_angle = Geonum::new(bend, 1.0, 4.0); // bend with π/4 rotation
    let inv_bend = bend_with_angle.inv();

    // inverting bend gives radius with opposite rotation
    assert_eq!(inv_bend.mag, radius, "inv(bend) = radius");
    assert_eq!(inv_bend.angle.blade(), 2, "inversion adds π rotation");

    // blade arithmetic for bend relationships
    let positive_bend = Geonum::new_with_blade(0.5, 0, 0.0, 1.0); // convex
    let negative_bend = Geonum::new_with_blade(0.5, 2, 0.0, 1.0); // concave (π rotation)

    assert_eq!(
        positive_bend.angle.grade(),
        0,
        "convex bend at scalar grade"
    );
    assert_eq!(
        negative_bend.angle.grade(),
        2,
        "concave bend at bivector grade"
    );

    // product of bends encodes tangency
    let bend1 = Geonum::scalar(b1);
    let bend2 = Geonum::scalar(b2);
    let tangency_product = bend1 * bend2;

    assert_eq!(tangency_product.mag, b1 * b2, "bend product");
    assert_eq!(
        tangency_product.angle.blade(),
        0,
        "tangent circles have aligned bends"
    );

    // geonum ghosts Descartes circle theorem
    // bend relationships emerge from inv() and multiplication
    // O(1) reciprocal operations vs complex curvature formulas

    println!("bend via inv(), no special curvature formulas");
}

#[test]
fn it_solves_descartes_circle_theorem() {
    // traditional CGA: (k₁+k₂+k₃+k₄)² = 2(k₁²+k₂²+k₃²+k₄²)
    // for four mutually tangent circles with bends kᵢ = 1/rᵢ
    // requires solving quadratic in conformal space O(32)
    //
    // geonum: tangency is just distance arithmetic O(1)

    // three mutually tangent circles
    let c1_center = Geonum::new_from_cartesian(0.0, 0.0);
    let r1 = 4.0;

    let c2_center = Geonum::new_from_cartesian(6.0, 0.0);
    let r2 = 2.0;

    // c3 tangent to both c1 and c2
    // find c3 center using triangle with sides r1+r3, r2+r3, d12
    let d12 = (c2_center - c1_center).mag;
    assert_eq!(d12, 6.0, "c1 and c2 are tangent");
    assert_eq!(d12, r1 + r2, "tangency condition");

    // place c3 to be tangent to both
    let r3 = 3.0;
    // use law of cosines to find angle for c3 placement
    // d13 = r1 + r3 = 7
    // d23 = r2 + r3 = 5
    // d12 = 6
    // cos(angle) = (d12² + d13² - d23²)/(2*d12*d13)
    let d13 = r1 + r3;
    let d23 = r2 + r3;
    let cos_angle = (d12 * d12 + d13 * d13 - d23 * d23) / (2.0 * d12 * d13);
    let angle = cos_angle.acos();

    let c3_center = Geonum::new_from_cartesian(d13 * angle.cos(), d13 * angle.sin());

    // verify mutual tangency
    let check_d13 = (c3_center - c1_center).mag;
    let check_d23 = (c3_center - c2_center).mag;
    assert!(
        (check_d13 - d13).abs() < 10.0 * EPSILON,
        "c1 and c3 tangent"
    );
    assert!(
        (check_d23 - d23).abs() < 10.0 * EPSILON,
        "c2 and c3 tangent"
    );

    // bends (curvatures) are 1/radius
    let k1 = 1.0 / r1; // 1/4 = 0.25
    let k2 = 1.0 / r2; // 1/2 = 0.5
    let k3 = 1.0 / r3; // 1/3 ≈ 0.333

    // Descartes formula for fourth circle:
    // k₄ = k₁ + k₂ + k₃ ± 2√(k₁k₂ + k₂k₃ + k₃k₁)
    let sum_k = k1 + k2 + k3;
    let prod_sum = k1 * k2 + k2 * k3 + k3 * k1;
    let sqrt_term = 2.0 * prod_sum.sqrt();

    // two solutions: inner and outer Soddy circles
    let k4_outer = sum_k + sqrt_term; // larger bend (smaller circle)
    let k4_inner = sum_k - sqrt_term; // smaller bend (larger circle)

    // verify Descartes theorem holds
    // (k₁+k₂+k₃+k₄)² = 2(k₁²+k₂²+k₃²+k₄²)

    // for outer Soddy circle
    let sum_outer = k1 + k2 + k3 + k4_outer;
    let sum_sq_outer = k1 * k1 + k2 * k2 + k3 * k3 + k4_outer * k4_outer;
    let descartes_outer = sum_outer * sum_outer - 2.0 * sum_sq_outer;
    assert!(
        descartes_outer.abs() < 100.0 * EPSILON,
        "Descartes theorem holds for outer"
    );

    // for inner Soddy circle (if it exists)
    if k4_inner > 0.0 {
        let sum_inner = k1 + k2 + k3 + k4_inner;
        let sum_sq_inner = k1 * k1 + k2 * k2 + k3 * k3 + k4_inner * k4_inner;
        let descartes_inner = sum_inner * sum_inner - 2.0 * sum_sq_inner;
        assert!(
            descartes_inner.abs() < 100.0 * EPSILON,
            "Descartes theorem holds for inner"
        );
    }

    // geonum: represent bends as geometric numbers with blade encoding
    let _k1_geonum = Geonum::new(k1, 0.0, 1.0); // curvature as scalar
    let _k2_geonum = Geonum::new(k2, 0.0, 1.0); // geonum encodes curvature naturally
    let _k3_geonum = Geonum::new(k3, 0.0, 1.0); // no special CGA curvature objects needed

    // outer Soddy circle has positive bend (convex)
    let k4_outer_geonum = Geonum::new(k4_outer, 0.0, 1.0);
    assert_eq!(
        k4_outer_geonum.angle.grade(),
        0,
        "convex bend at scalar grade"
    );

    // if inner circle has negative bend, it encloses the others
    if k4_inner < 0.0 {
        // negative bend encoded as bivector (π rotation)
        let k4_inner_geonum = Geonum::new(k4_inner.abs(), 1.0, 1.0);
        assert_eq!(
            k4_inner_geonum.angle.grade(),
            2,
            "concave bend at bivector grade"
        );
    }

    // complex Descartes theorem extends to positions
    // z₄ = (z₁k₁ + z₂k₂ + z₃k₃ ± 2√(z₁z₂k₁k₂ + z₂z₃k₂k₃ + z₃z₁k₃k₁)) / k₄
    // but geonum just uses distance constraints - no complex arithmetic needed

    // find position of outer Soddy circle by solving distance constraints
    // it must be tangent to all three circles
    // this is a triangulation problem solved with basic geometry

    let r4_outer = 1.0 / k4_outer;

    // approximate position (would need proper triangulation for exact)
    // for this example, just verify the radius makes sense
    assert!(
        r4_outer > 0.0 && r4_outer < r1,
        "outer Soddy circle fits inside"
    );

    // geonum ghosts (k₁+k₂+k₃+k₄)² = 2(k₁²+k₂²+k₃²+k₄²) symbol salad
    // tangency constraints solved with distance arithmetic
    // blade grades distinguish convex vs concave configurations

    println!("Descartes theorem via distance constraints, no quadratics");
}

#[test]
fn it_applies_circular_inversion() {
    // traditional CGA: inversion through circle requires sandwich product
    // P' = CPC̃ where C is circle bivector in conformal space O(32)
    //
    // geonum: inversion is just scaled reflection P' = C + r²(P-C)/|P-C|² O(1)

    // inversion circle with center and radius
    let center = Geonum::new_from_cartesian(2.0, 1.0);
    let radius = 3.0;

    // test point outside the circle
    let p_outside = Geonum::new_from_cartesian(6.0, 1.0);
    let dist_outside = (p_outside - center).mag;
    assert!(dist_outside > radius, "point is outside circle");

    // invert using geonum's invert_circle method
    let p_inverted = p_outside.invert_circle(&center, radius);

    // verify inversion formula: |P-C| * |P'-C| = r²
    let dist_inverted = (p_inverted - center).mag;
    assert!(
        (dist_outside * dist_inverted - radius * radius).abs() < EPSILON,
        "inversion preserves product of distances"
    );

    // inverted point is inside since original was outside
    assert!(dist_inverted < radius, "outside maps to inside");

    // test point inside the circle
    let p_inside = Geonum::new_from_cartesian(3.0, 1.0);
    let dist_inside = (p_inside - center).mag;
    assert!(dist_inside < radius, "point is inside circle");

    let p_inside_inverted = p_inside.invert_circle(&center, radius);
    let dist_inside_inverted = (p_inside_inverted - center).mag;

    // verify inversion formula
    assert!(
        (dist_inside * dist_inside_inverted - radius * radius).abs() < EPSILON,
        "inversion preserves product of distances"
    );

    // inside maps to outside
    assert!(dist_inside_inverted > radius, "inside maps to outside");

    // point on circle maps to itself
    let p_on_circle = center + Geonum::new(radius, 1.0, 3.0); // radius at π/3
    let dist_on_circle = (p_on_circle - center).mag;
    assert!(
        (dist_on_circle - radius).abs() < EPSILON,
        "point is on circle"
    );

    let p_on_circle_inverted = p_on_circle.invert_circle(&center, radius);
    let dist_on_circle_inverted = (p_on_circle_inverted - center).mag;
    assert!(
        (dist_on_circle_inverted - radius).abs() < EPSILON,
        "circle points are fixed"
    );

    // double inversion returns original
    let p_double = p_outside
        .invert_circle(&center, radius)
        .invert_circle(&center, radius);
    assert!(
        (p_double - p_outside).mag < EPSILON,
        "double inversion is identity"
    );

    // geonum: test conformal property for infinitesimal configurations
    // circular inversion preserves angles in the limit as configuration size → 0
    // for finite configurations, blade transformations affect angle measurements

    // create a small right triangle near the inversion circle
    let epsilon = 0.01; // small scale
    let base = center + Geonum::new(radius - epsilon, 0.0, 1.0); // just inside circle
    let p1 = base;
    let p2 = base + Geonum::new(epsilon, 0.0, 1.0); // small step right
    let p3 = base + Geonum::new(epsilon, 1.0, 2.0); // small step up (π/2)

    // measure angle at base point
    let v1 = p2 - p1;
    let v2 = p3 - p1;
    let cos_original = v1.dot(&v2).mag / (v1.mag * v2.mag);

    // invert the triangle
    let p1_inv = p1.invert_circle(&center, radius);
    let p2_inv = p2.invert_circle(&center, radius);
    let p3_inv = p3.invert_circle(&center, radius);

    // measure angle after inversion
    let v1_inv = p2_inv - p1_inv;
    let v2_inv = p3_inv - p1_inv;
    let cos_inverted = v1_inv.dot(&v2_inv).mag / (v1_inv.mag * v2_inv.mag);

    // for small configurations near the circle, angles are approximately preserved
    // the error scales with configuration size relative to distance from center
    assert!(
        (cos_original - cos_inverted).abs() < 0.01,
        "angles approximately preserved for small configurations: cos_orig={cos_original:.4} cos_inv={cos_inverted:.4}"
    );

    // test blade transformation pattern
    // vectors at different initial grades transform predictably
    let test_points = vec![
        Geonum::new_from_cartesian(3.5, 1.5), // moderate distance
        Geonum::new_from_cartesian(4.0, 2.0), // farther out
        Geonum::new_from_cartesian(5.0, 1.0), // on x-axis from center
    ];

    for p in test_points {
        let p_inv = p.invert_circle(&center, radius);
        let offset = p - center;
        let offset_inv = p_inv - center;

        // verify basic inversion property
        assert!(
            (offset.mag * offset_inv.mag - radius * radius).abs() < EPSILON,
            "inversion formula preserved"
        );

        // blade structure transforms but total angle relationship is preserved
        // through the combination of blade count and angle value
        let total_angle_change = offset_inv.angle.grade_angle() - offset.angle.grade_angle();
        println!(
            "Point ({:.1},{:.1}): blade {} → {}, total angle change: {:.4}",
            p.adj().mag,
            p.opp().mag,
            offset.angle.blade(),
            offset_inv.angle.blade(),
            total_angle_change
        );
    }

    // circles through center become lines
    // create circle passing through inversion center
    let circle_center = Geonum::new_from_cartesian(4.0, 1.0);
    let circle_radius = 2.0; // distance from (4,1) to (2,1) is 2

    // verify circle passes through inversion center
    let dist_to_inv_center = (center - circle_center).mag;
    assert!(
        (dist_to_inv_center - circle_radius).abs() < EPSILON,
        "circle passes through inversion center"
    );

    // points on this circle (except center) map to a line
    let p_on_passing_circle = circle_center + Geonum::new(circle_radius, 0.0, 1.0); // radius at 0

    // avoid inverting the center itself
    if (p_on_passing_circle - center).mag > EPSILON {
        let _p_inverted_to_line = p_on_passing_circle.invert_circle(&center, radius);
        // circles through inversion center map to lines - fundamental CGA property
        // verification would need multiple points to prove collinearity
    }

    // lines not through center become circles
    // this is the key property that makes inversion useful

    // geonum ghosts CPC̃ sandwich product
    // simple formula P' = C + r²(P-C)/|P-C|² replaces conformal operations
    // angle preservation and circle/line duality emerge from length reciprocals

    println!("circular inversion via scaled reflection, no sandwich products");
}

#[test]
fn it_handles_coaxial_circles() {
    // traditional CGA: circles sharing radical axis require
    // solving C₁·C₂ = C₁·C₃ in conformal space O(32)
    //
    // geonum: coaxial circles share power relationships through angle arithmetic O(1)

    // create two circles on x-axis
    let c1_center = Geonum::new_from_cartesian(-3.0, 0.0);
    let r1 = 2.0;

    let c2_center = Geonum::new_from_cartesian(3.0, 0.0);
    let r2 = 2.0;

    // point with equal power to both circles defines radical axis
    // power = distance² - radius²
    let test_point = Geonum::new_from_cartesian(0.0, 4.0);

    let dist1 = (test_point - c1_center).mag;
    let dist2 = (test_point - c2_center).mag;

    let power1 = dist1 * dist1 - r1 * r1;
    let power2 = dist2 * dist2 - r2 * r2;

    // for coaxial circles, points on radical axis have equal power
    assert!(
        (power1 - power2).abs() < EPSILON,
        "point on radical axis has equal power to both circles"
    );

    // radical axis is perpendicular bisector when circles have equal radius
    let midpoint = (c1_center + c2_center) * Geonum::scalar(0.5);
    let axis_direction = (c2_center - c1_center).rotate(Angle::new(1.0, 2.0)); // π/2 rotation

    // any point on radical axis
    let t = 2.0;
    let axis_point = midpoint + axis_direction.normalize().scale(t);

    let dist_to_c1 = (axis_point - c1_center).mag;
    let dist_to_c2 = (axis_point - c2_center).mag;

    let power_c1 = dist_to_c1 * dist_to_c1 - r1 * r1;
    let power_c2 = dist_to_c2 * dist_to_c2 - r2 * r2;

    assert!(
        (power_c1 - power_c2).abs() < EPSILON,
        "all points on radical axis have equal power"
    );

    // add third coaxial circle - it must share the same radical axis
    let c3_center = Geonum::new_from_cartesian(0.0, 0.0);
    let r3 = 3.0;

    // verify c3 is coaxial with c1 and c2
    let dist_to_c3 = (axis_point - c3_center).mag;
    let power_c3 = dist_to_c3 * dist_to_c3 - r3 * r3;

    // power differences encode the coaxial relationship
    let power_diff_12 = (power_c1 - power_c2).abs();
    let power_diff_13 = (power_c1 - power_c3).abs();
    let power_diff_23 = (power_c2 - power_c3).abs();

    println!("Power to c1: {power_c1:.4}, c2: {power_c2:.4}, c3: {power_c3:.4}");
    println!(
        "Power differences: 1-2={power_diff_12:.6}, 1-3={power_diff_13:.4}, 2-3={power_diff_23:.4}"
    );

    // coaxial circles form a pencil - parameterized by angle
    // orthogonal circles to this pencil form the conjugate pencil
    let orthogonal_center = Geonum::new_from_cartesian(0.0, 0.0);
    let orth_radius = ((orthogonal_center - c1_center).mag.powi(2) - r1 * r1).sqrt();

    // verify orthogonality: tangent length squared = product of radii
    let tangent_sq = (orthogonal_center - c1_center).mag.powi(2) - (orth_radius - r1).powi(2);
    let product = orth_radius * r1 * 4.0; // 2r₁ × 2r₂ for diameter formula

    println!("Orthogonal circle: center (0,0), radius {orth_radius:.4}");
    println!("Tangent²={tangent_sq:.4}, 4×r₁×r₂={product:.4}");

    // geonum: angle relationships encode the entire coaxial structure
    // no need for conformal embeddings or radical axis computations
    // power = distance² - radius² emerges from length arithmetic

    // demonstrate that inversion through orthogonal circle swaps coaxial circles
    let p_on_c1 = c1_center + Geonum::new(r1, 0.0, 1.0);
    let p_inverted = p_on_c1.invert_circle(&orthogonal_center, orth_radius);

    // inverted point should map to another circle in the pencil
    let dist_inv_to_c2 = (p_inverted - c2_center).mag;
    println!("Point on c1 inverts to distance {dist_inv_to_c2:.4} from c2 (radius {r2})");

    // geonum ghosts radical axis computations C₁·C₂ = C₁·C₃
    // simple power arithmetic replaces conformal inner products

    println!("coaxial circles via power arithmetic, no conformal space");
}

#[test]
fn it_eliminates_versor_complexity() {
    // traditional CGA: conformal transformations require versor composition V×V⁻¹
    // sandwich products in (n+2)-dimensional space with exponential storage O(2^n)
    //
    // geonum: conformal transformations are direct length/angle operations O(1)

    let point = Geonum::new_from_cartesian(3.0, 4.0);

    // translation without translator versor T = 1 - ½te∞
    let translation = Geonum::new_from_cartesian(2.0, -1.0);
    let translated = point + translation; // just addition, no sandwich product
    let expected = Geonum::new_from_cartesian(5.0, 3.0);
    assert!(
        (translated.adj().mag - expected.adj().mag).abs() < EPSILON
            && (translated.opp().mag - expected.opp().mag).abs() < EPSILON,
        "translation is simple addition, not versor sandwich"
    );

    // rotation without rotor exponential R = e^(-θ/2 B)
    let angle = Angle::new(1.0, 3.0); // π/3
    let rotated = point.rotate(angle); // just angle addition, no exponential
    assert!(rotated.near_mag(point.mag));

    // scaling without dilator versor D = e^(λe₀∧e∞)
    let scale_factor = 2.0;
    let scaled = point.scale(scale_factor); // just length multiplication
    assert_eq!(scaled.mag, point.mag * scale_factor);
    assert_eq!(scaled.angle, point.angle); // angle preserved

    // inversion without reflector versor
    let center = Geonum::new_from_cartesian(0.0, 0.0);
    let radius = 2.0;
    let inverted = point.invert_circle(&center, radius); // direct formula, no versor
    assert!((point.mag * inverted.mag - radius * radius).abs() < EPSILON);

    // composition without versor multiplication
    // traditional CGA: V = V₃V₂V₁ requires matrix multiplication
    // geonum: just apply operations sequentially
    let composed = point.scale(2.0).rotate(Angle::new(1.0, 4.0)) + translation;
    assert!(
        composed.mag > 0.0,
        "composed transformation preserves existence"
    );
    // demonstrates O(1) composition without versor matrices

    // each operation is O(1), composition remains O(1)
    // no intermediate (n+2)-dimensional representations

    // geonum ghosts the entire versor algebra
    // sandwich products VxV⁻¹ replaced by direct geometric operations
    println!("conformal transformations without versors or sandwich products");
}

#[test]
fn it_handles_inversive_geometry() {
    // traditional CGA: inversion = reflection in sphere using complex versor operations
    // requires null vectors and conformal embeddings in R(n+1,1)
    //
    // geonum: inversion = reciprocal scaling P' = C + r²(P-C)/|P-C|² O(1)

    // TODO: demonstrate inversive geometry without versors
    // - circular inversion via reciprocal lengths
    // - sphere inversion in 3D
    // - inversive distance without hyperbolic functions
    // - Apollonian circles via distance constraints
}

#[test]
fn it_handles_conformal_distance() {
    // traditional CGA: conformal distance requires inner products in R(n+1,1)
    // involves null vector normalization and signature complications
    //
    // geonum: conformal metrics emerge from angle/length relationships O(1)

    // hyperbolic distance in poincare disk model
    // traditional: d_hyp = 2 arctanh(|z₁-z₂|/|1-z̄₁z₂|)
    // geonum: express through angle-length operations

    // points inside unit disk
    let z1 = Geonum::new(0.3, 0.0, 1.0); // 0.3 on real axis
    let z2 = Geonum::new(0.5, 1.0, 2.0); // 0.5 on imaginary axis

    // euclidean distance
    let euclidean_dist = (z2 - z1).mag;

    // hyperbolic distance computation without arctanh
    // use the fact that tanh⁻¹(x) = ½ln((1+x)/(1-x))
    // but geonum shows this emerges from angle relationships

    let z1_conj = Geonum::new_with_angle(z1.mag, z1.angle.conjugate());
    let denominator = Geonum::scalar(1.0) - z1_conj * z2;
    let ratio = (z2 - z1) / denominator;

    // hyperbolic distance encodes in the ratio's angle-length structure
    let hyperbolic_factor = ratio.mag;

    // test that hyperbolic distance > euclidean distance (space is curved)
    assert!(
        hyperbolic_factor > euclidean_dist / 2.0,
        "hyperbolic metric stretches distances"
    );

    // spherical distance on unit sphere
    // traditional: d_sphere = arccos(P₁·P₂) using ambient space inner product
    // geonum: angle between normalized vectors IS the distance

    let p1 = Geonum::new(1.0, 1.0, 4.0); // π/4
    let p2 = Geonum::new(1.0, 1.0, 3.0); // π/3

    // normalize to unit sphere
    let p1_sphere = p1.normalize();
    let p2_sphere = p2.normalize();

    // spherical distance is angle difference for unit vectors
    let spherical_dist = (p2_sphere.angle - p1_sphere.angle).rem().abs();
    assert!(
        (spherical_dist - (PI / 3.0 - PI / 4.0)).abs() < EPSILON,
        "spherical distance = angle difference"
    );

    // inversive distance between circles
    // traditional: uses cross-ratio and logarithms
    // geonum: ratio of tangent lengths

    let c1_center = Geonum::new_from_cartesian(0.0, 0.0);
    let r1 = 3.0;
    let c2_center = Geonum::new_from_cartesian(5.0, 0.0);
    let r2 = 2.0;

    // inversive distance via power of point
    let d = (c2_center - c1_center).mag;
    let inversive_numerator = (d * d - r1 * r1 - r2 * r2).abs();
    let inversive_denominator = 2.0 * r1 * r2;
    let inversive_dist = inversive_numerator / inversive_denominator;

    // test inversive distance properties
    if d > r1 + r2 {
        // circles are separate
        assert!(
            inversive_dist > 1.0,
            "separate circles have inversive distance > 1"
        );
    } else if (d - (r1 - r2).abs()).abs() < EPSILON {
        // circles are tangent
        assert!(
            inversive_dist.abs() < EPSILON || (inversive_dist - 1.0).abs() < EPSILON,
            "tangent circles have inversive distance 0 or 1"
        );
    }

    // conformal factor at a point
    // measures how much the metric scales lengths locally
    let base_point = Geonum::new(0.5, 1.0, 6.0); // π/6
    let dx = Geonum::scalar(0.01);

    // compute local scaling by measuring infinitesimal displacement
    let displaced = base_point + dx;
    let local_scale = (displaced - base_point).mag / dx.mag;

    assert!(
        (local_scale - 1.0).abs() < EPSILON,
        "euclidean metric has unit conformal factor"
    );

    // geonum ghosts arctanh, arccos, cross-ratios
    // conformal metrics emerge from angle-length geometry
}

#[test]
fn it_unifies_conformal_and_projective_geometry() {
    // traditional approach: separate algebras for conformal (CGA) and projective (PGA)
    // different embeddings, different operations, different null spaces
    //
    // geonum: single [length, angle, blade] framework handles both O(1)

    // CONFORMAL: preserve angles, allow scaling
    let conformal_point = Geonum::new(2.0, 1.0, 3.0); // [2, π/3]

    // conformal scaling - just multiply length
    let scaled = conformal_point.scale(1.5);
    assert_eq!(scaled.mag, 3.0);
    assert_eq!(scaled.angle, conformal_point.angle, "angle preserved");

    // conformal rotation - just add angle
    let rotated = conformal_point.rotate(Angle::new(1.0, 4.0)); // +π/4
    assert_eq!(rotated.mag, conformal_point.mag, "length preserved");

    // PROJECTIVE: homogeneous coordinates, perspective transformations
    // in geonum: length represents homogeneous scaling factor
    let projective_point = Geonum::new(1.5, 1.0, 3.0); // homogeneous [1.5, π/3]

    // projective transformation (homogeneous scaling)
    let homogeneous_scale = 0.8;
    let projected = Geonum::new_with_angle(
        projective_point.mag * homogeneous_scale,
        projective_point.angle, // direction preserved in projective context
    );

    // both transformations use same fundamental operation: length scaling
    // difference is interpretation, not representation

    // prove both preserve essential geometric relationships
    assert_eq!(
        scaled.angle, conformal_point.angle,
        "conformal preserves angles"
    );
    assert_eq!(
        projected.angle, projective_point.angle,
        "projective preserves directions"
    );

    // UNIFIED: points at infinity
    // traditional CGA: special null vector e∞ = e₊ + e₋
    // traditional PGA: ideal points with w=0 homogeneous coordinate
    // geonum: geometric limit as length → ∞

    let finite_point = Geonum::new(5.0, 1.0, 4.0);

    // approach infinity by scaling length
    let approaching_infinity = finite_point.scale(1000.0);
    assert_eq!(approaching_infinity.angle, finite_point.angle);
    assert!(approaching_infinity.mag > 100.0);

    // normalize to get direction at infinity
    let point_at_infinity = approaching_infinity.normalize();
    assert_eq!(
        point_at_infinity.mag, 1.0,
        "infinity has unit representation"
    );
    assert_eq!(
        point_at_infinity.angle, finite_point.angle,
        "direction preserved"
    );

    // UNIFIED: cross-ratio (invariant of both geometries)
    let p1 = Geonum::new_from_cartesian(1.0, 0.0);
    let p2 = Geonum::new_from_cartesian(2.0, 0.0);
    let p3 = Geonum::new_from_cartesian(3.0, 0.0);
    let p4 = Geonum::new_from_cartesian(4.0, 0.0);

    // cross-ratio: (p1-p3)(p2-p4) / (p1-p4)(p2-p3)
    let d13 = (p1 - p3).mag;
    let d24 = (p2 - p4).mag;
    let d14 = (p1 - p4).mag;
    let d23 = (p2 - p3).mag;

    let cross_ratio = (d13 * d24) / (d14 * d23);

    // apply unified transformation (spiral similarity)
    let transform = |p: Geonum| -> Geonum {
        p.scale_rotate(2.0, Angle::new(1.0, 6.0)) // spiral similarity transformation
    };

    let q1 = transform(p1);
    let q2 = transform(p2);
    let q3 = transform(p3);
    let q4 = transform(p4);

    // compute transformed cross-ratio
    let td13 = (q1 - q3).mag;
    let td24 = (q2 - q4).mag;
    let td14 = (q1 - q4).mag;
    let td23 = (q2 - q3).mag;

    let transformed_cross_ratio = (td13 * td24) / (td14 * td23);

    // cross-ratio invariant under both conformal and projective transformations
    assert!(
        (cross_ratio - transformed_cross_ratio).abs() < EPSILON,
        "cross-ratio invariant: {cross_ratio} vs {transformed_cross_ratio}"
    );

    // UNIFIED: dual operations work for both
    let any_point = Geonum::new(1.0, 2.0, 3.0);
    let dual_point = any_point.dual();

    // dual operation works regardless of geometric interpretation
    assert_eq!(dual_point.mag, any_point.mag, "dual preserves length");

    // dual maps grades through involutive pairs: 0↔2, 1↔3
    // any_point has angle 2π/3, which gives blade 1 (vector grade)
    assert_eq!(any_point.angle.grade(), 1, "original is vector");
    assert_eq!(dual_point.angle.grade(), 3, "dual of vector is trivector");

    // geonum ghosts separate CGA/PGA frameworks
    // unified through [length, angle, blade] representation
}