E = EllipticCurve('20.a3')
if E.has_cm():
print("E has Complex Multiplication!")
P = Primes()
prime_max = 500000
theta = []
final_p = 0
for p in P:
if p >= prime_max:
break
if not E.is_good(p):
continue
a_p = 1 + p - E.Np(p)
theta_p = arccos(a_p/2/sqrt(p))
theta.append(theta_p)
final_p = p
HH = histogram(theta, bins=41, alpha = 0.7, color = 'skyblue', edgecolor = 'skyblue', density = True, axes_labels=[r'$\theta$', ''])
ymax = HH.get_minmax_data()['ymax']
SS = plot(2/pi*sin(x)^2, (x,0,pi), color = 'purple', legend_label = r'$\dfrac{2}{\pi}\sin^2(x)$')
SS.set_legend_options(shadow=False)
TT = text(r'$E: '+latex(E)+'$\n$p: 2, 3, ..., ' + final_p + '$', (0.1, ymax), horizontal_alignment='left', vertical_alignment = 'top')
YLABEL = text('relative frequency', (-0.3, ymax/3), horizontal_alignment='right', vertical_alignment = 'bottom', rotation=90, fontsize = 'large', color='black')
(HH+SS+TT+YLABEL).save('distribution of theta_p.svg')