#OK to post homework
#Yuxuan Yang, April 9th, Assignment 21
with(combinat):

#1
JTseq:=proc(N) local n,La,Lb,L:
La:=aSeq(N):
Lb:=bSeq(N):
L:=[]:
 for n from 1 to N do
   if La[n]=0 then
    if n mod 2=1 or Lb[n/2]=0 then
     if IsSquareFree(n) then
      L:=[op(L),n]:
     fi:
    fi:
   fi:
 od:
L:  
end:



#2
IsSquareFree:=proc(n) local f,i:
f:=ifactors(n):
for i from 1 to nops(f[2]) do
 if f[2][i][2]<>1 then
  RETURN(false):
 fi:
od:
true:
end:
#JTseqW(N): All the terms <=N of the sequence of congruent numbers, those that are areas of right-angled triangle
#with rational sides. Shuold be called the Tunnell numbers; Using the Wikipedia version
JTseqW:=proc(N) local n,L:
#print(`Needs more work`):
L:=[]:
for n from 1 to N do
if IsCon(n) and IsSquareFree(n) then
   L:=[op(L),n]:
  fi:
od:
L:
end:



#4
#The main theorem of Sasa has a typo, sigma_i should be (3-i,4-i,...n+1-i) instead of (3-i,4-i,...n+2-i), 
#which gave me a hard time when debugging 

#Completeseq(m) implements the clever construction of Jerrold Tunnell's brilliant academic son Sasa Radomirovic
#and outputs the complete sequence with all permuations of [0,1,2,...,m-1] as a subsequence.
#Completeseq(13) will give the Example 9 in the paper.
Completeseq:=proc(m) local n,sig,i,j,k,x,y,med,tau,tau1,tau2,k0:
if m<7 then 
print(`m needs to be at least 7.`):
RETURN(FAIL):
fi:
n:=3*ceil((m-1)/3)-1:
sig:=[[seq(0..n-1)],[seq(0..n-1)],seq([seq(j-i mod n,j=3..n+1)],i=3..n-1),[seq(3..n-1),0,1,2],[seq(3..n-1),0,1,2]]:
x:=n:
y:=n+1:
med:=[seq(
op([   op(ins(sig[3*j-2],n-1,y)),   x,   op(sig[3*j-1]),   x,   op(ins(sig[3*j],2,y)),   x])
,j=2..(n+1)/3-1)]:
tau:=[x,op(sig[1]),y,x,op(sig[2]),x,y,op(sig[3]),x,op(med),op(sig[n-1]),y,x,op(sig[n]),x,y,op(sig[n+1]),x]:
#main theorem case (m=3k+1)
if m mod 3=1 then RETURN(tau): fi:
#one less element case (m=3k)
tau1:=[]:k0:=1:
for k from nops(tau) by -1 to 1 do
if tau[k]>n-1 then
 tau1:=[tau[k]-1,op(tau1)]:
elif tau[k]<n-1 then 
 tau1:=[tau[k],op(tau1)]:
else
 k0:=k:
fi:
od:
tau1:=[op(k0..-1,tau1)]:
if m mod 3=0 then RETURN(tau1): fi:
#two less elements case (m=3k-1)
tau2:=[]:k0:=1:
for k from nops(tau1) by -1 to 1 do
if tau1[k]>n-2 then
 tau2:=[tau1[k]-1,op(tau2)]:
elif tau1[k]<n-2 then 
 tau2:=[tau1[k],op(tau2)]:
else
 k0:=k:
fi:
od:
tau2:=[op(k0..-1,tau2)]:
if m mod 3=2 then RETURN(tau2): fi:
end:
#ins(sq,k,n):input a sequence sq, a position k and a number n, output the sequence obtained by inserting n to the position k
#of the sequence sq. e.g. ins([1,3,5],2,6)=[1,6,3,5]:
ins:=proc(sq,k,n) local i: [op(1..k-1,sq),n,op(k..nops(sq),sq)]:end:

#http://oeis.org/A062714 has more information

