feanor_math/algorithms/

discrete_log.rs

use std::alloc::Global;
use std::collections::HashMap;
use std::marker::PhantomData;
use std::mem::replace;
use std::rc::Rc;

use feanor_serde::dependent_tuple::DeserializeSeedDependentTuple;
use feanor_serde::impl_deserialize_seed_for_dependent_struct;
use feanor_serde::map::DeserializeSeedMapped;
use feanor_serde::newtype_struct::{DeserializeSeedNewtypeStruct, SerializableNewtypeStruct};
use feanor_serde::seq::{DeserializeSeedSeq, SerializableSeq};
use serde::{Deserialize, Serialize};

use crate::algorithms::eea::{inv_crt, signed_gcd};
use crate::algorithms::int_bisect::root_floor;
use crate::algorithms::int_factor::factor;
use crate::algorithms::linsolve::LinSolveRingStore;
use crate::algorithms::linsolve::smith::{determinant_using_pre_smith, pre_smith};
use crate::algorithms::lll::exact::lll;
use crate::algorithms::matmul::{MatmulAlgorithm, STANDARD_MATMUL};
use crate::divisibility::{DivisibilityRing, DivisibilityRingStore};
use crate::field::FieldStore;
use crate::group::{HashableGroupEl, MultGroup, *};
use crate::homomorphism::Homomorphism;
use crate::integer::{BigIntRing, int_cast};
use crate::iters::multi_cartesian_product;
use crate::matrix::transform::TransformTarget;
use crate::matrix::*;
use crate::ordered::OrderedRingStore;
use crate::pid::PrincipalIdealRingStore;
use crate::primitive_int::StaticRing;
use crate::ring::*;
use crate::rings::finite::FiniteRingStore;
use crate::rings::rational::RationalField;
use crate::rings::zn::{ZnRing, ZnRingStore, zn_big};
use crate::serialization::{DeserializeWithRing, SerializeWithRing};

const ZZ: StaticRing<i64> = StaticRing::<i64>::RING;
const ZZbig: BigIntRing = BigIntRing::RING;

/// Represents a subgroup of an [`AbelianGroupBase`] by a set of generators.
/// Supports computing discrete logarithms, i.e. representing a given element
/// as a combination of the generators.
///
/// Note that the used algorithms have a worst case complexity of `O(sqrt(ord^n))`
/// where `ord` is the given multiple of the orders of each generator, and `n`
/// is the number of generators. However, if `ord` is smooth, much faster algorithms
/// are used.
#[stability::unstable(feature = "enable")]
pub struct SubgroupBase<G: AbelianGroupStore> {
    parent: G,
    generators: Vec<GroupEl<G>>,
    order_multiple: El<BigIntRing>,
    order_factorization: Vec<(i64, usize)>,
    /// the `(i, j)`-th entry has rows that form a basis of the relation lattice of
    /// the set `n/pi^j g1, ..., n/pi^j gk` (where `n` is the order of the group,
    /// and the `pi^ei` are its prime power factors)
    scaled_relation_lattices: Vec<Vec<OwnedMatrix<i64>>>,
    /// the `(i, j, k)`-th entry contains `sum_l row[l] n/pi^(j + 1) gl`, where
    /// `row` is the `k`-th row of `scaled_relation_lattice[i, j]`; These values
    /// are important, since they form a basis of the `p`-torsion subgroup of
    /// `< n/pi^(j + 1) g1, ..., n/pi^(j + 1) gk >`
    scaled_generating_sets: Vec<Vec<Vec<GroupEl<G>>>>,
}

/// [`GroupStore`] of [`SubgroupBase`]
#[stability::unstable(feature = "enable")]
#[allow(type_alias_bounds)]
pub type Subgroup<G: AbelianGroupStore> = GroupValue<SubgroupBase<G>>;

impl<G: AbelianGroupStore> Subgroup<G> {
    /// Creates a new [`GeneratingSet`] representing the subgroup generated
    /// by the given generators.
    ///
    /// The value `order_multiple` should be a multiple of the order of every
    /// generator, including generators that will be added later on via
    /// [`GeneratingSet::add_generator()`].
    #[stability::unstable(feature = "enable")]
    pub fn new(group: G, order_multiple: El<BigIntRing>, generators: Vec<GroupEl<G>>) -> Self {
        let n = generators.len();
        if n == 0 {
            return GroupValue::from(SubgroupBase {
                parent: group,
                generators: Vec::new(),
                order_multiple: ZZbig.clone_el(&order_multiple),
                order_factorization: factor(ZZbig, order_multiple)
                    .into_iter()
                    .map(|(p, e)| (int_cast(p, ZZ, ZZbig), e))
                    .collect(),
                scaled_generating_sets: Vec::new(),
                scaled_relation_lattices: Vec::new(),
            });
        } else {
            let mut result = Self::new(group, order_multiple, Vec::new());
            for g in generators {
                result = result.add_generator(g);
            }
            return result;
        }
    }

    /// Returns the group that this group is a subgroup of.
    #[stability::unstable(feature = "enable")]
    pub fn parent(&self) -> &G { self.get_group().parent() }

    /// Returns the order of the subgroup, i.e. the number of elements.
    #[stability::unstable(feature = "enable")]
    pub fn subgroup_order(&self) -> El<BigIntRing> { self.get_group().subgroup_order() }

    /// Returns the stored generating set of the subgroup.
    #[stability::unstable(feature = "enable")]
    pub fn generators(&self) -> &[GroupEl<G>] { self.get_group().generators() }

    /// Adds a generator to this subgroup, returning a new, larger subgroup.
    #[stability::unstable(feature = "enable")]
    pub fn add_generator(self, new_gen_base: GroupEl<G>) -> Self { Self::from(self.into().add_generator(new_gen_base)) }

    /// Checks whether the given element of the parent group is contained
    /// in the subgroup.
    #[stability::unstable(feature = "enable")]
    pub fn contains(&self, element: &GroupEl<G>) -> bool { self.get_group().contains(element) }

    /// Writes the given element of the parent group as a combination of the
    /// subgroup generators, if this exists.
    #[stability::unstable(feature = "enable")]
    pub fn dlog(&self, target: &GroupEl<G>) -> Option<Vec<i64>> { self.get_group().dlog(target) }

    /// Returns an iterator over all elements of the subgroup.
    #[stability::unstable(feature = "enable")]
    pub fn enumerate_elements<'a>(&'a self) -> impl use<'a, G> + Clone + Iterator<Item = GroupEl<G>> {
        self.get_group().enumerate_elements()
    }
}

impl<G: AbelianGroupStore> SubgroupBase<G> {
    #[stability::unstable(feature = "enable")]
    pub fn parent(&self) -> &G { &self.parent }

    /// The number of elements in the subgroup generated by this generating set.
    #[stability::unstable(feature = "enable")]
    pub fn subgroup_order(&self) -> El<BigIntRing> {
        let mut result = ZZbig.one();
        let n = self.generators.len();
        if n == 0 {
            return result;
        }
        for i in 0..self.order_factorization.len() {
            let (p, e) = self.order_factorization[i];
            let relation_lattice = self.scaled_relation_lattices[i][e].data();
            let Zpne = zn_big::Zn::new(ZZbig, ZZbig.pow(int_cast(p, ZZbig, StaticRing::<i64>::RING), e * n));
            let mod_pne = Zpne.can_hom(&StaticRing::<i64>::RING).unwrap();
            let relation_lattice_det = determinant_using_pre_smith(
                &Zpne,
                OwnedMatrix::from_fn(relation_lattice.row_count(), relation_lattice.col_count(), |k, l| {
                    mod_pne.map(*relation_lattice.at(k, l))
                })
                .data_mut(),
                Global,
            );
            ZZbig.mul_assign(
                &mut result,
                signed_gcd(
                    ZZbig.clone_el(Zpne.modulus()),
                    Zpne.smallest_positive_lift(relation_lattice_det),
                    ZZbig,
                ),
            );
        }
        return result;
    }

    /// Returns a multiple of the order of each element in the subgroup
    /// generated by this generating set.
    #[stability::unstable(feature = "enable")]
    pub fn order_multiple(&self) -> &El<BigIntRing> { &self.order_multiple }

    /// Returns a set of generators of this subgroup.
    #[stability::unstable(feature = "enable")]
    pub fn generators(&self) -> &[GroupEl<G>] { &self.generators }

    /// Extends the generating set by an additional generator, which is likely
    /// to grow the represented subgroup.
    ///
    /// The new generator must be of order dividing [`GeneratingSet::order_multiple()`].
    #[stability::unstable(feature = "enable")]
    pub fn add_generator(self, new_gen_base: GroupEl<G>) -> Self {
        let group = &self.parent;
        assert!(group.is_identity(&group.pow(&new_gen_base, &self.order_multiple)));
        let ZZ_to_ZZbig = ZZbig.can_hom(&ZZ).unwrap();

        let mut scaled_relation_lattices = Vec::new();
        let mut scaled_generating_sets = Vec::new();
        for p_idx in 0..self.order_factorization.len() {
            let (p, e) = self.order_factorization[p_idx];
            let p_bigint = ZZ_to_ZZbig.map(p);
            let power = ZZbig
                .checked_div(&self.order_multiple, &ZZbig.pow(ZZbig.clone_el(&p_bigint), e))
                .unwrap();
            let gens = self.generators.iter().map(|g| group.pow(g, &power)).collect::<Vec<_>>();
            let new_gen = group.pow(&new_gen_base, &power);

            let n = self.generators.len();

            let mut main_relation_matrix = OwnedMatrix::zero(n + 1, n + 1, ZZ);
            for i in 0..n {
                for j in 0..n {
                    *main_relation_matrix.at_mut(i, j) = *self.scaled_relation_lattices[p_idx][e].at(i, j);
                }
            }
            *main_relation_matrix.at_mut(n, n) = -ZZ.pow(p, e);
            for k in 0..e {
                if let Some(dlog) =
                    self.padic_dlog(p_idx, e, &group.pow(&new_gen, &ZZbig.pow(ZZbig.clone_el(&p_bigint), k)))
                {
                    *main_relation_matrix.at_mut(n, n) = -ZZ.pow(p, k);
                    for j in 0..n {
                        *main_relation_matrix.at_mut(n, j) = dlog[j];
                    }
                    break;
                }
            }
            debug_assert!(main_relation_matrix.data().row_iter().all(|row| group.is_identity(
                &(0..n).fold(group.pow(&new_gen, &ZZ_to_ZZbig.map(row[n])), |current, i| {
                    group.op(current, group.pow(&gens[i], &ZZ_to_ZZbig.map(row[i])))
                })
            )));

            let mut result = Vec::with_capacity(e + 1);
            result.push(main_relation_matrix);
            for _ in 0..e {
                result.push(Self::relation_lattice_basis_downscale_p(
                    result.last().unwrap().data(),
                    p,
                ));
            }
            result.reverse();
            scaled_relation_lattices.push(result);

            let mut generating_sets = Vec::new();
            for i in 0..e {
                let generating_set = scaled_relation_lattices.last().unwrap()[i]
                    .data()
                    .row_iter()
                    .map(|row| {
                        let scale = ZZbig.pow(ZZbig.clone_el(&p_bigint), e - i - 1);
                        let result = (0..n).fold(
                            group.pow(&new_gen, &ZZ_to_ZZbig.mul_ref_map(&scale, &row[n])),
                            |current, j| {
                                group.op(current, group.pow(&gens[j], &ZZ_to_ZZbig.mul_ref_map(&scale, &row[j])))
                            },
                        );
                        debug_assert!(group.is_identity(&group.pow(&result, &p_bigint)));
                        result
                    })
                    .collect::<Vec<_>>();
                generating_sets.push(generating_set);
            }
            scaled_generating_sets.push(generating_sets);
        }

        return Self {
            generators: self
                .generators
                .iter()
                .map(|g| group.clone_el(g))
                .chain([new_gen_base])
                .collect(),
            order_multiple: ZZbig.clone_el(&self.order_multiple),
            order_factorization: self.order_factorization.clone(),
            scaled_generating_sets,
            scaled_relation_lattices,
            parent: self.parent,
        };
    }

    /// # Algorithm
    ///  
    /// We are working over `G = ord/p^e global_group`, in which every element
    /// has order dividing `p^e`. Clearly, it is generated by the global generators,
    /// scaled by `ord/p^e`.
    ///
    /// We want to compute a dlog of `x` w.r.t. `g1, ..., gn`. For this, we use the exact sequence
    /// ```text
    ///   0  ->  H  ->  G  ->  G/H  ->  0
    /// ```
    /// where `H = { a in G | p a = 0 }` is the `p`-torsion subgroup. Note that the
    /// power-of-`p` map gives an isomorphism `G/H -> pG`, which allows us to recursively
    /// solve dlog in `G/H`. Hence, we want to solve dlog in `H`, which we can do using
    /// the baby-giant step method - if we can find a generating set of `H`. We find it
    /// using the already provided basis of the relation modulo of the generators.
    fn padic_dlog(&self, p_idx: usize, e: usize, target: &GroupEl<G>) -> Option<Vec<i64>> {
        let group = &self.parent;
        let ZZ_to_ZZbig = ZZbig.can_hom(&ZZ).unwrap();

        let n = self.generators.len();
        if n == 0 {
            return if group.is_identity(target) {
                Some(Vec::new())
            } else {
                None
            };
        } else if e == 0 {
            debug_assert!(group.is_identity(target));
            return Some((0..n).map(|_| 0).collect());
        }

        let p = self.order_factorization[p_idx].0;
        debug_assert!(group.is_identity(&group.pow(target, &ZZbig.pow(ZZ_to_ZZbig.map(p), e))));

        let power = ZZbig
            .checked_div(&self.order_multiple, &ZZbig.pow(int_cast(p, ZZbig, ZZ), e))
            .unwrap();
        let gens = self.generators.iter().map(|g| group.pow(g, &power)).collect::<Vec<_>>();

        // here we use the power-of-`p` map and the fact that `G/H ~ pG` to compute the dlog in
        // `G/H`
        let G_mod_H_dlog = self.padic_dlog(p_idx, e - 1, &group.pow(target, &ZZ_to_ZZbig.map(p)))?;
        debug_assert!(group.eq_el(
            &group.pow(target, &ZZ_to_ZZbig.map(p)),
            &(0..n).fold(group.identity(), |current, i| {
                group.op(current, group.pow(&gens[i], &ZZ_to_ZZbig.map(p * G_mod_H_dlog[i])))
            })
        ));

        // delta is now in H, i.e. is a p-torsion element
        let delta = (0..n).fold(group.clone_el(target), |current, i| {
            group.op(current, group.pow(&gens[i], &ZZ_to_ZZbig.map(-G_mod_H_dlog[i])))
        });
        debug_assert!(group.is_identity(&group.pow(&delta, &ZZ_to_ZZbig.map(p))));

        let H_generators = &self.scaled_generating_sets[p_idx][e - 1];

        let H_dlog_wrt_H_gens = baby_giant_step(
            group,
            delta,
            H_generators,
            &(0..n).map(|_| int_cast(p, ZZbig, ZZ)).collect::<Vec<_>>(),
        )?;
        let H_dlog = {
            let mut result = (0..n).map(|_| 0).collect::<Vec<_>>();
            STANDARD_MATMUL.matmul(
                TransposableSubmatrix::from(Submatrix::from_1d(&H_dlog_wrt_H_gens, 1, n)),
                TransposableSubmatrix::from(self.scaled_relation_lattices[p_idx][e - 1].data()),
                TransposableSubmatrixMut::from(SubmatrixMut::from_1d(&mut result, 1, n)),
                ZZ,
            );
            result
        };

        let result = G_mod_H_dlog
            .into_iter()
            .zip(H_dlog)
            .map(|(x, y)| x + y)
            .collect::<Vec<_>>();
        debug_assert!(group.eq_el(
            target,
            &(0..n).fold(group.identity(), |current, i| {
                group.op(current, group.pow(&gens[i], &ZZ_to_ZZbig.map(result[i])))
            })
        ));

        return Some(result);
    }

    /// Returns `true` if the given element is contained in this subgroup.
    #[stability::unstable(feature = "enable")]
    pub fn contains(&self, element: &GroupEl<G>) -> bool { self.dlog(element).is_some() }

    /// Computes a discrete logarithm of `target` w.r.t. the stored set
    /// if generators, or `None` if `target` is not in the subgroup generated by
    /// these generators
    #[stability::unstable(feature = "enable")]
    pub fn dlog(&self, target: &GroupEl<G>) -> Option<Vec<i64>> {
        let group = &self.parent;
        let ZZ_to_ZZbig = ZZbig.can_hom(&ZZ).unwrap();

        let n = self.generators.len();
        if n == 0 {
            return if group.is_identity(target) {
                Some(Vec::new())
            } else {
                None
            };
        }

        let mut current_dlog = (0..n).map(|_| 0).collect::<Vec<_>>();
        let mut current_order = (0..n).map(|_| 1).collect::<Vec<_>>();

        for p_idx in 0..self.order_factorization.len() {
            let (p, e) = self.order_factorization[p_idx];
            let power = ZZbig
                .checked_div(&self.order_multiple, &ZZbig.pow(int_cast(p, ZZbig, ZZ), e))
                .unwrap();
            let padic_dlog = self.padic_dlog(p_idx, e, &group.pow(target, &power))?;
            for j in 0..n {
                current_dlog[j] = inv_crt(current_dlog[j], padic_dlog[j], &current_order[j], &ZZ.pow(p, e), ZZ);
                current_order[j] *= ZZ.pow(p, e);
            }
        }
        debug_assert!(group.eq_el(
            target,
            &(0..n).fold(group.identity(), |current, i| group.op(
                current,
                group.pow(&self.generators[i], &ZZ_to_ZZbig.map(current_dlog[i]))
            ))
        ));

        return Some(current_dlog);
    }

    fn padic_rectangular_form(&self, p_idx: usize) -> Vec<(GroupEl<G>, usize)> {
        let group = &self.parent;
        let (p, e) = self.order_factorization[p_idx];
        let power = ZZbig
            .checked_div(&self.order_multiple, &ZZbig.pow(int_cast(p, ZZbig, ZZ), e))
            .unwrap();
        let n = self.generators.len();

        if n == 0 {
            return Vec::new();
        }

        let Zpne = zn_big::Zn::new(ZZbig, ZZbig.pow(int_cast(p, ZZbig, StaticRing::<i64>::RING), e * n));
        let mod_pne = Zpne.can_hom(&StaticRing::<i64>::RING).unwrap();
        let relation_lattice = self.scaled_relation_lattices[p_idx][e].data();
        let mut relation_lattice_mod_pne =
            OwnedMatrix::from_fn(relation_lattice.row_count(), relation_lattice.col_count(), |k, l| {
                mod_pne.map(*relation_lattice.at(k, l))
            });
        let mut generators = self.generators.iter().map(|g| group.pow(g, &power)).collect::<Vec<_>>();

        struct TransformGenerators<'a, G: AbelianGroupStore> {
            group: &'a G,
            generators: &'a mut [GroupEl<G>],
        }
        impl<'a, G: AbelianGroupStore> TransformTarget<zn_big::ZnBase<BigIntRing>> for TransformGenerators<'a, G> {
            fn transform<S: Copy + RingStore<Type = zn_big::ZnBase<BigIntRing>>>(
                &mut self,
                ring: S,
                i: usize,
                j: usize,
                transform: &[El<zn_big::Zn<BigIntRing>>; 4],
            ) {
                let transform_inv_det = ring
                    .invert(&ring.sub(
                        ring.mul_ref(&transform[0], &transform[3]),
                        ring.mul_ref(&transform[1], &transform[2]),
                    ))
                    .unwrap();
                let inv_transform = [
                    ring.smallest_positive_lift(ring.mul_ref(&transform[3], &transform_inv_det)),
                    ring.smallest_positive_lift(ring.negate(ring.mul_ref(&transform[1], &transform_inv_det))),
                    ring.smallest_positive_lift(ring.negate(ring.mul_ref(&transform[2], &transform_inv_det))),
                    ring.smallest_positive_lift(ring.mul_ref(&transform[0], &transform_inv_det)),
                ];
                let new_gens = (
                    self.group.op(
                        self.group.pow(&self.generators[i], &inv_transform[0]),
                        self.group.pow(&self.generators[j], &inv_transform[1]),
                    ),
                    self.group.op(
                        self.group.pow(&self.generators[i], &inv_transform[2]),
                        self.group.pow(&self.generators[j], &inv_transform[3]),
                    ),
                );
                self.generators[i] = new_gens.0;
                self.generators[j] = new_gens.1;
            }
        }

        pre_smith(
            &Zpne,
            &mut (),
            &mut TransformGenerators {
                group,
                generators: &mut generators,
            },
            relation_lattice_mod_pne.data_mut(),
        );

        return generators
            .into_iter()
            .enumerate()
            .map(|(i, g)| {
                (
                    g,
                    int_cast(
                        ZZbig.ideal_gen(
                            Zpne.modulus(),
                            &Zpne.smallest_positive_lift(Zpne.clone_el(relation_lattice_mod_pne.at(i, i))),
                        ),
                        ZZ,
                        ZZbig,
                    ) as usize,
                )
            })
            .collect();
    }

    /// Returns a list (g[i], l[i]) such that every element of the subgroup
    /// can be uniquely written as `prod_i g[i]^k[i]` with `0 <= k[i] < l[i]`.
    #[stability::unstable(feature = "enable")]
    pub fn rectangular_form<'a>(&'a self) -> Vec<(GroupEl<G>, usize)> {
        (0..self.order_factorization.len())
            .flat_map(|p_idx| self.padic_rectangular_form(p_idx))
            .collect()
    }

    /// Returns an iterator that yields every element contained in the subgroup
    /// exactly once.
    #[stability::unstable(feature = "enable")]
    pub fn enumerate_elements<'a>(&'a self) -> impl use<'a, G> + Clone + Iterator<Item = GroupEl<G>> {
        let rectangular_form = Rc::new(self.rectangular_form());
        multi_cartesian_product(
            rectangular_form
                .iter()
                .map(|(_, l)| 0..*l)
                .collect::<Vec<_>>()
                .into_iter(),
            move |pows| {
                pows.iter()
                    .enumerate()
                    .fold(self.parent().identity(), |current, (i, e)| {
                        self.parent().op(
                            current,
                            self.parent()
                                .pow(&rectangular_form[i].0, &int_cast(*e as i64, ZZbig, ZZ)),
                        )
                    })
            },
            |_, x| *x,
        )
    }

    /// Takes a matrix whose rows form a basis of the relation lattice w.r.t. group
    /// elements `g1, ..., gn` and computes a matrix whose rows form a basis of the relation
    /// lattice of `p g1, ..., p gn`.
    fn relation_lattice_basis_downscale_p<V>(basis: Submatrix<V, i64>, p: i64) -> OwnedMatrix<i64>
    where
        V: AsPointerToSlice<i64>,
    {
        let n = basis.row_count();
        assert_eq!(n, basis.col_count());

        let QQ = RationalField::new(ZZbig);
        let ZZ_to_QQ = QQ.inclusion().compose(QQ.base_ring().can_hom(&ZZ).unwrap());
        let as_ZZ = |x| int_cast(ZZbig.checked_div(QQ.num(x), QQ.den(x)).unwrap(), ZZ, ZZbig);

        let mut dual_basis = OwnedMatrix::identity(n, 2 * n, &QQ);
        let mut Binv = dual_basis.data_mut().submatrix(0..n, n..(2 * n));
        let mut rhs = OwnedMatrix::identity(n, n, &QQ);
        QQ.solve_right(
            OwnedMatrix::from_fn(n, n, |i, j| ZZ_to_QQ.map_ref(basis.at(i, j))).data_mut(),
            rhs.data_mut(),
            Binv.reborrow(),
        )
        .assert_solved();
        Binv.reborrow()
            .row_iter()
            .flat_map(|row| row.iter_mut())
            .for_each(|x| ZZ_to_QQ.mul_assign_map(x, p));

        let mut identity = OwnedMatrix::identity(n, n, &QQ);
        lll(
            &QQ,
            identity.data(),
            dual_basis.data_mut(),
            &QQ.div(&ZZ_to_QQ.map(9), &ZZ_to_QQ.map(10)),
            Global,
            false,
        );

        let mut result_QQ = rhs;
        QQ.solve_right(
            dual_basis.data_mut().submatrix(0..n, n..(2 * n)),
            identity.data_mut(),
            result_QQ.data_mut(),
        )
        .assert_solved();

        let result = OwnedMatrix::from_fn(n, n, |i, j| as_ZZ(result_QQ.at(i, j)));

        return result;
    }
}

impl<G: AbelianGroupStore> PartialEq for SubgroupBase<G> {
    fn eq(&self, other: &Self) -> bool {
        self.parent().get_group() == other.parent().get_group()
            && other.generators().iter().all(|g| self.contains(g))
            && self.generators().iter().all(|g| other.contains(g))
    }
}

impl<G: AbelianGroupStore> AbelianGroupBase for SubgroupBase<G> {
    type Element = GroupEl<G>;

    fn clone_el(&self, x: &Self::Element) -> Self::Element { self.parent().clone_el(x) }

    fn eq_el(&self, lhs: &Self::Element, rhs: &Self::Element) -> bool { self.parent().eq_el(lhs, rhs) }

    fn hash<H: std::hash::Hasher>(&self, x: &Self::Element, hasher: &mut H) { self.parent().hash(x, hasher) }

    fn identity(&self) -> Self::Element { self.parent().identity() }

    fn inv(&self, x: &Self::Element) -> Self::Element { self.parent().inv(x) }

    fn is_identity(&self, x: &Self::Element) -> bool { self.parent().is_identity(x) }

    fn op(&self, lhs: Self::Element, rhs: Self::Element) -> Self::Element { self.parent().op(lhs, rhs) }

    fn pow(&self, x: &Self::Element, e: &El<BigIntRing>) -> Self::Element { self.parent().pow(x, e) }

    fn op_ref(&self, lhs: &Self::Element, rhs: &Self::Element) -> Self::Element { self.parent().op_ref(lhs, rhs) }

    fn op_ref_snd(&self, lhs: Self::Element, rhs: &Self::Element) -> Self::Element {
        self.parent().op_ref_snd(lhs, rhs)
    }
}

impl<G: AbelianGroupStore + Serialize> Serialize for SubgroupBase<G>
where
    G::Type: SerializableElementGroup,
{
    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
    where
        S: serde::Serializer,
    {
        #[derive(Serialize)]
        struct SubgroupData<'a, Gens: Serialize> {
            order_multiple: SerializeWithRing<'a, BigIntRing>,
            generators: Gens,
            group: (),
        }
        SerializableNewtypeStruct::new(
            "Subgroup",
            (
                self.parent(),
                SubgroupData {
                    order_multiple: SerializeWithRing::new(&self.order_multiple, ZZbig),
                    generators: SerializableSeq::new(
                        self.generators
                            .iter()
                            .map(|g| SerializeWithGroup::new(g, self.parent())),
                    ),
                    group: (),
                },
            ),
        )
        .serialize(serializer)
    }
}

impl<'de, G: AbelianGroupStore + Clone + Deserialize<'de>> Deserialize<'de> for SubgroupBase<G>
where
    G::Type: SerializableElementGroup,
{
    fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
    where
        D: serde::Deserializer<'de>,
    {
        use serde::de::DeserializeSeed;

        struct DeserializeSeedSubgroupData<G: AbelianGroupStore>
        where
            G::Type: SerializableElementGroup,
        {
            group: G,
        }

        impl_deserialize_seed_for_dependent_struct! {
            <{ 'de, G }> pub struct SubgroupData<{'de, G}> using DeserializeSeedSubgroupData<G> {
                order_multiple: El<BigIntRing>: |_| DeserializeWithRing::new(ZZbig),
                generators: Vec<GroupEl<G>>: |master: &DeserializeSeedSubgroupData<G>| {
                    let group_clone = master.group.clone();
                    DeserializeSeedSeq::new((0..).map(move |_| DeserializeWithGroup::new(group_clone.clone())), Vec::new(), |mut current, next| { current.push(next); current })
                },
                group: G: |master: &DeserializeSeedSubgroupData<G>| {
                    let group_clone = master.group.clone();
                    DeserializeSeedMapped::new(PhantomData::<()>, move |()| group_clone)
                }
            } where G: AbelianGroupStore + Clone, G::Type: SerializableElementGroup
        }

        DeserializeSeedNewtypeStruct::new(
            "Subgroup",
            DeserializeSeedDependentTuple::new(PhantomData::<G>, |group| DeserializeSeedSubgroupData { group }),
        )
        .deserialize(deserializer)
        .map(|data| Subgroup::new(data.group, data.order_multiple, data.generators).into())
    }
}

impl<G: AbelianGroupStore> Clone for SubgroupBase<G>
where
    G: Clone,
{
    fn clone(&self) -> Self {
        Self {
            parent: self.parent.clone(),
            generators: self.generators.iter().map(|g| self.parent.clone_el(g)).collect(),
            order_factorization: self.order_factorization.clone(),
            order_multiple: self.order_multiple.clone(),
            scaled_generating_sets: self
                .scaled_generating_sets
                .iter()
                .map(|sets| {
                    sets.iter()
                        .map(|set| set.iter().map(|g| self.parent.clone_el(g)).collect())
                        .collect()
                })
                .collect(),
            scaled_relation_lattices: self
                .scaled_relation_lattices
                .iter()
                .map(|x| x.iter().map(|x| x.clone_matrix(StaticRing::<i64>::RING)).collect())
                .collect(),
        }
    }
}

impl<R> Subgroup<MultGroup<R>>
where
    R: RingStore,
    R::Type: ZnRing + HashableElRing + DivisibilityRing,
{
    #[stability::unstable(feature = "enable")]
    pub fn for_zn(group: MultGroup<R>, generators: Vec<GroupEl<MultGroup<R>>>) -> Self {
        let n = generators.len();
        if n == 0 {
            let n_factorization = factor(ZZbig, group.underlying_ring().size(ZZbig).unwrap());
            let mut order_factorization = n_factorization
                .into_iter()
                .flat_map(|(p, e)| {
                    factor(ZZbig, ZZbig.sub_ref_fst(&p, ZZbig.one()))
                        .into_iter()
                        .chain([(p, e - 1)])
                })
                .collect::<Vec<_>>();
            order_factorization.sort_unstable_by(|(pl, _), (pr, _)| ZZbig.cmp(pl, pr));
            order_factorization.dedup_by(|(p1, e1), (p2, e2)| {
                if ZZbig.eq_el(p1, p2) {
                    *e2 += *e1;
                    true
                } else {
                    false
                }
            });
            order_factorization.retain(|(_, e)| *e > 0);
            let order = ZZbig.prod(
                order_factorization
                    .iter()
                    .map(|(p, e)| ZZbig.pow(ZZbig.clone_el(p), *e)),
            );
            let order_factorization = order_factorization
                .into_iter()
                .map(|(p, e)| (int_cast(p, ZZ, ZZbig), e))
                .collect();
            return Self::from(SubgroupBase {
                parent: group,
                generators: Vec::new(),
                order_multiple: order,
                order_factorization,
                scaled_generating_sets: Vec::new(),
                scaled_relation_lattices: Vec::new(),
            });
        } else {
            let mut result = Self::for_zn(group, Vec::new());
            for g in generators {
                result = result.add_generator(g);
            }
            return result;
        }
    }
}

/// Computes the a vector `k` with entries `1 <= k[i] <= dlog_bounds[i]` such that
/// `generators^k = value` (`generators` is a list of elements of an abelian group).
///
/// If there is no such vector, then `None` is returned. If there are multiple such
/// vectors, any one of them is returned. In the 1d-case, it is guaranteed that this
/// is the smallest one, but in the multidimensional case, no such guarantee can be made
/// (in particular, the vector in general won't be the shortest one w.r.t. any natural
/// ordering like lex or degrevlex).
///
/// Note: The vector `k` is required to have positive entries. In particular, this
/// function won't return the zero vector if the given element is the identity.
/// This can have unexpected consequences, like
/// ```
/// # use feanor_math::algorithms::discrete_log::*;
/// # use feanor_math::integer::*;
/// # use feanor_math::group::*;
/// # use feanor_math::primitive_int::*;
/// let group = AddGroup::new(StaticRing::<i64>::RING);
/// assert_eq!(
///     Some(vec![1]),
///     baby_giant_step(&group, 0, &[0], &[BigIntRing::RING.power_of_two(10)])
/// );
/// ```
///
/// # Implementation notes
///
/// The complexity of the algorithm is `O(sqrt(prod_i dlog_bounds[i]))`.
/// Thus, when possible, `order_bound[i]` should be the order of `generators[i]`
/// in the group.
///
/// Why do we need a group? Because we search for collisions `ab = ac`, and assume
/// that this implies `b = c`. So actually,  a cancelable abelian monoid would be sufficient...
///
/// Why don't we use Pollard's rhos? Because Pollard's rho cannot deterministically
/// detect the case that `value` is not in the subgroup generated by `generators`.
/// It can do so with high probability, but only if the used hash function satisfies
/// certain properties. With BSGS, the correctness does not depend on the used hash
/// function (although performance does, of course).
///
/// # Example
///
/// ```rust
/// # use feanor_math::ring::*;
/// # use feanor_math::group::*;
/// # use feanor_math::homomorphism::*;
/// # use feanor_math::rings::zn::*;
/// # use feanor_math::integer::*;
/// # use feanor_math::primitive_int::StaticRing;
/// # use feanor_math::rings::zn::zn_64::*;
/// # use feanor_math::wrapper::*;
/// # use feanor_math::algorithms::discrete_log::*;
/// let ring = Zn::new(17);
/// let group = MultGroup::new(ring);
/// let x = group.from_ring_el(ring.int_hom().map(9)).unwrap();
/// assert_eq!(
///     Some(vec![3]),
///     baby_giant_step(
///         &group,
///         group.pow(&x, &int_cast(3, BigIntRing::RING, StaticRing::<i64>::RING)),
///         &[x],
///         &[BigIntRing::RING.power_of_two(4)]
///     )
/// );
/// ```
#[stability::unstable(feature = "enable")]
pub fn baby_giant_step<G>(
    group: G,
    value: GroupEl<G>,
    generators: &[GroupEl<G>],
    dlog_bounds: &[El<BigIntRing>],
) -> Option<Vec<i64>>
where
    G: AbelianGroupStore,
{
    let n = generators.len();
    assert_eq!(n, dlog_bounds.len());
    if generators.is_empty() {
        if group.is_identity(&value) {
            return Some(Vec::new());
        } else {
            return None;
        }
    }
    let ns = dlog_bounds
        .iter()
        .map(|n| int_cast(root_floor(ZZbig, n.clone(), 2), ZZ, ZZbig) + 1)
        .collect::<Vec<_>>();
    let count = int_cast(ZZbig.prod(ns.iter().map(|n| int_cast(*n, ZZbig, ZZ))), ZZ, ZZbig);
    let mut baby_step_table: HashMap<HashableGroupEl<_>, i64> = HashMap::with_capacity(count as usize);

    // fill baby step table
    {
        let mut current_els = (0..n).map(|_| group.clone_el(&value)).collect::<Vec<_>>();
        let mut current_idxs = (0..n).map(|_| 0).collect::<Vec<_>>();
        for idx in 0..count {
            _ = baby_step_table.insert(HashableGroupEl::new(&group, group.clone_el(&current_els[n - 1])), idx);

            let mut i = n - 1;
            while current_idxs[i] == ns[i] - 1 {
                if i == 0 {
                    assert!(idx + 1 == count);
                    break;
                }
                current_idxs[i] = 0;
                i -= 1;
            }
            current_idxs[i] += 1;
            current_els[i] = group.op_ref_snd(replace(&mut current_els[i], group.identity()), &generators[i]);
            for j in (i + 1)..n {
                current_els[j] = group.clone_el(&current_els[i]);
            }
        }
    }

    let giant_steps = generators
        .iter()
        .zip(ns.iter())
        .map(|(g, n)| group.pow(g, &int_cast(*n, ZZbig, ZZ)))
        .collect::<Vec<_>>();
    // iterate through giant steps
    {
        let start_el = giant_steps.iter().fold(group.identity(), |x, y| group.op_ref_snd(x, y));
        let mut current_els = (0..n).map(|_| group.clone_el(&start_el)).collect::<Vec<_>>();
        let mut current_idxs = (0..n).map(|_| 1).collect::<Vec<_>>();
        for idx in 0..count {
            if let Some(bs_idx) =
                baby_step_table.get(&HashableGroupEl::new(&group, group.clone_el(&current_els[n - 1])))
            {
                let mut bs_idx = *bs_idx;
                let mut result = current_idxs.clone();
                for j in (0..n).rev() {
                    let bs_idxs_j = bs_idx % ns[j];
                    bs_idx /= ns[j];
                    result[j] = result[j] * ns[j] - bs_idxs_j;
                }
                if (0..dlog_bounds.len()).all(|j| ZZbig.is_leq(&int_cast(result[j], ZZbig, ZZ), &dlog_bounds[j])) {
                    debug_assert_eq!(n, result.len());
                    return Some(result);
                }
            }

            let mut i = n - 1;
            while current_idxs[i] == ns[i] {
                if i == 0 {
                    assert!(idx + 1 == count);
                    break;
                }
                current_idxs[i] = 1;
                i -= 1;
            }
            current_idxs[i] += 1;
            current_els[i] = group.op_ref_snd(replace(&mut current_els[i], group.identity()), &giant_steps[i]);
            for j in (i + 1)..n {
                current_els[j] = group.clone_el(&current_els[i]);
            }
        }
    }

    return None;
}

/// Computes the discrete log in the group `(Z/nZ)*`.
///
/// This only called `finite_field_discrete_log` for backwards compatibility,
/// it now completely supports `(Z/nZ)*` for any `n`.
pub fn finite_field_discrete_log<R: RingStore>(value: El<R>, base: El<R>, Zn: R) -> Option<i64>
where
    R::Type: ZnRing + HashableElRing,
{
    let group = MultGroup::new(Zn);
    let generators = vec![group.from_ring_el(base).unwrap()];
    let subgroup = Subgroup::for_zn(group, generators);
    return subgroup
        .dlog(&subgroup.parent().from_ring_el(value).unwrap())
        .map(|res| res[0]);
}

/// Computes the multiplicative order in the group `(Z/nZ)*`.
pub fn multiplicative_order<R: RingStore>(x: El<R>, Zn: R) -> i64
where
    R::Type: ZnRing + HashableElRing,
{
    let group = MultGroup::new(Zn);
    let gen_set = Subgroup::for_zn(group, Vec::new());
    let Zn = gen_set.parent().underlying_ring();

    let mut result = ZZbig.one();
    for (p, e) in &gen_set.get_group().order_factorization {
        let mut current = Zn.pow_gen(
            Zn.clone_el(&x),
            &ZZbig
                .checked_div(
                    &gen_set.get_group().order_multiple,
                    &ZZbig.pow(int_cast(*p, ZZbig, ZZ), *e),
                )
                .unwrap(),
            ZZbig,
        );
        while !Zn.is_one(&current) {
            current = Zn.pow(current, *p as usize);
            ZZbig.mul_assign(&mut result, int_cast(*p, ZZbig, ZZ));
        }
    }
    return int_cast(result, ZZ, ZZbig);
}

#[cfg(test)]
struct ProdGroupBase<G: AbelianGroupStore, const N: usize>(G);

#[cfg(test)]
impl<G: AbelianGroupStore, const N: usize> PartialEq for ProdGroupBase<G, N> {
    fn eq(&self, other: &Self) -> bool { self.0.get_group() == other.0.get_group() }
}

#[cfg(test)]
impl<G: AbelianGroupStore, const N: usize> AbelianGroupBase for ProdGroupBase<G, N> {
    type Element = [GroupEl<G>; N];

    fn clone_el(&self, x: &Self::Element) -> Self::Element { from_fn(|i| self.0.clone_el(&x[i])) }

    fn eq_el(&self, lhs: &Self::Element, rhs: &Self::Element) -> bool { (0..N).all(|i| self.0.eq_el(&lhs[i], &rhs[i])) }

    fn op(&self, lhs: Self::Element, rhs: Self::Element) -> Self::Element {
        from_fn(|i| self.0.op_ref(&lhs[i], &rhs[i]))
    }

    fn hash<H: std::hash::Hasher>(&self, x: &Self::Element, hasher: &mut H) {
        for i in 0..N {
            self.0.hash(&x[i], hasher)
        }
    }

    fn inv(&self, x: &Self::Element) -> Self::Element { from_fn(|i| self.0.inv(&x[i])) }

    fn identity(&self) -> Self::Element { from_fn(|_| self.0.identity()) }
}

#[cfg(test)]
use std::array::from_fn;

#[cfg(test)]
use oorandom::Rand64;

#[cfg(test)]
use crate::RANDOM_TEST_INSTANCE_COUNT;
#[cfg(test)]
use crate::algorithms::matmul::ComputeInnerProduct;
#[cfg(test)]
use crate::assert_matrix_eq;
#[cfg(test)]
use crate::group::AddGroup;
#[cfg(test)]
use crate::rings::zn::zn_static::Zn;

#[test]
fn test_baby_giant_step() {
    for base_bound in [21, 26, 31, 37] {
        let dlog_bound = [int_cast(base_bound, ZZbig, ZZ)];
        let G = AddGroup::new(ZZ);
        assert_eq!(Some(vec![6]), baby_giant_step(&G, 6, &[1], &dlog_bound));
        assert_eq!(None, baby_giant_step(&G, 0, &[1], &dlog_bound));

        let G = AddGroup::new(Zn::<20>::RING);
        assert_eq!(Some(vec![20]), baby_giant_step(&G, 0, &[1], &dlog_bound));
        assert_eq!(Some(vec![10]), baby_giant_step(&G, 10, &[1], &dlog_bound));
        assert_eq!(Some(vec![5]), baby_giant_step(&G, 0, &[16], &dlog_bound));
    }

    let G = AddGroup::new(ZZ);

    // the collision point is at 96
    assert_eq!(
        Some(vec![9 - 1, 6 - 1]),
        baby_giant_step(&G, 85, &[10, 1], &from_fn::<_, 2, _>(|_| int_cast(8, ZZbig, ZZ)))
    );
    // the collision point is at 105
    assert_eq!(
        Some(vec![10 - 2, 5 - 0]),
        baby_giant_step(&G, 85, &[10, 1], &from_fn::<_, 2, _>(|_| int_cast(21, ZZbig, ZZ)))
    );
    // the collision point is at 90
    assert_eq!(
        Some(vec![6 - 0, 30 - 5]),
        baby_giant_step(&G, 85, &[10, 1], &from_fn::<_, 2, _>(|_| int_cast(31, ZZbig, ZZ)))
    );
}

#[test]
fn test_padic_relation_lattice() {
    let G = AddGroup::new(Zn::<81>::RING);

    let subgroup = Subgroup::new(&G, int_cast(81, ZZbig, ZZ), vec![1]);
    assert_matrix_eq!(ZZ, [[-81]], subgroup.get_group().scaled_relation_lattices[0][4]);
    assert_matrix_eq!(ZZ, [[-27]], subgroup.get_group().scaled_relation_lattices[0][3]);
    assert_matrix_eq!(ZZ, [[-9]], subgroup.get_group().scaled_relation_lattices[0][2]);
    assert_matrix_eq!(ZZ, [[-3]], subgroup.get_group().scaled_relation_lattices[0][1]);
    assert_matrix_eq!(ZZ, [[1]], subgroup.get_group().scaled_relation_lattices[0][0]);

    let subgroup = Subgroup::new(&G, int_cast(81, ZZbig, ZZ), vec![3, 6]);
    let matrix = &subgroup.get_group().scaled_relation_lattices[0][4];
    assert_eq!(-27, *matrix.at(0, 0));
    assert_eq!(-1, *matrix.at(1, 1));
    assert_eq!(
        0,
        ZZ.get_ring()
            .inner_product_ref(matrix.data().row_at(1).iter().zip([3, 6].iter()))
            % 81
    );

    let subgroup = Subgroup::new(&G, int_cast(81, ZZbig, ZZ), vec![3, 9]);
    let matrix = &subgroup.get_group().scaled_relation_lattices[0][4];
    assert_eq!(-27, *matrix.at(0, 0));
    assert_eq!(-1, *matrix.at(1, 1));
    assert_eq!(
        0,
        ZZ.get_ring()
            .inner_product_ref(matrix.data().row_at(1).iter().zip([3, 9].iter()))
            % 81
    );

    let subgroup = Subgroup::new(&G, int_cast(81, ZZbig, ZZ), vec![6, 18, 9]);
    let matrix = &subgroup.get_group().scaled_relation_lattices[0][4];
    assert_eq!(-27, *matrix.at(0, 0));
    assert_eq!(-1, *matrix.at(1, 1));
    assert_eq!(-1, *matrix.at(2, 2));
    assert_eq!(
        0,
        ZZ.get_ring()
            .inner_product_ref(matrix.data().row_at(1).iter().zip([6, 18, 9].iter()))
            % 81
    );
    assert_eq!(
        0,
        ZZ.get_ring()
            .inner_product_ref(matrix.data().row_at(2).iter().zip([6, 18, 9].iter()))
            % 81
    );

    let G = GroupValue::from(ProdGroupBase(AddGroup::new(Zn::<81>::RING)));

    let subgroup = Subgroup::new(&G, int_cast(81, ZZbig, ZZ), vec![[1, 4], [1, 1]]);
    let matrix = &subgroup.get_group().scaled_relation_lattices[0][4];
    assert_eq!(-81, *matrix.at(0, 0));
    assert_eq!(-27, *matrix.at(1, 1));
    assert_eq!(
        0,
        ZZ.get_ring()
            .inner_product_ref(matrix.data().row_at(1).iter().zip([1, 1].iter()))
            % 81
    );
    assert_eq!(
        0,
        ZZ.get_ring()
            .inner_product_ref(matrix.data().row_at(1).iter().zip([4, 1].iter()))
            % 81
    );

    let G = GroupValue::from(ProdGroupBase(AddGroup::new(Zn::<8>::RING)));

    let subgroup = Subgroup::new(&G, int_cast(8, ZZbig, ZZ), vec![[6, 3, 5], [6, 2, 6], [4, 5, 7]]);
    let matrix = &subgroup.get_group().scaled_relation_lattices[0][3];
    assert_eq!(-8, *matrix.at(0, 0));
    assert_eq!(-4, *matrix.at(1, 1));
    assert_eq!(-2, *matrix.at(2, 2));
    assert_eq!(
        0,
        ZZ.get_ring()
            .inner_product_ref(matrix.data().row_at(1).iter().zip([6, 6, 4].iter()))
            % 8
    );
    assert_eq!(
        0,
        ZZ.get_ring()
            .inner_product_ref(matrix.data().row_at(1).iter().zip([3, 2, 5].iter()))
            % 8
    );
    assert_eq!(
        0,
        ZZ.get_ring()
            .inner_product_ref(matrix.data().row_at(1).iter().zip([5, 6, 7].iter()))
            % 8
    );
    assert_eq!(
        0,
        ZZ.get_ring()
            .inner_product_ref(matrix.data().row_at(2).iter().zip([6, 6, 4].iter()))
            % 8
    );
    assert_eq!(
        0,
        ZZ.get_ring()
            .inner_product_ref(matrix.data().row_at(2).iter().zip([3, 2, 5].iter()))
            % 8
    );
    assert_eq!(
        0,
        ZZ.get_ring()
            .inner_product_ref(matrix.data().row_at(2).iter().zip([5, 6, 7].iter()))
            % 8
    );
}

#[test]
fn random_test_dlog() {
    let ring = Zn::<1400>::RING;
    let i = ring.can_hom(&ZZ).unwrap();
    let mut rng = Rand64::new(0);
    let G = GroupValue::from(ProdGroupBase(AddGroup::new(ring)));
    let rand_gs = |rng: &mut Rand64| from_fn::<_, 3, _>(|_| ring.random_element(|| rng.rand_u64()));

    for _ in 0..RANDOM_TEST_INSTANCE_COUNT {
        let gs = from_fn::<_, 3, _>(|_| rand_gs(&mut rng));
        let subgroup = Subgroup::new(&G, int_cast(1400, ZZbig, ZZ), gs.into());

        let coeffs = rand_gs(&mut rng);
        let val = (0..3).fold(G.identity(), |current, i| {
            G.op(current, G.pow(&gs[i], &int_cast(coeffs[i] as i64, ZZbig, ZZ)))
        });
        let dlog = subgroup.dlog(&val);
        println!(
            "{:?} * x + {:?} * y + {:?} * z = {:?} mod 1400",
            gs[0], gs[1], gs[2], val
        );
        if let Some(dlog) = dlog {
            for k in 0..3 {
                assert_el_eq!(
                    ring,
                    val[k],
                    ring.sum([
                        i.mul_map(gs[0][k], dlog[0]),
                        i.mul_map(gs[1][k], dlog[1]),
                        i.mul_map(gs[2][k], dlog[2])
                    ])
                );
            }
            println!("checked solution");
        }
    }

    for _ in 0..RANDOM_TEST_INSTANCE_COUNT {
        let gs = from_fn::<_, 3, _>(|_| rand_gs(&mut rng));
        let subgroup = Subgroup::new(&G, int_cast(1400, ZZbig, ZZ), gs.into());

        let val = rand_gs(&mut rng);
        let dlog = subgroup.dlog(&val);
        println!(
            "{:?} * x + {:?} * y + {:?} * z = {:?} mod 1400",
            gs[0], gs[1], gs[2], val
        );
        if let Some(dlog) = dlog {
            for k in 0..3 {
                assert_el_eq!(
                    ring,
                    val[k],
                    ring.sum([
                        i.mul_map(gs[0][k], dlog[0]),
                        i.mul_map(gs[1][k], dlog[1]),
                        i.mul_map(gs[2][k], dlog[2])
                    ])
                );
            }
            println!("checked solution");
        } else {
            let mut gen_matrix = OwnedMatrix::from_fn(3, 3, |i, j| gs[j][i]);
            let mut value = OwnedMatrix::from_fn(3, 1, |i, _| val[i]);
            let mut res = OwnedMatrix::zero(3, 1, ring);
            let solved = ring.solve_right(gen_matrix.data_mut(), value.data_mut(), res.data_mut());
            println!("[{}, {}, {}]", res.at(0, 0), res.at(1, 0), res.at(2, 0));
            if solved.is_solved() {
                for k in 0..3 {
                    assert_el_eq!(
                        ring,
                        val[k],
                        ring.sum([
                            ring.mul(gs[0][k], *res.at(0, 0)),
                            ring.mul(gs[1][k], *res.at(1, 0)),
                            ring.mul(gs[2][k], *res.at(2, 0))
                        ])
                    );
                }
                assert!(solved == crate::algorithms::linsolve::SolveResult::NoSolution);
            }
            println!("has no solution");
        }
    }
}

#[test]
fn test_zn_dlog() {
    let ring = Zn::<51>::RING;
    let g1 = ring.int_hom().map(37);
    let g2 = ring.int_hom().map(35);

    assert_eq!(0, finite_field_discrete_log(ring.one(), g1, ring).unwrap() % 16);
    assert_eq!(0, finite_field_discrete_log(ring.one(), g2, ring).unwrap() % 2);
    assert_eq!(1, finite_field_discrete_log(g2, g2, ring).unwrap() % 2);
    for i in 0..16 {
        assert_eq!(
            i,
            finite_field_discrete_log(ring.pow(g1, i as usize), g1, ring).unwrap() % 16
        );
    }
    for i in 0..16 {
        assert_eq!(
            None,
            finite_field_discrete_log(ring.mul(ring.pow(g1, i as usize), g2), g1, ring)
        );
    }

    assert_eq!(16, multiplicative_order(g1, ring));
    assert_eq!(8, multiplicative_order(ring.pow(g1, 2), ring));
    assert_eq!(4, multiplicative_order(ring.pow(g1, 4), ring));
    assert_eq!(2, multiplicative_order(ring.pow(g1, 8), ring));
    assert_eq!(2, multiplicative_order(g2, ring));
}

#[test]
fn test_zn_subgroup_size() {
    let ring = Zn::<153>::RING;
    let group = MultGroup::new(ring);
    let g1 = group.from_ring_el(ring.int_hom().map(2)).unwrap();
    let g2 = group.from_ring_el(ring.int_hom().map(37)).unwrap();

    let mut subgroup = Subgroup::for_zn(group.clone(), Vec::new());
    assert_el_eq!(ZZbig, ZZbig.int_hom().map(1), subgroup.subgroup_order());

    let next_gen = subgroup.parent().clone_el(&g1);
    subgroup = subgroup.add_generator(next_gen);
    assert_el_eq!(ZZbig, ZZbig.int_hom().map(24), subgroup.subgroup_order());

    let next_gen = subgroup.parent().pow(&g1, &ZZbig.int_hom().map(2));
    subgroup = subgroup.add_generator(next_gen);
    assert_el_eq!(ZZbig, ZZbig.int_hom().map(24), subgroup.subgroup_order());

    let next_gen = subgroup.parent().clone_el(&g2);
    subgroup = subgroup.add_generator(next_gen);
    assert_el_eq!(ZZbig, ZZbig.int_hom().map(96), subgroup.subgroup_order());

    let generating_set = Subgroup::for_zn(group, vec![g2]);
    assert_el_eq!(ZZbig, ZZbig.int_hom().map(16), generating_set.subgroup_order());
}

#[test]
fn test_enumerate_elements() {
    let ring = Zn::<45>::RING;
    let group = AddGroup::new(ring);

    assert_eq!(
        vec![ring.zero()],
        Subgroup::new(group.clone(), int_cast(45, ZZbig, ZZ), Vec::new())
            .enumerate_elements()
            .collect::<Vec<_>>()
    );

    let subgroup = Subgroup::new(group, int_cast(45, ZZbig, ZZ), vec![9, 15]);
    let mut elements = subgroup.enumerate_elements().collect::<Vec<_>>();
    elements.sort_unstable();
    assert_eq!((0..45).step_by(3).collect::<Vec<_>>(), elements);
}