q := 2;
d := 2;
r := 5;
u := 2;

//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; 


// Count how many pairs of A,Z gives an expansion Z,ZA,..., ZA^(r(d-1)) of full dimension r(d-1).
// A is genrated according to the parameters, q,d,r given with a structure described in equation (24).
// It can be easily modified to see the whole "spectrum", i.e. how many tests gave us some space of dimension: n,n-1,...,1.
Test := procedure(q,r,u,d,trials)
// Set parameters
n := (d-1)*r;
Fq := GF(q);
VecFqn := VectorSpace(Fq, n);

D := DiagonalMatrix(Fq,[1: i  in [1..n-r]]);
D := HorizontalJoin(D, ZeroMatrix(Fq,n-r,r));


counter := [0:i in [0..n]];

//SetSeed(11);
for trial in [1..trials] do
	// Generate the first r rows of A as random.
	A := Matrix(Fq, r,n, [Random(Fq): i in [1..r*n]]);
	// Add the subsequent rows (I|0).
	A := VerticalJoin(A,D);
	
	// Generate a random matrix z having u rows.
	z := Matrix(Fq, u,n, [Random(Fq): i in [1..u*n]]);
	
	// The space generated by rows of Z.
	Row_Sample := sub<VecFqn | Rows(z)>;
	
	// Expanding the space RowSpan(Z).
	for i in [1..n] do
		Row_Sample := Row_Sample + sub<VecFqn | Rows(z*(A^i))>;
	end for;
	
	// Increase the counter at the corresponding dimension achieved by the above process.
	counter[Dimension(Row_Sample) + 1] +:= 1; 
end for;	

//Compare with theoretical expected values.
// Expected value := (Number of Attempts)*(Probability of full Expansion).
print("\n\n\n");
printf "Parameters: q = %o, r = %o, u = %o, d = %o \n", q,r,u,d;
// Lower bound for the probability in Corollary 1.
print("\n\nLower bound Uniform: ");
trials*Real(1 - (q^(-u+1))*(q-1)^-1);

// Exact probability in Corollary 1.
print("\nExpected Uniform: ");
trials*Real(Hq(q,r*(d-1)+u-1)/Hq(q,u-1));

// Number of succesful expansions during the test.
print("\nMeasured: ");
counter[n+1];

// Probability if A was sampled from companion matrices of same size.
print("\nExpected bound Companion: ");
trials*Real((1 - (1/q)^u));
print("\n\n\n");

end procedure;


// An example of a test.
q:=2;
r:=3;
u:=3;
d:=9;

Test(q,r,u,d,1000);