Trac #16719: replace gap.eval with libgap calls in combinat/combinat.py

As first reported in #15625 (comment:7) with `lucas_number1`, the same
problem is encountered with gap evaluation of Stirling numbers:
{{{
sage: stirling_number1(1000,2)
------------------------------------------------------------------------
---
TypeError                                 Traceback (most recent call
last)
<ipython-input-22-b29f3b3e2003> in <module>()
----> 1 stirling_number1(Integer(1000),Integer(2))

/home/ralf/sage/local/lib/python2.7/site-
packages/sage/combinat/combinat.pyc in stirling_number1(n, k)
    650     """
    651     return
Integer(gap.eval("Stirling1({0},{1})".format(Integer(n),
--> 652
Integer(k))))
    653
    654

/home/ralf/sage/local/lib/python2.7/site-packages/sage/rings/integer.so
in sage.rings.integer.Integer.__init__
(build/cythonized/sage/rings/integer.c:7902)()

TypeError: unable to convert x (=<integer 301...000 (2566 digits)>) to
an integer
}}}
The way to go would be to replace `gap.eval` with `libgap.eval` and it's
faster, too:
{{{
sage: timeit('Integer(gap.eval("Stirling1(100,2)"))')
625 loops, best of 3: 419 µs per loop
sage: timeit('libgap.eval("Stirling1(100,2)").sage()')
625 loops, best of 3: 125 µs per loop
sage: timeit('libgap.eval("Stirling1(1000,2)").sage()')
125 loops, best of 3: 6.45 ms per loop
}}}
Addendum: even faster would be to use `(n-1)!*H(n-1)` for
`stirling(n,2)`, dependent on #16671:
{{{
sage: harmonic_number(99)*factorial(99)==stirling_number1(100,2)
True
sage: timeit('harmonic_number(99)*factorial(99)')
625 loops, best of 3: 56.9 µs per loop
sage: timeit('harmonic_number(999)*factorial(999)')
625 loops, best of 3: 119 µs per loop
}}}

URL: http://trac.sagemath.org/16719
Reported by: rws
Ticket author(s): Ralf Stephan
Reviewer(s): Jeroen Demeyer, Volker Braun

Author: Release Manager

src/sage/combinat/combinat.py

@@ -147,11 +147,10 @@
 #
 #                  http://www.gnu.org/licenses/
 #*****************************************************************************
-from sage.interfaces.all import gap, maxima
-from sage.rings.all import QQ, ZZ, Integer
+from sage.interfaces.all import maxima
+from sage.rings.all import ZZ, QQ, Integer, infinity
 from sage.rings.arith import bernoulli, binomial
 from sage.rings.polynomial.polynomial_element import Polynomial
-from sage.misc.sage_eval import sage_eval
 from sage.libs.all import pari
 from sage.misc.prandom import randint
 from sage.misc.misc import prod
@@ -379,7 +378,7 @@ def bell_number(n, algorithm='flint', **options):
 
     elif algorithm == 'gap':
         from sage.libs.gap.libgap import libgap
-        return libgap.eval("Bell(%s)" % n).sage()
+        return libgap.Bell(n).sage()
 
     elif algorithm == 'dobinski':
         # Hardcode small cases. We only proved the algorithm below
@@ -544,7 +543,8 @@ def fibonacci(n, algorithm="pari"):
     if algorithm == 'pari':
         return ZZ(pari(n).fibonacci())
     elif algorithm == 'gap':
-        return ZZ(gap.eval("Fibonacci({})".format(n)))
+        from sage.libs.gap.libgap import libgap
+        return libgap.Fibonacci(n).sage()
     else:
         raise ValueError("no algorithm {}".format(algorithm))
 
@@ -578,22 +578,15 @@ def lucas_number1(n,P,Q):
         13
         sage: lucas_number1(7,1,-2)
         43
-
-    ::
-
         sage: lucas_number1(5,2,3/5)
         229/25
         sage: lucas_number1(5,2,1.5)
-        Traceback (most recent call last):
-        ...
-        TypeError: no canonical coercion from Real Field with 53 bits of precision to Rational Field
+        1/4
 
     There was a conjecture that the sequence `L_n` defined by
     `L_{n+2} = L_{n+1} + L_n`, `L_1=1`,
     `L_2=3`, has the property that `n` prime implies
-    that `L_n` is prime.
-
-    ::
+    that `L_n` is prime. ::
 
         sage: lucas = lambda n : Integer((5/2)*lucas_number1(n,1,-1)+(1/2)*lucas_number2(n,1,-1))
         sage: [[lucas(n),is_prime(lucas(n)),n+1,is_prime(n+1)] for n in range(15)]
@@ -615,8 +608,9 @@ def lucas_number1(n,P,Q):
 
     Can you use Sage to find a counterexample to the conjecture?
     """
-    ans=gap.eval("Lucas(%s,%s,%s)[1]"%(QQ._coerce_(P),QQ._coerce_(Q),ZZ(n)))
-    return sage_eval(ans)
+    n = ZZ(n);  P = QQ(P);  Q = QQ(Q)
+    from sage.libs.gap.libgap import libgap
+    return libgap.Lucas(P, Q, n)[0].sage()
 
 def lucas_number2(n,P,Q):
     r"""
@@ -662,8 +656,9 @@ def lucas_number2(n,P,Q):
         sage: [lucas_number2(n,1,-1) for n in range(10)]
         [2, 1, 3, 4, 7, 11, 18, 29, 47, 76]
     """
-    ans=gap.eval("Lucas(%s,%s,%s)[2]"%(QQ._coerce_(P),QQ._coerce_(Q),ZZ(n)))
-    return sage_eval(ans)
+    n = ZZ(n);  P = QQ(P);  Q = QQ(Q)
+    from sage.libs.gap.libgap import libgap
+    return libgap.Lucas(P, Q, n)[1].sage()
 
 
 def stirling_number1(n, k):
@@ -687,8 +682,9 @@ def stirling_number1(n, k):
 
     Indeed, `S_1(n,k) = S_1(n-1,k-1) + (n-1)S_1(n-1,k)`.
     """
-    return Integer(gap.eval("Stirling1({0},{1})".format(Integer(n),
-                                                        Integer(k))))
+    n = ZZ(n);  k = ZZ(k)
+    from sage.libs.gap.libgap import libgap
+    return libgap.Stirling1(n, k).sage()
 
 
 def stirling_number2(n, k, algorithm=None):
@@ -809,13 +805,15 @@ def stirling_number2(n, k, algorithm=None):
          Traceback (most recent call last):
          ...
          ValueError: unknown algorithm: CloudReading
-     """
+    """
+    n = ZZ(n);  k = ZZ(k)
     if algorithm is None:
         return _stirling_number2(n, k)
     elif algorithm == 'gap':
-        return ZZ(gap.eval("Stirling2(%s,%s)"%(ZZ(n),ZZ(k))))
+        from sage.libs.gap.libgap import libgap
+        return libgap.Stirling2(n, k).sage()
     elif algorithm == 'maxima':
-        return ZZ(maxima.eval("stirling2(%s,%s)"%(ZZ(n),ZZ(k))))
+        return ZZ(maxima.eval("stirling2(%s,%s)"%(n, k)))
     else:
         raise ValueError("unknown algorithm: %s" % algorithm)
 
@@ -2129,7 +2127,6 @@ def an_element(self):
         return self.f(self.cc.an_element())
 
 ##############################################################################
-from sage.rings.all import infinity
 class InfiniteAbstractCombinatorialClass(CombinatorialClass):
     r"""
     This is an internal class that should not be used directly.  A class which
@@ -2249,7 +2246,7 @@ def tuples(S,k):
             ans.append(y)
     return ans
 
-def number_of_tuples(S,k):
+def number_of_tuples(S, k):
     """
     Return the size of ``tuples(S,k)``. Wraps GAP's ``NrTuples``.
 
@@ -2262,10 +2259,12 @@ def number_of_tuples(S,k):
         sage: number_of_tuples(S,2)
         25
     """
-    ans=gap.eval("NrTuples(%s,%s)"%(S,ZZ(k)))
-    return ZZ(ans)
+    k = ZZ(k)
+    from sage.libs.gap.libgap import libgap
+    S = libgap.eval(str(S))
+    return libgap.NrTuples(S, k).sage()
 
-def unordered_tuples(S,k):
+def unordered_tuples(S, k):
     """
     Return the set of all unordered tuples of length ``k`` of the
     set ``S``. Wraps GAP's ``UnorderedTuples``.
@@ -2296,8 +2295,10 @@ def unordered_tuples(S,k):
         sage: unordered_tuples(["a","b","c"],2)
         ['aa', 'ab', 'ac', 'bb', 'bc', 'cc']
     """
-    ans=gap.eval("UnorderedTuples(%s,%s)"%(S,ZZ(k)))
-    return eval(ans)
+    k = ZZ(k)
+    from sage.libs.gap.libgap import libgap
+    S = libgap.eval(str(S))
+    return libgap.UnorderedTuples(S, k).sage()
 
 def number_of_unordered_tuples(S,k):
     """
@@ -2310,8 +2311,9 @@ def number_of_unordered_tuples(S,k):
         sage: number_of_unordered_tuples(S,2)
         15
     """
-    ans=gap.eval("NrUnorderedTuples(%s,%s)"%(S,ZZ(k)))
-    return ZZ(ans)
+    from sage.libs.gap.libgap import libgap
+    S = libgap.eval(str(S))
+    return libgap.NrUnorderedTuples(S, k).sage()
 
 def unshuffle_iterator(a, one=1):
     r"""

src/sage/libs/gap/element.pyx

@@ -167,7 +167,9 @@ cdef GapElement make_any_gap_element(parent, libGAP_Obj obj):
         sage: type(T_CHAR)
         <type 'sage.libs.gap.element.GapElement_String'>
 
-        sage: t = libgap.eval("UnorderedTuples(['a', 'b', 'c'], 2)"); t
+        sage: libgap.eval("['a', 'b', 'c']")   # gap strings are also lists of chars
+        "abc"
+        sage: t = libgap.UnorderedTuples('abc', 2);  t
         [ "aa", "ab", "ac", "bb", "bc", "cc" ]
         sage: t[1]
         "ab"

src/sage/libs/gap/gap_functions.py

@@ -40,6 +40,7 @@
   'BaseOfGroup',
   'Basis',
   'BasisVectors',
+  'Bell', 
   'Binomial',
   'BlockMatrix',
   'Blocks',