geometric-traits 0.2.2

//! [`PermutationGroup`], a group of permutations of `0..degree` stored as a
//! stabilizer chain.

use alloc::{borrow::Cow, collections::BTreeMap, vec::Vec};

/// Schreier tree label of a point outside a level's orbit.
const UNREACHED: usize = usize::MAX;
/// Schreier tree label of a level's base point.
const ROOT: usize = usize::MAX - 1;

/// Error returned by [`PermutationGroup`] for input outside its domain.
#[derive(Clone, Copy, Debug, PartialEq, Eq, thiserror::Error)]
pub enum PermutationGroupError {
    /// A generator does not list exactly one image per point.
    #[error("generator {generator} has {length} images, but the group acts on {degree} points")]
    GeneratorLength {
        /// Position of the generator in the input.
        generator: usize,
        /// Number of images the generator lists.
        length: usize,
        /// Number of points the group acts on.
        degree: usize,
    },
    /// A generator repeats an image or maps a point outside `0..degree`.
    #[error("generator {generator} is not a permutation of 0..{degree}")]
    NotAPermutation {
        /// Position of the generator in the input.
        generator: usize,
        /// Number of points the group acts on.
        degree: usize,
    },
    /// A point lies outside `0..degree`.
    #[error("point {point} lies outside 0..{degree}")]
    PointOutOfRange {
        /// The offending point.
        point: usize,
        /// Number of points the group acts on.
        degree: usize,
    },
    /// A set lists the same point more than once.
    #[error("point {point} appears more than once in the set")]
    DuplicatePoint {
        /// The repeated point.
        point: usize,
    },
}

/// Smallest image of a set, returned by [`PermutationGroup::smallest_image`].
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct SmallestImage {
    /// The lexicographically smallest image, sorted ascending.
    pub image: Vec<usize>,
    /// A group element mapping the input set onto `image`, where `element[p]`
    /// is the image of point `p`.
    pub element: Vec<usize>,
}

/// A chain level fixing earlier active points and moving `base_point`.
#[derive(Clone, Debug)]
struct Level {
    base_point: usize,
    /// Indices into the strong generators that generate this level's group.
    generators: Vec<usize>,
    /// Orbit of `base_point` in discovery order.
    orbit: Vec<usize>,
    /// Strong generator labelling the Schreier tree edge into each point, or
    /// [`ROOT`] and [`UNREACHED`].
    tree: Vec<usize>,
}

impl Level {
    fn reaches(&self, point: usize) -> bool {
        self.tree[point] != UNREACHED
    }

    fn divide(&self, inverses: &[Vec<usize>], permutation: &mut [usize], point: usize) {
        divide(&self.tree, inverses, permutation, point);
    }
}

/// A group of permutations of the points `0..degree`, held as a base and
/// strong generating set.
///
/// A permutation is a slice where `permutation[p]` is the image of point `p`,
/// the convention of `CanonicalLabelingResult::automorphism_generators`.
/// Construction runs the deterministic Schreier–Sims algorithm ([Sims, 1970](https://doi.org/10.1016/B978-0-08-012975-4.50020-5))
/// over the base `0, 1, …, degree - 1`, storing only the levels whose group
/// moves its base point, so every pointwise stabilizer of `0..p` is a stored
/// level.
///
/// # Examples
///
/// ```
/// use geometric_traits::prelude::*;
///
/// // (0 1)(2 3)(4 5): eight elements, and {1, 3, 5} has smallest image {0, 2, 4}.
/// let group =
///     PermutationGroup::new(6, [[1, 0, 2, 3, 4, 5], [0, 1, 3, 2, 4, 5], [0, 1, 2, 3, 5, 4]])
///         .unwrap();
/// assert_eq!(group.order(), Some(8));
/// assert!(group.contains(&[1, 0, 3, 2, 4, 5]));
/// assert!(!group.contains(&[2, 1, 0, 3, 4, 5]));
/// assert_eq!(group.orbit(4), Ok(vec![4, 5]));
/// let smallest = group.smallest_image(&[1, 3, 5]).unwrap();
/// assert_eq!(smallest.image, vec![0, 2, 4]);
/// assert_eq!(smallest.element, vec![1, 0, 3, 2, 5, 4]);
/// ```
#[derive(Clone, Debug)]
pub struct PermutationGroup {
    degree: usize,
    strong_generators: Vec<Vec<usize>>,
    inverses: Vec<Vec<usize>>,
    /// Levels by ascending base point.
    levels: Vec<Level>,
}

impl PermutationGroup {
    /// Builds the group on `0..degree` generated by `generators`.
    ///
    /// # Errors
    ///
    /// Returns [`PermutationGroupError::GeneratorLength`] or
    /// [`PermutationGroupError::NotAPermutation`] for the first generator that
    /// is not a permutation of `0..degree`.
    pub fn new<I>(degree: usize, generators: I) -> Result<Self, PermutationGroupError>
    where
        I: IntoIterator,
        I::Item: AsRef<[usize]>,
    {
        let mut seen = vec![false; degree];
        let mut pending = Vec::new();
        for (index, generator) in generators.into_iter().enumerate() {
            let generator = generator.as_ref();
            if generator.len() != degree {
                return Err(PermutationGroupError::GeneratorLength {
                    generator: index,
                    length: generator.len(),
                    degree,
                });
            }
            if !is_permutation(generator, &mut seen) {
                return Err(PermutationGroupError::NotAPermutation { generator: index, degree });
            }
            pending.push((0, generator.to_vec()));
        }
        Ok(Self::build(degree, pending, None))
    }

    /// Builds on `active`, discarding residues that fix every active point.
    fn build(
        degree: usize,
        mut pending: Vec<(usize, Vec<usize>)>,
        active: Option<&[bool]>,
    ) -> Self {
        pending.reverse();
        let mut group = Self {
            degree,
            strong_generators: Vec::new(),
            inverses: Vec::new(),
            levels: Vec::new(),
        };
        // Each residue fixes the active points below `start`.
        while let Some((start, mut permutation)) = pending.pop() {
            if let Some(position) = group.sift(&mut permutation, start, active) {
                group.add_strong_generator(start, position, permutation, &mut pending);
            }
        }
        group
    }

    /// Returns the number of points the group acts on.
    #[inline]
    #[must_use]
    pub fn degree(&self) -> usize {
        self.degree
    }

    /// Returns the number of elements, or `None` when it exceeds `u128::MAX`.
    #[must_use]
    pub fn order(&self) -> Option<u128> {
        self.levels.iter().try_fold(1_u128, |order, level| {
            order.checked_mul(u128::try_from(level.orbit.len()).ok()?)
        })
    }

    /// Returns whether `permutation` is an element of the group.
    ///
    /// A slice that is not a permutation of `0..degree` is not an element.
    #[must_use]
    pub fn contains(&self, permutation: &[usize]) -> bool {
        if permutation.len() != self.degree
            || !is_permutation(permutation, &mut vec![false; self.degree])
        {
            return false;
        }
        self.sift(&mut permutation.to_vec(), 0, None).is_none()
    }

    /// Returns the orbit of `point`, sorted ascending.
    ///
    /// # Errors
    ///
    /// Returns [`PermutationGroupError::PointOutOfRange`] when `point` is not
    /// below the degree.
    pub fn orbit(&self, point: usize) -> Result<Vec<usize>, PermutationGroupError> {
        if point >= self.degree {
            return Err(PermutationGroupError::PointOutOfRange { point, degree: self.degree });
        }
        let generators = self.levels.first().map_or(&[][..], |level| &level.generators[..]);
        let mut reached = vec![false; self.degree];
        reached[point] = true;
        let mut orbit = vec![point];
        let mut next = 0;
        while let Some(&current) = orbit.get(next) {
            for &label in generators {
                let image = self.strong_generators[label][current];
                if !core::mem::replace(&mut reached[image], true) {
                    orbit.push(image);
                }
            }
            next += 1;
        }
        orbit.sort_unstable();
        Ok(orbit)
    }

