#OK to post homework
#Lucy Martinez, 03-26-2024, Assignment 18


# Problem 1: Write a procedure MakeCubic(L,n) that inputs a binary string
# of length 20 and outputs the cubic polynomial in n
# BinToIn([op(1..5,L)]+BinToIn([op(6..10,L)]*n+BinToIn([op(11..15,L)]*n^2
# +BinToIn([op(16..20,L)]*n^3FindPer(SR(INI,L,n))]
MakeCubic:=proc(L,n) local i,B,p1,p2,p3:

if nops(L)<>20 then
 return(FAIL):
fi:

if not (convert(L,set)={0} or convert(L,set)={1} or convert(L,set)={0,1}) then
 return(FAIL):
fi:
BinToIn([op(1..5,L)])+BinToIn([op(6..10,L)])*n+BinToIn([op(11..15,L)])*n^2+BinToIn([op(16..20,L)])*n^3:
end:


# Problem 2: Write a procedure EncodeDES(M,K,r) that inputs a message M in binary
# (M has to be of even length) (the usal DES has M is of lenght 64),
# a key of length 20*r (a random binary string only known to you and the recipient),
# uses the key to generate r cubic functions
# (the usual DES has r=16, but the key-length is always 56, and it does not use cubics)
# uses the key to generate a list of r cubic polynomials [F1,...,Fr],
# and then applies the Feistel procedure r times, (composing F1, F2, ..., Fr)
# to the message M
EncodeDES:=proc(M,K,r) local i,j,F,f,L,newM:
if nops(M) mod 2 <>0 then
 return(FAIL):
fi:

if nops(K)<>20*r then
 return(FAIL):
fi:

F:=[]:
for i from 1 to r do
 L:=K[20*(i-1)+1..20*i]:
 F:=[op(F),MakeCubic(L,n)]:
od:

newM:=M:
for f in F do
 newM:=Feistel(newM,unapply(f,n)):
od:
newM:
end:



# Problem 3: Write a procedure DecodeDES(M,K,r) such that for ANY key of length 20*r, 
# and any message M we have DecodeDES(EncodeDES(M,K,r),K,r)=M
DecodeDES:=proc(M,K,r) local i,j,F,f,L,newM:
if nops(M) mod 2 <>0 then
 return(FAIL):
fi:

if nops(K)<>20*r then
 return(FAIL):
fi:

F:=[]:
for i from 1 to r do
 L:=K[20*(i-1)+1..20*i]:
 #to use revFeistel, we need to apply the last cubic function that was applied in 
 #the encoding process
 F:=[MakeCubic(L,n),op(F)]:
od:

newM:=M:
for f in F do
 newM:=revFeistel(newM,unapply(f,n)):
od:
newM:
end:


####################################### C18.txt
Help18:=proc():
 print(`RW(k), BinToIn(L),InToBin(n,k),BinFun(F,L),Feistel(LR,F),revFeistel(LR,F)`):
end:


#q=2
#RW(k): a random binary word of length k
RW:=proc(k) local ra,i:
ra:=rand(0..1):
[seq(ra(),i=1..k)]:
end:


#BinToIn(L): inputs a binary word and outputs the integer whose binary
# representation it is
#convert from binary to base 10
BinToIn:=proc(L) local k,i:
k:=nops(L):
#ex: [1,0,1]=1*2^2+0*2^1+1*2^0
add(L[i]*2^(k-i),i=1..k):
end:


#InToBin(n,k): inputs a non-neg integer n<2^k into its binary representation with
#0's in front if necessary so that it is of length k
InToBin:=proc(n,k) local i:
#option remember:
if (n>=2^k or n<0) or not (type(n,integer)) then
 return(FAIL):
fi:

if k=1 then
 if n=0 then
  return([0]):
 else
  return([1]):
 fi:
fi:

if n<2^(k-1) then
 return([0,op(InToBin(n,k-1))]):
else
#here we need to remove 2^(k-1)
 return([1,op(InToBin(n-2^(k-1),k-1))]):
fi:
end:


#BinFun(F,L): converts a function on the integers (mod 2^k) to a function
# on binary words
BinFun:=proc(F,L) local k:
k:=nops(L):
InToBin(F(BinToIn(L)) mod 2^k,k):
end:

#Feistel(LR,F): The Feistel transform that takes [L,R]->[R+F(L) mod 2, L]
#For example Feistel([1,1,0,1],n->n^5+n);
Feistel:=proc(LR,F) local k,L,R:
k:=nops(LR):
if k mod 2<>0 then
 return(FAIL):
fi:
L:=[op(1..k/2,LR)]:
R:=[op(k/2+1..k,LR)]:

[op(R + BinFun(F,L) mod 2) , op(L) ]:
end:


#revFeistel(LR,F): The reverse Feistel transform that takes [L,R]->[R,L+F(R) mod 2]
#For example revFeistel([0,0,1,0],n->n^5+n);
revFeistel:=proc(LR,F) local k,L,R:
k:=nops(LR):
if k mod 2<>0 then
 return(FAIL):
fi:
L:=[op(1..k/2,LR)]:
R:=[op(k/2+1..k,LR)]:

[op(R),op(L + BinFun(F,R) mod 2)]:
end:

####################################### C17.txt
Help17:=proc():
 print(`Pols(x,d),PolsE(x,d),IsIr(P),IsPr(P,x),WisW(d,x)`):
end:


#Pols(x,d): The set of all polynomials of degree <=d in x mod 2
# that have 1 (constant term) in it
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)} mod 2:
end:

#PolsE(x,d): The set of polynomials of exactly degree d
# i.e. we only take polynomials that are 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 primitive?
IsPr:=proc(P,x) local d,m,i:
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' article is wrong`):
 return(FAIL):
fi:
true:
end:


#WisW(d,x): inputs a pos. integer d and outputs the list of sets
# of polynomials of degree exactly d such that
# (i) neither
# (ii) irreducible but not primitive
# (iii) primitive but not irreducible
# (iv) both 
WisW:=proc(d,x) local S,s,Si,Sp,Sip,Sne:
#Si:=set of polys of degree d (mod 2) that are irreducible but not primitive
#Sp:=set of polys of degree d (mod 2) that are NOT irreducible but primitive
#Sip:=set of polys of degree d (mod 2) that are both irreducible and primitive
#Sne:=set of polys of degree d (mod 2) that are neither irreducible nor primitive
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 IsPr(s,x) and not IsIr(s) 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:


##################################################### from class:
Help16:=proc():
 print(`SR(INI,L,n),IsPer(L,t),FindPer(L)`):
end:

#SR=Shift Register from Clifford Cocks summary
#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 0s and 1s 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:
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(L): finds the period of the list L
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: