Beyond the Sato-Tate Conjecture

William Stein and Chris Swierczewski (joint work with Barry Mazur)

24 April, 2007

1 Some Functions

• X_{a}^{b}(T) = \frac{∫_{a}^{T} \sqrt{1 − x^{2}}dx}{∫_{a}^{b} \sqrt{1 − x^{2}}dx} = area under arc of semicircle

• Y_{a,C}^{b}(T) = \frac{#\{primes p < C : a < \frac{a_{p}}{2\sqrt{p}} < T \}}{#\{primes p < C : a < \frac{a_{p}}{2\sqrt{p}} < b\}}.

• Δ_{a}^{b}(C) = ∫_{a}^{b}(X_{a}^{b}(T) − Y_{a,C}^{b}(T))^{2}dT = the L_{2} norm of the difference of X_{a}^{b} and Y_{a,C}^{b}.

• e_{a}^{b}(C) = − log_{C}(Δ_{a}^{b}(C)) = − log(Δ_{a}^{b}(C))/ log(C) - Measures behavior of Δ_{a}^{b}(C) as a function of C.

• ϵ_{a}^{b}(E) = lim sup_{C→+∞} e_{a}^{b}(C) - The lim sup of e functions.

2 The Code

We compute the data on the convergence of the Δ_{a}^{b}(C) function using the following code.

#auto
from math import asin, log, sqrt
def line1(xmin,xmax):
    return line([(xmin,1),(xmax,1)], rgbcolor=(1,0,0))
def Xab(a,b):
    bb = (asin(b)/2r + b*sqrt(1r-b^2r)/2r)
    aa = (asin(a)/2r + a*sqrt(1r-a^2r)/2r)
    def X(T):
        return (asin(T)/2r + T*sqrt(1r-T^2r)/2r - aa)/(bb - aa)
    return X

The Sato-Tate conjecture implies that the values following function, Y_{a,C}^{b}(T) as C approaches infinity, is equal to X_{a}^{b}(T). X_{a}^{b}(T) has a closed-form definite integral shown in the code above. In SAGE, we define a SatoTate class determined by a chosen elliptic curve E.

import bisect
class SatoTate:
    def __init__(self, E):
        self._E = E
    def __repr__(self):
        return "Sato-Tate data for ..."

def anlist(self, n):
return self._E.anlist(n)
def normalized_aplist(self, n):
anlist = self.anlist(n)
two = float(2)
v = [float(anlist[p])/(two*sqrt(p)) for p in prime_range(n)]
return v
def sorted_aplist(self, n):
v = self.normalized_aplist(n)
v.sort()
return v
def YCab(self, Cmax, a=-1, b=1):
v = self.sorted_aplist(Cmax)
denom = bisect.bisect_right(v, float(b)) - bisect.bisect_left(v, float(a))
try:
   normalize = float(1)/denom
except:
   def Y(T):
       return 1.0r
   return Y
start_pos = bisect.bisect_left(v, float(a))
def Y(T):
    # find position that T would go in if it were inserted
    # in the sorted list v.
    n = bisect.bisect_right(v, float(T)) - start_pos
    return n * normalize
return Y

In the definition of the function Y_{a,C}^{b}(T), we create a sorted a_{p} list from the attached elliptic curve where a_{p} is defined by

a_{p} = p − 1 − #E(𝔽_{p})  (1)

where #E(𝔽_{p}) is the number of points on an elliptic curve E over the finite field 𝔽_{p}. After normalization, a function is returned that computes the number of elements in the list between a and T. The bisect left and bisect right commands are a fast way of finding these upper and lower bound elements.

def xyplot(self, C, a=-1, b=1):
"""
Return the quantile-quantile plot for given a,b, up to C.
"""
Y = self.YCab(C,a=a,b=b)
X = Xab(a=a,b=b)
pX = plot(X, a, b, rgbcolor=(1,0,0))
pY = plot(Y, a, b, rgbcolor=(0,0,1))
return pX + pY
def qqplot(self, C, a=-1, b=1):
"""
Return the quantile-quantile plot for given a,b, up to C.
"""
Y = self.YCab(C,a=a,b=b)

X = Xab(a=a,b=b)
pl = parametric_plot((X, Y), a,b)
ll = line([(0,0), (1.1,1.1)], rgbcolor=(1,0,0))
return pl+ll
def Delta(self, C, a, b, max_points=300):
"""
Delta_{a}^{b} function:
INPUT: C - cutoff
     a,b - evaluate over the interval (a,b)
     max_points - number of points used in numerical integral
"""
key = (C,a,b,max_points)
try:
   return self._delta[key]
except AttributeError:
   self._delta = {}
except KeyError:
   pass
X = Xab(a,b)
Y = self.YCab(C,a,b)
def h(T):
    return (X(T) - Y(T))^2r
val, err = integral_numerical(h, a, b, max_points=max_points, algorithm='qag', rule=1,
eps_abs=1e-10, eps_rel=1e-10)
self._delta[key] = (val, err)
return val, err
def theta(self, C, a=-1, b=1, max_points=300):
val, err = self.Delta(C, a, b, max_points=max_points)
return -log(val)/log(C), val, err
def theta_interval(self, C, a=-1, b=1, max_points=300):
val, err = self.Delta(C, a, b, max_points=max_points)
return -log(val-abs(err))/log(C), -log(val+abs(err))/log(C)

Here we define the Δ_{a}^{b}(C) function as an L_{2}-norm on the difference of the X_{a}^{b}(T) and Y_{a,C}^{b}(T) functions:

Δ_{a}^{b}(C) = ∫_{a}^{b}(X_{a}^{b}(T) − Y_{a,C}^{b}(T))^{2}dT  (2)

We use an adaptive numerical integration scheme (which is part of the GSL library, and which Josh Kantor made easy to use in SAGE) which allows us to (hopefully) put a bound on the error that is likely to crop up. Another function called theta returns the value and error of the delta function along with e_{a}^{b}(C) = − log(Δ_{a}^{b}(C))/ log(C). These e_{a}^{b}(C) values are what we compute and plot in the following code.

def compute_theta(self, Cmax, plot_points=30, a=-1, b=1, max_points=300, verbose=False):
a,b = (float(a), float(b))
def f(C):
   z = self.theta(C, a, b, max_points=max_points)
   if verbose: print C, z
   return z[0]
return [(x,f(x)) for x in range(100, Cmax, int(Cmax/plot_points))]
def compute_theta_interval(self, Cmax, plot_points=30, a=-1, b=1, max_points=300, verbose=False):