//Hq function.
// Compute (1-q^-1)(1-q^-2)...(1-q^-n)
Hq := function(q,n)
	res := 1;
	for i in [1..n] do
		res := res*(1-q^(-i));
	end for;
	return res;
end function; 

// Initial parameters.
q := 2;
d := 4;
r := 5;
u := 3;
// Dimension of the code
k := 10;

// n chosen such that n-k = r+u.
n := r+u+k; 
// How many interation we will perform in the expanding phase, corresonds to t in the paper.
iter := r*(d-1);
// Value of m must be big enough.
m := 2*(d + r*(d-1))*r + 20; 

// Number of trials.
trials:=100;


counter := 0;
counter_xp := 0;

// Field Fqm and Fq.
Fqm<x> := GF(q^m);
Fq := GF(q);

// Fqm seen as Fq^m.
VecFqm := VectorSpace(Fq, m);


// Rank Support of dimension d of the parity check H generated by some elemnt alpha. (V_{alpha,d} in the paper).
alpha := Random(Fqm);
SuppH := sub<VecFqm | [VecFqm ! ElementToSequence(alpha^i): i in [0..d-1]]>;
// Check V_a,d has indeed dimension d, if not choose a different alpha.
while( Dimension(SuppH) lt d) do
	alpha := Random(Fq^m);
    SuppH := sub<VecFqm | [VecFqm ! ElementToSequence(alpha^i): i in [0..d-1]]>;
end while;

// Parity check matrix.
H := Matrix(Fqm, n-k, n, [Fqm ! ElementToSequence(Random(SuppH)): i in [1..(n-k)*n]]);
// Check H hs indeed rank n-k, so the code has indeed dimension k.
// If not generate another H.
while( Rank(H) lt n-k) do
	H := Matrix(Fqm, n-k, n, [Fqm ! ElementToSequence(Random(SuppH)): i in [1..(n-k)*n]]); 
end while;

// Test the code with different errors and keep count of how many succesful support recovery we have.
for trial in [1..trials] do

// Create a subspace  SuppE of (Fq)^m of dimension r. The support of the error.
SuppE := sub<VecFqm | [Random(VecFqm): i in [1..r]]>;
// Check column support is really of dimension r.
while( Dimension(SuppE) lt r) do
	SuppE := sub<VecFqm | [Random(VecFqm): i in [1..r]]>;
end while;
// Basis SuppE as elements of Fqm.
BasisSuppE := [Fqm ! ElementToSequence(e): e in Basis(SuppE)];

// Create an error having support in SuppE.
e := Matrix(Fqm, 1, n, [Fqm ! ElementToSequence(Random(SuppE)): i in [1..n]]);

// Compute the relative syndrome.
s := e*Transpose(H);

// Syndrome Support.
SuppS := sub<VecFqm | [VecFqm ! ElementToSequence(s[1][i]): i in [1..n-k]]>;
// Basis SuppS as elements of Fqm.
BasisSuppS := [Fqm ! ElementToSequence(v): v in Basis(SuppS)];

// Expand S.
SuppS_xp := SuppS + sub<VecFqm | [VecFqm ! ElementToSequence((alpha^i)*s): s in BasisSuppS, i in [1..iter]]>;
// Basis of the expanded syndrome subspace as elements of Fqm.
BasisSuppS_xp := [Fqm ! ElementToSequence(v): v in Basis(SuppS_xp)];

// Recover the support of E from the expanded syndrome space.
// Intersect SuppS_xp and a^{-t-d+1}SuppS_exp.
SuppE_rec := SuppS_xp meet sub<VecFqm | [VecFqm ! ElementToSequence((alpha^(-iter-d+1))*s): s in BasisSuppS_xp]>;

// Check if the recovered space is indeed equal to the original error support.
if SuppE eq SuppE_rec then
	counter := counter + 1;
end if;

//Keep trace of how many succesful attempts we got so far out of the total number of attempts.
// Comment out to speed a bit up.
printf "%o out of %o\n", counter, trial;
end for;



// Prob that the  M_t(A,Z) has full rank 
print("\nExpected Uniform: ");
trials*Real(Hq(q,r*(d-1)+u-1)/Hq(q,u-1));

// Prob that both M_t(A,Z) has full rank  and X_1 is full rank.
// That should be the value of succesful recoveries we expect to observe.
print("\nExpected recoveries in the uniform case (d=2): ");
trials*Real(Hq(q,r*(d-1)+u-1)/Hq(q,u-1) * (1-q^(-u)*(q-1)^-1));

// Observed succesful recoveries.
print("\nObserved succesful recoveries: ");
counter;