geoit 0.0.2

//! Exact trigonometry for constructible angles.
//!
//! Computes `cos(nπ/d)` and `sin(nπ/d)` as exact [`RadicalElement`] values
//! for rational multiples of π whose denominators yield constructible regular
//! polygons. Supported base angles after quadrant reduction:
//!
//! | Angle | cos | sin |
//! |-------|-----|-----|
//! | 0 | 1 | 0 |
//! | π/6 | √3/2 | 1/2 |
//! | π/4 | √2/2 | √2/2 |
//! | π/3 | 1/2 | √3/2 |
//! | π/2 | 0 | 1 |
//! | 5π/12 | via addition formula | |
//! | π/12 | via addition formula | |
//!
//! Quadrant reduction, negation, and supplementary angle identities
//! extend these to all rational multiples of π with the above denominators.

use crate::scalar::radical::{RadicalElement, RadicalLayer, RadicalTower};
use crate::scalar::rat::Rat;
use std::sync::Arc;

/// Error from exact trigonometric computation.
#[derive(Clone, Debug)]
pub enum TrigError {
    /// Denominator was zero.
    ZeroDenominator,
    /// The angle nπ/d is not constructible with the supported base angles.
    /// The values are the reduced fraction.
    NonConstructible { num: i64, den: i64 },
}

impl std::fmt::Display for TrigError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            TrigError::ZeroDenominator => write!(f, "denominator is zero"),
            TrigError::NonConstructible { num, den } => write!(
                f,
                "angle {}π/{} is not constructible with supported base angles",
                num, den
            ),
        }
    }
}

impl std::error::Error for TrigError {}

fn gcd_u64(mut a: u64, mut b: u64) -> u64 {
    while b != 0 {
        let t = b;
        b = a % b;
        a = t;
    }
    a
}

/// Compute `(cos(n·π/d), sin(n·π/d))` as exact radical elements.
///
/// The angle is specified as a rational multiple of π: `n/d` times π.
/// Returns an error if the denominator is zero or the angle's base
/// component is not in the supported constructible set.
///
/// # Examples
///
/// ```ignore
/// // cos(π/3) = 1/2, sin(π/3) = √3/2
/// let (c, s) = exact_cos_sin_pi(1, 3).unwrap();
/// assert_eq!(c.to_rat(), Some(Rat::new(1, 2)));
///
/// // cos(π) = -1, sin(π) = 0
/// let (c, s) = exact_cos_sin_pi(1, 1).unwrap();
/// assert_eq!(c.to_rat(), Some(Rat::from(-1)));
/// ```
pub fn exact_cos_sin_pi(num: i64, den: i64) -> Result<(RadicalElement, RadicalElement), TrigError> {
    if den == 0 {
        return Err(TrigError::ZeroDenominator);
    }

    // Normalize: reduce fraction, ensure den > 0.
    let g = gcd_u64(num.unsigned_abs(), den.unsigned_abs()) as i64;
    let (mut n, d) = if den > 0 {
        (num / g, den / g)
    } else {
        (-num / g, -den / g)
    };

    // Map n into [0, 2d) so the angle is in [0, 2π).
    n = ((n % (2 * d)) + 2 * d) % (2 * d);

    // Quadrant reduction. Let α = n·π/d ∈ [0, 2π).
    //   Q1: 0 ≤ α < π/2   →  2n < d   → cos(α) = +C, sin(α) = +S
    //   Q2: π/2 ≤ α < π   →  d ≤ 2n < 2d → cos(π−α) = −C, sin(π−α) = +S
    //   Q3: π ≤ α < 3π/2  →  2d ≤ 2n < 3d → cos(α−π) = −C, sin(α−π) = −S
    //   Q4: 3π/2 ≤ α < 2π →  3d ≤ 2n < 4d → cos(2π−α) = +C, sin(2π−α) = −S
    let two_n = 2 * n;
    if two_n < d {
        // Q1
        base_cos_sin(n, d)
    } else if two_n < 2 * d {
        // Q2: use cos(π − α) = −cos(α), sin(π − α) = sin(α)
        let (c, s) = base_cos_sin(d - n, d)?;
        Ok((c.neg(), s))
    } else if two_n < 3 * d {
        // Q3: use cos(α − π) = −cos(α), sin(α − π) = −sin(α)
        let (c, s) = base_cos_sin(n - d, d)?;
        Ok((c.neg(), s.neg()))
    } else {
        // Q4: use cos(2π − α) = cos(α), sin(2π − α) = −sin(α)
        let (c, s) = base_cos_sin(2 * d - n, d)?;
        Ok((c, s.neg()))
    }
}

