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.)
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).
Specify hypergeometric data using tuples α=(α1,…,αn),β=(β1,…,βn),
plus a prime p. In this demonstration, we insist that β contain no repeated elements nor 0 or 1, and p−1 be divisible by the least common denominator of β. (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 sortedp=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.)
Optional: compute the connection matrix
associated to the hypergeometric differential equation
where D=tdtd. (Normalization as per Beukers-Heckman, Monodromy for the hypergeometric function nFn−1.)
The Clausen-Thomae hypergeometric seriesnFn−1(β1,…,βn−1α1,…,αn∣∣∣∣t)=k=0∑∞(β1)k⋯(βn−1)k(α1)k⋯(αn)kk!zk,(α)k=α(α+1)⋯(α+k−1)
has the property that for i=1,…,n, the Puiseux series
(where βi−βi+1 is omitted; it contributes the k!) is a solution of the hypergeometric equation. Using this, we construct a formal solution matrixU (and its inverse) for which
namely, the i-th column of U consists of y,(D+1−βi)(y),…,(D+1−βi)n−1(y) where y is the i-th solution from above with the factor of z1−βi removed (because I don't want to deal with fractional exponents).
fromsage.rings.padics.miscimportgauss_sumaspadic_gauss_sumdefD(x,t):returnt*x.derivative()defformal_solution_matrix(H,R,tprec,K=QQ):alpha,beta=H.alpha_beta()n=H.degree()l=forbinbeta:alpha2=[a2-b+1fora2inalpha]beta2=[b2-b+1forb2inbetaifb2!=b]tmp=[K(1)]foriinrange(1,tprec):tmp.append(tmp[-1]*(prod(a2+i-1fora2inalpha2)/prod(b2+i-1forb2inbeta2)/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()foriinrange(n-1):matlist.append([(1-beta[i])*matlist[-1][i]+D(matlist[-1][i],t)foriinrange(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
Compute the (conjectural) Frobenius matrix F0 on the formal solution basis. Note that this matrix must satisfy the commutation relation
Because the βj are distinct and p−1 is divisible by all denominators, F0 is forced to be diagonal, and the (j,j)-entry must be a scalar multiple of t−p+1+⌊pβj⌋. 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)
where ωp is the Teichmüller character, M and v↦γv are invariants of the hypergeometric data, and gp(a)=u∈Fp∗∑ωp(u)−aζpu
is a Gauss sum. The latter is computed using the Gross-Koblitz formula:
for π the root of πp−1=−p for which ζp≡1+π(modπ2),
Here we rely on the fact that (p−1)βj∈Z for all j. If this fails, but the βj are still distinct, then F0 is not scalar but rather a weighted permutation matrix. Things get much more complicated if the βj are not all distinct: for example, for the Dwork quintic family
a calculation of Shapiro shows that
where ∗ is a particular nonzero multiple of ζ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.
definitial_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)forrinrange(p-1)]gamma=H.gamma_array()gamma_prods=K=Qp(p,prec=pprec)forjinrange(n):s=1-beta[j]sg=sum(gauss_list[frac(v*s)*(p-1)]*gamma[v]forvingamma)//(p-1)gamma_prods.append((-1)^sg*K.teichmuller(H.M_value())^(Integer(s*(p-1)))*prod(gauss_list[frac(v*s)*(p-1)]^gamma[v]forvingamma))m=min(H.zigzag(b)forbinbeta)perm=[beta.index(frac(p*b))forbinbeta]# 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,lambdai,j:0ifi!=perm[j]elsegamma_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
on the formal solution basis. Note that σ(N0)=N0 because N0 is a scalar matrix.
Convert the coefficients of the entries of F from Qp to Z/peZ. 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 F0 was off by a conjugation which I am currently unable to predict.)
... 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.)
Beware that we got a bit lucky here: the reduction of F modulo pe is not actually supposed to be in (Z/peZ)[t±], but rather in (Z/peZ)[t±,(tp−1+⋯+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.)
Theoretical and computational evidence both suggest that the factor of (tp−1+⋯+1)−1 required to clear denominators modulo pe 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 ∞ behave differently from 1 because we are using the Frobenius lift t↦tp which branches over both 0 and ∞.)
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.)
defeuler_factor(H,t0,p,pprec=20,tprec=None,e=7):iftprec==None:tprec=7*pK=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(lambdax:sigma(x,p,tprec))F=U*F0*UinvsigF=F.apply_map(lambdax:x.add_bigoh(tprec-p))pe_ring=IntegerModRing(p^e)F_mod_pe=F.apply_map(lambdax:x.change_ring(pe_ring).truncate(tprec-p))arg=pe_ring(K.teichmuller(t0))l=[i(arg).lift_centered()foriinF_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)