#Code for the odd prime characteristic function
def chi_p(n):
x=0
if is_prime(n):
x=1
if is_even(n):
x=0
return x
|
|
#Code that produces the polynomial F_N(z) for fixed integer N
z=var('z')
def F(N,z):
sum=0
for k in range(N):
innersum=0
for n in range(1,N):
innersum=innersum+chi_p(n)*z^(k*n)
sum=sum+innersum^2
return expand(sum)
|
|
#Listing the first few factorizations to notice cyclotomic factors
x = PolynomialRing(RationalField(), 'x').gen()
for i in range(6,10):
factor(F(i,x))
(x^4 - x^2 + 1) * (x^46 + x^44 - x^40 - x^38 + 3*x^36 + 4*x^34 + x^32 - 3*x^30 - 2*x^28 + 3*x^26 + 5*x^24 + 2*x^22 - 2*x^18 - x^16 + 2*x^14 + 5*x^12 + 3*x^10 - x^8 - 3*x^6 + 4*x^2 + 4) (x^6 - x^5 + x^4 - x^3 + x^2 - x + 1) * (x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) * (x^48 - x^46 + x^38 + x^36 - x^34 - x^32 + 3*x^28 - 3*x^26 + 2*x^24 + x^20 - x^18 - 2*x^16 + 3*x^14 - x^10 + x^8 + x^6 - 4*x^2 + 4) (x^8 + 1) * (x^90 - x^82 + 3*x^76 + x^74 - 3*x^68 - x^66 + 2*x^64 + 4*x^62 + 3*x^60 + x^58 - 2*x^56 - 4*x^54 + 2*x^52 - x^50 + 5*x^48 + 4*x^46 - 2*x^44 + 4*x^42 - x^40 - 4*x^38 + 2*x^36 - 2*x^34 + 6*x^32 + 4*x^30 + x^28 + 2*x^26 - 4*x^24 - 2*x^18 + 9*x^16 + 3*x^12 + 3*x^10 - 7*x^8 + x^6 + 9) (x^6 - x^3 + 1) * (x^6 + x^3 + 1) * (x^100 - x^94 + x^86 + 2*x^84 + x^82 - x^80 - 2*x^78 - x^76 + 3*x^72 + 4*x^68 - x^66 + x^64 - 4*x^62 + 3*x^58 + x^54 - 2*x^52 + 4*x^50 + 4*x^48 - x^46 - x^44 - 5*x^42 + 3*x^40 + 6*x^36 - 2*x^34 + x^32 + x^30 + 4*x^28 - x^26 - 4*x^24 - 2*x^22 + 2*x^20 + 7*x^18 - x^16 - x^14 + 2*x^12 + 3*x^10 + 2*x^8 - 8*x^6 + 9) (x^4 - x^2 + 1) * (x^46 + x^44 - x^40 - x^38 + 3*x^36 + 4*x^34 + x^32 - 3*x^30 - 2*x^28 + 3*x^26 + 5*x^24 + 2*x^22 - 2*x^18 - x^16 + 2*x^14 + 5*x^12 + 3*x^10 - x^8 - 3*x^6 + 4*x^2 + 4) (x^6 - x^5 + x^4 - x^3 + x^2 - x + 1) * (x^6 + x^5 + x^4 + x^3 + x^2 + x + 1) * (x^48 - x^46 + x^38 + x^36 - x^34 - x^32 + 3*x^28 - 3*x^26 + 2*x^24 + x^20 - x^18 - 2*x^16 + 3*x^14 - x^10 + x^8 + x^6 - 4*x^2 + 4) (x^8 + 1) * (x^90 - x^82 + 3*x^76 + x^74 - 3*x^68 - x^66 + 2*x^64 + 4*x^62 + 3*x^60 + x^58 - 2*x^56 - 4*x^54 + 2*x^52 - x^50 + 5*x^48 + 4*x^46 - 2*x^44 + 4*x^42 - x^40 - 4*x^38 + 2*x^36 - 2*x^34 + 6*x^32 + 4*x^30 + x^28 + 2*x^26 - 4*x^24 - 2*x^18 + 9*x^16 + 3*x^12 + 3*x^10 - 7*x^8 + x^6 + 9) (x^6 - x^3 + 1) * (x^6 + x^3 + 1) * (x^100 - x^94 + x^86 + 2*x^84 + x^82 - x^80 - 2*x^78 - x^76 + 3*x^72 + 4*x^68 - x^66 + x^64 - 4*x^62 + 3*x^58 + x^54 - 2*x^52 + 4*x^50 + 4*x^48 - x^46 - x^44 - 5*x^42 + 3*x^40 + 6*x^36 - 2*x^34 + x^32 + x^30 + 4*x^28 - x^26 - 4*x^24 - 2*x^22 + 2*x^20 + 7*x^18 - x^16 - x^14 + 2*x^12 + 3*x^10 + 2*x^8 - 8*x^6 + 9) |
#function PlotRootsF gives a point plot of the roots for a given N
#Note: If the polynomial were not squarefree, this would loop forever
from sage.rings.polynomial.complex_roots import complex_roots
def PlotRootsF(n):
L=len(complex_roots(F(n,x),skip_squarefree=True))
l=[]
for j in range(L):
l.append((real(complex_roots(F(n,x),skip_squarefree=True)[j][0]).center(), imaginary(complex_roots(F(n,x),skip_squarefree=True)[j][0]).center()))
return point(l,rgbcolor=hue(1),size=10)
|
|
#The next three lines show the zeros in relation to the unit circle
#Note: Many of the roots are only slightly inside or outside the circle. Hence a numeric check is required for an accurate answer to part b.
p1=PlotRootsF(5)
p2=circle((0,0),1)
p1+p2
|
|
p1=PlotRootsF(6)
p2=circle((0,0),1)
p1+p2
|
|
p1=PlotRootsF(7)
p2=circle((0,0),1)
p1+p2
|
|
|
|