1650.sagews

"""
Author: Stefan Keil, 2010-02-21.

This is an extended Version of
http://sage.math.washington.edu/sato-tate/
by William Stein.

Let E/K be an elliptic curve over a number field K, which is not the rational numbers QQ. For each finite place v of K, for which E has good reduction, one can calculate the number
a_v := q_v+1-|E(FF_v)|,
where q_v=Nv is the usual norm of v and |E(FF_v)| is the number of points of the mod v reduced curve over the finite residuo field FF_v := O_K / v, which is the ring of integers of K modulo v.
Hasse says |a_v| <= 2 \sqrt q_v,
so one can normalize a_v/(2 \sqrt q_v) = cos theta_v,
with an unique angle theta_v in [0,pi].
The Sato-Tate-Conjecture says, that if one runs over all finite places v of good reduction, then the theta_v (ordered by increasing norm Nv) are equidistributed with respect to the so called Sato-Tate-Measure 2/pi sin^2 theta d theta.

This Worksheet calculates the distribution of the angles theta_v of a given elliptic curve defined over a number field K up to all finite places v of good reduction of E, with norm Nv=q_v <= N, whereas the intervall [0,pi] will be divided into B bins. Additionally it draws the curve of the Sato-Tate-Measure for comparison.
"""

# returns all theta_v of a given elliptic curve E over a number field K up to Nv=q_v<N
def sato_tate(E, K, N):
    all_theta_v = []
    for v in primes_of_K_up_to_norm_N(K,N):
        if E.has_good_reduction(v):
            all_theta_v.append(theta_of(E, v))
    return all_theta_v

# calculates theta_v of a given elliptic curve E and prime v
def theta_of(E, v):
    E_red=E.reduction(v)
    q_v=v.norm()
    cardinality_of_E_red=E_red.cardinality()
    a_v=q_v+1-cardinality_of_E_red
    theta_v = acos(a_v/(2*sqrt(q_v)))
    return theta_v

# returns all primes v of a given number field K with Nv=q_v<N
def primes_of_K_up_to_norm_N(K, N):
    primes =[]
    for p in prime_range(N+1):
        principle_ideal_p = K.ideal(p)
        primfactors_of_principle_ideal_p = principle_ideal_p.prime_factors()
        # take only primes with norm < N+1
        for v in primfactors_of_principle_ideal_p:
            if v.norm()<N+1: 
                primes.append(v)
    return primes

# this function divides the intervall [0,pi] into B bins and then counts how many of the theta_v
# are in which bin, and then returns this distribution and the number of all theta_v
def distribution_of_all_theta_v_in_B_bins(all_theta_v, B):
    normalize = float(B/pi)
    vals = {}
    distribution = dict([(i,0) for i in range(B)])
    for theta_v in all_theta_v:
        # n is the number of the bin
        n = int(normalize*float(theta_v))
        # in case of theta_v=pi
        if n > B-1:
            n=B-1
        distribution[n] += 1
    return distribution, len(all_theta_v)

# this will make the blue recangles according to the given distribution, the number of bins B,
# and the total number of thetas
def graph(distribution, B, number_of_thetas):
    s = Graphics()
    left = float(0); right = float(pi)
    length = right - left
    w = length/B
    k = 0
    for i, n in distribution.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/(number_of_thetas*w)), (w*i+left, n/(number_of_thetas*w))], \
            rgbcolor=(0,0,0.5))
    return s

# plots curve of Sato-Tate measure
def sato_tate_measure():
    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)

# returns the plotted distribution of all theta_v with Nv=q_v<N of an elliptic curve E over K
# divided in B bins and adds the curve of the Sato-Tate measure for comparison
def graph_ellcurve(E, K, B, N):
    all_theta_v = sato_tate(E, K, N)
    distribution, number_of_thetas = distribution_of_all_theta_v_in_B_bins(all_theta_v,B)
    return graph(distribution, B, number_of_thetas) + sato_tate_measure()

# defines a number field K and an elliptic curve E over K
# CAUTION:
# K cannot be QQ
K.<i> = NumberField(x^2 + 1)
E = EllipticCurve(K,[2*i,3])
E

Elliptic Curve defined by y^2 = x^3 + 2*i*x + 3 over Number Field in i with defining polynomial x^2 + 1

# starts the program and shows the distribution of all theta_v with Nv=q_v<N 
# of an elliptic curve E over a number field K divided in B bins
g = graph_ellcurve(E,K,B=40,N=10000)
show(g)

# starts the program and shows the distribution of all theta_v with Nv=q_v<N 
# of an elliptic curve E over a number field K divided in B bins
g = graph_ellcurve(E,K,B=50,N=10000)
show(g)

K.<t> = NumberField(x^2 + 3)
E = EllipticCurve(K,[-5,-3*t])
E

Elliptic Curve defined by y^2 = x^3 + (-5)*x + (-3*t) over Number Field in t with defining polynomial x^2 + 3

E.global_integral_model()

Elliptic Curve defined by y^2 = x^3 + (-5)*x + (-3*t) over Number Field in t with defining polynomial x^2 + 3

# starts the program and shows the distribution of all theta_v with Nv=q_v<N 
# of an elliptic curve E over a number field K divided in B bins
g = graph_ellcurve(E,K,B=50,N=10000)
show(g)

K.<t> = NumberField(x^2 + 167)
E = EllipticCurve(K,[t,2])
E

Elliptic Curve defined by y^2 = x^3 + t*x + 2 over Number Field in t with defining polynomial x^2 + 167

E.global_integral_model()

Elliptic Curve defined by y^2 = x^3 + t*x + 2 over Number Field in t with defining polynomial x^2 + 167

# starts the program and shows the distribution of all theta_v with Nv=q_v<N 
# of an elliptic curve E over a number field K divided in B bins
g = graph_ellcurve(E,K,B=50,N=10000)
show(g)

K.<t> = NumberField(x^5 + 4*x^3 -17)
E = EllipticCurve(K,[t,1])
E

Elliptic Curve defined by y^2 = x^3 + t*x + 1 over Number Field in t with defining polynomial x^5 + 4*x^3 - 17

E.global_integral_model()

Elliptic Curve defined by y^2 = x^3 + t*x + 1 over Number Field in t with defining polynomial x^5 + 4*x^3 - 17

# starts the program and shows the distribution of all theta_v with Nv=q_v<N 
# of an elliptic curve E over a number field K divided in B bins
g = graph_ellcurve(E,K,B=50,N=10000)
show(g)