gh-39891: Speed up binomial
    
As first observed in #39379, when
`n` and `k` are small sage integers, `binomial(n, k)` are much slower
than `n.binomial(k)`.

This pull request fixes it. Now `%timeit binomial(20, 10)` is only twice
slower than `%timeit 20.binomial(10)` — previously it was **30** times
slower.

There's also a minor patch where some code is skipped when all arguments
are instances of `Element`. It is certainly a more common case for the
input to be `Element` than otherwise, so we optimize for it.

Another minor patch involves storing imported modules as C global
variables instead of re-import every time.

I think using not-overriding `_method_arguments` as an indicator of
"don't call the method"
(then silently catch `TypeError` and throw it away) is rather bad
design,
but this is not in scope here.

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->

- [x] The title is concise and informative.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [ ] I have updated the documentation and checked the documentation
preview.
    
URL: #39891
Reported by: user202729
Reviewer(s): Travis Scrimshaw

Author: Release Manager

src/sage/arith/misc.py

@@ -3959,7 +3959,7 @@ def binomial(x, m, **kwds):
     P = parent(x)
     x = py_scalar_to_element(x)
 
-    # case 1: native binomial implemented on x
+    # case 1: native binomial implemented on x (see also dont_call_method_on_arg)
     try:
         return P(x.binomial(m, **kwds))
     except (AttributeError, TypeError):

src/sage/functions/other.py

@@ -1715,6 +1715,17 @@ def _binomial_sym(self, n, k):
         from sage.misc.misc_c import prod
         return prod(n - i for i in range(k)) / factorial(k)
 
+    def _method_arguments(self, n, k):
+        """
+        See :meth:`sage.symbolic.function.BuiltinFunction._method_arguments`.
+
+        TESTS::
+
+            sage: binomial._method_arguments(10, 5)
+            (10, 5)
+        """
+        return (n, k)
+
     def _eval_(self, n, k):
         """
         EXAMPLES::

src/sage/functions/transcendental.py

@@ -390,7 +390,7 @@ def _evalf_(self, n, x, parent=None, algorithm=None):
         """
         return _mpmath_utils_call(_mpmath_zeta, x, 1, n, parent=parent)
 
-    def _method_arguments(self, k, x, **args):
+    def _method_arguments(self, k, x):
         r"""
         TESTS::
 

src/sage/symbolic/function.pxd

@@ -3,20 +3,20 @@ from sage.structure.sage_object cimport SageObject
 cdef class Function(SageObject):
     cdef unsigned int _serial
     cdef int _nargs
-    cdef object _name
-    cdef object _alt_name
-    cdef object _latex_name
+    cdef str _name
+    cdef str _alt_name
+    cdef str _latex_name
     cdef object _conversions
     cdef object _evalf_params_first
     cdef _is_registered(self)
     cdef _register_function(self)
 
 cdef class BuiltinFunction(Function):
-    cdef object _preserved_arg
+    cdef int _preserved_arg  # 0 if none, otherwise in [1.._nargs], see function.pyx
     cdef _is_registered(self)
 
 cdef class GinacFunction(BuiltinFunction):
-    cdef object _ginac_name
+    cdef str _ginac_name
     cdef _is_registered(self)
     cdef _register_function(self)
 

src/sage/symbolic/function.pyx

@@ -141,6 +141,7 @@ from sage.structure.sage_object cimport SageObject
 from sage.structure.element cimport Element, parent, Expression
 from sage.misc.lazy_attribute import lazy_attribute
 
+from sage.structure.parent cimport Parent
 from sage.structure.coerce cimport (coercion_model,
         py_scalar_to_element, is_numpy_type, is_mpmath_type)
 from sage.structure.richcmp cimport richcmp
@@ -157,6 +158,7 @@ try:
 except ImportError:
     register_or_update_function = None
 
+cdef object SR = None, PolynomialRing_commutative = None, MPolynomialRing_polydict_domain = None
 
 # List of functions which ginac allows us to define custom behavior for.
 # Changing the order of this list could cause problems unpickling old pickles.
@@ -886,7 +888,7 @@ cdef class BuiltinFunction(Function):
             sage: c(pi/2)                                                               # needs sage.symbolic
             0
         """
-        self._preserved_arg = preserved_arg
+        self._preserved_arg = 0 if preserved_arg is None else preserved_arg
         if preserved_arg and (preserved_arg < 1 or preserved_arg > nargs):
             raise ValueError("preserved_arg must be between 1 and nargs")
 
@@ -913,6 +915,13 @@ cdef class BuiltinFunction(Function):
         univariate functions. Multivariate symbolic functions should override
         it as appropriate.
 
+        Note that it is mandatory for multivariate symbolic functions to
+        override this function in order to allow delegating to the method
+        implemented on the element.
+        For example, ``binomial(n, k)`` tries to call ``n.binomial(k)``
+        because :meth:`sage.functions.other.Function_binomial._method_arguments`
+        is implemented.
+
         EXAMPLES::
 
             sage: zeta._method_arguments(1)
@@ -983,7 +992,7 @@ cdef class BuiltinFunction(Function):
             'foo'
         """
         res = None
-        if args and not hold:
+        if args and not hold and not all(isinstance(arg, Element) for arg in args):
             # try calling the relevant math, cmath, mpmath or numpy function.
             # And as a fallback try the custom self._eval_numpy_ or
             # self._eval_mpmath_
@@ -1046,17 +1055,20 @@ cdef class BuiltinFunction(Function):
                 res = super().__call__(
                         *args, coerce=coerce, hold=hold)
 
-        # Convert the output back to the corresponding
-        # Python type if possible.
+        cdef Parent arg_parent
         if any(isinstance(x, Element) for x in args):
             if (self._preserved_arg
                     and isinstance(args[self._preserved_arg-1], Element)):
                 arg_parent = parent(args[self._preserved_arg-1])
-                from sage.symbolic.ring import SR
+                global SR
+                if SR is None:
+                    from sage.symbolic.ring import SR
                 if arg_parent is SR:
                     return res
-                from sage.rings.polynomial.polynomial_ring import PolynomialRing_commutative
-                from sage.rings.polynomial.multi_polynomial_ring import MPolynomialRing_polydict_domain
+                global PolynomialRing_commutative, MPolynomialRing_polydict_domain
+                if PolynomialRing_commutative is None:
+                    from sage.rings.polynomial.polynomial_ring import PolynomialRing_commutative
+                    from sage.rings.polynomial.multi_polynomial_ring import MPolynomialRing_polydict_domain
                 if isinstance(arg_parent, (PolynomialRing_commutative,
                                            MPolynomialRing_polydict_domain)):
                     try:
@@ -1072,6 +1084,8 @@ cdef class BuiltinFunction(Function):
         if not isinstance(res, Element):
             return res
 
+        # Convert the output back to the corresponding
+        # Python type if possible.
         p = res.parent()
         from sage.rings.complex_double import CDF
         from sage.rings.integer_ring import ZZ