q:=256;
r:=2;
s:=r+3;
u:=r+4;
v:=r+4;
m:=s*(u+v);
n:=r*(r+8);
F<w>:=GF(q);

Pol<[x]>:=PolynomialRing(F,n);
Poly<[y]>:=PolynomialRing(Pol,n);
Poll<[z]>:=PolynomialRing(Poly,m);

//Linear map T --------------------------------------------
T:=[];
repeat
	MT:=RandomMatrix(F,n,n);
until IsInvertible(MT) eq true;

for loop:=1 to n do
	T[loop]:=Pol!0;
	for i:=1 to n do
		T[loop]:=T[loop]+MT[loop][i]*x[i];
	end for;
end for;

//Linear map S -------------------------------------------------
S:=[];
repeat MS:=RandomMatrix(F,m,m);
until IsInvertible(MS) eq true;

for loop:=1 to m do
	S[loop]:=Poll!0;
	for i:=1 to m do
		S[loop]:=S[loop]+MS[loop][i]*z[i];
	end for;
end for;

// Central map F ----------------------------------------------------
A:=ZeroMatrix(Poly,s,r);
B:=ZeroMatrix(Poly,r,u);
C:=ZeroMatrix(Poly,r,v);
BP:=ZeroMatrix(Pol,r,u);
CP:=ZeroMatrix(Pol,r,v);

/*
for i:=1 to s do
	for j:=1 to r do
		A[i][j]:=y[(i-1)*r+j];
	end for;
end for;
*/

/*for i:=1 to s do
	for j:=1 to r do
		A[i][j]:=y[Random(1,n)];
	end for;
end for;
*/

for i:=1 to s do
	for j:=1 to r do
		for k:=1 to n do
			A[i][j]:=A[i][j]+Random(F)*y[k];
		end for;
	end for;
end for;

for i:=1 to r do
	for j:=1 to u do
		for k:=1 to n do
			a:=Random(F);
			B[i][j]:=B[i][j]+ a*y[k];
			BP[i][j]:=BP[i][j]+a*x[k];
		end for;
	end for;
end for;

for i:=1 to r do
	for j:=1 to v do
		for k:=1 to n do
			b:=Random(F);
			C[i][j]:=C[i][j]+b*y[k];
			CP[i][j]:=CP[i][j]+b*x[k];
		end for;
	end for;
end for;

E1:=A*B;
E2:=A*C;

Q:=[];
for i:=1 to s do
	for j:=1 to u do
		Q:=Q cat [E1[i][j]];
	end for;
end for;
for i:=1 to s do
	for j:=1 to v do
		Q:=Q cat [E2[i][j]];
	end for;
end for;

MQ:=[];
for loop:=1 to m do
	MQ[loop]:=ZeroMatrix(F,n,n);
	for i:=1 to n do
		for j:=1 to n do
			if i ne j then
				MQ[loop][i][j]:=MonomialCoefficient(Q[loop],y[i]*y[j]);
			end if;
		end for;
	end for;
end for;

// Public key computation
for loop:=1 to m do
	for i:=1 to n do
		Q[loop]:=Evaluate(Q[loop],y[i],T[i]);
	end for;
end for;

P:=S;
Pk:=[];
for loop:=1 to m do
	for i:=1 to m do
		P[loop]:=Evaluate(P[loop],z[i],Q[i]);
	end for;
	Pk[loop]:=MonomialCoefficient(MonomialCoefficient(P[loop],1),1);
end for;

// Output ----------------------------------------------------------------------------

printf "Write public_key.txt \n \n";
SetOutputFile("public_key.txt":Overwrite:=true);
printf "q:= %o; \n \n", q;
printf "r:= %o; \n \n", r; 
printf "s:= %o ; \n \n", s;
printf "u:= %o ; \n \n", u;
printf"v:= %o ; \n \n", v;
printf "n:= %o; \n \n", n;
printf "m:=s*(u+v); \n \n";
printf "F<w>:=GaloisField(q); \n \n";
printf "Pol<[x]>:=PolynomialRing(F,n+r*s); \n \n";
printf "Pk:= %o ; \n \n",Pk;
UnsetOutputFile();
 
printf "Write private_key.txt \n \n";
SetOutputFile("private_key.txt":Overwrite:=true);
printf "q:= %o; \n \n", q;
printf "r:= %o; \n \n", r; 
printf "s:= %o ; \n \n", s;
printf "u:= %o ; \n \n", u;
printf"v:= %o ; \n \n", v;
printf "n:= %o; \n \n", n;
printf "m:=s*(u+v); \n \n";
printf "F<w>:=GaloisField(q); \n \n";
printf "Pol<[x]>:=PolynomialRing(F,n+r*s); \n \n";
printf "Poly<[y]>:=PolynomialRing(Pol,n); \n \n";
printf "S:= %o ; \n \n", Eltseq(MS);
printf "T:= %o ; \n \n", Eltseq(MT);
printf "A:= %o ; \n \n", Eltseq(A);
printf "B:= %o ; \n \n", Eltseq(BP);
printf "C:= %o ; \n \n", Eltseq(CP);
UnsetOutputFile();