Sharedfrobenius.ipynbOpen in CoCalc
Author: Kiran Kedlaya
Views : 41

### $L$-functions via deformations: from hyperelliptic curves to hypergeometric motives

##### Arithmetic of Hyperelliptic Curves, Abdus Salam International Centre for Theoretical Physics (ICTP), September 8, 2017

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.

In [1]:
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).

In [2]:
optional = True


Specify hypergeometric data using tuples $\alpha = (\alpha_1,\dots,\alpha_n), \beta = (\beta_1,\dots,\beta_n)$, plus a prime $p$. In this demonstration, we insist that $\beta$ contain no repeated elements nor 0 or 1, and $p-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.)

In [3]:
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.)

In [4]:
pprec = 25
tprec = 7*p
e = 7


Define a few global variables.

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


Optional: compute the connection matrix $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 = \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 = t \frac{d}{dt}$. (Normalization as per Beukers-Heckman, Monodromy for the hypergeometric function ${}_nF_{n-1}$.)

In [6]:
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 ${}_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,\dots,n$, the Puiseux series $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 $\beta_i - \beta_i+1$ is omitted; it contributes the $k!$) is a solution of the hypergeometric equation. Using this, we construct a formal solution matrix $U$ (and its inverse) for which $U^{-1} N U + U^{-1} D(U) = N_0, \qquad N_0 = \mathrm{Diag}(\beta_1-1, \dots, \beta_n - 1);$ namely, the $i$-th column of $U$ consists of $y,(D+1-\beta_i)(y),\dots,(D+1-\beta_i)^{n-1}(y)$ where $y$ is the $i$-th solution from above with the factor of $z^{1-\beta_i}$ removed (because I don't want to deal with fractional exponents).

In [7]:
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))
#
# 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 $U^{-1} N U + U^{-1} D(U) = N_0.$

In [8]:
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 $F_0$ on the formal solution basis. Note that this matrix must satisfy the commutation relation $N_0F_0 - pF_0 \sigma(N_0) + D(F_0) = 0.$ Because the $\beta_j$ are distinct and $p-1$ is divisible by all denominators, $F_0$ is forced to be diagonal, and the $(j,j)$-entry must be a scalar multiple of $t^{-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) $\omega_p(M^{(p-1)(1-\beta_i)}) \prod_{v} g_p((p-1)(1-\beta_i) v)^{\gamma_v}$ where $\omega_p$ is the Teichmüller character, $M$ and $v \mapsto \gamma_v$ are invariants of the hypergeometric data, and $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 $\pi^{p-1} = -p$ for which $\zeta_p \equiv 1 + \pi \pmod{\pi^2}$, $g_p(a) = \pi^{-a} \Gamma_p(a/(p-1)) \qquad (a=0,\dots,p-2).$ Here we rely on the fact that $(p-1) \beta_j \in \mathbb{Z}$ for all $j$. If this fails, but the $\beta_j$ are still distinct, then $F_0$ is not scalar but rather a weighted permutation matrix. Things get much more complicated if the $\beta_j$ are not all distinct: for example, for the Dwork quintic family $\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 $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 $\zeta_3$. It may be possible to guess the off-diagonal terms in general using ideas of Vologodsky; they should involve $p$-adic multiple zeta values.

In [9]:
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 $N_0F_0 - pF_0 \sigma(N_0) + D(F_0) = 0$ on the formal solution basis. Note that $\sigma(N_0) = N_0$ because $N_0$ is a scalar matrix.

In [10]:
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 $F$ on the original basis from the identity $F_0 = U^{-1} F \sigma(U)$, where $\sigma$ denotes the substitution $t \mapsto t^p$.

In [11]:
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))]

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: $NF - pF \sigma(N) + D(F) = 0$ on the original basis.

In [12]:
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 $F$ from $\QQ_p$ to $\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 $p$-adic denominators. (Beware that this might fail for legitimate reasons; for $p=13$, I get a matrix which is not $p$-adically integral but still has the right characteristic polynomial. This means that the matrix $F_0$ was off by a conjugation which I am currently unable to predict.)

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


Evaluate $F$ at the Teichmüller lift of some value of $t$, then take the characteristic polynomial to get the Euler factor at $p$ for the hypergeometric motive in this pencil at parameter $t$.

In [14]:
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...

In [15]:
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 $t$ corresponds to an input value of $1/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 $F$ modulo $p^e$ is not actually supposed to be in $(\ZZ/p^e \ZZ)[t^{\pm}]$, but rather in $(\ZZ/p^e\ZZ)[t^{\pm}, (t^{p-1} + \cdots + 1)^{-1}]$. We can see this by increasing $e$. (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.)

In [16]:
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
In [17]:
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 $(t^{p-1} + \cdots + 1)^{-1}$ required to clear denominators modulo $p^e$ is $e$ plus a constant. (This $e$ should be thought of as $(d-2)e$ where $d$ is the number of singularities of the connection, which is 3 here. The singularities at $0$ and $\infty$ behave differently from $1$ because we are using the Frobenius lift $t \mapsto t^p$ which branches over both $0$ and $\infty$.)

If this is confirmed (and we make precise the required working precision in the $p$-adic and $t$-adic directions), this indeed gives an algorithm for computing full Euler factors at $p$ with a linear dependence on $p$. (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.)

In [18]:
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.)

In [19]:
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
In [20]:
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.

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


> 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...)

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