Mazur - Finding Meaning in Error Terms (SAGE)

g = prime_pi.plot(0,25, rgbcolor=(1,0,0))
g += text(r"$\pi(N)$", (10,8), fontsize=30, rgbcolor=(0,0,0))
g.save('pi25.png')
g.save('pi25.eps')
g.save('pi25.pdf')

pi25.eps pi25.pdf

g = prime_pi.plot(0,100000, rgbcolor=(1,0,0))
g += text(r"$\pi(N)$", (50000,8000), fontsize=30, rgbcolor=(0,0,0))
g.save('pi100000.png')
g.save('pi100000.eps')
g.save('pi100000.pdf')

pi100000.eps pi100000.pdf

PREC  = 50
mu = moebius   # the Moebius function
PI = float(pi)
C = CDF
I = C.0
rho = [None] + zeta_zeros()
# Precompute coefficients of Riemann's function.
coeffs = [1]

def init(prec):
    global coeffs
    coeffs = [1]
    n_factorial = 1.0
    for n in xrange(1, PREC):
        n_factorial *= n
        zeta_value = float(abs(zeta(n+1)))
        coeffs.append(float(1.0/(n_factorial*n*zeta_value)))

init(PREC)

class Riemann:
    def __init__(self, kmax, xmax, prec=50):
        self.xmax = xmax
        self.kmax = kmax
        self.prec = prec
        self.N = int(math.log(xmax)/math.log(2))
        self.rho_k = [0] + [CDF(0.5, zeta_zeros()[k-1]) for k in range(1,kmax+1)]
        self.rho = [[0]+[rho_k / n for n in range(1, self.N+1)]  for rho_k in self.rho_k]
        self.mu = [float(x) for x in moebius.range(0,self.N+2)]
        self._init_coeffs()

    def _init_coeffs(self):
        self.coeffs = [1]
        n_factorial = 1.0
        for n in xrange(1, self.prec):
            n_factorial *= n
            zeta_value = float(abs(zeta(n+1)))
            self.coeffs.append(float(1.0/(n_factorial*n*zeta_value)))
    
    def _check(self, x, k):
        if x > self.xmax:
             raise ValueError, "x (=%s) must be at most %s"%(x, self.xmax)
        if k > self.kmax:
            raise ValueError, "k (=%s) must be at most %s"%(k, self.kmax)
         
    def R(self, x):
        y = math.log(x)
        z = y
        a = float(1)
        for n in xrange(1,self.prec):
            a += self.coeffs[n]*z
            z *= y
        return a

    def Rk(self, x, k):
        return self.R(x) + self.Sk(x, k)

    def plot_Rk(self, k, xmax=None, thickness=0.3, **kwds):
        """
        Plot pi(x) and R_k up to xmax.
        """
        if not xmax:
            xmax = self.xmax
        else:
            xmax = min(self.xmax, xmax)
        f = plot(lambda x: self.Rk(x, k), 10, xmax, thickness=thickness, **kwds)
        g = prime_pi.plot(0, xmax, hue=0)
        return f+g

    def plot_Rk_all(self, kmax, step=1, xmin=10, xmax=None, **kwds):
        """
        Plot pi(x) and every R_k, for k<=kmax, up to xmax.

        Each successive R_k is darker and darker blue (ending at black).
        """
        if not xmax:
            xmax = self.xmax
        else:
            xmax = min(self.xmax, xmax)
        return [plot(lambda x: self.Rk(x, k), xmin=xmin, xmax=xmax,   \
                     rgbcolor=(0,0,float(kmax-k)/kmax), **kwds) for \
                     k in range(1, kmax+1, step)]

    def animate_Rk(self, kmax, step=1, 
                   xmin=10, xmax=None, **kwds):
        if not xmax:
            xmax = self.xmax
        else:
            xmax = min(self.xmax, xmax)
        v = self.plot_Rk_all(kmax, step=step, xmin=xmin, xmax=xmax, **kwds)
        g = plot(prime_pi, 0, xmax)
        G =[z+g for z in v]    
        return animate(G)

    def Sk(self, x, k):
        """
        Compute approximate correction term, so Rk(x,k) = R(x) + Sk(x,k)
        """
        self._check(x, k)

        log_x = math.log(x)  # base e
    
        # This is from equation 32 on page 978 of Riesel-Gohl.
        term1 = sum([self.mu[n] for n in xrange(1,self.N+1)]) / (2*log_x) + \
                   (1/PI) * math.atan(PI/log_x)
    
        # This is from equation 19 on page 975
        term2 = sum(self.Tk(x, v) for v in xrange(1,k+1))
        return term1 + term2

    def Tk(self, x, k):
        """
        Compute sum from 1 to N of
           mu(n)/n * ( -2*sqrt(x) * cos(im(rho_k/n)*math.log(x) \
                - arg(rho_k/n)) / ( pi_over_2 * math.log(x) )
        """
        self._check(x, k)
        x = float(x)
        log_x = math.log(x)
        val = float(0)
        rho_k = self.rho_k[k]
        rho = self.rho[k]
        for n in xrange(1, self.N+1):
            rho_k_over_n = rho[n]
            mu_n = self.mu[n]
            if mu_n != 0:
                z = Ei( rho_k_over_n * log_x)
                val += (mu_n/float(n)) * (2*z).real()
        return -val

%time
Rm = Riemann(30,xmax=100, prec=50)
g = Rm.plot_Rk(29)
g.save('riemann.png')

CPU time: 9.70 s,  Wall time: 9.70 s

CPU time: 9.70 s,  Wall time: 9.70 s

g.save('riemann.pdf')
g.save('riemann.eps')

riemann.eps riemann.pdf

stop = 100
T = theta_qexp(stop)^24
for p in primes(stop):
    print p, '\t', T[p]

2 	1104
3 	16192
5 	1362336
7 	44981376
11 	6631997376
13 	41469483552
17 	793229226336
19 	2697825744960
23 	22063059606912
29 	282507110257440
31 	588326886375936
37 	4119646755044256
41 	12742799887509216
43 	21517654506205632
47 	57242599902057216
53 	214623041906680992
59 	698254765677746880
61 	1007558483942335776
67 	2827903926520931136
71 	5351602023957373056
73 	7264293802635839712
79 	17319684851070915840
83 	29819539398107307072
89 	64258709626203556320
97 	165626956557080594016

2 	1104
3 	16192
5 	1362336
7 	44981376
11 	6631997376
13 	41469483552
17 	793229226336
19 	2697825744960
23 	22063059606912
29 	282507110257440
31 	588326886375936
37 	4119646755044256
41 	12742799887509216
43 	21517654506205632
47 	57242599902057216
53 	214623041906680992
59 	698254765677746880
61 	1007558483942335776
67 	2827903926520931136
71 	5351602023957373056
73 	7264293802635839712
79 	17319684851070915840
83 	29819539398107307072
89 	64258709626203556320
97 	165626956557080594016

#auto
    
def sato_tate(E, N):
    return [acos(E.ap(p)/(2*sqrt(p))) for p in prime_range(N+1) if N%p != 0]

def dist(v, b, left=float(0), right=float(pi)):
    """
    We divide the interval between left (default: 0) and 
    right (default: pi) up into b bins.
   
    For each number in v (which must left and right), 
    we find which bin it lies in and add this to a counter.
    This function then returns the bins and the number of
    elements of v that lie in each one. 

    ALGORITHM: To find the index of the bin that a given 
    number x lies in, we multiply x by b/length and take the 
    floor. 
    """
    length = right - left
    normalize = float(b/length)
    vals = {}
    d = dict([(i,0) for i in range(b)])
    for x in v:
        n = int(normalize*(float(x)-left))
        d[n] += 1
    return d, len(v)

def graph(d, b, num=5000, left=float(0), right=float(pi)):
    s = Graphics()
    left = float(left); right = float(right)
    length = right - left
    w = length/b
    k = 0
    for i, n in d.iteritems():
        k += n
        # ith bin has n objects in it. 
        s += polygon([(w*i+left,0), (w*(i+1)+left,0), \
                     (w*(i+1)+left, n/(num*w)), (w*i+left, n/(num*w))],\ 
                     rgbcolor=(0,0,0.5))
    return s

def sin2():
    PI = float(pi)
    return plot(lambda x: (2/PI)*math.sin(x)^2, 0,pi, plot_points=200, 
              rgbcolor=(0.3,0.1,0.1), thickness=2)

def graph_ellcurve(E, b=10, num=5000):
    v = sato_tate(E, num)
    d, total_number_of_points = dist(v,b)
    return graph(d, b, total_number_of_points) + sin2()

def sato_tate_noacos(E, N):
    return [E.ap(p)/(2*math.sqrt(p)) for p in prime_range(N+1) if N%p != 0]

def graph_ellcurve_noacos(E, b=10, num=5000):
    v = sato_tate_noacos(E, num)
    d, total_number_of_points = dist(v,b,-1,1)
    return graph(d, b, total_number_of_points,-1,1)

def delta_dist(b=10, num=5000): 
    D = delta_qexp(num)
    print "got delta"
    w = [float(D[p])/(2*float(p)^(5.5)) for p in prime_range(num+1)]
    print "normalized"
    v = [acos(x) for x in w]
    print "arccos"
    d, total_number_of_points = dist(v,b)
    print "distributed"
    return graph(d, b, total_number_of_points)

# First we compute the normalization so 
# that the resulting function has integral 1.

x = var('x')
f = -sin(acos(x))^2 / sqrt(1-x^2)
show(f)
print f.integrate(x,-1,1)

# It turns out to be 4/3
def f(x):
    if abs(x) == 1 or x < -1:
        return 0
    return (2/math.pi) * sqrt(1-x^2)

def sin2acos():
    PI = float(pi)
    return plot(f, -1.01,1, plot_points=200, \
              rgbcolor=(1,0,0), thickness=4,alpha=0.7)

x2�1p1�x2
				        pi
        			     - ---
					2

x2�1p1�x2
				        pi
        			     - ---
					2

%time

n = 10^6
D = delta_qexp(n)
v = [D[p]/(2*float(p)^(5.5)) for p in prime_range(n)]
print "normalized"
b = 200
d, total_number_of_points = dist(v,b,-1,1)
print "distributed"
g = graph(d, b, total_number_of_points, -1, 1)
print "saving..."
h = plot(sin2acos())

normalized
distributed
saving...
CPU time: 34.22 s,  Wall time: 34.22 s

normalized
distributed
saving...
CPU time: 34.22 s,  Wall time: 34.22 s

(g+h).show(ymin=0, ymax=1,figsize=[8,4])

(g+h).save('squares.eps', ymin=0, ymax=1, figsize=[8,4])

squares.eps

(g+h).save('squares.pdf', ymin=0, ymax=1, figsize=[8,4])

squares.pdf

%time

pmax = 10^6  # use primes up to this bound
g = graph_ellcurve(EllipticCurve('11a'),200, pmax)
show(g, xmin=0, ymin=0, ymax=0.8)

CPU time: 44.07 s,  Wall time: 44.07 s

CPU time: 44.07 s,  Wall time: 44.07 s

g.save('ellst.eps',xmin=0, ymin=0, ymax=0.8)
g.save('ellst.pdf',xmin=0, ymin=0, ymax=0.8)

ellst.eps ellst.pdf

#auto
def line1(xmin,xmax): return line([(xmin,1),(xmax,1)],
rgbcolor=(1,0,0))
def Xab(a,b): bb = (math.asin(b)/2r + b*math.sqrt(1r-b^2r)/2r) aa =
(math.asin(a)/2r + a*math.sqrt(1r-a^2r)/2r) def X(T): return
(math.asin(T)/2r + T*math.sqrt(1r-T^2r)/2r - aa)/(bb - aa) return X
import bisect
class SatoTate: def __init__(self, E): self._E = E def __repr__(self):
return "Sato-Tate data for %s"%self._E 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*math.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 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 -math.log(val)/math.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
-math.log(val-abs(err))/math.log(C),
-math.log(val+abs(err))/math.log(C) 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): a,b = (float(a), float(b)) vmin = [];
vmax = [] for C in range(100, Cmax, int(Cmax/plot_points)): zmin,zmax =
self.theta_interval(C, a, b, max_points=max_points) vmin.append((C,
zmin)) vmax.append((C, zmax)) if verbose: print C, zmin, zmax return
vmin, vmax def plot_theta_interval(self, Cmax, clr=(0,0,0), *args,
**kwds): vmin, vmax = self.compute_theta_interval(Cmax, *args, **kwds)
v = self.compute_theta(Cmax, *args, **kwds) grey = (0.7,0.7,0.7) return
line(vmin,rgbcolor=grey)+line(vmax,rgbcolor=grey) +
point(v,rgbcolor=clr) + line(v,rgbcolor=clr) + line1(0, Cmax) def
histogram(self, Cmax, num_bins): v = self.normalized_aplist(Cmax) d,
total_number_of_points = dist(v, num_bins) return
frequency_histogram(d, num_bins, total_number_of_points) + semicircle
def x_times_Delta(self, x): return x*self.Delta(x, -1,1,
max_points=500)[0]
S0 = SatoTate(EllipticCurve('11a'))
S1 = SatoTate(EllipticCurve('37a'))
S2 = SatoTate(EllipticCurve('389a'))
S3 = SatoTate(EllipticCurve('5077a'))
S4 = SatoTate(EllipticCurve([1,-1,0,-79,289]))
S5 = SatoTate(EllipticCurve([0, 0, 1, -79, 342]))
S6 = SatoTate(EllipticCurve([1, 1, 0, -2582, 48720]))
S7 = SatoTate(EllipticCurve([0, 0, 0, -10012, 346900]))
S8 = SatoTate(EllipticCurve([0, 0, 1, -23737, 960366]))
S28 =
SatoTate(EllipticCurve([1,-1,1,-20067762415575526585033208209338542750930230312178956502,34481611795030556467032985690390720374855944359319180361266008296291939448732243429
]))
S = [S0,S1,S2,S3,S4,S5,S6,S7,S8]
CM = SatoTate(EllipticCurve('32a'))
Equest = SatoTate(EllipticCurve('72a1')) # non-CM but doesn't satisfy
hypo of Taylor's theorem
# Our curves S satisfy Taylor's theorem
[min([e for p,e in X._E.conductor().factor()]) for X in S]

[1, 1, 1, 1, 1, 1, 1, 1, 1]

[1, 1, 1, 1, 1, 1, 1, 1, 1]

g = S0.qqplot(10^2)
g.show(xmin=0,ymin=0)
g.save('qq2.eps',xmin=0,ymin=0)
g.save('qq2.pdf',xmin=0,ymin=0)

qq2.eps qq2.pdf

g = S0.qqplot(10^3)
g.show(xmin=0,ymin=0)
g.save('qq3.eps',xmin=0,ymin=0)
g.save('qq3.pdf',xmin=0,ymin=0)

qq3.eps qq3.pdf

g = S0.qqplot(10^4)
g.show(xmin=0,ymin=0)
g.save('qq4.eps',xmin=0,ymin=0)
g.save('qq4.pdf',xmin=0,ymin=0)

qq4.eps qq4.pdf

# Find expression four the 24th power of theta in terms of Eisenstein series
# and cusp forms.

M = ModularForms(1,12)
theta24 = theta_qexp(30)^24
delta = delta_qexp(30); eis = eisenstein_series_qexp(12,30)
q = delta.parent().gen()
M.find_in_space(theta24, [delta, delta(q^2), delta(q^4), eis, eis(q^2), eis(q^4)])

[33152/691, 1525760/691, 135790592/691, 16/691, -32/691, 65536/691]

[33152/691, 1525760/691, 135790592/691, 16/691, -32/691, 65536/691]

def diff_count(v, **kwds):
    """
    INPUT:    
       v -- a list of [a2, a3, a5, a7, ...] of error terms, one for each prime.
    """
    over = 0
    under = 0
    diff = []
    P = primes_first_n(len(v)+1)
    i = 0
    for ap in v:
        p = P[i]; i += 1 
        if ap < 0:
            under += 1
        elif ap > 0:
            over += 1
        diff.append((p,over - under))
    return line(diff, **kwds)

%time
n = 10^6
t = delta_qexp(n)
v = [t[p] for p in primes(n)]
g = diff_count(v, rgbcolor=(0,0,0.5))
g.show()

CPU time: 44.75 s,  Wall time: 44.75 s

CPU time: 44.75 s,  Wall time: 44.75 s

g.save('deltarace.eps')
g.save('deltarace.pdf')

deltarace.eps deltarace.pdf

E = EllipticCurve('11a')

def units(p, ap):
    ap = CDF(ap)
    s = sqrt(ap^2 - 4*p)
    sqrtp = math.sqrt(p)
    return [ (ap + s)/(2*sqrtp), (ap - s)/(2*sqrtp) ]

def E_many_units(pmax):
    return sum([units(p,E.ap(p)) for p in prime_range(pmax+1)],[])

# Same graph in 3d

t = Tachyon(xres=1000, yres=800, camera_center=(1,1,2), look_at=(0,0,0), raydepth=4)
t.light((10,3,2), 1, (1,1,1))
t.light((10,-3,2), 1, (1,1,1))
t.texture('black', color=(0,0,0))
t.texture('red', color=(1,0,0))
t.texture('grey', color=(.95,.95,.95))
t.plane((0,0,0),(0,0,1),'grey')
t.cylinder((0,0,0),(1,0,0),.003,'black')
t.cylinder((0,0,0),(0,1,0),.003,'black')

for x, y in E_many_units(10000):
    t.sphere((x, -y, .01), .005, 'red')
          
t.save()

# distributed version, i.e., where the graph indicates something about the
# number of points at a given position.

E = EllipticCurve('11a')

def units(p, ap):
    ap = CDF(ap)
    s = sqrt(ap^2 - 4*p)
    sqrtp = math.sqrt(p)
    return [ (ap + s)/(2*sqrtp), (ap - s)/(2*sqrtp) ]

def E_units(p):
    return units(p, E.ap(p))

@CachedFunction
def E_many_units(pmax):
    return sum([units(p,E.ap(p)) for p in prime_range(pmax+1)], [])