#Checkcomplete(sq,m):This is n! algorithm!!! Check sq is a complete sequence with all permuations of [0,1,2,...,m-1] as a subsequence
Checkcomplete:=proc(sq,m) local i,per:
for per in permute([seq(0..m-1)]) do
if not Checkperm(sq,per) then
RETURN(false):
fi:
od:
true:
end:
#Checkperm(sq,per):Check sq has per as a subsequence
Checkperm:=proc(sq,per) local i,j:
j:=1:
for i from 1 to nops(sq) do
 if sq[i]=per[j] then
  j:=j+1:
 fi:
 if j>nops(per) then
  RETURN(true):
 fi:
od:
false:
end:


#Maple code for Lecture 21
#C21: April 7, 2022. In fond memory of Jerrold B. Tunnell (1950-2022)
Help21:=proc(): print(`aSeq(N), bSeq(N), Tn(n,A,B,C) , JTseq(N), IsCon(n), JTseqW(N) `): end:

#aSeq(N): The first N terms of the sequence of a(n) defined by Jerrold Tunnell in Inv. Math. 72 (1983), 323-334
aSeq:=proc(N) local q,n,g,f,theta2,theta4:
g:=q*mul(1-q^(8*n),n=1..N)*mul(1-q^(16*n),n=1..N):
theta2:=1+2*add(q^(2*n^2)   , n=1..trunc(sqrt(N/2))+1):
#theta4:=1+2*add(q^(4*n^2)   , n=1..trunc(sqrt(N/2))+1):
f:=g*theta2:
f:=taylor(f,q=0,N+1):
[seq(coeff(f,q,n),n=1..N)]:
end:

#bSeq(N): The first N terms of the sequence of b(n) defined by Jerrold Tunnell in Inv. Math. 72 (1983), 323-334
bSeq:=proc(N) local q,n,g,f,theta2,theta4:
g:=q*mul(1-q^(8*n),n=1..N)*mul(1-q^(16*n),n=1..N):
#theta2:=1+2*add(q^(2*n^2)   , n=1..trunc(sqrt(N/2))+1):
theta4:=1+2*add(q^(4*n^2)   , n=1..trunc(sqrt(N/2))+1):
f:=g*theta4:
f:=taylor(f,q=0,N+1):
[seq(coeff(f,q,n),n=1..N)]:
end:

#JTseq(N): The first N terms of the sequence of congruent numbers, those that are areas of right-angled triangle
#with rational sides. Shuold be called the Tunnell numbers
#JTseq:=proc(N) local n,La,Lb,L:
#print(`Needs more work!`):

#La:=aSeq(N):
#Lb:=bSeq(N):
#L:=[]:
# for n from 1 to N do
#   if La[n]=0 then
#     L:=[op(L),n]:
#   fi:
# od:

#L:   
#  end:

#Tn(n,A,B,C): the set of triples [x,y,z] in Z^3 such that n=A*x^2+B*y^2+C*z^2

Tn:=proc(n,A,B,C) local x,y,z,S:
 S:={}:
 for z from -trunc(sqrt(n/C)) to trunc(sqrt(n/C)) do
 
    for x from -trunc(sqrt((n-C*z^2)/A)) to trunc(sqrt((n-C*z^2)/A)) do
 
        y:=sqrt((n-A*x^2-C*z^2)/B):
  
        if type(y,integer) then
           S:=S union {[x,y,z],[x,-y,z]}:
        fi:
    od:
 od:
S:
end:


#IsCon(n): Is n such that (mod. Birch-S-D) is the area of right-angled triangle with RATIONAL sides
IsCon:=proc(n):
#print(`Needs more work`):
if n mod 2=1 and 2*nops(Tn(n,2,1,32))=nops(Tn(n,2,1,8)) then
   RETURN(true):
 elif n mod 2=0 and  2*nops(Tn(n,8,2,64))=nops(Tn(n,8,2,16)) then
  RETURN(true):
fi:
false:
end: