@vv Is it possible to ask you to run your code for n=8 and s=10 and  f= x^8+x^4+x^3+x^2+1. I could not obtain no solutions. Thanks 

@vv I appreciate you taking the time to post this nice solution. I'm a math teacher and always I tell to my student that working with Maple is not only a skill, but also is an art of math thought. Thanks again 

@dharr I mean, because of we assumed that the polynomial f is an irreducible polynomial, then firstly we should check the singularity of A and after that its characteristic polynomial. In fact, as you said we could do the determinant test only when the polynomials match to improve efficiency (maybe the cost of determinant is less than the cost of characteristic polynomial ).

@dharr Thank you so much for your nice solutions. In fact, due to f is an irreducible polynomial f of degree n the integer number s satisfies s>=n and also the determinant of A should be non-zero. Unfortunately, I need the case n=8 and s=10 for some irreducible polynomials f of degree 8 over GF(2) :)

@Carl Love You answered not only this question, but also you solved my another question. You know, from your valuable answer, I realized that I did not know anything about Maple. 

@Carl Love Thanks for your nice solution. The factors are power of irreducible polynomials over F_2. I am waiting for your revise answer.

@dharr An interesting suggestion you shared. I try to check that whether the proposed method in the mentioned article, can be applied over finite fields with the characteristic two.  I appreciate you for your suggestion. 

@dharr Thank you for the response. I really appreciate you taking the time to write the recursive code in Maple. I have been working on MDS matrices over finite fields for two years and I know what you mean.  Thanks again professor Harrington. 

@dharr  Moreover, we can apply the simple command evalb(0 in convert(A^{-1}, list)) to check that whether in the matrix A^{-1} there is the entry zero. But the main problem is that the numbers of square sub-matrices are binomial(32, 16)-1 which can not be checked the singularity or adjoint of these sub-matrices with the proposed code. By the way, thanks for your comment. 

@Carl Love 

I would like to explain more about the application of MDS matrices. Assume that A is an n*n MDS matrix over F_q, the finite field with the q element where q is a prime power. Also suppose that I_n is an nxn identity matrix. Now consider an nx2n matrix G = (I_n | A) that is constructed by concatenating two nxn matrices I_n and A. Then, using G as a generator of a linear code, G produces an MDS code over  F_q. Moreover, MDS matrices are one of the most applicable matrices in the design of diffusion layers in block ciphers which can provide maximum diffusion.

I agree with you that the construction of  MDS matrices over F_{2^k} is desirable, but I suspect that the characteristic 2, substantially simplifies the computation. For example, suppose that A is an 8x8 matrix. It takes 3 seconds to verify that A is an MDS matrix over R(the field of real numbers ) and it takes 15 seconds to check that A is an  MDS matrix over F_{2^k} for some k. As I mentioned in reply to Kitonum, maybe the Algorithm 1 of the following paper can reduces the search space.

https://hal.inria.fr/hal-01506562/document 

Thanks for your notations. 

@Kitonum 

Thank you so much for your proposed clean procedure. I have a suggestion and want to know your idea about the proposed suggestion. 

Suppose that A is an n*n matrix. First we obtain the inverse of A, denoted with A^-1. If all entries of A^-1 are non-zero, then it means all determinant of square sub-matrices of A in size (n-1)*(n-1) are non-zero which results in there is no need to calculate these determinants and so on.

For details please see Algorithm 1 of the following paper
https://hal.inria.fr/hal-01506562/document

@Kitonum Ok, you right. Thanks again for your an interesting procedure. A strong point about your code is that the problem of Frobenius Coin problem can be solved with simplicity from your nice procedure.

@Kitonum Thanks for your valuable answer. I have a main question. Is it possible to ask you to explain with simplicity that why we should use the command "option remember" in our procedure. In fact, I would like to use the mentioned command in my programming but It is not clear to me that how to apply it. Thanks 

@vv I really appreciate your help. 

@vv 

I appreciate you taking the time to answer my questions with an interesting answer. 
I have one question. 
If [r1,r2,...,rn] is a list of non-negative integer numbers, then
how changes we shroud do in your second nice code until that we get
r1*i1+r2*i2+...+rn*in
In fact, we want to modify your code when  we have fixed coefficients. 
In other words, we have 

For i1 from 0 to k1 do
for i2 from 0 to k2 do
......................
for in from 0 to kn do
r1*i1+r2*i2+...+rn*in
end do; end do;...end do;
Thanks in advance