Trac #34854: molien_series() should not use GAP's VirtualCharacter

One can use GAP's PermutationCharacter() to avoid the troubles.

Currently we have a problem here:
{{{
sage: PG = PermutationGroup(["(1,2,3,4,5,6,7)","(5,6,7)"])
sage: PG.molien_series()
... weird error...
}}}

see the report in https://groups.google.com/d/msgid/sage-devel/42c6be87
-614e-4a21-9273-91948f0c6216n%40googlegroups.com

URL: https://trac.sagemath.org/34854
Reported by: dimpase
Ticket author(s): Dima Pasechnik
Reviewer(s): Travis Scrimshaw

Author: Release Manager

src/sage/groups/perm_gps/permgroup.py

@@ -4665,17 +4665,23 @@ def molien_series(self):
             sage: G = PermutationGroup([[(2,)]])
             sage: G.molien_series()
             1/(x^2 - 2*x + 1)
-        """
-        pi = self._libgap_().NaturalCharacter()
-        # because NaturalCharacter forgets about fixed points:
-        pi += self._libgap_().TrivialCharacter() * len(self.fixed_points())
 
-        # TODO: pi is a Character from a CharacterTable on self, however libgap
-        # does not know about this type and when adding two Characters just
-        # returns a plain List; this needs to be fixed on the libgap side but
-        # in the meantime we can fix by converting pi back to the right type
-        pi = libgap.VirtualCharacter(self._libgap_().CharacterTable(), pi)
+        TESTS:
+
+        Check that :trac:`34854` is fixed::
+
+            sage: PG = PermutationGroup(["(1,2,3,4,5,6,7)","(5,6,7)"])
+            sage: PG.molien_series()
+            (-x^18 + x^15 - x^12 + x^9 - x^6 + x^3 - 1)/(x^25 - x^24 - x^23 - x^22 + x^21 + 2*x^20 + x^19 - x^17 - x^16 - x^15 - x^13 + x^12 + x^10 + x^9 + x^8 - x^6 - 2*x^5 - x^4 + x^3 + x^2 + x - 1)
 
+        and 2 extra fixed points are correctly accounted for::
+
+            sage: PG1 = PermutationGroup(["(9,2,3,4,5,6,7)","(5,6,7)"])
+            sage: R.<x> = QQ[]
+            sage: PG.molien_series() == PG1.molien_series()*(1-x)^2
+            True
+        """
+        pi = self._libgap_().PermutationCharacter(list(self.domain()),libgap.OnPoints)
         M = pi.MolienSeries()
 
         R = QQ['x']