Avoid overflow in passing large exponents to singular

src/sage/libs/singular/polynomial.pxd

@@ -25,7 +25,7 @@ cdef int singular_polynomial_rmul (poly **ret, poly *p, RingElement q, ring *r)
 cdef int singular_polynomial_mul (poly **ret, poly *p, poly *q, ring *r) except -1
 cdef int singular_polynomial_sub (poly **ret, poly *p, poly *q, ring *r) noexcept
 cdef int singular_polynomial_div_coeff (poly **ret, poly *p, poly *q, ring *r) except -1
-cdef int singular_polynomial_pow (poly **ret, poly *p, unsigned long exp, ring *r) except -1
+cdef int singular_polynomial_pow (poly **ret, poly *p, int exp, ring *r) except -1
 cdef int singular_polynomial_neg(poly **ret, poly *p, ring *r) noexcept
 
 cdef object singular_polynomial_latex(poly *p, ring *r, object base, object latex_gens)

src/sage/libs/singular/polynomial.pyx

@@ -271,7 +271,7 @@
     ::
 
         sage: R.<x,y> = Integers(10)[]
         sage: l = [i*x+j*y+k for i in range(10) for j in range(10) for k in range(10)]
         sage: l.sort()
         sage: for i in range(len(l)):
         ....:     for b in l[:i]:
@@ -374,7 +374,7 @@
     sig_off()
     return 0
 
