diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst
index 09b78514e1c..bbb5332a2ec 100644
--- a/src/doc/en/reference/references/index.rst
+++ b/src/doc/en/reference/references/index.rst
@@ -5471,6 +5471,10 @@ REFERENCES:
             Hecke algebras of type A*. J. Algebraic Combin.
             **6** (1997), 59-87.
 
+.. [RRV2022] David P. Roberts and Fernando Rodriguez Villegas,
+     *Hypergeometric motives*, Notices Amer. Math. Soc., **69** vol. 6 (2022).
+     :doi:`10.1090/noti2491`.
+
 .. [RSS] :wikipedia:`Residual_sum_of_squares`, accessed 13th
          October 2009.
 
@@ -5520,6 +5524,9 @@ REFERENCES:
 .. [RV2007] Fernando Rodriguez Villegas. Experimental Number Theory.
             Oxford Graduate Texts in Mathematics 13, 2007.
 
+.. [RV2019] Fernando Rodriguez Villegas. *Mixed Hodge numbers and
+            factorial ratios*. :arxiv:`1907.02722`, 2019.
+
 .. [RW2008] Alexander Raichev and Mark C. Wilson. *Asymptotics of
             coefficients of multivariate generating functions:
             improvements for smooth points*, Electronic Journal of
diff --git a/src/sage/modular/hypergeometric_motive.py b/src/sage/modular/hypergeometric_motive.py
index eb7dfaafe75..49e13394129 100644
--- a/src/sage/modular/hypergeometric_motive.py
+++ b/src/sage/modular/hypergeometric_motive.py
@@ -42,6 +42,8 @@
 
 - [Roberts2015]_
 
+- [RRV2022]_
+
 - [BeCoMe]_
 
 - [Watkins]_
@@ -60,27 +62,31 @@
 
 from collections import defaultdict
 from itertools import combinations
+
 from sage.arith.misc import divisors, gcd, euler_phi, moebius, is_prime
 from sage.arith.misc import gauss_sum, kronecker_symbol
 from sage.combinat.integer_vector_weighted import WeightedIntegerVectors
 from sage.functions.generalized import sgn
 from sage.functions.log import log
 from sage.functions.other import floor, ceil, frac
+from sage.geometry.lattice_polytope import LatticePolytope
+from sage.matrix.constructor import matrix
 from sage.misc.cachefunc import cached_method
 from sage.misc.functional import cyclotomic_polynomial
 from sage.misc.misc_c import prod
 from sage.modular.hypergeometric_misc import hgm_coeffs
-from sage.rings.fraction_field import FractionField
+from sage.modules.free_module_element import vector
+from sage.rings.finite_rings.finite_field_constructor import GF
 from sage.rings.finite_rings.integer_mod_ring import IntegerModRing
+from sage.rings.fraction_field import FractionField
 from sage.rings.integer_ring import ZZ
 from sage.rings.padics.padic_generic_element import gauss_table
 from sage.rings.polynomial.polynomial_ring import polygen, polygens
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from sage.rings.power_series_ring import PowerSeriesRing
 from sage.rings.rational_field import QQ
-from sage.schemes.generic.spec import Spec
-from sage.rings.finite_rings.finite_field_constructor import GF
 from sage.rings.universal_cyclotomic_field import UniversalCyclotomicField
+from sage.schemes.generic.spec import Spec
 
 
 def characteristic_polynomial_from_traces(traces, d, q, i, sign):
@@ -201,7 +207,7 @@ def formule(u):
                 yield H
 
 
-def possible_hypergeometric_data(d, weight=None):
+def possible_hypergeometric_data(d, weight=None) -> list:
     """
     Return the list of possible parameters of hypergeometric motives (up to swapping).
 
@@ -220,7 +226,7 @@ def possible_hypergeometric_data(d, weight=None):
     return list(enumerate_hypergeometric_data(d, weight))
 
 
-def cyclotomic_to_alpha(cyclo):
+def cyclotomic_to_alpha(cyclo) -> list:
     """
     Convert a list of indices of cyclotomic polynomials
     to a list of rational numbers.
@@ -246,14 +252,12 @@ def cyclotomic_to_alpha(cyclo):
         sage: cyclotomic_to_alpha([2,3])
         [1/3, 1/2, 2/3]
     """
-    alpha = []
-    for d in cyclo:
-        for k in ZZ(d).coprime_integers(d):
-            alpha.append(QQ((k, d)))
+    alpha = [QQ((k, d)) for d in cyclo
+             for k in ZZ(d).coprime_integers(d)]
     return sorted(alpha)
 
 
-def alpha_to_cyclotomic(alpha):
+def alpha_to_cyclotomic(alpha) -> list:
     """
     Convert from a list of rationals arguments to a list of integers.
 
@@ -314,10 +318,10 @@ def capital_M(n):
         [1, 4, 27, 64, 3125, 432, 823543]
     """
     n = ZZ(n)
-    return QQ.prod(d ** (d * moebius(n / d)) for d in divisors(n))
+    return QQ.prod(d ** (d * moebius(n // d)) for d in divisors(n))
 
 
-def cyclotomic_to_gamma(cyclo_up, cyclo_down):
+def cyclotomic_to_gamma(cyclo_up, cyclo_down) -> dict:
     """
     Convert a quotient of products of cyclotomic polynomials
     to a quotient of products of polynomials `x^n - 1`.
@@ -392,7 +396,7 @@ def gamma_list_to_cyclotomic(galist):
             sorted(d for d in resu for k in range(-resu[d])))
 
 
-class HypergeometricData():
+class HypergeometricData:
     _gauss_table = {}
 
     def __init__(self, cyclotomic=None, alpha_beta=None, gamma_list=None):
@@ -444,7 +448,7 @@ def __init__(self, cyclotomic=None, alpha_beta=None, gamma_list=None):
             deg = sum(euler_phi(x) for x in cyclo_down)
             up_deg = sum(euler_phi(x) for x in cyclo_up)
             if up_deg != deg:
-                msg = 'not the same degree: {} != {}'.format(up_deg, deg)
+                msg = f'not the same degree: {up_deg} != {deg}'
                 raise ValueError(msg)
             cyclo_up.sort()
             cyclo_down.sort()
@@ -483,7 +487,7 @@ def __init__(self, cyclotomic=None, alpha_beta=None, gamma_list=None):
                                         for v in cyclo_up)
 
     # --- Internals ---
-    def __repr__(self):
+    def __repr__(self) -> str:
         """
         Return the string representation.
 
@@ -498,7 +502,7 @@ def __repr__(self):
         txt = "Hypergeometric data for {} and {}"
         return txt.format(list(self._alpha), list(self._beta))
 
-    def __eq__(self, other):
+    def __eq__(self, other) -> bool:
         """
         Return whether two data are equal.
 
@@ -515,7 +519,7 @@ def __eq__(self, other):
         return (self._alpha == other._alpha and
                 self._beta == other._beta)
 
-    def __ne__(self, other):
+    def __ne__(self, other) -> bool:
         """
         Return whether two data are unequal.
 
@@ -544,7 +548,7 @@ def __hash__(self):
         return hash((self._alpha, self._beta))
 
     # --- Parameters and invariants ---
-    def cyclotomic_data(self):
+    def cyclotomic_data(self) -> tuple:
         """
         Return the pair of tuples of indices of cyclotomic polynomials.
 
@@ -556,7 +560,7 @@ def cyclotomic_data(self):
         """
         return (list(self._cyclo_up), list(self._cyclo_down))
 
-    def alpha_beta(self):
+    def alpha_beta(self) -> tuple:
         """
         Return the pair of lists of rational arguments.
 
@@ -568,7 +572,7 @@ def alpha_beta(self):
         """
         return (list(self._alpha), list(self._beta))
 
-    def alpha(self):
+    def alpha(self) -> list:
         """
         Return the first tuple of rational arguments.
 
@@ -580,7 +584,7 @@ def alpha(self):
         """
         return list(self._alpha)
 
-    def beta(self):
+    def beta(self) -> list:
         """
         Return the second tuple of rational arguments.
 
@@ -592,7 +596,7 @@ def beta(self):
         """
         return list(self._beta)
 
-    def defining_polynomials(self):
+    def defining_polynomials(self) -> tuple:
         """
         Return the pair of products of cyclotomic polynomials.
 
@@ -606,7 +610,7 @@ def defining_polynomials(self):
         down = prod(cyclotomic_polynomial(d) for d in self._cyclo_down)
         return (up, down)
 
-    def gamma_array(self):
+    def gamma_array(self) -> dict:
         r"""
         Return the dictionary `\{v: \gamma_v\}` for the expression
 
@@ -624,7 +628,7 @@ def gamma_array(self):
         """
         return dict(self._gamma_array)
 
-    def gamma_list(self):
+    def gamma_list(self) -> list:
         r"""
         Return a list of integers describing the `x^n - 1` factors.
 
@@ -648,7 +652,7 @@ def gamma_list(self):
             resu += [sgn(n) * v] * abs(n)
         return resu
 
-    def wild_primes(self):
+    def wild_primes(self) -> list:
         r"""
         Return the wild primes.
 
@@ -661,14 +665,15 @@ def wild_primes(self):
             [2, 3, 5]
         """
         gamma = self.gamma_array()
-        return sorted({p for n in gamma.keys() for (p, _) in n.factor()})
+        return sorted({p for n in gamma for p, _ in n.factor()})
 
     def zigzag(self, x, flip_beta=False):
         r"""
         Count ``alpha``'s at most ``x`` minus ``beta``'s at most ``x``.
 
         This function is used to compute the weight and the Hodge numbers.
-        With `flip_beta` set to True, replace each `b` in `\beta` with `1-b`.
+        With `flip_beta` set to ``True``, replace each `b` in `\beta`
+        with `1-b`.
 
         .. SEEALSO::
 
@@ -765,7 +770,7 @@ def degree(self):
         """
         return self._deg
 
-    def hodge_numbers(self):
+    def hodge_numbers(self) -> list:
         """
         Return the Hodge numbers.
 
@@ -876,7 +881,7 @@ def hodge_function(self, x):
             return j
         return j - (i - x) * (k - 1)
 
-    def hodge_polygon_vertices(self):
+    def hodge_polygon_vertices(self) -> list:
         """
         Return the vertices of the Hodge polygon.
 
@@ -942,7 +947,7 @@ def E_polynomial(self, vars=None):
         gamma_plus = [g for g in gamma if g > 0]
         gamma_minus = [g for g in gamma if g < 0]
 
-        domain = set(d for g in gamma for d in divisors(g.abs()))
+        domain = {d for g in gamma for d in divisors(g.abs())}
 
         m_plus = {d: len([1 for g in gamma_plus if not g % d])
                   for d in domain}
@@ -997,7 +1002,7 @@ def M_value(self):
         """
         return self._M_value
 
-    def is_primitive(self):
+    def is_primitive(self) -> bool:
         """
         Return whether this data is primitive.
 
@@ -1037,7 +1042,7 @@ def primitive_index(self):
         """
         return gcd(self.gamma_list())
 
-    def has_symmetry_at_one(self):
+    def has_symmetry_at_one(self) -> bool:
         """
         If ``True``, the motive H(t=1) is a direct sum of two motives.
 
@@ -1117,8 +1122,8 @@ def canonical_scheme(self, t=None):
         gamma_neg = [u for u in self.gamma_list() if u < 0]
         N_pos = len(gamma_pos)
         N_neg = len(gamma_neg)
-        varX = ['X{}'.format(i) for i in range(N_pos)]
-        varY = ['Y{}'.format(i) for i in range(N_neg)]
+        varX = [f'X{i}' for i in range(N_pos)]
+        varY = [f'Y{i}' for i in range(N_neg)]
         ring = PolynomialRing(basering, varX + varY)
         gens = ring.gens()
         X = gens[:N_pos]
@@ -1131,6 +1136,44 @@ def canonical_scheme(self, t=None):
         ideal = ring.ideal([eq0, eq1, self.M_value() * eq2_neg - t * eq2_pos])
         return Spec(ring.quotient(ideal))
 
+    def lattice_polytope(self):
+        """
+        Return the associated lattice polytope.
+
+        This uses the matrix defined in section 3 of [RRV2022]_ and
+        section 3 of [RV2019]_.
+
+        EXAMPLES::
+
+            sage: from sage.modular.hypergeometric_motive import HypergeometricData as Hyp
+            sage: H = Hyp(gamma_list=[-5, -2, 3, 4])
+            sage: P = H.lattice_polytope(); P
+            2-d lattice polytope in 2-d lattice M
+            sage: P.polyhedron().f_vector()
+            (1, 4, 4, 1)
+            sage: len(P.points())
+            7
+
+        The Chebyshev example from [RV2019]_::
+
+            sage: H = Hyp(gamma_list=[-30, -1, 6, 10, 15])
+            sage: P = H.lattice_polytope(); P
+            3-d lattice polytope in 3-d lattice M
+            sage: len(P.points())
+            19
+            sage: P.polyhedron().f_vector()
+            (1, 5, 9, 6, 1)
+        """
+        l = len(self.gamma_list())
+        m = matrix(ZZ, l, 1, self.gamma_list())
+        ext_ker = m.kernel().basis_matrix().insert_row(0, vector(ZZ, [1] * l))
+        unique_relation = ext_ker.kernel().basis()[0]
+        removed = next(i for i, ci in enumerate(unique_relation)
+                       if i and abs(ci) == 1)
+        mat = matrix(ZZ, [v for i, v in enumerate(ext_ker)
+                          if i and i != removed])
+        return LatticePolytope(mat.transpose())
+
     # --- Operations on data ---
     def twist(self):
         r"""
@@ -1264,7 +1307,7 @@ def padic_H_value(self, p, f, t, prec=None, cache_p=False):
         If left unspecified, `prec` is set to the minimum `p`-adic precision
         needed to recover the Euler factor.
 
-        If `cache_p` is True, then the function caches an intermediate
+        If `cache_p` is ``True``, then the function caches an intermediate
         result which depends only on `p` and `f`. This leads to a significant
         speedup when iterating over `t`.
 
@@ -1357,7 +1400,7 @@ def padic_H_value(self, p, f, t, prec=None, cache_p=False):
             raise ValueError('p not prime')
         if not all(x.denominator() % p for x in self._alpha + self._beta):
             raise NotImplementedError('p is wild')
-        if (t.numerator()*t.denominator() % p == 0 or (t-1) % p == 0):
+        if (t.numerator() * t.denominator() % p == 0 or (t - 1) % p == 0):
             raise NotImplementedError('p is tame')
 
         if 0 in alpha:
@@ -1376,7 +1419,7 @@ def padic_H_value(self, p, f, t, prec=None, cache_p=False):
         # also: D = (self.weight() + 1 - m[0]) // 2
 
         if prec is None:
-            prec = ceil((self.weight() * f) / 2 + log(2*self.degree()+1, p))
+            prec = ceil((self.weight() * f) / 2 + log(2 * self.degree() + 1, p))
         use_longs = (p ** prec < 2 ** 31)
 
         gamma = self._gamma_array
@@ -1390,10 +1433,10 @@ def padic_H_value(self, p, f, t, prec=None, cache_p=False):
         else:
             gtab = gauss_table(p, f, prec, use_longs)
             trcoeffs = hgm_coeffs(p, f, prec, gamma, m, D, gtab, prec, use_longs)
-        sigma = trcoeffs[p-2]
+        sigma = trcoeffs[p - 2]
         p_ring = sigma.parent()
-        teich = p_ring.teichmuller(M/t)
-        for i in range(p-3, -1, -1):
+        teich = p_ring.teichmuller(M / t)
+        for i in range(p - 3, -1, -1):
             sigma = sigma * teich + trcoeffs[i]
         resu = ZZ(-1) ** m[0] * sigma / (1 - q)
         return IntegerModRing(p**prec)(resu).lift_centered()
@@ -1498,7 +1541,7 @@ def H_value(self, p, f, t, ring=None):
             raise ValueError('p not prime')
         if not all(x.denominator() % p for x in self._alpha + self._beta):
             raise NotImplementedError('p is wild')
-        if (t.numerator()*t.denominator() % p == 0 or (t-1) % p == 0):
+        if (t.numerator() * t.denominator() % p == 0 or (t - 1) % p == 0):
             raise NotImplementedError('p is tame')
 
         if 0 in alpha:
@@ -1710,7 +1753,7 @@ def euler_factor(self, t, p, cache_p=False):
             raise ValueError('p not prime')
         if not all(x.denominator() % p for x in self._alpha + self._beta):
             raise NotImplementedError('p is wild')
-        if (t.numerator()*t.denominator() % p == 0 or (t-1) % p == 0):
+        if (t.numerator() * t.denominator() % p == 0 or (t - 1) % p == 0):
             raise NotImplementedError('p is tame')
         # now p is good
         d = self.degree()