gh-37911: Fix overflow in hypergeometric trace

This PR fixes an overflow bug in the computation of the hypergeometric trace formula (#37910). URL: #37911 Reported by: kedlaya Reviewer(s):
sagemath · May 1, 2024 · 422c59e · 422c59e
2 parents 67f2675 + e3426d7
commit 422c59e
Show file tree

Hide file tree

Showing 2 changed files with 26 additions and 11 deletions.
diff --git a/src/sage/modular/hypergeometric_misc.pyx b/src/sage/modular/hypergeometric_misc.pyx
@@ -4,7 +4,7 @@ significantly from Cythonization.
 """
 from cpython cimport array
 from cysignals.signals cimport sig_check
-
+from sage.rings.integer cimport Integer
 
 cpdef hgm_coeffs(long long p, int f, int prec, gamma, m, int D,
                  gtable, int gtable_prec, bint use_longs):
@@ -37,7 +37,7 @@ cpdef hgm_coeffs(long long p, int f, int prec, gamma, m, int D,
 
     cdef int gl, j, k, l, v, gv
     cdef long long i, q1, w, w1, w2, q2, r, r1
-    cdef bint flip
+    cdef bint flip, use_longlongs
 
     q1 = p ** f - 1
     gl = len(gamma)
@@ -55,6 +55,7 @@ cpdef hgm_coeffs(long long p, int f, int prec, gamma, m, int D,
     # for efficiency.
     flip = (f == 1 and prec == 1 and gtable_prec == 1)
     ans = []
+    Rz = R.zero()
     if use_longs:
         q2 = p ** prec
         try:
@@ -64,13 +65,15 @@ cpdef hgm_coeffs(long long p, int f, int prec, gamma, m, int D,
             for r in range(q1):
                 gtab2[r] = gtable[r].lift() % q2
     else:
-        gtab2 = array.array('q', [0]) * q1
-        try:
-            for r in range(q1):
-                gtab2[r] = gtable[r]
-        except TypeError:
-            for r in range(q1):
-                gtab2[r] = gtable[r].lift()
+        use_longlongs = (Integer(p) ** prec < 2 ** 63)
+        if use_longlongs:
+            gtab2 = array.array('q', [0]) * q1
+            try:
+                for r in range(q1):
+                    gtab2[r] = gtable[r]
+            except TypeError:
+                for r in range(q1):
+                    gtab2[r] = gtable[r].lift()
     if f == 1:
         for r in range(q1):
             digit_count[r] = r
@@ -82,7 +85,6 @@ cpdef hgm_coeffs(long long p, int f, int prec, gamma, m, int D,
                 w += r1 % p
                 r1 //= p
             digit_count[r] = w
-    Rz = R.zero()
     ans = [None] * q1
     for r in range(q1):
         sig_check()
@@ -123,14 +125,21 @@ cpdef hgm_coeffs(long long p, int f, int prec, gamma, m, int D,
                     else:
                         for j in range(-gv):
                             w1 = w1 * w2 % q2
-                else:
+                elif use_longlongs:
                     w2 = gtab2[r1]
                     if gv > 0:
                         for j in range(gv):
                             u *= w2
                     else:
                         for j in range(-gv):
                             u1 *= w2
+                else:
+                    if gv > 0:
+                        for j in range(gv):
+                            u *= gtable[r1]
+                    else:
+                        for j in range(-gv):
+                            u1 *= gtable[r1]
             if use_longs:
                 u = R(w)
                 u1 = R(w1)

diff --git a/src/sage/modular/hypergeometric_motive.py b/src/sage/modular/hypergeometric_motive.py
@@ -1389,6 +1389,12 @@ def padic_H_value(self, p, f, t, prec=None, cache_p=False):
             ....:     print(s)
             p is tame
 
+        Check that :issue:`37910` is resolved::
+
+            sage: H = Hyp(alpha_beta=[[1/2,1/2,1/2,1/2,1/2,1/3,2/3,1/6,5/6], [0,0,0,0,0,0,0,0,0]])
+            sage: H.padic_H_value(151, 2, -512000)
+            50178940126155881
+
         REFERENCES:
 
         - [MagmaHGM]_