gh-35270: construct AdditiveAbelianGroupWrapper from (not necessarily…

… independent) generating set

    
Currently, constructing an `AdditiveAbelianGroupWrapper` requires prior
knowledge of an *independent* set of generators for the group. This
patch adds an algorithm (implemented following Sutherland's thesis) to
construct an `AdditiveAbelianGroupWrapper` from an arbitrary generating
set, whose members need not be independent. This is useful for all kinds
of things, including structure computations for finite abelian groups.
    
URL: #35270
Reported by: Lorenz Panny
Reviewer(s): Lorenz Panny, Travis Scrimshaw

src/doc/en/reference/references/index.rst

@@ -5889,6 +5889,11 @@ REFERENCES:
              isogenies of prime degree. Journal de Théorie des Nombres de Bordeaux,
              2012.
 
+.. [Suth2007] Andrew V. Sutherland, *Order Computations in Generic Groups*.
+              PhD Thesis, Massachusetts Institute of Technology,
+              June 2007.
+              https://math.mit.edu/~drew/sutherland-phd.pdf
+
 .. [Suth2008] Andrew V. Sutherland, *Structure computation and discrete
               logarithms in finite abelian p-groups*.
               Mathematics of Computation **80** (2011), pp. 477-500.

src/sage/groups/additive_abelian/additive_abelian_wrapper.py

@@ -50,6 +50,7 @@
 
 - David Loeffler (2010)
 - Lorenz Panny (2017): :meth:`AdditiveAbelianGroupWrapper.discrete_log`
+- Lorenz Panny (2023): :meth:`AdditiveAbelianGroupWrapper.from_generators`
 """
 
 # ****************************************************************************
@@ -71,6 +72,7 @@
 from sage.rings.integer_ring import ZZ
 from sage.categories.morphism import Morphism
 from sage.structure.element import parent
+from sage.structure.sequence import Sequence
 from sage.modules.free_module_element import vector
 
 from sage.misc.superseded import deprecated_function_alias
@@ -250,6 +252,29 @@ def _repr_(self):
         """
         return addgp.AdditiveAbelianGroup_fixed_gens._repr_(self) + " embedded in " + self.universe()._repr_()
 
+    def _element_constructor_(self, x, check=False):
+        r"""
+        Create an element from ``x``.
+
+        This may be either an element of self, an element of the ambient
+        group, or an iterable (in which case the result is the corresponding
+        product of the generators of self).
+
+        EXAMPLES::
+
+            sage: V = Zmod(8)**2
+            sage: G = AdditiveAbelianGroupWrapper(V, [[2,2],[4,0]], [4, 2])
+            sage: G(V([6,2]))
+            (6, 2)
+            sage: G([1,1])
+            (6, 2)
+            sage: G(G([1,1]))
+            (6, 2)
+        """
+        if parent(x) is self.universe():
+            return self.element_class(self, self.discrete_log(x), element=x)
+        return addgp.AdditiveAbelianGroup_fixed_gens._element_constructor_(self, x, check)
+
     def discrete_exp(self, v):
         r"""
         Given a list (or other iterable) of length equal to the number of
@@ -450,26 +475,53 @@ def torsion_subgroup(self, n=None):
                 ords.append(d)
         return AdditiveAbelianGroupWrapper(self.universe(), gens, ords)
 
-    def _element_constructor_(self, x, check=False):
+    @staticmethod
+    def from_generators(gens, universe=None):
         r"""
-        Create an element from x. This may be either an element of self, an element of the
-        ambient group, or an iterable (in which case the result is the corresponding
-        product of the generators of self).
+        This method constructs the subgroup generated by a sequence
+        of *finite-order* elements in an additive abelian group.
+
+        The elements need not be independent, hence this can be used
+        to perform tasks such as finding relations between some given
+        elements of an abelian group, computing the structure of the
+        generated subgroup, enumerating all elements of the subgroup,
+        and solving discrete-logarithm problems.
 
         EXAMPLES::
 
-            sage: V = Zmod(8)**2
-            sage: G = AdditiveAbelianGroupWrapper(V, [[2,2],[4,0]], [4, 2])
-            sage: G(V([6,2]))
-            (6, 2)
-            sage: G([1,1])
-            (6, 2)
-            sage: G(G([1,1]))
-            (6, 2)
+            sage: G = AdditiveAbelianGroup([15, 30, 45])
+            sage: gs = [G((1,2,3)), G((4,5,6)), G((7,7,7)), G((3,2,1))]
+            sage: H = AdditiveAbelianGroupWrapper.from_generators(gs); H
+            Additive abelian group isomorphic to Z/90 + Z/15 embedded in Additive abelian group isomorphic to Z/15 + Z/30 + Z/45
+            sage: H.gens()
+            ((12, 13, 14), (1, 26, 21))
+
+        TESTS:
+
+        Random testing::
+
+            sage: invs = []
+            sage: while not 1 < prod(invs) < 10^4:
+            ....:     invs = [randrange(1,100) for _ in range(randrange(1,20))]
+            sage: G = AdditiveAbelianGroup(invs)
+            sage: gs = [G.random_element() for _ in range(randrange(1,10))]
+            sage: H = AdditiveAbelianGroupWrapper.from_generators(gs)
+            sage: os = H.generator_orders()
+            sage: vecs = cartesian_product_iterator(list(map(range, os)))
+            sage: els = {sum(i*g for i,g in zip(vec, H.gens())) for vec in vecs}
+            sage: len(els) == prod(os)
+            True
         """
-        if parent(x) is self.universe():
-            return self.element_class(self, self.discrete_log(x), element=x)
-        return addgp.AdditiveAbelianGroup_fixed_gens._element_constructor_(self, x, check)
+        if not gens:
+            if universe is None:
+                raise ValueError('need universe if no generators are given')
+            return AdditiveAbelianGroupWrapper(universe, [], [])
+
+        if universe is None:
+            universe = Sequence(gens).universe()
+
+        basis, ords = basis_from_generators(gens)
+        return AdditiveAbelianGroupWrapper(universe, basis, ords)
 
 
 def _discrete_log_pgroup(p, vals, aa, b):
@@ -566,3 +618,209 @@ def _rec(j, k, c):
         return x
 
     return _rec(0, max(vals), b)
+
+
+def _expand_basis_pgroup(p, alphas, vals, beta, h, rel):
+    r"""
+    Given a basis of a `p`-subgroup of a finite abelian group
+    and an element lying outside the subgroup, extend the basis
+    to the subgroup spanned jointly by the original subgroup and
+    the new element.
+
+    Used as a subroutine in :func:`basis_from_generators`.
+
+    This function modifies ``alphas`` and ``vals`` in place.
+
+    ALGORITHM: [Suth2007]_, Algorithm 9.2
+
+    INPUT:
+
+    - ``p`` -- prime integer `p`
+    - ``alphas`` -- list; basis for a `p`-subgroup of an abelian group
+    - ``vals`` -- list; valuation at `p` of the orders of the ``alphas``
+    - ``beta`` -- element of the same abelian group as the ``alphas``
+    - ``h`` -- integer; valuation at `p` of the order of ``beta``
+    - ``rel`` -- list of integers; relation on ``alphas + [beta]``
+
+    OUTPUT: basis of the subgroup generated by ``alphas + [beta]``
+
+    EXAMPLES::
+
+        sage: from sage.groups.additive_abelian.additive_abelian_wrapper import _expand_basis_pgroup
+        sage: A = AdditiveAbelianGroup([9,3])
+        sage: alphas = [A((5,2))]
+        sage: beta = A((1,0))
+        sage: vals = [2]
+        sage: rel = next([ZZ(r),ZZ(s)] for s in range(9) for r in range(9) if s > 1 and not r*alphas[0] + s*beta)
+        sage: _expand_basis_pgroup(3, alphas, vals, beta, 2, rel)
+        sage: alphas
+        [(5, 2), (6, 2)]
+        sage: vals
+        [2, 1]
+        sage: len({i*alphas[0] + j*alphas[1] for i in range(3^2) for j in range(3^1)})
+        27
+    """
+    # The given assertions should hold, but were commented out for speed.
+
+    k = len(rel)
+    if not (isinstance(alphas, list) and isinstance(vals, list)):
+        raise TypeError('alphas and vals must be lists for mutability')
+    if not len(alphas) == len(vals) == k - 1:
+        raise ValueError(f'alphas and/or vals have incorrect length')
+#    assert not sum(r*a for r,a in zip(rel, alphas+[beta]))
+#    assert all(a.order() == p**v for a,v in zip(alphas,vals))
+
+    if rel[-1] < 0:
+        raise ValueError('rel must have nonnegative entries')
+
+    # step 1
+    min_r = rel[-1] or float('inf')
+    for i in range(k-1):
+        if not rel[i]:
+            continue
+        if rel[i] < 0:
+            raise ValueError('rel must have nonnegative entries')
+        q = rel[i].p_primary_part(p)
+        alphas[i] *= rel[i] // q
+        rel[i] = q
+        if q < min_r:
+            min_r = q
+    if min_r == float('inf'):
+        raise ValueError('rel must have at least one nonzero entry')
+    val_rlast = rel[-1].valuation(p)
+#    assert rel[-1] == p ** val_rlast
+#    assert not sum(r*a for r,a in zip(rel, alphas+[beta]))
+
+    # step 2
+    if rel[-1] == min_r:
+        for i in range(k-1):
+            beta += alphas[i] * (rel[i]//rel[-1])
+        alphas.append(beta)
+        vals.append(val_rlast)
+#        assert alphas[-1].order() == p**vals[-1]
+        return
+
+    # step 3
+    j = next(j for j,r in enumerate(rel) if r == min_r)
+    alphas[j] = sum(a * (r//rel[j]) for a,r in zip(alphas+[beta], rel))
+
+    # step 4
+    if not alphas[j]:
+        del alphas[j], vals[j]
+        if not alphas:
+            alphas.append(beta)
+            vals.append(val_rlast)
+#            assert alphas[-1].order() == p**vals[-1]
+            return
+
+    # step 5
+    beta_q = beta
+    for v in range(1, h):
+        beta_q *= p
+        try:
+            e = _discrete_log_pgroup(p, vals, alphas, -beta_q)
+        except TypeError:
+            continue
+        # step 6
+        _expand_basis_pgroup(p, alphas, vals, beta, h, list(e) + [p**v])
+        break
+    else:
+        alphas.append(beta)
+        vals.append(h)
+#    assert alphas[-1].order() == p**vals[-1]
+
+def basis_from_generators(gens, ords=None):
+    r"""
+    Given a generating set of some finite abelian group
+    (additively written), compute and return a basis of
+    the group.
+
+    .. NOTE::
+
+        A *basis* of a finite abelian group is a generating
+        set `\{g_1, \ldots, g_n\}` such that each element of the
+        group can be written as a unique linear combination
+        `\alpha_1 g_1 + \cdots + \alpha_n g_n` with each
+        `\alpha_i \in \{0, \ldots, \mathrm{ord}(g_i)-1\}`.
+
+    ALGORITHM: [Suth2007]_, Algorithm 9.1 & Remark 9.1
+
+    EXAMPLES::
+
+        sage: from sage.groups.additive_abelian.additive_abelian_wrapper import basis_from_generators
+        sage: E = EllipticCurve(GF(31337^6,'a'), j=37)
+        sage: E.order()
+        946988065073788930380545280
+        sage: (R,S), (ordR,ordS) = basis_from_generators(E.gens())
+        sage: ordR, ordS
+        (313157428926517503432720, 3024)
+        sage: R.order() == ordR
+        True
+        sage: S.order() == ordS
+        True
+        sage: ordR * ordS == E.order()
+        True
+        sage: R.weil_pairing(S, ordR).multiplicative_order() == ordS
+        True
+        sage: E.abelian_group().invariants()
+        (3024, 313157428926517503432720)
+    """
+    if not gens:
+        return [], []
+    if ords is None:
+        ords = [g.order() for g in gens]
+
+    from sage.arith.functions import lcm
+    lam = lcm(ords)
+    ps = sorted(lam.prime_factors(), key=lam.valuation)
+
+    gammas = []
+    ms = []
+    for p in ps:
+        pgens = [(o.prime_to_m_part(p) * g, o.p_primary_part(p)) for g, o in zip(gens, ords) if not o % p]
+        assert pgens
+        pgens.sort(key=lambda tup: tup[1])
+
+        alpha, ord_alpha = pgens.pop()
+        vals = [ord_alpha.valuation(p)]
+        alphas = [alpha]
+
+        while pgens:
+            beta, ord_beta = pgens.pop()
+            try:
+                dlog = _discrete_log_pgroup(p, vals, alphas, beta)
+            except TypeError:
+                pass
+            else:
+                continue
+
+            # step 4
+            val_beta = ord_beta.valuation(p)
+            beta_q = beta
+            for v in range(1, val_beta):
+                beta_q *= p
+#                assert beta_q == beta * p**v
+                try:
+                    e = _discrete_log_pgroup(p, vals, alphas, -beta_q)
+                except TypeError:
+                    continue
+                _expand_basis_pgroup(p, alphas, vals, beta, val_beta, list(e) + [p**v])
+#                assert all(a.order() == p**v for a,v in zip(alphas, vals))
+                break
+            else:
+                alphas.append(beta)
+                vals.append(val_beta)
+
+        for i, (v, a) in enumerate(sorted(zip(vals, alphas), reverse=True)):
+            if i < len(gammas):
+                gammas[i] += a
+                ms[i] *= p ** v
+            else:
+                gammas.append(a)
+                ms.append(p ** v)
+
+##    assert len({sum(i*g for i,g in zip(vec,gammas))
+##                for vec in __import__('itertools').product(*map(range,ms))}) \
+##               == __import__('sage').misc.misc_c.prod(ms)
+
+    return gammas, ms

src/sage/schemes/elliptic_curves/ell_finite_field.py

@@ -932,6 +932,10 @@ def abelian_group(self):
         - John Cremona: original implementation
         - Lorenz Panny (2021): current implementation
 
+        .. SEEALSO::
+
+            :meth:`AdditiveAbelianGroupWrapper.from_generators()<sage.groups.additive_abelian.additive_abelian_wrapper.AdditiveAbelianGroupWrapper.from_generators>`
+
         EXAMPLES::
 
             sage: E = EllipticCurve(GF(11),[2,5])