#Homework hw17 by  George Spahn

with(NumberTheory):
#CheckMordell(N) verifies that τ(p)*τ(q)=τ(pq) for all primes p,q, between 2 and N 

CheckMordell:=proc(N) local rs,i,j,a :
a:=5000:
if N*N < a then:
 a:=N*N:
fi:
rs:=RseqFF(24, a):

for i from 2 to N do
 for j from 2 to N do
  if i*j >= a then
   break:
  fi:
  if AreCoprime(i,j)  then
   if rs[i]*rs[j]<> rs[i*j] then
    print("something went wrong"):
   fi:
  fi:
 od:
od:  
print("all good"):
end: 

#CheckDeligne(N) outputs the prime between 2 and N such that |tau(p)|/(2 p^11/2)
# is as small as possible, followed by its valu

CheckDeligne:=proc(N) local rs,a,j,bestp,bestv,v :
a:=min(N,5000):
rs:=RseqFF(24, a):
bestp:=1:
bestv:=1:
for j from 1 to a do
 if isprime(j) then 
  v := evalf(abs(rs[j])/(2*j^(11/2))):
  if v<bestv then
   bestp:=j:
   bestv:=v:
  fi:
  
 fi:
od:
bestp,bestv:
end:

# For CheckDeligne(5000) I got 2927 was the best prime with a value of 
# 0.0004943967857
# The fact that 2927 is on the order of 5000 suggests that increasing N
# will yeild new champions.



#seq(FZ(r, 5000), r = 1 .. 30);
# 3, 7, 2, 9, 1560, 5, infinity, 3, infinity, 6, infinity, 

#   infinity, infinity, 4, 53, infinity, infinity, infinity, 

#   infinity, infinity, infinity, infinity, infinity, infinity, 

#   infinity, 9, infinity, infinity, infinity, infinity

# the sequence of non-infinities does not appear to be in the oeis
# if we have proofs for all the infinities up to say N, i would say the
# sequence up to N deserves to be in the OEIS





#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:

#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: