"""
AUTHOR: William Stein, 2006-08-20

Create graph of

  (1) the frequently histogram (distribution) of rounded values
        acos(a_p/2*sqrt(p))
      for all good primes p up to the second argument (e.g., 50000).

"""

[5]

#%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]

[0]

#%auto
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)

[1]

#%auto

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

[2]

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)

[52]

#%auto
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()

[7]

g = graph_ellcurve(EllipticCurve('11a'),200,150000)
show(g)

[3]

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

def graph_ellcurve(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)

[58]

show(graph_ellcurve(EllipticCurve('11a'),50,num=100000))

[59]

show(g,dpi=200)

[4]

f.integrate('x',0,pi)

[8]	%pi/2 %pi/2

#%auto

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)

[9]

g = delta_dist(30,100000)
show(g)

[11]	got delta normalized arccos distributed got delta normalized arccos distributed

show(g,xmin=0,ymin=0,xmax=pi,ymax=0.75,dpi=200)

[22]

save g

[23]

%time 
g20000 = delta_dist(40,20000)
print ""
show(g20000)

[24]	CPU time: 242.95 s, Wall time: 244.06 s CPU time: 242.95 s, Wall time: 244.06 s

%time 
g30000 = delta_dist(40,30000)
show(g30000)

[35]	CPU time: 559.59 s, Wall time: 560.95 s CPU time: 559.59 s, Wall time: 560.95 s

show(g30000,xmin=0,ymin=0,xmax=pi,ymax=0.75,dpi=200)

[40]

%time 
g40000 = delta_dist(40,40000)
print ""
show(g40000)

[25]

Traceback (most recent call last):
    g40000 = delta_dist(40,40000)
  File
"/home/was/papers/riemann_hypothesis/sage_notebook/worksheets/sato_tate/cells/code-13\
3.py", line 4, in delta_dist
    D = pari('q')*pari('prod(n=1,%s,1-q^n+O(q^%s))'%(num,num))**Integer(24)
  File "gen.pyx", line 5649, in gen._pari_trap
gen.PariError: valuation (valp) overflow (28)

Traceback (most recent call last):
    g40000 = delta_dist(40,40000)
  File "/home/was/papers/riemann_hypothesis/sage_notebook/worksheets/sato_tate/cells/code-133.py", line 4, in delta_dist
    D = pari('q')*pari('prod(n=1,%s,1-q^n+O(q^%s))'%(num,num))**Integer(24)
  File "gen.pyx", line 5649, in gen._pari_trap
gen.PariError: valuation (valp) overflow (28)

time D = delta_qexp(500)

[12]	CPU time: 3.76 s, Wall time: 3.78 s CPU time: 3.76 s, Wall time: 3.78 s

time n=10000; D = pari('q')*pari('prod(n=1,%s,1-q^n+O(q^%s))'%(n,n))^24

[17]

D[2]

[19]

-24

-24

delta_qexp(30)

[18]

q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 - 6048*q^6 - 16744*q^7 + 84480*q^8 -
113643*q^9 - 115920*q^10 + 534612*q^11 - 370944*q^12 - 577738*q^13 + 401856*q^14 +
1217160*q^15 + 987136*q^16 - 6905934*q^17 + 2727432*q^18 + 10661420*q^19 -
7109760*q^20 - 4219488*q^21 - 12830688*q^22 + 18643272*q^23 + 21288960*q^24 -
25499225*q^25 + 13865712*q^26 - 73279080*q^27 + 24647168*q^28 + 128406630*q^29 +
O(q^30)

q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 - 6048*q^6 - 16744*q^7 + 84480*q^8 - 113643*q^9 - 115920*q^10 + 534612*q^11 - 370944*q^12 - 577738*q^13 + 401856*q^14 + 1217160*q^15 + 987136*q^16 - 6905934*q^17 + 2727432*q^18 + 10661420*q^19 - 7109760*q^20 - 4219488*q^21 - 12830688*q^22 + 18643272*q^23 + 21288960*q^24 - 25499225*q^25 + 13865712*q^26 - 73279080*q^27 + 24647168*q^28 + 128406630*q^29 + O(q^30)

E = EllipticCurve('11')

[20]

[26]	Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10x - 20 over Rational Field Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10x - 20 over Rational Field

E.Lseries_zeros_in_interval(0,20, 0.1)

[27]

[(6.3626138940, -0.24609225610), (8.6035396196, 0.13661787362), (10.035509098,
0.0068255913393), (11.451258611, -0.23773823824), (13.568639059, -0.21768247004),
(15.914072603, 0.18324124877), (17.033610322, 0.34609712393), (17.941433571,
0.053793259038), (19.185724974, -0.17452048988)]

[(6.3626138940, -0.24609225610), (8.6035396196, 0.13661787362), (10.035509098, 0.0068255913393), (11.451258611, -0.23773823824), (13.568639059, -0.21768247004), (15.914072603, 0.18324124877), (17.033610322, 0.34609712393), (17.941433571, 0.053793259038), (19.185724974, -0.17452048988)]

L = E.Lseries_dokchitser()

[28]

[29]

Dokchitser L-function associated to Elliptic Curve defined by y^2 + y = x^3 - x^2 -
10*x - 20 over Rational Field

Dokchitser L-function associated to Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field

L(1+I)

[30]	0.25232984431224526 + 0.34591234236210699I 0.25232984431224526 + 0.34591234236210699I

L(1)

[31]	0.25384186085591071 0.25384186085591071

L(2+7*I)

[32]

verbose -1 (281: dokchitser.py, __call__) Warning: Loss of 4 decimal digits due to
cancellation
0.90048519563780138 - 0.25383923230303973*I

verbose -1 (281: dokchitser.py, __call__) Warning: Loss of 4 decimal digits due to cancellation
0.90048519563780138 - 0.25383923230303973*I

L(1+6.362613*I)

[33]

verbose -1 (281: dokchitser.py, __call__) Warning: Loss of 4 decimal digits due to
cancellation
-0.0000011925971238401708 - 0.0000019572423741573705*I

verbose -1 (281: dokchitser.py, __call__) Warning: Loss of 4 decimal digits due to cancellation
-0.0000011925971238401708 - 0.0000019572423741573705*I

n=200000
w = [de[p]/(2*float(p)^(5.5)) for p in prime_range(n)]
v = [acos(x) for x in w]
b = 100
d, total_number_of_points = dist(v,b)
g = graph(d, b, total_number_of_points)
g.save('/home/was/people/mazur/sato-tate-delta/sato_tate_delta-%s.png'%n,0,pi,0,0.75)
g.save('/home/was/people/mazur/sato-tate-delta/sato_tate_delta-%s.eps'%n,0,pi,0,0.75)

[38]

g.show()

[45]

b = 100
d, total_number_of_points = dist(v,b)
g = graph(d, b, total_number_of_points)
show(g,0,pi,0,0.75)

[39]

g.save('sato_tate_delta.eps',0,pi,0,0.75)

[41]	sato_tate_delta.eps

show(g,0,pi,0,0.75)

[42]

binomial?

[43]

n=10^6
D=load('/home/was/db/misc/delta-million.sobj')
print "loaded"
w = [D[p]/(2*float(p)^(5.5)) for p in prime_range(n)]
print "normalized"
v = [acos(x) for x in w]
print "arccos"
b = 100
d, total_number_of_points = dist(v,b)
print "distributed"
g = graph(d, b, total_number_of_points)
print "saving..."
g.save('/home/was/people/mazur/sato-tate-delta/sato_tate_delta-%s.png'%n,0,pi,0,0.75)
g.save('/home/was/people/mazur/sato-tate-delta/sato_tate_delta-%s.eps'%n,0,pi,0,0.75)

[44]	loaded normalized arccos distributed saving... loaded normalized arccos distributed saving...

g.show(xmin=0,xmax=pi,ymin=0, ymax=.7,dpi=350)

[46]

%time
n=10^6
Dmillion=load('http://sage.math.washington.edu/home/was/db/misc/delta-million.sobj')

[48]

Attempting to load remote file:
http://sage.math.washington.edu/home/was/db/misc/delta-million.sobj
Loading: [..................................................]
CPU time: 21.49 s,  Wall time: 21.81 s

Attempting to load remote file: http://sage.math.washington.edu/home/was/db/misc/delta-million.sobj
Loading: [..................................................]
CPU time: 21.49 s,  Wall time: 21.81 s

[53]

Sata-Tate for Delta Without Acos

Next we compute the distribution without first computing arc cos, i.e., the actual error distribution, instead of the distribution of angles.

Sata-Tate for Delta Without Acos

Next we compute the distribution without first computing arc cos, i.e., the actual error distribution, instead of the distribution of angles.

n=10^5
D = delta_qexp(n)

[49]

D=Dmillion
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..."
g.show()

[47]	normalized distributed saving... normalized distributed saving...

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

f = maxima('-sin(acos(x))^2/sqrt(1-x^2)')
show(f)
print f.integrate('x',-1,1)
PI = float(pi)

# It turns out to be 4/3
def f(x):
    if abs(x) == 1 or x < -1:
        return 0
    return (2/PI) * sqrt(1-x^2)
    #return (2/PI)*math.sin(math.acos(x))^2/(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)

[50]	\frac{x^2-1}{\sqrt{1-x^2}} -%pi/2 \frac{x^2-1}{\sqrt{1-x^2}} -%pi/2

(g+plot(sin2acos())).show(ymin=0,ymax=1,figsize=[10,5],fontsize=12)

[54]

(g+plot(sin2acos())).save('nature-sato-tate-circle.png',ymin=0,ymax=1,figsize=[10,5])

[55]

maxima('acos(x)').derivative('x')

[56]	-1/sqrt(1 - x^2) -1/sqrt(1 - x^2)

maxima('sin(acos(t))')

[57]	sqrt(1 - t^2) sqrt(1 - t^2)

f = maxima('-sin(acos(x))^2/sqrt(1-x^2)')

[60]

[62]

You can also download the live interactive SAGE worksheet [sato-tate.sws]. You'll need the free program SAGE to use it. Alternatively, check out this live worksheet.

Sata-Tate for Delta Without Acos

Sata-Tate for Delta Without Acos