    /// Returns the lexicographically smallest image of `set` under the group,
    /// comparing images as ascending vectors, with an element reaching it.
    ///
    /// The backtrack over the stabilizer chain follows Linton ([2004](https://doi.org/10.1145/1005285.1005319)), greedily moving candidates toward small points and merging equal images.
    /// Orbits disjoint from every candidate image are excluded from the base.
    /// The search is exponential in the worst case.
    ///
    /// # Errors
    ///
    /// Returns [`PermutationGroupError::PointOutOfRange`] for a point not below
    /// the degree and [`PermutationGroupError::DuplicatePoint`] for a point
    /// listed twice.
    pub fn smallest_image(&self, set: &[usize]) -> Result<SmallestImage, PermutationGroupError> {
        let mut seen = vec![false; self.degree];
        for &point in set {
            if point >= self.degree {
                return Err(PermutationGroupError::PointOutOfRange { point, degree: self.degree });
            }
            if core::mem::replace(&mut seen[point], true) {
                return Err(PermutationGroupError::DuplicatePoint { point });
            }
        }
        // Every candidate maps `set` into the image recorded as its key, and
        // all keys agree on the points below `position`.
        let mut candidates = BTreeMap::new();
        let (image, element) = self.greedy_image(0, set, (0..self.degree).collect());
        candidates.insert(image, element);
        let mut group = Cow::Borrowed(self);
        let mut orbit = Vec::new();
        let mut next_level = 0;
        for position in 0..self.degree {
            let placed = candidates
                .keys()
                .next()
                .is_none_or(|image: &Vec<usize>| image.last().is_none_or(|&last| last < position));
            if placed {
                break;
            }
            match group.levels.get(next_level) {
                Some(level) if level.base_point == position => {
                    let hit = candidates
                        .keys()
                        .any(|image| image.iter().any(|&point| level.reaches(point)));
                    if !hit {
                        seen.fill(false);
                        orbit.clear();
                        for &point in candidates.keys().flatten() {
                            if !core::mem::replace(&mut seen[point], true) {
                                orbit.push(point);
                            }
                        }
                        let mut next = 0;
                        while let Some(&point) = orbit.get(next) {
                            for &label in &level.generators {
                                let image = group.strong_generators[label][point];
                                if !core::mem::replace(&mut seen[image], true) {
                                    orbit.push(image);
                                }
                            }
                            next += 1;
                        }
                        let pending = level
                            .generators
                            .iter()
                            .map(|&label| (0, group.strong_generators[label].clone()))
                            .collect();
                        // Its kernel fixes every candidate image.
                        group = Cow::Owned(Self::build(self.degree, pending, Some(&seen)));
                        next_level = 0;
                        continue;
                    }
                    next_level += 1;
                    let mut children = BTreeMap::new();
                    for (image, element) in &candidates {
                        for &target in &level.orbit {
                            if image.binary_search(&target).is_err() {
                                continue;
                            }
                            let mut child = element.clone();
                            level.divide(&group.inverses, &mut child, target);
                            let (image, child) = group.greedy_image(next_level, set, child);
                            children.entry(image).or_insert(child);
                        }
                    }
                    candidates = children;
                }
                _ => {
                    let hit = candidates.keys().any(|image| image.binary_search(&position).is_ok());
                    candidates.retain(|image, _| image.binary_search(&position).is_ok() == hit);
                }
            }
        }
        let (image, element) =
            candidates.pop_first().expect("the identity candidate or one of its children survives");
        Ok(SmallestImage { image, element })
    }

    /// Moves the image of `set` under `element` by the levels from
    /// `first_level` on, mapping at each level the smallest image point in its
    /// orbit onto its base point.
    fn greedy_image(
        &self,
        first_level: usize,
        set: &[usize],
        mut element: Vec<usize>,
    ) -> (Vec<usize>, Vec<usize>) {
        let mut image = Vec::with_capacity(set.len());
        image_of(set, &element, &mut image);
        for level in &self.levels[first_level..] {
            let above = image.partition_point(|&point| point < level.base_point);
            if image.get(above) == Some(&level.base_point) {
                continue;
            }
            if let Some(&target) = image[above..].iter().find(|&&point| level.reaches(point)) {
                level.divide(&self.inverses, &mut element, target);
                image_of(set, &element, &mut image);
            }
        }
        (image, element)
    }

    /// Returns the first active point where sifting leaves a residue.
    fn sift(
        &self,
        permutation: &mut [usize],
        start: usize,
        active: Option<&[bool]>,
    ) -> Option<usize> {
        let first = self.levels.partition_point(|level| level.base_point < start);
        let mut position = start;
        for level in &self.levels[first..] {
            if let Some(moved) = first_moved(permutation, position, level.base_point, active) {
                return Some(moved);
            }
            let image = permutation[level.base_point];
            if !level.reaches(image) {
                return Some(level.base_point);
            }
            level.divide(&self.inverses, permutation, image);
            position = level.base_point + 1;
        }
        first_moved(permutation, position, self.degree, active)
    }

    /// Records a residue fixing every earlier active base point.
    fn add_strong_generator(
        &mut self,
        start: usize,
        position: usize,
        residue: Vec<usize>,
        pending: &mut Vec<(usize, Vec<usize>)>,
    ) {
        let label = self.strong_generators.len();
        self.inverses.push(invert(&residue));
        self.strong_generators.push(residue);
        let first = self.levels.partition_point(|level| level.base_point < start);
        let last = self.levels.partition_point(|level| level.base_point < position);
        if self.levels.get(last).is_none_or(|level| level.base_point != position) {
            // Until now `position` was fixed, so its group was the next
            // level's.
            let generators =
                self.levels.get(last).map(|next| next.generators.clone()).unwrap_or_default();
            let mut tree = vec![UNREACHED; self.degree];
            tree[position] = ROOT;
            self.levels.insert(
                last,
                Level { base_point: position, generators, orbit: vec![position], tree },
            );
        }
        for index in first..=last {
            self.extend_level(index, label, pending);
        }
    }

    /// Adds strong generator `label` to level `index`, extends its orbit, and
    /// queues the Schreier generators of every new orbit point and generator
    /// pair for the next level.
    fn extend_level(&mut self, index: usize, label: usize, pending: &mut Vec<(usize, Vec<usize>)>) {
        let Self { degree, strong_generators, inverses, levels } = self;
        let Level { base_point, generators, orbit, tree } = &mut levels[index];
        generators.push(label);
        let mut scratch = vec![0; *degree];
        let known = orbit.len();
        let mut next = known;
        for position in 0..known {
            let point = orbit[position];
            if let Some(schreier) =
                follow_edge(tree, orbit, strong_generators, inverses, (point, label), &mut scratch)
            {
                pending.push((*base_point + 1, schreier));
            }
        }
        while let Some(&point) = orbit.get(next) {
            for &label in generators.iter() {
                if let Some(schreier) = follow_edge(
                    tree,
                    orbit,
                    strong_generators,
                    inverses,
                    (point, label),
                    &mut scratch,
                ) {
                    pending.push((*base_point + 1, schreier));
                }
            }
            next += 1;
        }
    }
}

/// Replaces `permutation` with `permutation` followed by the inverse of the
/// transversal element of Schreier tree `tree` that reaches `point`.
fn divide(tree: &[usize], inverses: &[Vec<usize>], permutation: &mut [usize], mut point: usize) {
    loop {
        let label = tree[point];
        if label == ROOT {
            return;
        }
        let inverse = &inverses[label];
        for image in permutation.iter_mut() {
            *image = inverse[*image];
        }
        point = inverse[point];
    }
}

/// Follows the edge from `point` along strong generator `label`, growing the
/// Schreier tree or returning the Schreier generator the edge closes when it
/// is not the identity.
fn follow_edge(
    tree: &mut [usize],
    orbit: &mut Vec<usize>,
    strong_generators: &[Vec<usize>],
    inverses: &[Vec<usize>],
    (point, label): (usize, usize),
    scratch: &mut [usize],
) -> Option<Vec<usize>> {
    let generator = &strong_generators[label];
    let image = generator[point];
    if tree[image] == UNREACHED {
        tree[image] = label;
        orbit.push(image);
        return None;
    }
    for (index, entry) in scratch.iter_mut().enumerate() {
        *entry = index;
    }
    divide(tree, inverses, scratch, point);
    // `scratch` is the inverse of the transversal element `u` reaching
    // `point`, so this writes `u` followed by `generator`.
    let mut schreier = vec![0; scratch.len()];
    for (source, &target) in scratch.iter().enumerate() {
        schreier[target] = generator[source];
    }
    divide(tree, inverses, &mut schreier, image);
    schreier.iter().enumerate().any(|(point, &image)| point != image).then_some(schreier)
}

fn is_permutation(permutation: &[usize], seen: &mut [bool]) -> bool {
    seen.fill(false);
    permutation
        .iter()
        .all(|&image| seen.get_mut(image).is_some_and(|seen| !core::mem::replace(seen, true)))
}

fn invert(permutation: &[usize]) -> Vec<usize> {
    let mut inverse = vec![0; permutation.len()];
    for (point, &image) in permutation.iter().enumerate() {
        inverse[image] = point;
    }
    inverse
}

fn first_moved(
    permutation: &[usize],
    start: usize,
    end: usize,
    active: Option<&[bool]>,
) -> Option<usize> {
    (start..end)
        .find(|&point| active.is_none_or(|active| active[point]) && permutation[point] != point)
}

fn image_of(set: &[usize], element: &[usize], image: &mut Vec<usize>) {
    image.clear();
    image.extend(set.iter().map(|&point| element[point]));
    image.sort_unstable();
}