/// Compute cos(n·π/d) and sin(n·π/d) for n/d in [0, 1/2].
/// This is the first-quadrant base table plus addition formulas.
fn base_cos_sin(n: i64, d: i64) -> Result<(RadicalElement, RadicalElement), TrigError> {
    let one = || RadicalElement::from_rat(Rat::ONE);
    let zero = || RadicalElement::from_rat(Rat::ZERO);
    let half = || RadicalElement::from_rat(Rat::new(1, 2));

    if n == 0 {
        return Ok((one(), zero()));
    }

    // Reduce to lowest terms for table lookup.
    let g = gcd_u64(n.unsigned_abs(), d.unsigned_abs()) as i64;
    let (rn, rd) = (n / g, d / g);

    match (rn, rd) {
        // π/2: boundary of first quadrant
        (1, 2) => Ok((zero(), one())),
        // π/6 = 30°
        (1, 6) => Ok((sqrt_over_2(3), half())),
        // π/4 = 45°
        (1, 4) => {
            let v = sqrt_over_2(2);
            Ok((v.clone(), v))
        }
        // π/3 = 60°
        (1, 3) => Ok((half(), sqrt_over_2(3))),
        // π/12 = 15°: cos = (√6+√2)/4, sin = (√6−√2)/4
        // Directly in tower ℚ(√2, √3) with basis {1, √2, √3, √6}.
        (1, 12) => {
            let (cos, sin) = pi_over_12_direct();
            Ok((cos, sin))
        }
        // 5π/12 = 75°: cos(75°) = sin(15°), sin(75°) = cos(15°)
        (5, 12) => {
            let (cos15, sin15) = pi_over_12_direct();
            Ok((sin15, cos15))
        }
        _ => Err(TrigError::NonConstructible { num: rn, den: rd }),
    }
}

/// cos(π/12) and sin(π/12) directly in tower ℚ(√2, √3).
///
/// Basis of the 4-dimensional extension: {1, √3, √2, √6}
/// cos(π/12) = (√6 + √2)/4 = [0, 1/4, 0, 1/4]
/// sin(π/12) = (√6 − √2)/4 = [0, -1/4, 0, 1/4]
fn pi_over_12_direct() -> (RadicalElement, RadicalElement) {
    use crate::scalar::coeff::Coeff;
    let layer_sqrt2 = RadicalLayer::sqrt(Rat::from(2));
    let layer_sqrt3 = RadicalLayer::sqrt(Rat::from(3));
    let tower = Arc::new(RadicalTower::two(layer_sqrt2, layer_sqrt3));
    let q = Rat::new(1, 4);
    let cos = RadicalElement::new(
        vec![Coeff::ZERO, Coeff::ZERO, Coeff::Rat(q), Coeff::Rat(q)],
        Arc::clone(&tower),
    );
    let sin = RadicalElement::new(
        vec![Coeff::ZERO, Coeff::ZERO, Coeff::Rat(-q), Coeff::Rat(q)],
        tower,
    );
    (cos, sin)
}

/// Construct √d / 2 as a RadicalElement in tower ℚ(√d).
fn sqrt_over_2(d: i64) -> RadicalElement {
    use crate::scalar::coeff::Coeff;
    let layer = RadicalLayer::sqrt(Rat::from(d));
    let tower = Arc::new(RadicalTower::single(layer));
    RadicalElement::new(vec![Coeff::ZERO, Coeff::Rat(Rat::new(1, 2))], tower)
}

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

    fn approx(r: &RadicalElement) -> f64 {
        if let Some(rat) = r.to_rat() {
            return rat.num() as f64 / rat.den() as f64;
        }
        let coeffs = &r.coeffs;
        let tower = r.tower();
        if tower.layers.len() == 1 {
            let d = tower.layers[0].min_poly[0]; // -d for x²-d
            let sqrt_d = ((-d.num() as f64) / d.den() as f64).sqrt();
            let c0 = coeffs[0].to_rat().unwrap_or(Rat::ZERO);
            let c1 = coeffs[1].to_rat().unwrap_or(Rat::ZERO);
            c0.num() as f64 / c0.den() as f64 + c1.num() as f64 / c1.den() as f64 * sqrt_d
        } else {
            let c0 = coeffs[0].to_rat().unwrap_or(Rat::ZERO);
            c0.num() as f64 / c0.den() as f64
        }
    }

    #[test]
    fn cos_sin_zero() {
        let (c, s) = exact_cos_sin_pi(0, 1).unwrap();
        assert_eq!(c.to_rat(), Some(Rat::ONE));
        assert_eq!(s.to_rat(), Some(Rat::ZERO));
    }

    #[test]
    fn cos_sin_pi() {
        // cos(π) = -1, sin(π) = 0
        let (c, s) = exact_cos_sin_pi(1, 1).unwrap();
        assert_eq!(c.to_rat(), Some(Rat::from(-1)));
        assert!(s.is_zero());
    }

    #[test]
    fn cos_sin_pi_over_2() {
        let (c, s) = exact_cos_sin_pi(1, 2).unwrap();
        assert!(c.is_zero());
        assert_eq!(s.to_rat(), Some(Rat::ONE));
    }

    #[test]
    fn cos_sin_pi_over_3() {
        // cos(π/3) = 1/2
        let (c, s) = exact_cos_sin_pi(1, 3).unwrap();
        assert_eq!(c.to_rat(), Some(Rat::new(1, 2)));
        // sin(π/3) = √3/2 ≈ 0.866
        let s_approx = approx(&s);
        assert!((s_approx - 0.866025).abs() < 0.001);
    }

    #[test]
    fn cos_sin_pi_over_4() {
        let (c, s) = exact_cos_sin_pi(1, 4).unwrap();
        let c_approx = approx(&c);
        let s_approx = approx(&s);
        assert!((c_approx - std::f64::consts::FRAC_1_SQRT_2).abs() < 0.0001);
        assert!((s_approx - std::f64::consts::FRAC_1_SQRT_2).abs() < 0.0001);
        // cos = sin for π/4
        assert_eq!(c, s);
    }

    #[test]
    fn cos_sin_pi_over_6() {
        // cos(π/6) = √3/2, sin(π/6) = 1/2
        let (c, s) = exact_cos_sin_pi(1, 6).unwrap();
        assert_eq!(s.to_rat(), Some(Rat::new(1, 2)));
        let c_approx = approx(&c);
        assert!((c_approx - 0.866025).abs() < 0.001);
    }

    #[test]
    fn cos_sin_pi_over_12() {
        // cos(15°) ≈ 0.9659, sin(15°) ≈ 0.2588
        let (c, s) = exact_cos_sin_pi(1, 12).unwrap();
        // These are in tower ℚ(√2, √3) — not simple to approx,
        // but we can verify cos² + sin² = 1
        let c2 = c.mul(&c);
        let s2 = s.mul(&s);
        let sum = c2.add(&s2);
        assert_eq!(sum.to_rat(), Some(Rat::ONE));
    }

    #[test]
    fn cos_sin_2pi_over_3() {
        // cos(2π/3) = -1/2, sin(2π/3) = √3/2
        let (c, s) = exact_cos_sin_pi(2, 3).unwrap();
        assert_eq!(c.to_rat(), Some(Rat::new(-1, 2)));
        assert!(s.is_positive());
    }

    #[test]
    fn cos_sin_3pi_over_2() {
        // cos(3π/2) = 0, sin(3π/2) = -1
        let (c, s) = exact_cos_sin_pi(3, 2).unwrap();
        assert!(c.is_zero());
        assert_eq!(s.to_rat(), Some(Rat::from(-1)));
    }

    #[test]
    fn cos_sin_negative_angle() {
        // cos(-π/3) = cos(5π/3) = 1/2
        // sin(-π/3) = sin(5π/3) = -√3/2
        let (c, s) = exact_cos_sin_pi(-1, 3).unwrap();
        assert_eq!(c.to_rat(), Some(Rat::new(1, 2)));
        assert!(s.is_negative());
    }

    #[test]
    fn pythagorean_identity_all_base_angles() {
        // cos²(θ) + sin²(θ) = 1 for all supported angles
        for (n, d) in [
            (0, 1),
            (1, 6),
            (1, 4),
            (1, 3),
            (1, 2),
            (2, 3),
            (3, 4),
            (5, 6),
            (1, 1),
            (1, 12),
            (5, 12),
        ] {
            let (c, s) = exact_cos_sin_pi(n, d).unwrap();
            let sum = c.mul(&c).add(&s.mul(&s));
            assert_eq!(
                sum.to_rat(),
                Some(Rat::ONE),
                "cos²({}π/{}) + sin²({}π/{}) ≠ 1",
                n,
                d,
                n,
                d
            );
        }
    }

    #[test]
    fn non_constructible_rejected() {
        // π/7 is not constructible with our base set
        let err = exact_cos_sin_pi(1, 7).unwrap_err();
        assert!(matches!(err, TrigError::NonConstructible { .. }));
    }

    #[test]
    fn zero_denominator_rejected() {
        let err = exact_cos_sin_pi(1, 0).unwrap_err();
        assert!(matches!(err, TrigError::ZeroDenominator));
    }
}