diff --git a/src/sage/groups/abelian_gps/dual_abelian_group_element.py b/src/sage/groups/abelian_gps/dual_abelian_group_element.py
index aaabe107781..18fac135c94 100644
--- a/src/sage/groups/abelian_gps/dual_abelian_group_element.py
+++ b/src/sage/groups/abelian_gps/dual_abelian_group_element.py
@@ -41,7 +41,6 @@
 - Volker Braun (2012-11) port to new Parent base. Use tuples for immutables.
   Default to cyclotomic base ring.
 """
-
 # ****************************************************************************
 #       Copyright (C) 2006 William Stein <wstein@gmail.com>
 #       Copyright (C) 2006 David Joyner<wdjoyner@gmail.com>
@@ -53,62 +52,21 @@
 # (at your option) any later version.
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
-
-import operator
-
 from sage.arith.all import LCM
-from sage.misc.misc_c import prod
 from sage.groups.abelian_gps.element_base import AbelianGroupElementBase
-from functools import reduce
 
 
-def add_strings(x, z=0):
-    """
-    This was in sage.misc.misc but commented out. Needed to add
-    lists of strings in the word_problem method below.
-
-    Return the sum of the elements of x.  If x is empty,
-    return z.
-
-    INPUT:
-
-    - ``x`` -- iterable
-
-    - ``z`` -- the ``0`` that will be returned if ``x`` is empty.
-
-    OUTPUT:
-
-    The sum of the elements of ``x``.
-
-    EXAMPLES::
-
-        sage: from sage.groups.abelian_gps.dual_abelian_group_element import add_strings
-        sage: add_strings([], z='empty')
-        'empty'
-        sage: add_strings(['a', 'b', 'c'])
-        'abc'
-    """
-    if len(x) == 0:
-        return z
-    if not isinstance(x, list):
-        m = iter(x)
-        y = next(m)
-        return reduce(operator.add, m, y)
-    else:
-        return reduce(operator.add, x[1:], x[0])
-
-
-def is_DualAbelianGroupElement(x):
+def is_DualAbelianGroupElement(x) -> bool:
     """
     Test whether ``x`` is a dual Abelian group element.
 
     INPUT:
 
-    - ``x`` -- anything.
+    - ``x`` -- anything
 
     OUTPUT:
 
-    Boolean.
+    Boolean
 
     EXAMPLES::
 
@@ -169,10 +127,10 @@ def __call__(self, g):
         N = LCM(order)
         order_not = [N / o for o in order]
         zeta = F.zeta(N)
-        return F.prod(zeta ** (expsX[i] * expsg[i] * order_not[i])
+        return F.prod(zeta**(expsX[i] * expsg[i] * order_not[i])
                       for i in range(len(expsX)))
 
-    def word_problem(self, words, display=True):
+    def word_problem(self, words):
         """
         This is a rather hackish method and is included for completeness.
 
@@ -196,43 +154,21 @@ def word_problem(self, words, display=True):
             sage: w = a^7*b^3*c^5*d^4*e^4
             sage: x = a^3*b^2*c^2*d^3*e^5
             sage: y = a^2*b^4*c^2*d^4*e^5
-            sage: e.word_problem([u,v,w,x,y],display=False)
+            sage: e.word_problem([u,v,w,x,y])
             [[b^2*c^2*d^3*e^5, 245]]
-
-        The command e.word_problem([u,v,w,x,y],display=True) returns
-        the same list but also prints ``e = (b^2*c^2*d^3*e^5)^245``.
         """
-        ## First convert the problem to one using AbelianGroups
-        import copy
-        from sage.groups.abelian_gps.abelian_group import AbelianGroup
-        from sage.interfaces.gap import gap
-        M = self.parent()
-        G = M.group()
-        gens = M.variable_names()
-        g = prod([G.gen(i)**(self.list()[i]) for i in range(G.ngens())])
-        gap.eval("l:=One(Rationals)")            ## trick needed for LL line below to keep Sage from parsing
-        s1 = "gens := GeneratorsOfGroup(%s)"%G._gap_init_()
-        gap.eval(s1)
-        for i in range(len(gens)):
-            cmd = ("%s := gens["+str(i+1)+"]") % gens[i]
-            gap.eval(cmd)
-        s2 = "g0:=%s; gensH:=%s" % (str(g), words)
-        gap.eval(s2)
-        s3 = 'G:=Group(gens); H:=Group(gensH)'
-        gap.eval(s3)
-        phi = gap.eval("hom:=EpimorphismFromFreeGroup(H)")
-        l1 = gap.eval("ans:=PreImagesRepresentative(hom,g0)")
-        l2 = copy.copy(l1)
-        l4 = []
-        l3 = l1.split("*")
-        for i in range(1,len(words)+1):
-            l2 = l2.replace("x"+str(i),"("+str(words[i-1])+")")
-        l3 = eval(gap.eval("L3:=ExtRepOfObj(ans)"))
-        nn = eval(gap.eval("n:=Int(Length(L3)/2)"))
-        LL1 = eval(gap.eval("L4:=List([l..n],i->L3[2*i])"))         ## note the l not 1
-        LL2 = eval(gap.eval("L5:=List([l..n],i->L3[2*i-1])"))       ## note the l not 1
-        if display:
-            s = str(g)+" = "+add_strings(["("+str(words[LL2[i]-1])+")^"+str(LL1[i])+"*" for i in range(nn)])
-            m = len(s)
-            print("      ", s[:m-1], "\n")
-        return [[words[LL2[i]-1],LL1[i]] for i in range(nn)]
+        from sage.libs.gap.libgap import libgap
+        A = libgap.AbelianGroup(self.parent().gens_orders())
+        gens = A.GeneratorsOfGroup()
+        gap_g = libgap.Product([gi**Li for gi, Li in zip(gens, self.list())])
+        gensH = [libgap.Product([gi**Li for gi, Li in zip(gens, w.list())])
+                 for w in words]
+        H = libgap.Group(gensH)
+
+        hom = H.EpimorphismFromFreeGroup()
+        ans = hom.PreImagesRepresentative(gap_g)
+
+        resu = ans.ExtRepOfObj().sage()  # (indice, power, indice, power, etc)
+        indices = resu[0::2]
+        powers = resu[1::2]
+        return [[words[indi - 1], powi] for indi, powi in zip(indices, powers)]