algos 0.6.8

//! A self-contained, demonstrative (but general) Schur–Zassenhaus (Zassenhaus) algorithm in Rust.
//! This example works with finite groups given as subgroups of the symmetric group S_n (i.e., as
//! sets of permutations over {0..n-1}). It shows how to:
//!   1. Construct a group from generating permutations.
//!   2. Find a Sylow p-subgroup (by naive subgroup search).
//!   3. Check normality of that Sylow subgroup.
//!   4. Attempt to find a complement H such that G = N ⋊ H (if the normal Sylow p-subgroup is indeed normal).
//!
//! A real production-grade Schur–Zassenhaus implementation is large and involves many advanced
//! optimizations. This code, while "complete" in logic, is heavily brute-forced and only feasible
//! for small groups. It is purely demonstrative and will not handle large groups efficiently.
//!
//! ## Overview
//!   - `Permutation`: A permutation over {0..n-1}.
//!   - `PermutationGroup`: A set of permutations closed under composition and inverses, plus helpful
//!     group-theoretic methods (order, membership, subgroups, normality checks).
//!   - `schur_zassenhaus`: Attempts to factor a finite group `G` into a normal Sylow p-subgroup `N` and
//!     a complement `H` of matching order (if it exists).
//!
//! ## Example
//! ```
//! use algos::cs::combinatorial::zassenhaus::{Permutation, PermutationGroup, schur_zassenhaus};
//!
//! // Example: group S3 with permutations of {0,1,2}.
//! // We'll generate it from the transpositions (01), (12) (typical S3 generators).
//! let p1 = Permutation::new(vec![1, 0, 2]); // (0 1)
//! let p2 = Permutation::new(vec![0, 2, 1]); // (1 2)
//! let s3 = PermutationGroup::from_generators(&[p1, p2]);
//!
//! // Let's pick p=3. S3 has a normal Sylow 3-subgroup N = {e, (0 1 2), (0 2 1)} of order 3.
//! // The complement H would be a subgroup of order 2 (generated by a transposition).
//! let factor = schur_zassenhaus(&s3, 3);
//! if let Some((n_subgroup, h_subgroup)) = factor {
//!     assert_eq!(n_subgroup.order(), 3);
//!     assert_eq!(h_subgroup.order(), 2);
//! }
//! ```

use std::collections::HashSet;

/// A permutation of `0..n-1`.
#[derive(Clone, Debug, PartialEq, Eq, Hash)]
pub struct Permutation {
    /// `mapping[i]` = image of `i` under this permutation.
    mapping: Vec<usize>,
}

impl Permutation {
    /// Creates a new permutation from its `mapping`.
    pub fn new(mapping: Vec<usize>) -> Self {
        // Basic check that `mapping` is a valid permutation over {0..mapping.len()-1}.
        assert!(
            is_valid_permutation(&mapping),
            "Provided vector is not a valid permutation."
        );
        Permutation { mapping }
    }

    /// Returns the size of the set on which this permutation acts.
    pub fn len(&self) -> usize {
        self.mapping.len()
    }

    /// Returns true if this permutation acts on an empty set.
    pub fn is_empty(&self) -> bool {
        self.mapping.is_empty()
    }

    /// Permute an element `x`.
    pub fn apply(&self, x: usize) -> usize {
        self.mapping[x]
    }

    /// Returns the composition `self * other` (apply `other` first, then `self`).
    pub fn compose(&self, other: &Self) -> Self {
        // Both permutations must have the same size.
        assert_eq!(
            self.len(),
            other.len(),
            "Cannot compose permutations of different sizes."
        );
        let n = self.len();
        let mut new_map = vec![0; n];
        for (i, item) in new_map.iter_mut().enumerate().take(n) {
            *item = self.mapping[other.mapping[i]];
        }
        Permutation::new(new_map)
    }

    /// Returns the inverse of this permutation.
    pub fn inverse(&self) -> Self {
        let n = self.len();
        let mut inv = vec![0; n];
        for i in 0..n {
            inv[self.mapping[i]] = i;
        }
        Permutation::new(inv)
    }

    /// The identity permutation of size `n`.
    pub fn identity(n: usize) -> Self {
        Permutation::new((0..n).collect())
    }
}

// Simple helper to validate a permutation's mapping is one-to-one over {0..n-1}.
fn is_valid_permutation(mapping: &[usize]) -> bool {
    let n = mapping.len();
    let mut seen = vec![false; n];
    for &m in mapping {
        if m >= n || seen[m] {
            return false;
        }
        seen[m] = true;
    }
    true
}

/// A finite permutation group acting on {0..n-1}, stored explicitly as a set of permutations closed under composition.
#[derive(Clone, Debug)]
pub struct PermutationGroup {
    /// All elements (permutations) of the group, including the identity.
    elements: Vec<Permutation>,
    /// The size of the set on which these permutations act.
    n: usize,
}

impl PermutationGroup {
    /// Builds a group by generating all elements from the given set of permutations
    /// (together with the identity).
    /// This is a naive closure under composition and inverses, suitable only for small groups.
    pub fn from_generators(generators: &[Permutation]) -> Self {
        assert!(
            !generators.is_empty(),
            "Must provide at least one generator."
        );
        let n = generators[0].len();
        for g in generators {
            assert_eq!(g.len(), n, "Generators must have uniform permutation size.");
        }

        // Build closure under composition. BFS approach:
        let mut group_elems = HashSet::new();
        let id = Permutation::identity(n);
        group_elems.insert(id);
        for gen in generators {
            group_elems.insert(gen.clone());
        }

        let mut queue: Vec<Permutation> = group_elems.iter().cloned().collect();
        let mut idx = 0;
        while idx < queue.len() {
            // Clone current element to avoid borrow checker issues
            let current = queue[idx].clone();
            idx += 1;

            // Compose current with each generator
            for g in generators {
                let comp = current.compose(g);
                if group_elems.insert(comp.clone()) {
                    queue.push(comp);
                }
            }

            // Compose current with each known element
            let known_elems: Vec<_> = group_elems.iter().cloned().collect();
            for e in known_elems {
                let comp = current.compose(&e);
                if group_elems.insert(comp.clone()) {
                    queue.push(comp);
                }
            }
        }

        PermutationGroup {
            elements: group_elems.into_iter().collect(),
            n,
        }
    }

    /// Returns the number of elements in the group.
    pub fn order(&self) -> usize {
        self.elements.len()
    }

    /// Checks if `p` is in this group (by direct membership in `elements`).
    pub fn contains(&self, p: &Permutation) -> bool {
        self.elements.contains(p)
    }

    /// Returns the identity element of this group (we assume there's always exactly one).
    pub fn identity(&self) -> Option<Permutation> {
        self.elements.iter().find(|perm| is_identity(perm)).cloned()
    }

    /// Enumerate all subgroups of `self` in the simplest (exponential) manner.
    /// This is extremely expensive for larger groups; purely illustrative.
    pub fn all_subgroups(&self) -> Vec<PermutationGroup> {
        // We'll find all subsets of `self.elements` that form a subgroup.
        // Then build a PermutationGroup from them. This is 2^(|G|) in worst case.
        // Not feasible for large |G|.
        let all_elems: Vec<Permutation> = self.elements.clone();
        let mut subgroups = Vec::new();

        // Use bitmasks to pick subsets
        let total = 1 << self.order();
        'outer: for mask in 0..total {
            // Collect a subset
            let mut subset = Vec::new();
            for (i, elem) in all_elems.iter().enumerate().take(self.order()) {
                if (mask & (1 << i)) != 0 {
                    subset.push(elem.clone());
                }
            }
            if subset.is_empty() {
                continue; // skip
            }
            // Must contain identity
            if !subset.iter().any(is_identity) {
                continue;
            }
            // Must be closed under composition and inverses
            for a in &subset {
                for b in &subset {
                    let ab = a.compose(b);
                    if !subset.contains(&ab) {
                        continue 'outer;
                    }
                    let a_inv = a.inverse();
                    if !subset.contains(&a_inv) {
                        continue 'outer;
                    }
                }
            }
            // If we reach here, it's a subgroup
            // Build a new PermutationGroup from it directly
            let sub = PermutationGroup {
                elements: subset,
                n: self.n,
            };
            subgroups.push(sub);
        }
        subgroups
    }

    /// Check normality of `self` in `parent`, i.e. for all g in parent, x in self: g x g^-1 in self.
    pub fn is_normal_in(&self, parent: &PermutationGroup) -> bool {
        if !self.is_subgroup_of(parent) {
            return false;
        }
        for g in &parent.elements {
            let g_inv = g.inverse();
            for x in &self.elements {
                let gxg = g.compose(x).compose(&g_inv);
                if !self.contains(&gxg) {
                    return false;
                }
            }
        }
        true
    }

    /// Check if `self` is a subgroup of `parent`.
    pub fn is_subgroup_of(&self, parent: &PermutationGroup) -> bool {
        // All elements must be in `parent`.
        for e in &self.elements {
            if !parent.contains(e) {
                return false;
            }
        }
        true
    }
}

/// Check if a permutation is the identity mapping.
fn is_identity(p: &Permutation) -> bool {
    for (i, &m) in p.mapping.iter().enumerate() {
        if i != m {
            return false;
        }
    }
    true
}

/// Find a Sylow p-subgroup of G (if one exists).
/// A Sylow p-subgroup is a subgroup of order p^k where p^k divides |G| but p^(k+1) does not.
pub fn find_sylow_p_subgroup(group: &PermutationGroup, p: usize) -> Option<PermutationGroup> {
    let order = group.order();
    let mut p_power = 1;
    while order % (p_power * p) == 0 {
        p_power *= p;
    }
    group
        .all_subgroups()
        .into_iter()
        .find(|sub| sub.order() == p_power)
}

/// Find a complement H to N in G (if one exists).
/// A complement H is a subgroup of G such that:
///   1. G = NH (every element of G is a product of elements from N and H)
///   2. N ∩ H = {e} (N and H intersect only at the identity)
///   3. |N| * |H| = |G| (orders multiply)
///
/// This naive approach enumerates subgroups of G of the correct order and checks these conditions.
pub fn find_complement(g: &PermutationGroup, n: &PermutationGroup) -> Option<PermutationGroup> {
    let target_order = g.order() / n.order();

    'subgroup_search: for h_sub in g.all_subgroups() {
        if h_sub.order() == target_order {
            // Check trivial intersection: N ∩ H = {identity}.
            // We skip any that share non-identity elements.
            for x in &h_sub.elements {
                if !is_identity(x) && n.contains(x) {
                    continue 'subgroup_search;
                }
            }

            // Check if every element g in G can be written as n*h for some n in N, h in H.
            // This is extremely naive: try all pairs.
            if covers_whole_group(g, n, &h_sub) {
                return Some(h_sub);
            }
        }
    }
    None
}

/// Check if group G = N * H (setwise) by verifying that for every g in G,
/// there exist n in N and h in H with g = n*h. Very brute force.
fn covers_whole_group(g: &PermutationGroup, n: &PermutationGroup, h: &PermutationGroup) -> bool {
    'outer: for x in &g.elements {
        for nn in &n.elements {
            for hh in &h.elements {
                let comp = nn.compose(hh);
                if &comp == x {
                    continue 'outer;
                }
            }
        }
        return false; // didn't find representation for x
    }
    true
}

/// Main Schur–Zassenhaus procedure: factor G into (N, H) where N is a normal Sylow p-subgroup,
/// and H is a complement. Returns `None` if no such factorization is found.
pub fn schur_zassenhaus(
    group: &PermutationGroup,
    p: usize,
) -> Option<(PermutationGroup, PermutationGroup)> {
    // For the trivial group, there's no non-trivial factorization
    if group.order() == 1 {
        return None;
    }

    // 1. Find a Sylow p-subgroup of G by naive enumeration.
    let n_sub = find_sylow_p_subgroup(group, p)?;

    // 2. Check normality.
    if !n_sub.is_normal_in(group) {
        return None;
    }

    // 3. Attempt to find a complement H of order |G| / |N|.
    if group.order() % n_sub.order() != 0 {
        return None;
    }
    let h_sub = find_complement(group, &n_sub)?;

    // If found, we consider that the factorization is G = N ⋊ H.
    Some((n_sub, h_sub))
}

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

    #[test]
    fn test_s3_schur_zassenhaus_p3() {
        // S3 = group of permutations on {0,1,2}. We'll build from typical transposition generators.
        let p1 = Permutation::new(vec![1, 0, 2]); // (0 1)
        let p2 = Permutation::new(vec![0, 2, 1]); // (1 2)
        let s3 = PermutationGroup::from_generators(&[p1, p2]);

        // |S3|=6. A Sylow 3-subgroup N has order 3 and is normal in S3.
        // Schur–Zassenhaus says S3 splits as N ⋊ H for p=3. H would be of order 2.
        let factor = schur_zassenhaus(&s3, 3);
        assert!(factor.is_some());
        let (n_sub, h_sub) = factor.unwrap();
        assert_eq!(n_sub.order(), 3);
        assert_eq!(h_sub.order(), 2);
        // Check trivial intersection
        let mut intersection_size = 0;
        for x in &n_sub.elements {
            if h_sub.contains(x) {
                intersection_size += 1;
            }
        }
        // Intersection should contain exactly the identity.
        assert_eq!(intersection_size, 1);
    }

    #[test]
    fn test_s3_schur_zassenhaus_p2() {
        // S3 again, but p=2 => a Sylow 2-subgroup has order 2 (generated by a single transposition).
        let p1 = Permutation::new(vec![1, 0, 2]); // (0 1)
        let p2 = Permutation::new(vec![0, 2, 1]); // (1 2)
        let s3 = PermutationGroup::from_generators(&[p1, p2]);

        // A Sylow 2-subgroup might be {e, (0 1)} for instance, but it's not normal in S3.
        // So the algorithm won't find a normal Sylow 2-subgroup. We expect None.
        let factor = schur_zassenhaus(&s3, 2);
        assert!(factor.is_none());
    }

    #[test]
    fn test_trivial() {
        // The trivial group (single identity on {0}).
        let trivial = Permutation::new(vec![0]);
        let g = PermutationGroup::from_generators(&[trivial]);
        // For any p, there's no interesting factorization.
        assert!(schur_zassenhaus(&g, 2).is_none());
        assert!(schur_zassenhaus(&g, 3).is_none());
    }
}