@CachedFunction
def E_dist(pmax, bin_size):
    v = E_many_units(pmax)
    w = [(  (z[1]/z[0] if z[0] != 0 else 999999) , z) for z in v]
    w.sort()
    # Create each bin by finding a value z in v, then counting the number
    # of values in in v that are within r of z. 
    B = []; i = 0
    while i < len(w):
        bin = [w[i][1]+bin_size/2.0]
        i += 1
        while i < len(w) and abs(w[i][1] - bin[0]) < bin_size:
            bin.append(w[i][1])
            i += 1
        B.append(bin)
    return B

%time

# Distribution graph in 3d
pmax = 50000

t = Tachyon(xres=1000, yres=800, camera_center=(1.5,1.5,3), look_at=(0,0,0), raydepth=4)
#t.light((10,3,2), 1, (1,1,1))
t.light((10,-3,2), 1, (1,1,1))
t.texture('black', color=(0,0,0))
t.texture('red', color=(1,0,0))
t.texture('grey', color=(.95,.95,.95))
t.plane((0,0,0),(0,0,1),'grey')
t.cylinder((0,0,0),(1,0,0),.003,'black')
t.cylinder((0,0,0),(0,1,0),.003,'black')

for v in E_dist(pmax,.01):
    x,y = v[0]
    for i in range(len(v)):
       t.sphere((x, y, .01+i/30.0), .02, 'red')
          
t.save()

CPU time: 2.01 s,  Wall time: 2.88 s

CPU time: 2.01 s,  Wall time: 2.88 s

%time
n=10^6
E = EllipticCurve('11a')
view(E)
v = E.aplist(n)
g = diff_count(v, rgbcolor=(0,0,0.5))
g.show()
g.save('11a.eps')
g.save('11a.pdf')

y2+y=x3�x2�10x�20 
CPU time: 32.35 s,  Wall time: 32.35 s

y2+y=x3�x2�10x�20 
CPU time: 32.35 s,  Wall time: 32.35 s

11a.eps 11a.pdf

%time
n=10^6
E = EllipticCurve('37a')
view(E)
v = E.aplist(n)
g = diff_count(v, rgbcolor=(0,0,0.5))
g.show()
g.save('37a.eps')
g.save('37a.pdf')

y2+y=x3�x 
CPU time: 33.37 s,  Wall time: 33.37 s

y2+y=x3�x 
CPU time: 33.37 s,  Wall time: 33.37 s

37a.eps 37a.pdf

%time
n=10^6
E = EllipticCurve('389a')
view(E)
v = E.aplist(n)
g = diff_count(v, rgbcolor=(0,0,0.5))
g.show()
g.save('389a.eps')
g.save('389a.pdf')

y2+y=x3+x2�2x 
CPU time: 29.39 s,  Wall time: 29.39 s

y2+y=x3+x2�2x 
CPU time: 29.39 s,  Wall time: 29.39 s

389a.eps 389a.pdf

%time
n=10^6
E = EllipticCurve('5077a')
view(E)
v = E.aplist(n)
g = diff_count(v, rgbcolor=(0,0,0.5))
g.show()
g.save('5077a.eps')
g.save('5077a.pdf')

y2+y=x3�7x+6 
CPU time: 31.23 s,  Wall time: 31.30 s

y2+y=x3�7x+6 
CPU time: 31.23 s,  Wall time: 31.30 s

5077a.eps 5077a.pdf

Mazur - Finding Meaning in Error Terms

0 seconds ago by was

Figure 1.1: The step function �(N) counts the number of primes up to N

Figure 1.2: The smooth function slithering up the staircase of primes up to 100 is Riemann’s approximation that uses the “first” 29 zeroes of the Riemann zeta function.

Table of number of ways to write p as a sum of 24 squares.

Figure 3.1: Sato-Tate for the curve 11a

Figure 1.4: q-q-plot on (0;+1) for the cutoff C=100

Figure 1.5: q-q-plot on (0;+1) for the cutoff C=1000

Figure 1.6: The difference between undercounts and overcounts

Figure 2.1: E is the curve y2+y=x3�x2 . Sate-Tate on the unit circle in 3d

Figure 2.2: E is the curve y2+y=x3�x2 . The race between NE(p)<p+1 and NE(p)>p+1

Figure 2.3: E is the curve 37a. The race between NE(p)<p+1 and NE(p)>p+1

Figure 2.4: E is the curve 389a. The race between NE(p)<p+1 and NE(p)>p+1

Figure 2.5: E is the curve 5077a. The race between NE(p)<p+1 and NE(p)>p+1