##########################################################################
## Lehmerconj.txt Save this file as  Lehmerconj.txt                     ##
## to use it, stay in the same directory, get into Maple                ##
## (by typing: maple <Enter> ) and then type:                           ##
## read Lehmerconj.txt <Enter>                                          ##
## Then follow the instructions given there                             ##
##                                                                      ##
## Written by Yuxuan Yang, Doron Zeilberger and Zidong Zhang,           ##
## Rutgers University                                                   ##
## yy458@math.rutgers.edu                                               ##
## doronzeil@gmail.com                                                  ##
## zz430@rutgers.edu                                                    ##
##########################################################################
Help:=proc():
print(`tauF(r,k),pmodseqzero(r,N),FZ(r,N),taup(N),taupcheck(N,k)`):
end:



#Compute tau_r(k) using J.C.P. Miller recurrence.
tauF:=proc(r,k) local s,j:
option remember:
if k<0 then
 RETURN(0):
fi:
if k=0 then
  RETURN(1):
fi:
s:=0:
for j from 1 while (3*j-1)*j/2<=k do
 s:=s+(-1)^j*((r+1)*(3*j-1)*j/2-k)*tauF(r,k-(3*j-1)*j/2):
od:
for j from 1 while (3*j+1)*j/2<=k do
 s:=s+(-1)^j*((r+1)*(3*j+1)*j/2-k)*tauF(r,k-(3*j+1)*j/2):
od:
s/k:
end:

#Count the number of prime p<=N such that tau_r(p) module p equal 0
pmodseqzero:=proc(r,N) local k,j,zeros:
for j from 1 while ithprime(j)<=N do
od:
zeros:=0:
for k from 1 to j do
if tauF(r,ithprime(k)-1) mod ithprime(k) = 0 then
zeros:=zeros+1:
fi:
od:
zeros:
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 i:
for i from 1 to N while tauF(r,i)<>0 do 
od:
if i=N+1 then
 infinity:
else
i:
fi:

end:



#Outputs the sequence tau(n) for n from 1 to N
taup:=proc(N) local i,tau,known,knowni,s,j,m,n,taui:
known:={}:
tau:=Array([0$N]):
tau[1]:=1:
for i from 2 to trunc(sqrt(sqrt(N))) do 
 knowni:={}:
 if tau[i]=0 then 
  s:=0:
  for j from 1 while (3*j-1)*j/2<=i-1 do
   s:=s+(-1)^j*((24+1)*(3*j-1)*j/2-i+1)*tau[i-(3*j-1)*j/2]:
  od:
  for j from 1 while (3*j+1)*j/2<=i-1 do
   s:=s+(-1)^j*((24+1)*(3*j+1)*j/2-i+1)*tau[i-(3*j+1)*j/2]:
  od:
  tau[i]:=(s/(i-1)):
  knowni:=knowni union {i}:
  for m from 2 while i^m<=N do
   tau[i^m]:=tau[i]*tau[i^(m-1)] - i^11 * tau[i^(m-2)]:
   knowni:=knowni union {i^m}
  od:
  for m from 1 while i^m<=N do
   for n in known while n * i^m <=N do
    tau[n * i^m]:=tau[n]*tau[i^m]:
    knowni:=knowni union {n * i^m}:
   od:
  od:
  known:=known union knowni:
 fi:
od:
for i from trunc(sqrt(sqrt(N)))+1 to N do
 if tau[i]=0 then 
  taui:=[1]: 
  s:=0:
  for j from 1 while (3*j-1)*j/2<=i-1 do
   s:=s+(-1)^j*((24+1)*(3*j-1)*j/2-i+1)*tau[i-(3*j-1)*j/2]:
  od:
  for j from 1 while (3*j+1)*j/2<=i-1 do
   s:=s+(-1)^j*((24+1)*(3*j+1)*j/2-i+1)*tau[i-(3*j+1)*j/2]:
  od:
  taui:=[op(taui),s/(i-1)]:
  for m from 2 while i^m<=N do
   taui:=[op(taui),taui[2]*taui[m]-i^11*taui[m-1]]:
  od:
  for m from 1 while i^m<=N do
   for n from 2 while n * i^m<=N do
    tau[n * i^m]:=tau[n]*taui[m+1]:
   od:
  od:
  for m from 1 to nops(taui)-1 do
   tau[i^m]:=taui[m+1]:
  od:
 fi:
od:
tau:
end:

#find all prime p<=N such that tau(p)(mod p) is k
taupcheck:=proc(N,k) local i,tau,known,knowni,s,j,m,n,taui:
known:={}:
tau:=Array([0$N]):
tau[1]:=1:
for i from 2 to trunc(sqrt(sqrt(N))) do 
 knowni:={}:
 if tau[i]=0 then 
  s:=0:
  for j from 1 while (3*j-1)*j/2<=i-1 do
   s:=s+(-1)^j*((24+1)*(3*j-1)*j/2-i+1)*tau[i-(3*j-1)*j/2]:
  od:
  for j from 1 while (3*j+1)*j/2<=i-1 do
   s:=s+(-1)^j*((24+1)*(3*j+1)*j/2-i+1)*tau[i-(3*j+1)*j/2]:
  od:
  tau[i]:=(s/(i-1)):
  if tau[i] mod i = k mod i then 
   print(i):
  fi:
  knowni:=knowni union {i}:
  for m from 2 while i^m<=N do
   tau[i^m]:=tau[i]*tau[i^(m-1)] - i^11 * tau[i^(m-2)]:
   knowni:=knowni union {i^m}
  od:
  for m from 1 while i^m<=N do
   for n in known while n * i^m <=N do
    tau[n * i^m]:=tau[n]*tau[i^m]:
    knowni:=knowni union {n * i^m}:
   od:
  od:
  known:=known union knowni:
 fi:
od:
for i from trunc(sqrt(sqrt(N)))+1 to N do
 if tau[i]=0 then 
  taui:=[1]: 
  s:=0:
  for j from 1 while (3*j-1)*j/2<=i-1 do
   s:=s+(-1)^j*((24+1)*(3*j-1)*j/2-i+1)*tau[i-(3*j-1)*j/2]:
  od:
  for j from 1 while (3*j+1)*j/2<=i-1 do
   s:=s+(-1)^j*((24+1)*(3*j+1)*j/2-i+1)*tau[i-(3*j+1)*j/2]:
  od:
  taui:=[op(taui),s/(i-1)]:
  if taui[2] mod i = k mod i then 
   print(i):
  fi:
  for m from 2 while i^m<=N do
   taui:=[op(taui),taui[2]*taui[m]-i^11*taui[m-1]]:
  od:
  for m from 1 while i^m<=N do
   for n from 2 while n * i^m<=N do
    tau[n * i^m]:=tau[n]*taui[m+1]:
   od:
  od:
  for m from 1 to nops(taui)-1 do
   tau[i^m]:=taui[m+1]:
  od:
 fi:
od:
end: