LL-functions via deformations: from hyperelliptic curves to hypergeometric motives

Kiran S. Kedlaya, Department of Mathematics, University of California, San Diego; [email protected]
Arithmetic of Hyperelliptic Curves, Abdus Salam International Centre for Theoretical Physics (ICTP), September 8, 2017

About this demo

This file is a Jupyter notebook hosted on CoCalc (formerly SageMathCloud). The demo is running a custom version of Sage including an in-development port of the Magma hypergeometric motives (joint work with Frédéric Chapoton); see Sage's trac server. Some features, such as tame and wild primes, are not yet implemented. (For technical reasons, I've also compiled in some other recent modifications; see tickets 23780 and 23784 on trac.)

Here are the associated slides.

from sage.modular.hypergeometric_motive import HypergeometricData as HGData

Change the following parameter to True/False to enable/disable optional consistency checks. The non-optional code can be found gathered into a single function at the bottom of the notebook (which depends on the user-defined functions appearing along the way).

optional = True

Specify hypergeometric data using tuples α=(α1,,αn),β=(β1,,βn)\alpha = (\alpha_1,\dots,\alpha_n), \beta = (\beta_1,\dots,\beta_n), plus a prime pp. In this demonstration, we insist that β\beta contain no repeated elements nor 0 or 1, and p1p-1 be divisible by the least common denominator of β\beta. (Loosening the second and third conditions should be relatively easy; the first condition may be somewhat harder.)

H = HGData(alpha_beta=([1/6, 5/6, 1/8, 3/8, 5/8, 7/8], [1/3, 2/3, 1/12, 5/12, 7/12, 11/12]))
alpha, beta = H.alpha_beta() ## these come back sorted
p = 37

Set some auxiliary parameters: working precision in p, working precision in t, and desired final precision in p. (I haven't yet checked the calculations to figure out how these parameters should be set in general.)

pprec = 25
tprec = 7*p
e = 7

Define a few global variables.

n = len(alpha)
K = Qp(p, pprec)
R.<t> = LaurentSeriesRing(K, tprec)
n = H.degree() # = len(alpha)

Optional: compute the connection matrix N=(010000000001a0a1an2an1) N = \begin{pmatrix} 0 & -1 & \cdots & 0 & 0 \\ 0 & 0 & & 0 & 0 \\ \vdots & & \ddots & & \vdots \\ 0 & 0 & & 0 & -1 \\ a_0 & a_1 & \cdots & a_{n-2} & a_{n-1} \end{pmatrix} associated to the hypergeometric differential equation 0=1t1(ti=1n(D+αi)i=1n(D+βi1))(y)=Dn+a0Dn1++an1 0 = \frac{1}{t-1}\left( t \prod_{i=1}^n (D + \alpha_i) - \prod_{i=1}^n (D + \beta_i - 1) \right) (y) = D^n + a_0 D^{n-1} + \cdots + a_{n-1} where D=tddtD = t \frac{d}{dt}. (Normalization as per Beukers-Heckman, Monodromy for the hypergeometric function nFn1{}_nF_{n-1}.)

def companion_matrix(l):
    n = len(l)
    matlist = [[0 for j in range(n)] for i in range(n)]
    for i in range(n-1): matlist[i][i+1] = -1
    for i in range(n): matlist[n-1][i] = l[i]
    return(matrix(matlist))

if optional:
    pol2.<D> = PolynomialRing(R)
    diffop = (t*prod(D+a for a in alpha) - prod(D+b-1 for b in beta))/(t-1 + O(t^tprec))
    N = companion_matrix(diffop.list()[:-1])

The Clausen-Thomae hypergeometric series nFn1(α1,,αnβ1,,βn1t)=k=0(α1)k(αn)k(β1)k(βn1)kzkk!,(α)k=α(α+1)(α+k1) {}_nF_{n-1}\left( \left. \genfrac{}{}{0pt}{}{\alpha_1,\dots,\alpha_n}{\beta_1,\dots,\beta_{n-1}} \right| t \right) = \sum_{k=0}^\infty \frac{(\alpha_1)_k \cdots (\alpha_n)_k}{(\beta_1)_k \cdots (\beta_{n-1})_k} \frac{z^k}{k!}, \qquad (\alpha)_k = \alpha(\alpha+1) \cdots (\alpha+k-1) has the property that for i=1,,ni=1,\dots,n, the Puiseux series z1βinFn1(α1βi+1,,αnβi+1β1βi+1,,βiβi+1^,,βn1βi+1t) z^{1-\beta_i} {}_nF_{n-1}\left( \left. \genfrac{}{}{0pt}{}{\alpha_1-\beta_i+1,\dots,\alpha_n-\beta_i+1}{\beta_1-\beta_i+1,\dots,\widehat{\beta_i - \beta_i + 1}, \dots, \beta_{n-1}-\beta_i+1} \right| t \right) (where βiβi+1\beta_i - \beta_i+1 is omitted; it contributes the k!k!) is a solution of the hypergeometric equation. Using this, we construct a formal solution matrix UU (and its inverse) for which U1NU+U1D(U)=N0,N0=Diag(β11,,βn1); U^{-1} N U + U^{-1} D(U) = N_0, \qquad N_0 = \mathrm{Diag}(\beta_1-1, \dots, \beta_n - 1); namely, the ii-th column of UU consists of y,(D+1βi)(y),,(D+1βi)n1(y)y,(D+1-\beta_i)(y),\dots,(D+1-\beta_i)^{n-1}(y) where yy is the ii-th solution from above with the factor of z1βiz^{1-\beta_i} removed (because I don't want to deal with fractional exponents).

from sage.rings.padics.misc import gauss_sum as padic_gauss_sum

def D(x,t):
    return t*x.derivative()

def formal_solution_matrix(H, R, tprec, K=QQ):
    alpha, beta = H.alpha_beta()
    n = H.degree()
    l = []
    for b in beta:
        alpha2 = [a2-b+1 for a2 in alpha]
        beta2 = [b2-b+1 for b2 in beta if b2 != b]
        tmp = [K(1)]
        for i in range(1,tprec):
            tmp.append(tmp[-1] * (prod(a2+i-1 for a2 in alpha2)/prod(b2+i-1 for b2 in beta2)/i))
        l.append(R(tmp).add_bigoh(tprec))
#
# Aside: I could have done
#    l.append(ser(hypergeometric(alpha2, beta2, var('x')).series(var('x'), tprec).list()).add_bigoh(tprec))
# but this would have computed the series over Q, which is suboptimal for large p.
#
    matlist = [l]
    t = R.gen()
    for i in range(n-1):
        matlist.append([(1-beta[i])*matlist[-1][i] + D(matlist[-1][i],t) for i in range(n)])
    return(matrix(matlist))

U = formal_solution_matrix(H,R,tprec,K)
Uinv = U.inverse()

Consistency check: test that we have correctly constructed the formal solution matrix by verifying the equation U1NU+U1D(U)=N0. U^{-1} N U + U^{-1} D(U) = N_0.

if optional:
    test = Uinv * (N*U+U.apply_map(lambda x: D(x,t)))
    print([test[i][i] == (beta[i] - 1) for i in range(n)])
[True, True, True, True, True, True]

Compute the (conjectural) Frobenius matrix F0F_0 on the formal solution basis. Note that this matrix must satisfy the commutation relation N0F0pF0σ(N0)+D(F0)=0. N_0F_0 - pF_0 \sigma(N_0) + D(F_0) = 0. Because the βj\beta_j are distinct and p1p-1 is divisible by all denominators, F0F_0 is forced to be diagonal, and the (j,j)(j,j)-entry must be a scalar multiple of tp+1+pβjt^{-p+1+\lfloor p\beta_j \rfloor}. My conjecture for the scalar, based on the (conjectural) formula for the hypergeometric Euler factor at a prime of bad reduction, is (up to an overall twist) ωp(M(p1)(1βi))vgp((p1)(1βi)v)γv \omega_p(M^{(p-1)(1-\beta_i)}) \prod_{v} g_p((p-1)(1-\beta_i) v)^{\gamma_v} where ωp\omega_p is the Teichmüller character, MM and vγvv \mapsto \gamma_v are invariants of the hypergeometric data, and gp(a)=uFpωp(u)aζpu g_p(a) = \sum_{u \in \FF_p^*} \omega_p(u)^{-a} \zeta_p^{u} is a Gauss sum. The latter is computed using the Gross-Koblitz formula: for π\pi the root of πp1=p\pi^{p-1} = -p for which ζp1+π(modπ2)\zeta_p \equiv 1 + \pi \pmod{\pi^2}, gp(a)=πaΓp(a/(p1))(a=0,,p2). g_p(a) = \pi^{-a} \Gamma_p(a/(p-1)) \qquad (a=0,\dots,p-2). Here we rely on the fact that (p1)βjZ(p-1) \beta_j \in \mathbb{Z} for all jj. If this fails, but the βj\beta_j are still distinct, then F0F_0 is not scalar but rather a weighted permutation matrix. Things get much more complicated if the βj\beta_j are not all distinct: for example, for the Dwork quintic family α=(15,25,35,45),β=(0,0,0,0), \alpha = \left( \frac{1}{5}, \frac{2}{5}, \frac{3}{5}, \frac{4}{5} \right), \beta = (0,0,0,0), a calculation of Shapiro shows that F0=(p3000p20000p00001) F_0 = \begin{pmatrix} p^3 & 0 & 0 & * \\ 0 & p^2 & 0 & 0 \\ 0 & 0 & p & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} where * is a particular nonzero multiple of ζ3\zeta_3. It may be possible to guess the off-diagonal terms in general using ideas of Vologodsky; they should involve pp-adic multiple zeta values.

def initial_Frobenius_matrix(H,t,p,pprec=20):
    alpha, beta = H.alpha_beta()
    n = H.degree()
    gauss_list = [padic_gauss_sum(r,p,1,prec=pprec,factored=True) for r in range(p-1)]
    gamma = H.gamma_array()
    gamma_prods = []
    K = Qp(p,prec=pprec)
    for j in range(n):
        s = 1 - beta[j]
        sg = sum(gauss_list[frac(v*s)*(p-1)][0]*gamma[v] for v in gamma)//(p-1)
        gamma_prods.append((-1)^sg * K.teichmuller(H.M_value())^(Integer(s*(p-1))) *
           prod(gauss_list[frac(v*s)*(p-1)][1]^gamma[v] for v in gamma))

    m = min(H.zigzag(b) for b in beta)
    perm = [beta.index(frac(p*b)) for b in beta]
# In fact perm is the identity matrix. It is here as a placeholder for dealing with the case
# where p-1 need not be divisible by all denominators.
    return(matrix(n,n,lambda i,j: 0 if i != perm[j] else
            gamma_prods[j] * p^(H.zigzag(beta[j])-m) * t^(-p+1+floor(p*beta[j]))))

F0 = initial_Frobenius_matrix(H,t,p,pprec)

Consistency check: test the commutation relation between the connection matrix and the Frobenius structure N0F0pF0σ(N0)+D(F0)=0 N_0F_0 - pF_0 \sigma(N_0) + D(F_0) = 0 on the formal solution basis. Note that σ(N0)=N0\sigma(N_0) = N_0 because N0N_0 is a scalar matrix.

if optional:
    N0 = diagonal_matrix([b-1 for b in beta])
    test = N0*F0 - p*F0*N0 + F0.apply_map(lambda x: D(x,t))
    print(test == 0)
True

We now compute the Frobenius structure FF on the original basis from the identity F0=U1Fσ(U)F_0 = U^{-1} F \sigma(U), where σ\sigma denotes the substitution ttpt \mapsto t^p.

def sigma(y,p,tprec):
    R = y.parent()
    v = y.valuation()
    l = y.list()
    l = [(l[i//p] if i%p==0 else 0) for i in range(p*len(l))]
    return(R(l, p*v).add_bigoh(tprec))

Uinvsig = Uinv.apply_map(lambda x: sigma(x,p,tprec))
F = U * F0 * Uinvsig
F = F.apply_map(lambda x: x.add_bigoh(tprec-p))

Consistency check: test the commutation relation between the connection matrix and the Frobenius structure: NFpFσ(N)+D(F)=0 NF - pF \sigma(N) + D(F) = 0 on the original basis.

if optional:
    test = N*F - p*F*N.apply_map(lambda x: sigma(x,p,tprec)) + F.apply_map(lambda x: D(x,t))
    test = test.apply_map(lambda x: x.truncate(tprec-p))
    print(test==0)
True

Convert the coefficients of the entries of FF from Qp\QQ_p to Z/peZ\ZZ/p^e \ZZ. The fact that this does not raise an error is itself a consistency check, as it confirms some "magical" cancellation of pp-adic denominators. (Beware that this might fail for legitimate reasons; for p=13p=13, I get a matrix which is not pp-adically integral but still has the right characteristic polynomial. This means that the matrix F0F_0 was off by a conjugation which I am currently unable to predict.)

pe_ring = IntegerModRing(p^e)
F_mod_pe = F.apply_map(lambda x: x.change_ring(pe_ring).truncate(tprec-p))

Evaluate FF at the Teichmüller lift of some value of tt, then take the characteristic polynomial to get the Euler factor at pp for the hypergeometric motive in this pencil at parameter tt.

t0 = -5
arg = pe_ring(K.teichmuller(t0))
[i(arg).lift_centered() for i in F_mod_pe.charpoly().list()]
[50653, 15059, 3589, 520, 97, 11, 1]

Compare with the hypergeometric trace formula of Cohen-Rodriguez Villegas, as implemented in Sage...

H = HGData(alpha_beta=(alpha,beta))
print(H.euler_factor(t0, p))
50653*T^6 + 15059*T^5 + 3589*T^4 + 520*T^3 + 97*T^2 + 11*T + 1

... and Magma. Note that Magma uses a sign convention whereby what we call tt corresponds to an input value of 1/t1/t. (The Magma examples are static because I didn't attempt to set up Magma on CoCalc; they were run elsewhere.)

> H := HypergeometricData([1/6, 5/6, 1/8, 3/8, 5/8, 7/8],
> [1/3, 2/3, 1/12, 5/12, 7/12, 11/12]);
> EulerFactor(H, -15, 37);
50653*$.1^6 + 15059*$.1^5 + 3589*$.1^4 + 520*$.1^3 + 97*$.1^2 + 11*$.1 + 1

Beware that we got a bit lucky here: the reduction of FF modulo pep^e is not actually supposed to be in (Z/peZ)[t±](\ZZ/p^e \ZZ)[t^{\pm}], but rather in (Z/peZ)[t±,(tp1++1)1](\ZZ/p^e\ZZ)[t^{\pm}, (t^{p-1} + \cdots + 1)^{-1}]. We can see this by increasing ee. (To extract a higher working precision, we would have to multiply in this factor to get a polynomial, evaluate, then take the extra factor back out.)

F_mod_pe2 = F.apply_map(lambda x: x.change_ring(IntegerModRing(p^(e+2))).truncate(tprec-p))
F_mod_pe2[0][0].degree()
221
te = F_mod_pe2[0][0].parent().gen()
u = sum(te^i for i in range(p))
print (F_mod_pe2[0][0]*u^2).truncate(tprec-p).degree(), tprec-p
104 222

Theoretical and computational evidence both suggest that the factor of (tp1++1)1(t^{p-1} + \cdots + 1)^{-1} required to clear denominators modulo pep^e is ee plus a constant. (This ee should be thought of as (d2)e(d-2)e where dd is the number of singularities of the connection, which is 3 here. The singularities at 00 and \infty behave differently from 11 because we are using the Frobenius lift ttpt \mapsto t^p which branches over both 00 and \infty.)

If this is confirmed (and we make precise the required working precision in the pp-adic and tt-adic directions), this indeed gives an algorithm for computing full Euler factors at pp with a linear dependence on pp. (I suspect there is also an average-polynomial-time variant in the style of Harvey and Harvey-Sutherland.)

Let's repack the preceding code as a function, so we can try more examples. (Beware: if you lift this code, you must also lift the user-defined functions D, formal_solution_matrix, initial_Frobenius_matrix, sigma from earlier.)

def euler_factor(H,t0,p,pprec=20,tprec=None,e=7):
    if tprec==None: tprec=7*p

    K = Qp(p, pprec)
    R.<t> = LaurentSeriesRing(K, tprec)

    U = formal_solution_matrix(H,R,tprec,K)
    Uinv = U.inverse()
    F0 = initial_Frobenius_matrix(H,t,p,pprec)
    Uinvsig = Uinv.apply_map(lambda x: sigma(x,p,tprec))
    F = U * F0 * Uinvsig
    F = F.apply_map(lambda x: x.add_bigoh(tprec-p))
    pe_ring = IntegerModRing(p^e)
    F_mod_pe = F.apply_map(lambda x: x.change_ring(pe_ring).truncate(tprec-p))
    arg = pe_ring(K.teichmuller(t0))
    l = [i(arg).lift_centered() for i in F_mod_pe.charpoly().list()]
    l.reverse()
    S.<T> = PolynomialRing(ZZ)
    return(S(l))

# Did I copy everything correctly? Does this answer look familiar?
euler_factor(H,t0,p)
50653*T^6 + 15059*T^5 + 3589*T^4 + 520*T^3 + 97*T^2 + 11*T + 1

Let's check some more examples against Sage. (I omit the cross-checks against Magma.)

alpha = [1/3, 2/3, 1/12, 5/12, 7/12, 11/12]
beta = [1/6, 5/6, 1/10, 3/10, 7/10, 9/10]
p = 31
H = HGData(alpha_beta=(alpha,beta))
print euler_factor(H,-5,p) # Only valid mod p^7
print H.euler_factor(-5,p)
6212525767*T^5 + 1838730311*T^4 - 11633866*T^3 + 61721*T^2 - 303*T + 1 26439622160671*T^6 - 268913615343*T^5 + 1838730311*T^4 - 11633866*T^3 + 61721*T^2 - 303*T + 1
alpha = [1/4, 3/4, 1/8, 3/8, 5/8, 7/8]
beta = [1/3, 2/3, 1/5, 2/5, 3/5, 4/5]
p = 31
H = HGData(alpha_beta=(alpha,beta))
print euler_factor(H,-5,p)
print H.euler_factor(-5,p)
29791*T^6 + 5766*T^5 + 155*T^4 - 82*T^3 + 5*T^2 + 6*T + 1 29791*T^6 + 5766*T^5 + 155*T^4 - 82*T^3 + 5*T^2 + 6*T + 1

Let's try a bigger example with the same parameters.

from time import time
timer = time()
pol = euler_factor(H,-5,331)
print pol # Again, only accurate modulo p^7.
print(time()-timer)

What about Magma?

> H := HypergeometricData([1/4,3/4,1/8,3/8,5/8,7/8],
> [1/3,2/3,1/5,2/5,3/5,4/5]);
> time EulerFactor(H, -1/5, 331);
Runtime error: p^f is too large

What about Sage? (Don't actually run this, I haven't seen it terminate...)

#from time import time
#timer = time()
#print H.euler_factor(-5,331)
#print(time()-timer)