-cdef int singular_polynomial_pow(poly **ret, poly *p, unsigned long exp, ring *r) except -1:
+cdef int singular_polynomial_pow(poly **ret, poly *p, int exp, ring *r) except -1:
     """
     ``ret[0] = p**exp`` where ``p`` in ``r`` and ``exp`` > 0.
 

src/sage/libs/singular/singular.pyx

@@ -1731,7 +1731,6 @@ cdef int overflow_check(unsigned long e, ring *_ring) except -1:
     - ``_ring`` -- a pointer to some ring
 
     Whether an overflow occurs or not partially depends
-
     on the number of variables in the ring. See github issue
     :issue:`11856`. With Singular 4, it is by default optimized
     for at least 4 variables on 64-bit and 2 variables on 32-bit,

src/sage/rings/polynomial/multi_polynomial_libsingular.pyx

@@ -172,6 +172,8 @@
 
 from warnings import warn
 
+from libc.limits cimport INT_MAX
+
 from cpython.object cimport Py_NE
 from cysignals.memory cimport sig_malloc, sig_free
 from cysignals.signals cimport sig_on, sig_off
@@ -2425,6 +2427,54 @@
             Traceback (most recent call last):
             ...
             NotImplementedError: pow() with a modulus is not implemented for this ring
+
+        Check handling of large exponents (:issue:`40442`)::
+
+            sage: P.<x> = PolynomialRing(ZZ, implementation="singular")
+            sage: x^(2^31-1)
+            x^2147483647
+            sage: x^(2^31)
+            x^2147483648
+            sage: x^(2^32-1)
+            x^4294967295
+            sage: x^(2^32)
+            x^4294967296
+            sage: x^(2^63-1)
+            x^9223372036854775807
+            sage: x^(2^63)
+            Traceback (most recent call last):
+            ...
+            OverflowError: exponent overflow (9223372036854775808)
+            sage: x^(2^64-1)
+            Traceback (most recent call last):
+            ...
+            OverflowError: exponent overflow (18446744056529682436)
+            sage: x^(2^64)
+            Traceback (most recent call last):
+            ...
+            OverflowError: exponent overflow (18446744056529682436)
+            sage: P(1)^(2^64)
+            1
+
+            sage: P.<x,y> = ZZ[]
+            sage: y^(2^32-1)
+            y^4294967295
+            sage: y^(2^32)
+            Traceback (most recent call last):
+            ...
+            OverflowError: exponent overflow (4294967296)
+            sage: P(1)^(2^32)
+            1
+
+            sage: P.<x,y,z,t> = ZZ[]
+            sage: y^(2^16-1)
+            y^65535
+            sage: y^(2^16)
+            Traceback (most recent call last):
+            ...
+            OverflowError: exponent overflow (65536)
+            sage: P(1)^(2^16)  # known bug
+            ...UserWarning: error in Singular ignored: exponent 65536 is too large, max. is 65535...
         """
         if mod is not None:
             raise NotImplementedError(
@@ -2451,6 +2501,8 @@
             return self._parent._one_element / (self**(-exp))
         elif exp == 0:
             return self._parent._one_element
+        elif exp > INT_MAX:
+            return (self ** INT_MAX) ** (exp // INT_MAX) * self ** (exp % INT_MAX)
 
         cdef ring *_ring = self._parent_ring
         cdef poly *_p
@@ -5738,7 +5790,7 @@
 
         As such we test the computation is interruptible (previously it wasn't)::
 
             sage: alarm(0.5); h = f.resultant(g, x)
             Traceback (most recent call last):
             ...
             AlarmInterrupt

src/sage/rings/polynomial/plural.pyx

@@ -111,6 +111,7 @@ TESTS::
     True
 """
 from cysignals.memory cimport sig_malloc, sig_free
+from libc.limits cimport INT_MAX
 
 from sage.categories.algebras import Algebras
 from sage.cpython.string cimport char_to_str
@@ -1618,10 +1619,12 @@ cdef class NCPolynomial_plural(RingElement):
             sage: P = A.g_algebra(relations={y*x:-x*y + z},  order='lex')
             sage: P.inject_variables()
             Defining x, z, y
-            sage: (x^2^31) * x^2^31
+            sage: f = x^2^19; f
+            x^524288
+            sage: f*f
             Traceback (most recent call last):
             ...
-            OverflowError: exponent overflow (2147483648)
+            OverflowError: exponent overflow (1048576)
         """
         # all currently implemented rings are commutative
         cdef poly *_p
@@ -1688,10 +1691,30 @@ cdef class NCPolynomial_plural(RingElement):
             sage: P = A.g_algebra(relations={y*x:-x*y + z},  order='lex')
             sage: P.inject_variables()
             Defining x, z, y
+            sage: x^(2^20-1)
+            x^1048575
+            sage: x^2^20
+            Traceback (most recent call last):
+            ....
+            OverflowError: exponent overflow (1048576)
             sage: (x+y^2^31)^10
             Traceback (most recent call last):
             ....
-            OverflowError: exponent overflow (2147483648)
+            OverflowError: exponent overflow (2147483647)
+
+        In the example above, the maximum exponent is `2^{20}-1`,
+        and in the last example, ``2147483647`` instead of ``2147483648`` appears in the
+        error message because of an internal implementation detail.
+        When there are only two variables, the maximum exponent is `2^{32}-1`,
+        however because of a Singular bug, raising a variable directly to the `2^{32}-1`-th
+        power is very slow.
+
+            sage: A.<x,y> = FreeAlgebra(QQ)
+            sage: B.<x,y> = A.g_algebra(relations={y*x: -x*y}, order='lex')
+            sage: (x^65537)^65535
+            x^4294967295
+            sage: x^(2^31)  # known bug (very slow)
+            x^2147483648
 
         Check that using third argument raises an error::
 
@@ -1714,10 +1737,13 @@ cdef class NCPolynomial_plural(RingElement):
             except TypeError:
                 raise TypeError("non-integral exponents not supported")
 
+        # TODO duplication with multi_polynomial_libsingular and sage.structure.element
         if exp < 0:
             return 1/(self**(-exp))
         elif exp == 0:
             return (<NCPolynomialRing_plural>self._parent)._one_element
+        elif exp > INT_MAX:
+            return (self ** INT_MAX) ** (exp // INT_MAX) * self ** (exp % INT_MAX)
 
         cdef ring *_ring = (<NCPolynomialRing_plural>self._parent)._ring
         cdef poly *_p