#OK to post homework
#Blair Seidler, 2021-03-28 Assignment 17

with(combinat):

Help:=proc():
print(`CheckMordell(N),CheckDeligne(N),LehmerHunt(R,N)`):
end:

# 2. 
#It seems that the builtin expand/bigpow is really efficient. I couldn't find any
#details of the implementation, but I also couldn't find a way to make JC1 take
#a measurable amount of time even with millions of terms.

# 3.
#CheckMordell(N): verifies that tau(p)*tau(q)=tau(pq) for all primes p,q, between 2 and N.
#Please use RseqFF(24,N), with N as large as your computer would agree to.
CheckMordell:=proc(N) local R,i,j:
R:=RseqFF(24,N^2): #Need almost N^2 terms to check the two largest primes <= N
for i from 1 while ithprime(i)<=N do
  for j from i+1 while ithprime(j)<=N do
    if R[ithprime(i)]*R[ithprime(j)]<>R[ithprime(i)*ithprime(j)] then
      print(`failed on`,i,j):
      RETURN(false):
    fi:
  od:
od:
true:
end:

#### Results: I was able to run CheckMordell(900), which generated the first 810000
#### coefficients of the Ramanujan Tau function. Every time I tried to generate the
#### first million coefficients with RseqFF, the Maple kernel crashed.

# 4.
#CheckDeligne(N): outputs the primes between 2 and N such that |tau[(p)|/2 p11/2 is
#as small as possible, followed by its value (that hopefully is less than 2). Who is
#the lucky prime? What is the value of |tau(p)|/2 p11/2 for that lucky prime?
CheckDeligne:=proc(N) local R,i,p,v,smallest,sprime,largest,lprime:
R:=RseqFF(24,N):
sprime:=2:
lprime:=2:
smallest:=evalf(abs(R[2])/(2*2^(11/2))):
largest:=smallest:
for i from 2 while ithprime(i)<=N do
  p:=ithprime(i):
  v:=evalf(abs(R[p])/(2*p^(11/2))):
  if v>1 then
    print(`failed on`,p,v):
    RETURN(FAIL):
  fi:

  if v>largest then
    largest:=v:
    lprime:=p:
  elif v<smallest then
    smallest:=v:
    sprime:=p:
  fi:
od:
printf(`Smallest value %f for prime %d\n`,smallest,sprime):
printf(`Largest value %f for prime %d\n`,largest,lprime):
end:

#If you make N larger, do you get new champions? Perhaps you can sharpen Deligne's theorem
#and replace 2 by something smaller?

#CheckDeligne(100);
#Smallest value 0.008883 for prime 43
#Largest value 0.854586 for prime 47

#CheckDeligne(1000);
#Smallest value 0.002282 for prime 907
#Largest value 0.959407 for prime 103

#CheckDeligne(10000);
#Smallest value 0.000320 for prime 7993
#Largest value 0.980267 for prime 7589

#CheckDeligne(100000);
#Smallest value 0.000134 for prime 17729
#Largest value 0.998084 for prime 74471

#CheckDeligne(500000);
#Smallest value 0.000010 for prime 130079
#Largest value 0.998583 for prime 119809

#CheckDeligne(800000);
#Smallest value 0.000010 for prime 130079
#Largest value 0.999193 for prime 792119

# I think that we would need to find the LARGEST possible value in order to improve the constant.
# At N=800000, the constant is already very close to 2 and seems to be growing. If I had to guess,
# I would guess that 2 cannot be improved.

# 5.
#LehmerHunt(R,N): By taking N fairly large (say 50000, or whatever your computer will allow)
#find the first 30 terms of FZ(r,50000), for r from 1 to 40. Discarding the (conjectural)
#infinities, is it in the OEIS? Should it be? Which r would give a (possible) Lehmer conjecture?
LehmerHunt:=proc(R,N) local r:
[seq(FZ(r,N),r=1..R)]:
end:

#### Results: This doesn't feel like it accomplished much ####
# LehmerHunt(40, 100000)
# [4, 8, 3, 10, 1561, 6, 28018, 4, infinity, 7,
#  infinity, infinity, infinity, 5, 54, infinity, infinity, infinity, infinity, infinity,
#  infinity, infinity, infinity, infinity, infinity, 10, infinity, infinity, infinity, infinity,
#  infinity, infinity, infinity, infinity, infinity, infinity, infinity, infinity, infinity, infinity]

#Not sure what "first 30 terms" for each r from 1..40 meant, since FZ returns the position of
#the first zero coefficient, not a sequence.

# [4, 8, 3, 10, 1561, 6, 28018, 4] does not seem to be in the OEIS. If this sequence is correct
# (of which I have doubts), it probably should be. I would like to be able to get more terms.

#### From C17.txt ####

#C17.txt: Maple code for Lecture 17

Help17:=proc(): print(`Rseq(r,N),RseqF(r,N),RseqFF(r,N),JC(P,x,m,K),JC1(P,x,m,K),FZ(r,N) `): 
end:

#Rseq(r,N): The first N coefficients of q*eta(q)^r
#where eta(q)=(1-q)*(1-q^2)*(1-q^3)*...
Rseq:=proc(r,N) local f,i:

f:=taylor(q*mul(1-q^i,i=1..N)^r,q=0,N+1):
[seq(coeff(f,q,i),i=1..N)]:

end:

#RseqF(r,N): a Fast version of Rseq(r,N) (see above)
RseqF:=proc(r,N) local f,i:
f:=1:
for i from 1 while (3*i+1)*i/2<=N do
 f:=f+(-1)^i*q^((3*i-1)*i/2) +(-1)^i*q^((3*i+1)*i/2):
od:

 
f:=taylor(q*f^r,q=0,N+1):
[seq(coeff(f,q,i),i=1..N)]:

end:

#JC(P,x,m,K): inputs a polynomial P in the variable x, a pos. integer m
#and a pos. integer K, outputs the coeff. of x^n from n=0 to n=K
#In other words taylor(P^m,x=0,K+1):
#P^m=a(m,0)+a(m,1)*x+a(m,2)*x^2+...
#P=p[0]+p[1]*x+...+p[L]*x^L
#P=3+2*x+5*x^2
#p[0]=3, p[1]=2, p[2]=5
#a(m,k)=1/(k*p0)*Sum(p[i]*((m+1)*i-k)*a(m,k-i),i=1..L):
JC:=proc(P,x,m,K) local L,p0,M,ng,k:
L:=degree(P,x):
p0:=coeff(P,x,0):
M:=[p0^m]:

for k from 1 to K do
#ng is a(m,k):
ng:= add(coeff(P,x,i)*((m+1)*i-k)*M[k-i+1],i=1..min(L,k))/(k*p0): 
M:=[op(M),ng]:
od:
  
end:

#added after class
#JC1(P,x,m,K): same input and output as JC(P,x,m,K), for checking purposes
#
JC:=proc(P,x,m,K) local f:
f:=expand(P^m):
[seq(coeff(f,x,i),i=0..K)]:
end:



#Added after class
#RseqFF(r,N): a yet Fast version of RseqF(r,N) (see above) using the J.C.P. Miller recurrence
RseqFF:=proc(r,N) local f,i,q:
f:=1:
for i from 1 while (3*i+1)*i/2<=N do
 f:=f+(-1)^i*q^((3*i-1)*i/2) +(-1)^i*q^((3*i+1)*i/2):
od:

 
JC(f,q,r,N-1):

end:

#FZ(r,N): The first place where the coefficient of q^i in q*eta(q)^r is zero looking at the first N coefficeints of q*eta(q)^r.
#if they are all non-zero, it returns infinity
FZ:=proc(r,N) local L,i:
L:=RseqFF(r,N):

for i from 1 to nops(L) while L[i]<>0 do od:

if i=nops(L)+1 then
 RETURN(infinity):
fi:

i:

end: