geoit 0.0.2

//! Rotor construction via the Euler identity.
//!
//! In Clifford algebra, a rotor that rotates by angle θ in the plane
//! of a unit bivector B̂ (where B̂² = −1) is:
//!
//!   R = cos(θ/2) + sin(θ/2) · B̂
//!
//! The sandwich product R·x·R̃ then rotates x by θ in the B̂ plane.
//!
//! This module provides [`rotor_expr`] which builds the Expr AST for
//! a rotor given a bivector expression and a constructible angle
//! (rational multiple of π). No new Expr variants are needed — the
//! rotor is expressed as `Literal(cos) + Literal(sin) * bivector`.
//!
//! # Example
//!
//! ```ignore
//! // 90° rotation in the e₁e₂ plane:
//! let bv = Expr::Outer(Box::new(Expr::Generator(0)), Box::new(Expr::Generator(1)));
//! let rotor = rotor_expr(bv, 1, 2)?;  // angle = π/2
//! // rotor = cos(π/4) + sin(π/4) · (e₁∧e₂) = √2/2 + √2/2 · e₁₂
//! ```

use crate::governance::expr::Expr;
use crate::scalar::promote::Scalar;
use crate::scalar::trig::{exact_cos_sin_pi, TrigError};

/// Build a rotor Expr that rotates by `(angle_num / angle_den) · π` radians
/// in the plane of the given bivector expression.
///
/// The bivector expression should evaluate to a **unit bivector** (B² = −1)
/// at eval time. If B is not unit, the result is a general versor, not a
/// pure rotation.
///
/// The angle must be a constructible rational multiple of π. Supported
/// denominators (after reduction): 1, 2, 3, 4, 6, 12.
///
/// # Errors
///
/// Returns `TrigError` if the half-angle is not in the constructible set.
pub fn rotor_expr(unit_bivector: Expr, angle_num: i64, angle_den: i64) -> Result<Expr, TrigError> {
    // Rotor uses half-angle: cos(θ/2) + sin(θ/2) · B̂
    // θ = (angle_num / angle_den) · π
    // θ/2 = (angle_num / (2 · angle_den)) · π
    let (cos_half, sin_half) = exact_cos_sin_pi(angle_num, 2 * angle_den)?;

    let cos_expr = Expr::Literal(Scalar::from(cos_half));

    if sin_half.is_zero() {
        // θ = 0 or 2π: rotor is just ±1 (scalar)
        return Ok(cos_expr);
    }

    let sin_expr = Expr::Literal(Scalar::from(sin_half));
    let sin_bv = Expr::Mul(Box::new(sin_expr), Box::new(unit_bivector));

    if cos_half_is_zero(angle_num, angle_den) {
        // θ = π: rotor is just sin(π/2) · B̂ = B̂
        return Ok(sin_bv);
    }

    Ok(Expr::Add(Box::new(cos_expr), Box::new(sin_bv)))
}

/// Check if cos(θ/2) is zero without recomputing.
fn cos_half_is_zero(angle_num: i64, angle_den: i64) -> bool {
    // cos(θ/2) = 0 when θ/2 = π/2 + kπ, i.e., θ = π + 2kπ
    // In our representation: angle_num / angle_den = 1 + 2k
    // So angle_num = angle_den * (1 + 2k) for some integer k.
    if angle_den == 0 {
        return false;
    }
    let rem = angle_num.rem_euclid(2 * angle_den);
    rem == angle_den
}

/// Build a rotor Expr for rotation in a coordinate plane.
///
/// Convenience for the common case where the rotation plane is
/// spanned by two basis generators. Constructs the unit bivector
/// e_i ∧ e_j (which is automatically unit in Euclidean signature
/// where e_i² = e_j² = +1).
///
/// # Arguments
///
/// * `gen_i`, `gen_j` — generator indices (the bivector e_i ∧ e_j)
/// * `angle_num`, `angle_den` — angle = (num/den) · π
pub fn coordinate_rotor_expr(
    gen_i: u8,
    gen_j: u8,
    angle_num: i64,
    angle_den: i64,
) -> Result<Expr, TrigError> {
    let bivector = Expr::Outer(
        Box::new(Expr::Generator(gen_i)),
        Box::new(Expr::Generator(gen_j)),
    );
    rotor_expr(bivector, angle_num, angle_den)
}

/// Build the full rotation Expr: sandwich(R, x) = R · x · R̃
///
/// Given a rotor expression R and a target expression x, produces
/// the sandwich product that applies the rotation.
pub fn apply_rotor(rotor: Expr, target: Expr) -> Expr {
    Expr::Sandwich(Box::new(rotor), Box::new(target))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::algebra::mv::Mv;
    use crate::algebra::ops;
    use crate::algebra::signature::Signature;
    use crate::governance::expr::EvalContext;

    fn eval_in_vga3(expr: &Expr, params: &[Scalar]) -> Mv {
        let sig = Signature::new(0, 0, 3).unwrap();
        let ctx = EvalContext {
            params,
            sig: &sig,
            derived_gens: &[],
            constructions: &[],
            mv_table: &[],
            governances: &[],
            mv_governance_indices: &[],
            embeddings: &[],
            morphisms: &[],
            probe_mv: None,
            object_mv: None,
        };
        expr.eval(&ctx)
    }

    #[test]
    fn rotor_90_in_e12_plane() {
        // Rotate e₁ by 90° in the e₁e₂ plane → should give e₂
        let r = coordinate_rotor_expr(0, 1, 1, 2).unwrap(); // π/2
        let target = Expr::Generator(0); // e₁
        let rotated = apply_rotor(r, target);
        let result = eval_in_vga3(&rotated, &[]);
        // R = cos(π/4) + sin(π/4)·e₁₂ = √2/2(1 + e₁₂)
        // R·e₁·R̃ should be e₂
        let c0 = result.coefficient(0b001); // e₁
        let c1 = result.coefficient(0b010); // e₂
        assert!(c0.is_zero(), "e₁ component should be zero, got {:?}", c0);
        assert!(!c1.is_zero(), "e₂ component should be nonzero");
    }

    #[test]
    fn rotor_180_in_e12_plane() {
        // Rotate e₁ by 180° in e₁e₂ plane → should give -e₁
        let r = coordinate_rotor_expr(0, 1, 1, 1).unwrap(); // π
        let target = Expr::Generator(0);
        let rotated = apply_rotor(r, target);
        let result = eval_in_vga3(&rotated, &[]);
        let c0 = result.coefficient(0b001);
        assert_eq!(
            c0,
            Scalar::from(-1i64),
            "180° rotation of e₁ should give -e₁"
        );
    }

    #[test]
    fn rotor_360_is_identity() {
        // Rotate by 2π → identity (R = -1, but R·x·R̃ = (-1)·x·(-1) = x)
        let r = coordinate_rotor_expr(0, 1, 2, 1).unwrap(); // 2π
        let target = Expr::Generator(0);
        let rotated = apply_rotor(r, target);
        let result = eval_in_vga3(&rotated, &[]);
        let c0 = result.coefficient(0b001);
        assert_eq!(c0, Scalar::from(1i64), "360° rotation should be identity");
    }

    #[test]
    fn rotor_60_preserves_norm() {
        // Rotation preserves norm: |R·v·R̃|² = |v|²
        let sig = Signature::new(0, 0, 3).unwrap();
        let r_expr = coordinate_rotor_expr(0, 1, 1, 3).unwrap(); // π/3 = 60°
        let v_expr = Expr::Add(
            Box::new(Expr::Mul(
                Box::new(Expr::Literal(Scalar::from(3i64))),
                Box::new(Expr::Generator(0)),
            )),
            Box::new(Expr::Mul(
                Box::new(Expr::Literal(Scalar::from(4i64))),
                Box::new(Expr::Generator(1)),
            )),
        );

        let ctx = EvalContext {
            params: &[],
            sig: &sig,
            derived_gens: &[],
            constructions: &[],
            mv_table: &[],
            governances: &[],
            mv_governance_indices: &[],
            embeddings: &[],
            morphisms: &[],
            probe_mv: None,
            object_mv: None,
        };
        let v = v_expr.eval(&ctx);
        let rotated_expr = apply_rotor(r_expr, v_expr);
        let rv = rotated_expr.eval(&ctx);

        let norm_v = ops::norm_squared(&v, &sig);
        let norm_rv = ops::norm_squared(&rv, &sig);
        assert_eq!(norm_v, norm_rv, "rotation should preserve norm²");
    }

    #[test]
    fn rotor_non_constructible_rejected() {
        let err = coordinate_rotor_expr(0, 1, 1, 7); // π/7
        assert!(err.is_err());
    }
}