# HW 17 - Pablo Blanco. Spring 2024.
# OK to post

# returns the positive degree sequence of a polynomial's monomial terms
PolToList:=proc(P,x) local f,L,d:
d:=degree(P,x):
f:=P:
L:=[]:

while d >= 1 do:
 L:=[d,op(L)]:
 f:=f-coeff(f,x^d)*x^d:
 d:=degree(f,x):
od:
L:
end:

#1. Yes, MaxPer(d,x)=WisW(d,x)4] for all d=3..10.

MaxPer:=proc(d,x) local S,O,L,s:
 S:=PolsE(x,d) mod 2:
 O:={}:
 for s in S do: 
  L:=PolToList(s,x):
  if FindPer(SR([1,0$(d-1)],L,2^(d+1)))=2^d-1 then:
   O:=O union {s}:
  fi:
 od:
 O:
end:

#2.
# it should be true that
# qMaxPer(3,d,x)=qWisW(3,d,x)[4] but I can't find the bug on my qSR so that qMaxPer runs properly!
# after some testing, though, I think my qMaxPer should be bug-free so it only requires the correct qSR implementation.

#qPols(q,x,d): The SET of all polynomials of degree <=d in x mod q such that the constant term is 1
qPols:=proc(q,x,d) local S,s,c:
option remember:
if d=0 then
 RETURN({seq(c,c=1..q-1)}):
fi:
S:=qPols(q,x,d-1):
S union {seq(seq(s+c*x^d,s in S),c=1..q-1)}:
end:

#qPolsE(q,x,d): The set of polynomials of EXACTLY degree d
qPolsE:=proc(q,x,d) local S,s,c:
S:=qPols(q,x,d-1):
{seq(seq(s+c*x^d,s in S),c=1..q-1)}:
end:

#qIsIr(q,P): Is P irreducible mod q?
qIsIr:=proc(q,P)
option remember:
evalb(Factor(P) mod q=P mod q):
end:

#qIsPr(q,P,x): is the polynomial P in x of degree d, say, primitive?
qIsPr:=proc(q,P,x) local d,m,i:
option remember:
d:=degree(P,x):
m:=q^d-1:
 for i from 1 to m-1 do
   if rem(x^i,P,x) mod q=1 then
     RETURN(false):
   fi:
od:
if rem(x^m,P,x) mod q<>1 then
   print(`Something bad happened, Cocks' article is wrong!, or more likely (according to George)`):
   print(`we messed up`):
    RETURN(FAIL):
fi:
true: 
end:

#qWisW(q,d,x): inputs a pos. integer d and outputs the list of sets

qWisW:=proc(q,d,x) local S,s, Si, Sp, Sip, Sne:

S:=qPolsE(q,x,d):
Si:={}: Sp:={}: Sip:={}: Sne:={}:

for s in S do
 if qIsIr(q,s) and not qIsPr(q,s,x) then
   Si:=Si union {s}:
  elif  not qIsIr(q,s) and  qIsPr(q,s,x) then
    Sp:=Sp union {s}:
  elif   qIsIr(q,s) and  qIsPr(q,s,x) then
    Sip:=Sip union {s}:
   else
     Sne:=Sne union {s}:
  fi:
 od:

[Sne,Si,Sp,Sip]:
end:

qMaxPer:=proc(q,d,x) local S,O,L,s,i:
 S:=qPolsE(q,x,d):
 O:={}:
 for s in S do: 
  L:=PolToList(s,x):
  L:=[seq([coeff(s,x^i),L[i]], i=1..nops(L))]:

  if FindPer(qSR(q,[1,0$(d-1)],L,q^(d+1)))=q^d-1 then:
   O:=O union {s}:
  fi:
 od:
 O:
end:



############################################################################################# HW END


Help:=proc(): print(`Pols(x,d), PolsE(x,d), IsIr(P), IsPr(P,x), WisW(d,x) `): end:

#old code
#C16.txt, March 18, 2024
Help16:=proc(): print(`SR(INI,L,n), IsPer1(L,t), FindPer(L) `):end:

#SR(INI,L,n): inputs a 0-1 list of length L[nops(L)]
# a list of increasing positive integers L, and outputs
#the sequence of 0-1 generated by them of length n
#outputs the list of length n generated by the recurrence
#x[t]=x[t-L[1]]+...+x[t-L[nops(L)]
SR:=proc(INI,L,n) local r,i,M,ng:
r:=nops(L):
if nops(INI)<>L[-1] then
 RETURN(FAIL):
fi:

if not (convert(INI,set)={0,1} or convert(INI,set)={1}) then
 RETURN(FAIL):
fi:

if not (type(L,list) and {seq(type(L[i],posint),i=1..nops(L))}={true} and sort(L)=L) then
   RETURN(FAIL):
fi:
if not type(n,posint) then
     RETURN(FAIL):
fi:

M:=INI:

while(nops(M))<n do
 ng:=add(M[-L[i]],i=1..r) mod 2:
 M:=[op(M),ng]:
od:
M:

end:

#IsPer(L,t): Is L periodic of period t?
IsPer:=proc(L,t) local i:
if {seq(L[i]-L[i+t],i=1..nops(L)-t)}={0} then
  true:
else
 false:
fi:
end:

FindPer:=proc(L) local t:
for t from 1 to trunc(nops(L)/2) while not IsPer(L,t) do od:
if t=trunc(nops(L)/2)+1 then
RETURN(FAIL):
else
 RETURN(t):
fi:
end:

#end old code

#start new code for C17.txt
#Pols(x,d): The SET of all polynomials of degree <=d in x mod 2 such that the constant term is 1
Pols:=proc(x,d) local S,s:
option remember:
if d=0 then
 RETURN({1}):
fi:
S:=Pols(x,d-1):
S union {seq(s+x^d,s in S)}:
end:

#PolsE(x,d): The set of polynomials of EXACTLY degree d
PolsE:=proc(x,d) local S,s:
S:=Pols(x,d-1):
{seq(s+x^d, s in S)}:
end:

#IsIr(P): Is P irreducible mod 2?
IsIr:=proc(P)
option remember:
evalb(Factor(P) mod 2=P mod 2):
end:

#IsPr(P,x): is the polynomial P in x of degree d, say, primitive?
IsPr:=proc(P,x) local d,m,i:
option remember:
d:=degree(P,x):
m:=2^d-1:
 for i from 1 to m-1 do
   if rem(x^i,P,x) mod 2=1 then
     RETURN(false):
   fi:
od:
if rem(x^m,P,x) mod 2<>1 then
   print(`Something bad happened, Cocks' aritcle is wrong!, or more likely (according to George)`):
   print(`we messed up`):
    RETURN(FAIL):
fi:
true: 
end:

#WisW(d,x): inputs a pos. integer d and outputs the list of sets
#of pol. of degree exactly d such that
#(i) neither (ii) Irred. but not Primitive (iii) Primitive but not irreducible (iv) both
WisW:=proc(d,x) local S,s, Si, Sp, Sip, Sne:
#Si:=set of pol. of degree d (mod 2) that are irreducible but NOT primitive
#Sp:=set of pol. of degree d (mod 2) that are NOT irreducible but  primitive
#Sip:=set of pol. of degree d (mod 2) that are BOTH irreducible and  primitive
#Sne= neither
S:=PolsE(x,d):
Si:={}: Sp:={}: Sip:={}: Sne:={}:

for s in S do
 if IsIr(s) and not IsPr(s,x) then
   Si:=Si union {s}:
  elif  not IsIr(s) and  IsPr(s,x) then
    Sp:=Sp union {s}:
  elif   IsIr(s) and  IsPr(s,x) then
    Sip:=Sip union {s}:
   else
     Sne:=Sne union {s}:
  fi:
 od:

[Sne,Si,Sp,Sip]:
end: