Code for Matrix Dynamics

Consider a set F={0,1}. Consider a set M of all possible 4 dimensional square matrices over the finite field F. Here we know 1+1=0, 0+1=1+0=1 and 0+0=0 and multiplications is straight forward. Clearly the cardinality of the set M is 2^16. Now define a map f defined on M as f(X)=((X*X)+C) Here '*' denotes the matrix modulo 2 multiplication and '+' denotes the matrix modulo 2 addition. It is good to name all the matrices by their decimal value. The matrix [0 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 0] would be called X_0. Similarly [ 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 1] as X_1 and so on. Also X_(2^16) is [1 1 1 1; 1 1 1 1; 1 1 1 1; 1 1 1 1].

We need to write a code which can takes each of these matrices and map to others as per definition of the map. Finally a list should be saved as an xls file.

Any help is highly is desired.