# OK to post

# PART 1: MakeCubic() PROCEDURE:
# -----------------------

MakeCubic:=proc(L,n) local func:
	func:=proc(n):
		return 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;

	return func;
end;

# PART 2: EncodeDES() PROCEDURE:
# -----------------------

EncodeDES:=proc(M,K,r) local i,func,C:
	if nops(M) mod 2 <> 0 then
		return FAIL;
	fi;

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

	C:=M;
	for i from 1 to r do:
		func:=MakeCubic(K[20*(i-1)+1..20*i], n);
		C:=Feistel(C,func);
	od;

	return C;
end;

# PART 3: DecodeDES() PROCEDURE:
# -----------------------

DecodeDES:=proc(M,K,r) local i,j,func,C:
	if nops(M) mod 2 <> 0 then
		return FAIL;
	fi;

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

	C:=M;
	for i from 1 to r do:
		j:=r-i+1;
		func:=MakeCubic(K[20*(j-1)+1..20*j], n);
		C:=revFeistel(C,func);
	od;

	return C;
end;

############################################################## OLD CODE ##############################################################

#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 L and outputs the integer n whose binary representation
#it is
BinToIn:=proc(L) local k,i:
k:=nops(L):
#[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 neccessary 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
RETURN([1,op(InToBin(n-2^(k-1),k-1))]):
fi:

end:

#BinFun(F,L): conmverts 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)]
#For example Feistel([1,1,0,1],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: